How to estimate the integral involving the distance function

For small $t$ the dominant contribution comes from the boundary; starting from a point $y_0$ on the boundary and integrating inwards in a perpendicular direction we have $\int_{0}^\infty e^{-r^2/t}dr=\tfrac{1}{2}\sqrt{\pi t}$. Multiplication by the area $S$ of the boundary gives the estimate $$I\rightarrow \tfrac{1}{2}S\sqrt{\pi t}.$$ As a check, for a disc of radius $R$ one has $$I=\pi ^{3/2} R \sqrt{t} \,\text{erf}(R/\sqrt{t})+\pi t \left(e^{-R^2/t}-1\right)$$ which for small $t$ tends to $\pi R\sqrt{\pi t}+{\cal O}(t)$.
Similarly, for an $L\times L$ square one has $$I=2 L \sqrt{\pi t} \,\text{erf}\left(\tfrac{1}{2}L/\sqrt{t}\right)+4 t \left(e^{-L^2/4 t}-1\right)\rightarrow 2L\sqrt{\pi t}+{\cal O}(t).$$

The estimate you seek is reminiscent of H. Weyl's tube formula. I will give you some pointers referring for more details to section 9.3.5. of these lectures.

Denote by $r$ the distance to $\newcommand{\pa}{\partial}$ $\pa \Omega$ $\newcommand{\bn}{\boldsymbol{n}}$ and by $\bn$ the innner pointing unit normal of $\pa \Omega$. $\newcommand{\bp}{\boldsymbol{p}}$ There exists $r_0>0$ such that the map $\newcommand{\bR}{\mathbb{R}}$ $$ \pa \Omega \times [0,r_0)\ni(\bp,r)\stackrel{\Phi}{\longmapsto} \bp+r\bn\in \bR^n $$ is a diffeomorphism onto the region $\DeclareMathOperator{\dist}{dist}$ $$A_{\rho_0}:=\big\{\; x\in \bar{\Omega};\;\;\dist(x,\pa \Omega)<\rho_0\;\big\}.$$ Then $$\int_{A_{\rho_0}} e^{-r^2/t} dx=\int_{\pa\Omega\times [0,r_0)}\Phi^*\big(e^{-\rho^2/t} dx\Big)=\int_{\pa\Omega\times [0,r_0)}e^{-\rho^2/t}\Phi^*(dx). $$ The pullback of the Euclidean volume form $dx$ via the map $\Phi$ is described explicitly in the above reference. It has the form $\newcommand{\eQ}{\mathscr{Q}}$ $\DeclareMathOperator{\tr}{tr}$ $$ \Phi^*(dx)= \eQ_\bp(r)dV_{\pa \Omega}(\bp) dr, $$ where, for each $\bp\in \pa \Omega$ $(r)$, $\eQ_\bp$ is a polynomial of degree $n-1$ in $r$ $$ \eQ_\bp(r)=\sum_{j=0}^{n-1}c_j(\bp) r^j=\det\big(1-r S_\bp\big), $$ where $S_{\bp}$ is the second fundamental form of the hypersurface $\pa \Omega$ at the point $\bp$ defined in terms of the inner normal $\bn$. More precisely if $(x^i)$ are local coordinates on $\pa \Omega$ near $\bp$, then $$ S_\bp(\pa_{x^i},\pa_{x^j})=\big(\; \bn(\bp),\pa^2_{x^ix^j}\bn(\bp)\;\big), $$ where $(-,-)$ denotes the canonical inner product on $\bR^n$.

Thus $$\int_{A_{\rho_0}} e^{-r^2/t} dx=\int_0^{r_0} e^{-r^2/t}\left(\int_{\pa\Omega} \eQ_{\bp}(r)dV_{\pa \Omega}(\bp)\right) dr. $$ The integral $$ K(r):=\int_{\pa\Omega} \eQ_{\bp}(r)dV_{\pa \Omega}(\bp) $$ appears in Weyl's tube formula and, more precisely $\newcommand{\bom}{\boldsymbol{\omega}}$ (see Eq. (9.3.18) in the above reference)

$$ K(r)=\sum_{k=0}^{n-1} (-1)^{n-1-k}\bom_{n-k}r^{n-k}\mu_k(\Omega), $$ $$ =\bom_0\mu_{n-1}(\Omega)r-\bom_1\mu_{m-2}(\Omega)r^2+\cdots +(-1)^{n-1}\bom_n\mu_0(\Omega)r^n, $$ where $\bom_m$ denotes the volume of the $m$-dimensional Euclidean unit ball and $\mu_k(\Omega)$ is the curvature measure of degree $k$. (You need to be careful about various sign conventions. In the above reference the second fundamental form is defined using the outer normal.)$\DeclareMathOperator{\vol}{vol}$ For example $$ \mu_{n-1}(\Omega)= \frac{1}{2} \vol_{n-1}(\pa\Omega), $$ $$ \mu_{n-2}(\Omega)=\frac{1}{2\pi}\int_{\pa \Omega} \tr S_\bp dV_{\pa\Omega}(\bp), $$ where $\tr S_\bp$ is the mean curvature of $\pa\Omega$ at $\bp$. Also, $\mu_0(\Omega)$ is the Euler characteristic of $\Omega$.

