As Leo Moos suggested in the comments, in any dimension $d$, this is a simple consequence of the implicit function theorem.
The implicit function theorem implies that every point $x\in\phi^{-1}(0)$ has a neighborhood $\Omega_x$ that is a diffeomorphic image $\varphi_x(Q_{\delta_x})$ of a box $Q_{\delta_x}=\{|y_i|\leq\delta_x,\;1\leq i \leq d\}$, such that $\phi(\varphi(y))=|\nabla \phi(x)|\cdot y_1$, and such that $\varphi_x(0)=x$ and $D\varphi_x(0)$ is a rotation. (Details: assume wlog that $x=0$, and, by rotating, that $\nabla \phi (0)/|\nabla \phi (0)|=e_1,$ the first basis vector. Consider the function $g:\mathbb{R}^d\to \mathbb{R}^d$ defined by $g(x)=(\phi(x)/|\nabla \phi (0)|,x_2,\dots,x_d)$. Then $Dg(0)$ is the identity, and $\varphi_x$ is the inverse $g^{-1}$ provided by the inverse function theorem.) Given $\varepsilon>0$, by choosing smaller $\delta_x$ if necessary, we can ensure that the restriction of $\varphi_x$ to the leaves $\{y_1=h\}$ distorts the $(d-1)$-area by no more than $1+\varepsilon$, i. e., in $Q_{\delta_x}$,
$$1-\varepsilon<\left(\det_{2\leq i,j\leq d}\left(\partial_{y_i}\varphi\cdot\partial_{y_j}\varphi\right)\right)^\frac12\leq 1+\varepsilon.$$
By compactness, we can choose a finite cover $\Omega_i=\Omega_{x_i}$, $i=1,\dots,n$ of $\phi^{-1}(0)$. Then, for $\epsilon_0>0$ small enough, we have $W_{\epsilon_0}:=\{x:|\phi(x)|\leq\epsilon_0\}\subset \cup_{i=1}^n\Omega_i$. Put $\Omega_0=\Omega\setminus W_{\epsilon_0}$, then $\Omega_0,\Omega_1,\dots,\Omega_n$ is a finite cover of $\Omega$, and we can pick a partition of unity $f_0,\dots,f_n$ subordinate to that cover. We have for $|h|<\epsilon_0$,
$$
\mathcal{H}^{d-1}(\phi^{-1}(h))=\sum_{i=1}^n\int f_i d\mathcal{H}^{d-1}(\phi^{-1}(h)),
$$
hence, by changing the variable,
$$
(1-\varepsilon)\sum_{i=1}^n\int_{y_1=h} f_i\circ\varphi_i\leq \mathcal{H}^{d-1}(\phi^{-1}(h))\leq (1+\varepsilon)\sum_{i=1}^n\int_{y_1=h} f_i\circ\varphi_i.
$$
As $h\to 0$, the terms in the right-hand side tend to
$$
(1+\varepsilon)\int_{y_1=0}f_i\circ\varphi_i\leq (1+\varepsilon)^2\int f_i d\mathcal{H}^{d-1}(\phi^{-1}(0)),
$$
and similarly for the LHS. This means that
$$
(1-\varepsilon)^2\mathcal{H}^{d-1}(\phi^{-1}(0))\leq\liminf_{h\to 0} \mathcal{H}^{d-1}(\phi^{-1}(h)),
$$
$$
\limsup_{h\to 0} \mathcal{H}^{d-1}(\phi^{-1}(h))\leq (1+\varepsilon)^2\mathcal{H}^{d-1}(\phi^{-1}(0)),
$$
and since $\varepsilon>0$ is arbitrary, we are done.