I somewhat doubt that a closed form exists for general $n$ even admitting familiar special functions. Even so, there are some observations that might be useful.
First, I'll generalize your notation a bit. It will be useful to look at "standard" simplices of varying sizes, which I'll refer to as follows:
$$
\Delta_n(r)=\left\{\vec x\in\mathbb R^n:\ x_i\ge 0,\ \sum_i x_i\le r\right\}
$$
Also, to keep prefactors to a minimum, I'll use a rescaling of the error function that I'll call $E$:
$$
E(x)=\int_0^xe^{-y^2/2}dy=\sqrt{\frac{\pi}{2}}\text{erf}(x/\sqrt{2})
$$
We can write the family of integrals you're interested in as:
$$
I_n(r)=\int_{\Delta_n(r)}e^{-(x_1^2+...+x_n^2)/2}dx_1...dx_n
$$
This integral can be factorized somewhat when written in the typical simplical format.
$$
I_n(r)=\int_0^re^{-x_1^2/2}dx_1\int_0^{r-x_1}e^{-x_2^2/2}dx_2\int_0^{r-x_1-x_2}e^{-x_3^2/2}dx_3...\int_0^{r-x_1-...-x_{n-1}}e^{-x_n^2/2}dx_n
$$
Here we can observe a recurrence relation, Notice that if we remove the leftmost integral, we get the same expression with one less integral and $r-x_1$ in place of $r$.
$$
I_{n+1}(r)=\int_0^re^{-x^2/2}I_n(r-x)dx=\int_0^re^{-(r-x)^2/2}I_n(x)dx
$$
Now we can enumerate some cases. $n=0$ isn't necessarily meaningful, but this choice satisfies the recurrence relation.
$$
I_0(r)=1
$$
$n=0$ is a of course just a Gaussian integral.
$$
I_1(r)=E(r)
$$
$n=2$ can be solved using Mike Earnest's insightful observation that the 4 such simplical regions can be rotated into a square, yielding a separable integral. This can also be shown (tediously) by differentiating with respect or $r$ under the integral sign in the recurrence relation for $n=1$.
$$
I_2(r)=\frac 14\int_{-r/\sqrt{2}}^{r/\sqrt{2}}\int_{-r/\sqrt{2}}^{r/\sqrt{2}}e^{-(x^2+y^2)/2}dy\ dx=E^2(r/\sqrt{2})
$$
We can use the same trick to derive a 2nd order recurrence relation with a single integral, but It's not clearly more useful than the first.
$$
I_{n+2}(r)=\sqrt{2}\int_0^r e^{-(r-x)^2/4}E((r-x)/2)I_n(x)dx
$$
If you want to compute these values numerically, it is best to transform the functional recurrence relation into a discrete one. We can do that by assuming that $I_n(r)$ is an analytic function of $r$, and writing it in terms of its Maclaurian expansion $\sum_p a_p^{(n)}r^p$. As a shorthand, let $g_k$ be the Maclaurian coefficients of the Gaussian.
$$
\sum_pa^{(n+1)}_pr^p=\int_0^r\sum_{k=0}^\infty g_k(r-x)^k\sum_{l=0}^\infty a^{(n)}_lx^ldr,\ \ \ g_k=\begin{cases}\frac{(-1)^{k/2}}{2^{k/2}(k/2)!} & k\ \text{even} \\ 0 & k\ \text{odd} \end{cases}
$$
$$
\sum_pa^{(n+1)}_pr^p=\int_0^r\sum_{k=0}^\infty \sum_{l=0}^\infty\sum_{m=0}^kg_ka^{(n)}_l(-1)^m\binom{k}{m}r^mx^{l+k-m}dr
$$
$$
\sum_pa^{(n+1)}_pr^p=\sum_{k=0}^\infty \sum_{l=0}^\infty\sum_{m=0}^kg_ka^{(n)}_l\frac{(-1)^m}{k+l-m+1}\binom{k}{m}r^{k+l+1}
$$
Now, we let $c(p,l)=\sum_{m=0}^{p-l-1}g_{p-l-1}\frac{(-1)^m}{p-m}\binom{p-l-1}{m}$
$$
\sum_pa^{(n+1)}_pr^p=\sum_{k=0}^\infty \sum_{l=0}^\infty c(k+l+1,l)a^{(n)}_lr^{k+l+1}
$$
Reindexing to collect like powers:
$$
\sum_pa^{(n+1)}_pr^p=\sum_{p=1}^\infty\sum_{l=0}^{p-1}c(p, l)a^{(n)}_lr^p
$$
We obtain a somewhat serviceable recurrence relation. We can start with the initial condition $a^{(0)}_p=\delta_{p0}$.
$$
a^{(n+1)}_p=\sum_{l=0}^{p-1}c(p, l)a^{(n)}_l
$$
To evaluate $I_n(1)$, these coefficients can be generated reasonably quickly, and the resulting series converges quite quickly. Here are the results of a simple python script, which agrees with the exact results we already have:
$$
\begin{array}{|c|c|} \hline
n & I_n(1) \\ \hline
0 & 1 \\ \hline
1 & 0.855624392 \\ \hline
2 & 0.425560334 \\ \hline
3 & 0.1438817189 \\ \hline
4 & 0.0365313603 \\ \hline
5 & 7.40668651\cdot 10^{-3} \\ \hline
6 & 1.24877137\cdot 10^{-3} \\ \hline
7 & 1.80133459\cdot 10^{-4} \\ \hline
8 & 2.27017316\cdot 10^{-5} \\ \hline
9 & 2.54005636e\cdot 10^{-6} \\ \hline
10 & 2.55531558\cdot 10^{-7} \\ \hline
\end{array}
$$
If $n$ is very large, there is also a useful approximation. The integral $\int_{\Delta_n(1)}x_1^2$ is perfectly doable, which allows us to write the mean value of the radius squared $\rho^2=x_1^2+...+x_n^2$ on the standard simplex.
$$
\langle\rho^2\rangle=\frac{2n}{(n+1)(n+2)}
$$
Using this, along with the fact that the distribution of $\rho^2$ becomes increasingly sharp for large $n$, allows us to approximate the Gaussian as a constant.
$$
I_n(1)\approx\frac{1}{n!}e^{\frac{-n}{(n+1)(n+2)}},\ \ \ n>>1
$$
Numerically, it appears the relative error of this approximation tends toward zero as $n$ increases.
EDIT
The natural first choice when confronted with an integral equation is to use an integral transform, especially since the relation is essentially an n-fold "convolution". Fourier Transforms are not so cooperative, since they have the "nice" convolution identity
$$
\mathcal F\left[\int_{-\infty}^{\infty}f(x-y)g(y)dy\right](k)=\mathcal F[f(x)](k)\mathcal F[g(x)](k)
$$
Which doesn't quite match the integral we have. Things can be forced by introducing some Heaviside functions, but the result is a mess.
Laplace transforms, however, have exactly the convolution identity we need:
$$
\mathcal L\left[\int_0^x f(x-y)g(y)dy\right](s)=\mathcal L[f(x)](s)\mathcal L[g(x)](s)
$$
The functions involved transform in a not to terrible way:
$$
\mathcal L[1](s)=\frac 1s,\ \ \ \mathcal L[e^{-x^2/2}](s)=e^{s^2/2}\text{erfc}(s/\sqrt{2})
$$
Putting everything together:
$$
\mathcal L[I_n(r)](s)=\frac 1s\left(e^{s^2/2}\text{erfc}(s/\sqrt{2})\right)^n
$$
Transforming back to real space is a single contour integral, which could be useful numerically. However, though the closed form remains elusive.