The asymptotics of $$ J(t)=\int_0^{r_0} e^{-r^2/t}K(r) dr, $$ can be determined easily by making the change in variables $s=r^2/t$, $r=\sqrt{st}$ so that $$ J(t)=\frac{\sqrt{t}}{2}\int_0^{r_0^2/t} e^{-s}K(\sqrt{st}) s^{-1/2} ds $$ $$ =\frac{1}{2}\sum_{k=0}^{n-1}(-1)^{n-1-k}t^{\frac{n-k}{2}}\mu_k(\Omega)\int_0^{r_0^2/t} e^{-s} s^{\frac{n-k}{2}-1} ds. $$ Observe that as $t\searrow 0$ $$ \int_0^{r_0^2/t} e^{-s} s^{\frac{n-k}{2}-1} ds=\;\underbrace{\int_0^{\infty} e^{-s} s^{\frac{n-k}{2}-1} ds}_{\Gamma\big( \frac{n-k}{2}\big)}\;+ O\big(t^{N}\big),\;\;\forall N>0. $$ Hence $$ J(t)=\frac{1}{2}\sum_{k=0}^{n-1}(-1)^{n-1-k}t^{\frac{n-k}{2}}\Gamma\Big(\;\frac{n-k}{2}\;\Big)\mu_k(\Omega) + O\big(t^{N}\big),\;\;\forall N>0. $$ Finally $$ |I(t)-J(t)|=\int_{\dist(x,\pa \Omega)>r_0} e^{-\dist(x,\pa\Omega)^2/t} dx $$ $$ \leq e^{-r_0^2/t}\vol(\Omega)=O(t^{N}),\;\;\forall N>0. $$ Hence

$$ I(t)=\frac{1}{2}\sum_{k=0}^{n-1}(-1)^{n-1-k}t^{\frac{n-k}{2}}\Gamma\Big(\;\frac{n-k}{2}\;\Big)\mu_k(\Omega) + O\big(t^{N}\big),\;\;\forall N>0.$$

The leading term of this yields the estimate indicated by Carlo Beenakker.

About the tube formula The tube formula for smooth domains or convex domains or more generally sets with positive reach are both special cases of the very general kinematic formulas. The Brunn-Minkowski formula is also a special case of the kinematic formula.

\begin{aligned} & \int_{0}^{+\infty} e^{-\frac{s^{2}}{t}} \mathcal{L}_{n-1}(\{x \in \Omega \mid d(x, \partial \Omega)=s\}) d s \\ =& \int_{0}^{+\infty} e^{-\frac{s^{2}}{t}} \mathcal{L}_{n-1}(\partial\Omega) d s+\int_{0}^{+\infty} e^{-\frac{s^{2}}{t}}\left(\mathcal{L}_{n-1}\left(\left\{x\in \Omega \mid d(x, \partial \Omega)=s\right)-\mathcal{L}_{n-1}(\partial \Omega)\right) d s\right. \end{aligned}

Let us estimate the error term. It split into the estimate far from the boundary and the distortion estimate of the level set near the boundary.

The estimate far from the boundary:
$$ \int_{0}^{+\infty} e^{-\frac{s^{2}}{t}} \mathcal{L}_{n-1}(\partial \Omega) d \delta=\frac{1}{2} \sqrt{\pi t} \mathcal{L}_{n-1}(\partial \Omega) $$

$$ \int_{0}^{A} e^{-\frac{s^{2}}{t}} \mathcal{L}_{n-1}(\partial \Omega) d s=\frac{1}{2} \sqrt{\pi t} \mathcal{L}_{n-1}(\partial \Omega)+O\left(\sqrt{t} \cdot e^{-\frac{A^{2}}{t}} \mathcal{L}_{N-1}{(\partial \Omega)}\right) $$ The distortion estimate of the level set near the boundary:
Then estimate the distortion of $\mathcal{L}_{n-1}\left(\{x \in \Omega \mid d(x, \partial \Omega)=s\})-\mathcal{L}_{n-1}(\partial \Omega)\right.$ $$ \left.\mathcal{L}_{n-1}(\{x \in \Omega \mid d(x, \partial \Omega)=s\}\right)-\mathcal{L}_{n-1}(\partial \Omega)=O\left(s^{n-1}\right)=\frac{2 \pi^{\frac{n}{2}}}{\Gamma\left(\frac{n}{2}\right)} s^{n-1}+o\left(s^{n-1}\right) $$ This estimate come from a similar way in proving Brunn-Minkowski expansion, which give the $\mathcal{L}_{n}(A+B_r)\sim \sum_{i=0}^n \lambda_i\mathcal{L}_{n-i}(A)\mathcal{L}_{i}(B_r)$ at least for $A$ is a convex domain. and drop the high order term $\mathcal{L}_{n-i}(A)\mathcal{L}_{i}(B_r), i\geq2$, which is minor when $t\to 0$.

Combine together:
Now a suitable choice of $s$ relies on $t$ will given a reasonable estimate for the error term.,

$$E(t)=\int_{0}^{+\infty} e^{-\frac{s^{2}}{t}} \mathcal{L}_{n-1}(\{x \in \Omega \mid d(x, \partial \Omega)=s\}) d s-\frac{1}{2}\sqrt{\pi t}\mathcal{L}_{n-1}(\partial \Omega)$$

In fact, a suitable $s_0$ is near the maximum of the minimum. $$ s_0=\max _{s} \min\left\{\frac{2 \pi^{\frac{1}{2}}}{\Gamma\left(\frac{n}{2}\right)} s^{n-1}+o\left(s^{n-1}\right), \quad \sqrt{t} \cdot e^{-\frac{s^{2}}{t}} \cdot \mathcal{L}_{n-1}(\partial \Omega)\right\} $$

I have not finished the calculation of the explicit order(the best choice of $s_0$), but I think this can give some good enough order of error term in the general case.

Then at least for $n\geq 3$, $\Omega\subset \mathbb{R}^n$, take $s_0\sim t^{\frac{1}{2}+\delta}, \delta>0$, we know $E(t)=o(t^{\frac{1}{2}+\delta^{'}}), \delta^{'}>0$.

In dimension 2 case it seems we even do not have $I \rightarrow \frac{1}{2} S \sqrt{\pi t}$ because the error term is out of control.