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:

1) 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$.

2) 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.

3) 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.

4) 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_1^+$, hence the function $\frac{\delta}{2n}(|x-e_n|^2 - \frac{1}{2})$ would lie above $u_R$ in $B_1^+$ for $R$ large, contradicting that $u \geq 0$.