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)$$\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)