Let $\Omega$ be a bounded or unbounded domain in $\mathbf R^{3}$ with a smooth boundary $S$ and a normal vector given by $n$. Now, we consider the following second-order elliptic problem with Neumann condition (in the variable $u$):

$$ \left \{\qquad \begin{aligned} \operatorname{div}(h \nabla u) & = f & \mbox{in}\, \Omega\\ h \tfrac{\partial u}{\partial n} & =g & \mbox{on} \, S \end{aligned}\right .$$

for some functions $f,g,h$ (it comes from studying a nonstationary problem for a nonhomogeneous incompressible fluid).

From a result of Agmon, Douglis and Nirenberg (in [$*$]) it is known that the problem is solvable in $W^2_p(\Omega)$ and the following estimate holds:

$$\Vert\nabla u(t)\Vert_{W_{p}^{1}(\Omega)}\leq K(\Vert h\Vert_{C^{1}\langle\overline{\Omega})})\bigl(\Vert f\Vert_{L_{p}(\Omega)}+\Vert g\Vert_{W_{p}^{1- 1/p}(S)}\bigr).$$

I do not understand the original proof in [$*$] for that inequality. Although it seems to be some sort of classical estimate for solutions of second-order elliptic partial differential equations, I couldn't find it in a modern book of elliptic equations, or derive it from other classical estimates. Is there any other recent reference for that estimate or does someone know how to deduce it?

**Reference**

[$*$] S. Agmon, A. Douglas and L. Nuremberg, "Estimates near the boundary for solution of elliptic partial differential equations satisfying general boundary conditions I", Communications on Pure and Applied Mathematics 12 (1959), 623-727, 10.1002/cpa.3160120405, MR0125307, Zbl 0093.10401.