Edit. Corrected a wrong statement about the structure of the functional $J$, thanks to Carlo's answer, and added more notes.
The major difficulty of the derivation, as noted by Deane Yang in his comment, is the structure of the second term of the functional \ref{1}, called $J_2$ in the following, i.e.
$$\DeclareMathOperator{\Dm}{d\!} \DeclareMathOperator{\Div}{\nabla\cdot}
J_2(u)=\int_\Omega \chi_{\{u>0\}}\Dm x
$$
Nevertheless we see that the following equation holds,
$$
J_2(u)=\int_\Omega \chi_{\{u>0\}}\Dm x=\int_\Omega H(u(x))\Dm x
$$
where $H(y)$, $y\in \Bbb R$ is the standard Heaviside function,
since we have that
$$
\chi_{\{u>0\}} = H(u(x))\quad \forall x\in \Omega,
$$
as it's easy to prove by elementary set theoretic considerations.
We can then try to use this relation and say that
$$
\begin{split}
\left.\frac{\Dm J}{\Dm \epsilon}(u+\epsilon \varphi)\right|_{\epsilon=0}
& = \left.\frac{\Dm}{\Dm \epsilon}\int_\Omega H\big(u(x)+\epsilon \varphi(x)\big)\Dm x \right|_{\epsilon=0} \\
& = -\int_\Omega H^\prime\big(u(x)+\epsilon\varphi(x)\big)\varphi(x)\Dm x \\
& = -\int_\Omega \varphi(x)\Dm \delta_{\partial\{u(x)>0\}} .
\end{split}
$$
where $\delta_{\partial\{u>0\}}$ is the Dirac measure supported on the boundary of subset of $\Omega$ where $u>0$. Here the derivative respect to the parameter $\epsilon$ as well as the gradient $\nabla H(u)$ (see below) must be intended in the sense of distributions.
Finally, let us approach the functional \ref{1}, or better its form considered by by Alt and Caffarelli (reference [1] of the question, p. 105):
$$
\begin{split}
J(u)=J_1(u)+J_2(u) & = \int_{\Omega\color{red}{\cap\{u>0\}}} |\nabla u|^2 \Dm x+ \int_{\Omega\color{black}{\cap{\{u > 0\}}}} \Dm x \\
& = \int_\Omega \color{red}{\chi_{\{u>0\}}}|\nabla u|^2\Dm x + \int_\Omega \chi_{\{u>0\}}\Dm x \\
& = \int_\Omega \color{red}{H(u(x))}|\nabla u(x)|^2\Dm x + \int_\Omega H(u(x))\Dm x
\end{split}
$$
and applying the procedure developed above and the standard du Bois-Reymond's lemma we get
$$
\begin{split}
\left.\frac{\Dm J}{\Dm \epsilon}(u+\epsilon \varphi)\right|_{\epsilon=0}
& = \frac{\Dm}{\Dm \epsilon}\left[\int_\Omega {H(u+\epsilon \varphi)}|\nabla (u+\epsilon \varphi)|^2\Dm x + \int_\Omega H(u+\epsilon \varphi)\Dm x\right]_{\epsilon=0 }\\
& = \left[-\int_\Omega {H^\prime(u+\epsilon \varphi)}|\nabla (u+\epsilon \varphi)|^2\varphi\Dm x \right. \\
& \qquad +\left.2 \int_\Omega H(u+\epsilon \varphi)\big(\nabla u\cdot\nabla\varphi+\epsilon|\nabla\varphi|^2\big)\Dm x\right]_{\epsilon=0 } \\
&\qquad\qquad - \int_\Omega \varphi(x)\Dm \delta_{\partial\{u(x)>0\}} \\
& = -\int_\Omega |\nabla u(x)|^2\varphi(x)\Dm \delta_{\partial\{u(x)>0\}} \\
& \qquad-\int_\Omega \varphi(x)\Dm \delta_{\partial\{u(x)>0\}} \\
& \qquad\qquad + 2\int_\Omega H(u)\nabla u\cdot\nabla\varphi\Dm x \\
& = -\int_\Omega \big[|\nabla u(x)|^2+ 1\big] \varphi(x)\Dm \delta_{\partial\{u(x)>0\}} + 2\int_\Omega H(u)\nabla u\cdot\nabla\varphi\Dm x \\
& = -\int_\Omega \big[|\nabla u(x)|^2+ 1\big] \varphi(x)\Dm \delta_{\partial\{u(x)>0\}} - 2\int_\Omega \big[\Div \big(H(u)\nabla u\big)\big]\varphi\Dm x \\
& = -\int_\Omega \big[|\nabla u(x)|^2+ 1\big] \varphi(x)\Dm \delta_{\partial\{u(x)>0\}} \\
&\qquad - 2\int_\Omega \big(H(u)\Delta u\big)\varphi\Dm x \\
&\qquad\qquad - 2\int_\Omega \big(\nabla H(u)\cdot\nabla u\big)\varphi\Dm x \\
& = -\int_\Omega \big[|\nabla u(x)|^2+ 1\big] \varphi(x)\Dm \delta_{\partial\{u(x)>0\}} \\
&\qquad - 2\int_\Omega \big(H(u)\Delta u\big)\varphi\Dm x \\
&\qquad\qquad + 2\int_\Omega H^\prime(u)\lvert\nabla u\rvert^2\varphi\Dm x \\
& = -\int_\Omega \big[|\nabla u(x)|^2+ 1\big] \varphi(x)\Dm \delta_{\partial\{u(x)>0\}} \\
& \qquad +2 \int_\Omega |\nabla u(x)|^2\varphi(x)\Dm \delta_{\partial\{u(x)>0\}} \\
& \qquad\qquad - 2\int_\Omega \big(H(u)\Delta u\big)\varphi\Dm x\\
& = \int_\Omega \big[|\nabla u(x)|^2 - 1\big] \varphi(x)\Dm \delta_{\partial\{u(x)>0\}} - 2\int_{\Omega\cap\{u(x)>0\}} (\Delta u)\varphi\Dm x\\
\\
& \implies u=\min J \iff
\begin{cases}
\Delta u (x)= 0 & x\in \Omega\cap{\{u(x)>0\}} \\
\left.\begin{split}
&u(x)=0 \\
&|\nabla u(x)| =1
\end{split}\right\} & x\in \Omega\cap{\partial\{u(x)>0\}} \\
\end{cases}
\end{split}
$$
Notes
- The result above coincides exactly with the one given by Alt and Caffarelli (reference [1] of the question, §0, p. 106 eq. (0.1)), when considering the functional $$ f(x, u,\nabla u)= \lvert\nabla u\rvert^2 + Q^2 $$ where $Q$ is a given function. In the case analysed here we have $Q^2\equiv 1$, but there absolutely no difference in the derivation respect to the general case.
- The only "modern" reference dealing with the work [1] of Alt and Caffarelli is the monograph [A1] written by their coAuthor Avner Friedman. Nevertheless he does not explicitly calculate the Euler-Lagrange equations of the functional involved as he use the direct method in the development of the analysis, as it is customary done for free boundary problems.
Reference
[A1] Avner Friedman, Variational principles and free-boundary problems, Pure and Applied Mathematics. A Wiley-Interscience Publication. New York: John Wiley & Sons, Inc. IX, 710 p. (1982)., MR679313, Zbl 0564.49002.