Here is a post that I made on sci.math a while ago, regarding a method I have used.
Analysis of $we^w$
For $w\gt0$, $we^w$ increases monotonically from $0$ to $\infty$. When $w\lt0$, $we^w$ is negative.
Thus, for $x\gt0$, $\mathrm{W}(x)$ is positive and well-defined and increases monotonically.
For $w\lt0$, $we^w$ reaches a minimum of $-1/e$ at $w=-1$. On $(-1,0)$, $w e^w$
increases monotonically from $-1/e$ to $0$. On $(-\infty,-1)$, $w e^w$ decreases
monotonically from $0$ to $-1/e$. Thus, on $(-1/e,0)$, $\mathrm{W}(x)$ can have one of two values, one in $(-1,0)$ and another in $(-\infty,-1)$. The value in $(-1,0)$ is called the principal value.
The iteration
Using Newton's method to solve $we^w$ yields the following iterative
step for finding $\mathrm{W}(x)$:
$$
w_{\text{new}}=\frac{xe^{-w}+w^2}{w+1}
$$
Initial values of $w$
For the principal value, when $-1/e\le x\lt0$, and when $0\le x\le10$, use
$w=0$. When $x\gt10$, use $w=\log(x)-\log(\log(x))$.
For the non-principal value, when $x$ is in $[-1/e,-.1]$, use $w=-2$; and
if $x$ is in $(-.1,0)$, use $w=\log(-x)-\log(-\log(-x))$.
This says that for the branch you want, use the iteration with an initial value of $w=0$.
We can use a procedure known as "bootstrapping" to determine an approximation for the Lambert $W$ function. Let's go back to its definition.
For $x > 0$ the equation
$$
we^w = x
$$
has exactly one positive solution $w = W(x)$ which increases with $x$. Note that $(w,x) = (1,e)$ is one such solution, so if $x > e$ then $w > 1$. By taking logarithms of both sides of the equation we get
$$
\log w + w = \log x
$$
or
$$
w = \log x - \log w. \tag{1}
$$
When $x > e$ we therefore have
$$
w = \log x - \log w < \log x.
$$
In other words, our first approximation is that
$$
1 < w < \log x \tag{2}
$$
when $x > e$. We then have
$$
0 < \log w < \log\log x,
$$
and plugging this into $(1)$ yields
$$
\log x - \log \log x < w < \log x, \tag{3}
$$
where the left side is positive for $x > 1$. Taking logarithms as before yields
$$
\log\log x + \log\left(1 - \frac{\log\log x}{\log x}\right) < \log w < \log\log x,
$$
and upon substituting this back into $(1)$ we get
$$
\log x - \log\log x < w < \log x - \log\log x - \log\left(1 - \frac{\log\log x}{\log x}\right).
$$
Since $w = W(x)$ we have shown that
$$
\log x - \log\log x < W(x) < \log x - \log\log x - \log\left(1 - \frac{\log\log x}{\log x}\right) \tag{4}
$$
for $x > e$.
In your particular case we're interested in $W(e^{x+a})$, for which we have
$$
x+a - \log(x+a) < W(e^{x+a}) < x+a - \log(x+a) - \log\left(1 - \frac{\log(x+a)}{x+a}\right)
$$
for $x+a > 1$. In this sense we have
$$
W(e^{x+a}) \approx x+a - \log(x+a) = x\left(1 - \frac{\log(x+a) - a}{x}\right) \tag{5}
$$
when $x+a$ is large. Now by applying Taylor series a couple times we see that, for $x$ large and $a \ll x$,
$$
\begin{align}
\frac{\log x - a}{x+1} &= \frac{\log x - a}{x} \cdot \frac{1}{1+\frac{1}{x}} \\
&\approx \frac{\log x - a}{x} \left(1-\frac{1}{x}\right) \\
&= \frac{\log x - a}{x} - \frac{\log x - a}{x^2} \\
&= \frac{\log(x+a-a) - a}{x} - \frac{\log x - a}{x^2} \\
&= \frac{\log(x+a) + \log\left(1-\frac{a}{x+a}\right) - a}{x} - \frac{\log x - a}{x^2} \\
&= \frac{\log(x+a) - a}{x} + \frac{\log\left(1-\frac{a}{x+a}\right)}{x} - \frac{\log x - a}{x^2} \\
&\approx \frac{\log(x+a) - a}{x} - \frac{a}{x(x+a)} - \frac{\log x - a}{x^2} \\
&\approx \frac{\log(x+a) - a}{x}.
\end{align}
$$
We may then conclude from $(5)$ that
$$
W(e^{x+a}) \approx x \left(1 - \frac{\log x - a}{x+1}\right)
$$
for $x$ large and $a \ll x$.
Best Answer
Considering the bounds
$$\log (x)-\log (\log (x))+\frac12\frac{\log (\log (x))}{ \log (x)} < W(x)$$ $$W(x)< \log (x)-\log (\log (x))+\frac e{e-1 }\frac{ \log (\log (x))}{ \log (x)}$$ for a specific range, we can numerically minimize $$\Phi(k)=\int_a^b \Big[\log (x)-\log (\log (x))+k\frac{\log (\log (x))}{ \log (x)}-W(x)\Big]^2\,dx$$
For $a=5$ and $b=105$, this gives $k\sim 0.881076$ but the formula is a bit more complex than your.
Trying something of the same shape as your $$W(x) \sim a+b \log_e(x+c)$$ a nonlinear regression gives $(R^2>0.999999)$ $$\begin{array}{clclclclc} \text{} & \text{Estimate} & \text{Standard Error} & \text{Confidence Interval} \\ a & -0.201141 & 0.004453 & \{-0.209981,-0.192302\} \\ b & +0.774451 & 0.000974 & \{+0.772518,+0.776385\} \\ c & +2.288360 & 0.034073 & \{+2.220730,+2.356000\} \\ \end{array}$$ which leads to a maximum absolute error of $0.005$.
Congratulations for your idea !
Edit
Since we are using totally empirical models, let us try
$$W(x) \sim a+b \Big[\log_e(x^d+c)\Big]^f$$
$$\begin{array}{clclclclc} \text{} & \text{Estimate} & \text{Standard Error} & \text{Confidence Interval} \\ a & 0.256813 & 0.001340 & \{0.254152,0.259474\} \\ b & 0.815983 & 0.001525 & \{0.812955,0.819012\} \\ c & 0.547202 & 0.001079 & \{0.545059,0.549345\} \\ d & 0.678026 & 0.000656 & \{0.676723,0.679329\} \\ f & 1.172580 & 0.000318 & \{1.171940,1.173210\} \\ \end{array}$$ which reduces the previous sum of squares by a factor close to $800,000$ (!!) and leads to a maximum absolute error equal to $5\times 10^{-6}$. Better, isn't it ?