A more general form of the problem is: $$\Delta u = g$$ with $g > 0$, $u \geq 0$ and $u(x',0) = 0$ for all $x' \in \mathbb{R}^{n-1}$. Here are some things we can say:
If $0 < g < K$ then by $W^{2,p}$ estimates and Sobolev embedding we have an interior $C^{1,\alpha}$ estimate (for any $\alpha$), giving in particular the (scaling-invariant) estimate $$|\nabla u|_{B_{r/2}(x)} \leq \frac{C}{r}\|u\|_{L^{\infty}(B_r(x))}$$ for all $B_r(x) \subset \Omega$, where $C$ depends only on $K$.
If $0 \leq u \leq M$ for some $M$ apriori, then using barriers of the form $-C(|x-(x_0',-1)|^{2-n}-1)$ for any $x_0'$ and $C$ large depending on $M$ we easily obtain $$0 \leq u \leq C(M)x_n.$$ Of course, we must have this apriori bound since quadratic growth like $x_n^2$ is allowed otherwise.
Thus, 2) gives a boundary gradient bound, and when coupled with the scaling invariant interior estimate of 1) we get a full gradient bound.
As a remark, if $g$ is bounded we in fact get that $u_n$ is $C^{\alpha}$ up to the boundary by a boundary Harnack inequality of Krylov, which involves using the Harnack inequality to "improve" the trapping planes from 2) at smaller scales. Philosophically, ellipticity gives that the nice boundary and boundary data have some influence when we step in.
As another remark, if $g \geq \delta > 0$ then $u$ must have quadratic growth. Indeed, if not, the rescalings $u_R(x) = \frac{1}{R^2}u(Rx)$ satisfy the same conditions but go to $0$ on $B_2^+$, hence the function $\frac{\delta}{2n}(|x-e_n|^2 - \frac{1}{2})$ would lie above $u_R$ in $B_2^+$ for $R$ large, contradicting that $u \geq 0$.