**Edit:** *An earlier version of this answer missed a factor of two, now corrected, thanks to Daniele Tampieri.*

We seek the variation of the functional
$$L[u]=\int_\Omega\left(|\nabla u|^2+1\right)\chi_{u>0}\,dx.$$
The indicator function restricts the integration to the volume $V\subset\Omega$ where $u\geq 0$.     
<sub>It is essential that the indicator multiplies all terms, that was my initial mistake.</sub>

Introduce an infinitesimal variation $u(x)\mapsto u(x)+\epsilon(x)$, and compute the variation in $L$ to first order in $\epsilon$. Let $S$ be the surface boundary of $V$ on which $u=0$. The variation is
$$L[u+\epsilon]=\int_\Omega\left(|\nabla (u+\epsilon)|^2+1\right)\chi_{u+\epsilon>0}\,dx$$
$$=L[u]+\int_V 2(\nabla u)\cdot(\nabla\epsilon) \,dx + \int_{S} \left(|\nabla u|^2+1\right)\frac{\epsilon}{|\nabla u|}\,ds+{\cal O}(\epsilon^2)$$
$$=L[u]-2\int_V \epsilon\nabla^2 u \,dx +\int_{S} \epsilon\left(2n\cdot\nabla  u+\frac{|\nabla u|^2+1}{|\nabla u|}\right)\,ds+{\cal O}(\epsilon^2)$$
$$=L[u]-2\int_V \epsilon\nabla^2 u \,dx +\int_{S} \epsilon\left(-|\nabla  u|+\frac{1}{|\nabla u|}\right)\,ds+{\cal O}(\epsilon^2).$$
On the second line I used the identity
$$\int f(u)\frac{d}{du}\chi_{u>0}\,dx=\int_S f(u)|n\cdot\nabla u|^{-1}\,ds=\int_S f(u)|\nabla u|^{-1}\,ds,$$
with $n$ a unit vector normal to the surface $S$ and pointing outward. (Note that $\nabla u$ has only a normal component on $S$.) On the third line I carried out a partial integration of the volume integral, on the fourth line I used that $n\cdot\nabla u=-|\nabla u|$ on $S$.    

The variation of $L$ vanishes for arbitrary $\epsilon$ if the integrand of the volume integral $\int_V$ and the integrand of the surface integral $\int_{S}$ each vanish identically, so 
$$\nabla^2 u=0\;\;\text{in}\;\;V,$$
$$|\nabla u|^2=1\;\;\text{on}\;\;S.$$

---

More generally, following the same steps to minimize
$$L[u]=\int_V\left( f(x)|\nabla u|^2+h(u)+g(x)\right)\chi_{u>0}\,dx$$
gives the variation
$$\delta L[u]=\int_V\epsilon\biggl(-2\nabla\cdot(f(x)\nabla u)+h'(u)\biggr)\,dx+\int_S\epsilon\left(2f(x)(n\cdot\nabla u)+\frac{f(x)|\nabla u|^2+h(u)+g(x)}{|\nabla u|}\right)\,ds$$
hence the Euler-Lagrange equations
$$2\nabla\cdot(f(x)\nabla u)=h'(u)\;\;\text{in}\;\;V,$$
$$f(x)|\nabla u|^2=h(u)+g(x)\;\;\text{on}\;\;S.$$
*I note that this generalization differs from Eq. (2.8) of [this source](https://hal.sorbonne-universite.fr/hal-00881760/file/ABP28.pdf).*