$\newcommand{\ep}{\varepsilon}$
This is indeed a matter of splitting the integral. In polar coordinates, for some real $R>0$,
\begin{equation}
I(\ep) = \int_0^{2\pi}(J_t(\ep)+K_t(\ep))\,dt,
\end{equation}
where
\begin{equation}
J_t(\ep):=\int_0^{r_*}\frac{\sqrt{1+|\nabla h|^2}}{\big((h+\ep)^2+r^2\big)^2}\,r\,dr,
\end{equation}
\begin{equation}
K_t(\ep):=\int_{r_*}^R\frac{\sqrt{1+|\nabla h|^2}}{\big((h+\ep)^2+r^2\big)^2}\,r\,dr,
\end{equation}
$h:=h(r\cos t,r\sin t)$, $\nabla h:=\nabla h(r\cos t,r\sin t)$,

$r_*=r_*(\ep)>0$ varies with $\ep\downarrow0$ so that
$$\ep<<r_*<<\ep^{1/2}$$
(e.g., one may take $r_*=\ep^{3/4}$),
$a\ll b$ means $|a|=O(b)$, and $a<<b$ means $|a|=o(b)$.

We have $h\ll r^2$ and $|\nabla h|\ll r$, so that (as $\ep\downarrow0$) uniformly in $r\in[0,r_*]$ we have $h\ll r_*^2<<\ep$, $h+\ep\sim\ep$, $(h+\ep)^2+r^2\sim\ep^2+r^2$, $|\nabla h|<<1$, and $\sqrt{1+|\nabla h|^2}\sim1$. So,
\begin{equation}
J_t(\ep)\sim\int_0^{r_*}\frac{r\,dr}{\big(\ep^2+r^2\big)^2}
=\frac1{2\ep^2}\,\int_0^{r_*^2/\ep^2}\frac{ds}{(1+s)^2}
\sim\frac1{2\ep^2}\,\int_0^\infty\frac{ds}{(1+s)^2}
=\frac1{2\ep^2}. \tag{1}
\end{equation}
On the other hand,
\begin{equation}
K_t(\ep)\ll\int_{r_*}^\infty\frac{r\,dr}{(r^2)^2}\le\frac1{r_*^2}<<\frac1{\ep^2}.
\end{equation}
Thus, indeed
\begin{equation}
I(\ep)\sim\frac\pi{\ep^2}.
\end{equation}