Integral on level sets Let $g_\epsilon : K \subset \mathbb{R}^d \rightarrow \mathbb{R}$ (more regularity can be assumed if necessary) be defined on a compact set (with regular boundary) $K \subset \mathbb{R}^d$, and the familly $g_\epsilon$ is indexed by $\epsilon \in [0,1]$. $\|g_\epsilon - g_0\|_{L^\infty(K)} \le c \epsilon$ where $c$ does not depend on $\epsilon$, and $\|\nabla(g_\epsilon - g_0)\|_{L^\infty(K)} \le c \epsilon$. How to prove that
\begin{align*}
    \left| \int_{g_\epsilon^{-1}(\mu)} \nabla g_0 \cdot d n - \int_{g_0^{-1}(\mu)} \nabla g_0 \cdot d n \right| \le c \epsilon
\end{align*}
where $c$ does not depend on $\epsilon$ ? Also we assume that for any $x \in g_0^{-1}(\mu)$, $\nabla g_0 (x)\neq 0$
 A: I'll give a picture, though not entering into all details. Writing $g(\epsilon,x):=g_\epsilon(x)$ and $I:=]-\delta,\delta[$, a natural assumption is that
$g:I\times K\to\mathbb R$ is $C^1$, and $\mu$ is a regular value for $g$.
Therefore  $\Sigma:=g^{-1}(\mu)=\{(\epsilon,x)\in I\times K: g(\epsilon,x)=\mu\}$ is a $C^1$ co-dimension $1$ submanifold of $\mathbb{R}^{d+1}$. To avoid technicalities, let's also assume $\Sigma$ is contained in the interior of  $I\times K$.
Note that for any $y=(\epsilon,x)\in\Sigma$,   $T_y\Sigma=\ker dg(y)=\{(t,u)\in \mathbb{R}\times \mathbb{R}^{d}:t\,\partial_1 g(\epsilon,x)+(u\cdot\nabla g_\epsilon(x))=0\}$. So one has  $\nabla g_\epsilon(x)=0$  iff  $T_y\Sigma=\{0\}\times\mathbb{R}^{d}$  (note that $\nabla g_\epsilon(x)=0$ implies $\partial_1 g(\epsilon,x)\neq0$ because one has $dg(y)\neq0$ by the regularity assumption).
Let $P:\Sigma\ni (\epsilon,x)\mapsto \epsilon\in \mathbb R$, the restriction to $\Sigma$ of the projection onto the first coordinate. Note that the level sets of $P_\Sigma$, denoted $\{P_\Sigma=\epsilon\}$, are your integration domains $g_\epsilon^{-1}(\mu)$ (or, more precisely, $\{\epsilon\}\times g_\epsilon^{-1}(\mu)$. Moreover, by the above remark, for any  $(\epsilon,x)\in\Sigma$, one has  $\nabla g_\epsilon(x)=0$ iff   $(\epsilon,x)$ is a critical point of $P_\Sigma$ at level $\epsilon$.
Now you are in a position to apply the usual deformation lemma of the  Lusternik-Schnirelmann theory (check for instance  Variational Methods by M.Struwe, Thm 3.4 , Chap. II, where you will also find the details for construction of a pseudo-gradient field, in case $g$ is only assumed to be $C^1$), here to the $C^1$ function $P_\Sigma$. For any nbd $N$ of  the critical set of $P_\Sigma$ at level $0$, (which is $\{0\}\times\text{crit}(g_0)$, as said), you may break the integration domain as $\big(g_0^{-1}(\mu)\setminus N\big)\cup\big( g_0^{-1}(\mu)\cap N\big)$. The first set can be deformed into the level set  $\{P_\Sigma=\epsilon\}=g_\epsilon^{-1}(\mu)$. The other set gives a small contribute, since $\nabla g_0$ is small on $N$. This way you should be able to compare  your integral  $\int_{g_\epsilon^{-1}(\mu)}\nabla g_0dn$ with $\int_{g_0^{-1}(\mu)}\nabla g_0dn$. Note that to get a difference $O(\epsilon)$ you need to have the gradient flow $\eta^\epsilon(x)=x+O(\epsilon)$. Also, some care is needed to compare $N$ and $\epsilon$. (Assuming more regularity on $g$ may help, and may be necessary)
