Let $c \in {\mathbb R}$ and $d \in {\mathbb R}$. Without loss of generality we can consider a slightly different function:
\begin{equation}
f(c,d):= \int\limits_{-\infty}^\infty \exp(-t^2) \mbox{erfc}(t-c) \mbox{erfc}(t-d) dt
\end{equation}
Now we differentiate with respect to $c$. We have:
\begin{eqnarray}
\partial_c f(c,d) &=& \lim\limits_{\epsilon \rightarrow 0} \frac{\sqrt{2}}{\sqrt{\pi}} \exp(-\frac{c^2}{2})
\left( \sqrt{\pi} + 2 \sqrt{\pi} T(\epsilon,\frac{1}{\sqrt{2}},\frac{2d-c}{\sqrt{2}}) - 2 \sqrt{\pi} T(\epsilon,\frac{1}{\sqrt{2}},\frac{-2d+c}{\sqrt{2}})\right)\\
&=&\sqrt{2} \exp(-\frac{c^2}{2}) \left( 1-erf(\frac{c-2 d}{\sqrt{6}})\right)
\end{eqnarray}
where $T(h,a,b)$ is the generalized Owen's T function Generalized Owen's T function .
Now, all we need to do is to integrate. Since $f(-\infty,d)=0$ we integrate over $c$ from minus infinity to $c$. We have:
\begin{eqnarray}
f(c,d) &=& \sqrt{2} \left( \sqrt{\frac{\pi}{2}}(1+erf(\frac{c}{\sqrt{2}}) + 2erf(\frac{d}{\sqrt{2}})) + \sqrt{2 \pi} 2 T(c,\frac{1}{\sqrt{3}},-\frac{2 d}{\sqrt{3}})\right) \\
&=& \frac{1}{3 \sqrt{\pi}} \left( 12 \pi T\left(c,\frac{c-2 d}{\sqrt{3} c}\right)+12 \pi T\left(d,\frac{d-2 c}{\sqrt{3} d}\right)-6 \arctan\left(\frac{c-2 d}{\sqrt{3} c}\right)-6 \arctan\left(\frac{d-2 c}{\sqrt{3} d}\right)+3
\pi \text{erf}\left(\frac{c}{\sqrt{2}}\right)+3 \pi \text{erf}\left(\frac{d}{\sqrt{2}}\right)+4 \pi\right)
\end{eqnarray}
where $T(h,a)$ is Owen's T function https://en.wikipedia.org/wiki/Owen%27s_T_function .
In[1060]:= {c, d} = RandomReal[{-2, 2}, 2, WorkingPrecision -> 50];
NIntegrate[Exp[-t^2] Erfc[t - c] Erfc[t - d], {t, -Infinity, Infinity}]
(4 \[Pi] - 6 ArcTan[(c - 2 d)/(Sqrt[3] c)] -
6 ArcTan[(-2 c + d)/(Sqrt[3] d)] + 3 \[Pi] Erf[c/Sqrt[2]] +
3 \[Pi] Erf[d/Sqrt[2]] + 12 \[Pi] OwenT[c, (c - 2 d)/(Sqrt[3] c)] +
12 \[Pi] OwenT[d, (-2 c + d)/(Sqrt[3] d)])/(3 Sqrt[\[Pi]])
Out[1061]= 0.483318
Out[1062]= 0.48331807775609703646923225386370751256977344715
Update: Now let us take four real numbers $a_1 \in {\mathbb R}$, $a_2 \in {\mathbb R}$, $c\in {\mathbb R}$ and $d \in {\mathbb R}$ and consider a more general integral:
\begin{equation}
f^{(a_1,a_2)}(c,d):= \int\limits_{-\infty}^\infty \exp(-t^2) \mbox{erfc}(a_1 t-c) \mbox{erfc}(a_2 t-d) dt
\end{equation}
Then by doing the same calculations as above we easily arrive at the following formula:
\begin{eqnarray}
&&f^{(a_1,a_2)}(c,d)= \frac{1}{\sqrt{\pi}} \left(\right.\\
&&
4 \pi T\left(\frac{\sqrt{2} c}{\sqrt{a_1^2+1}},\frac{a_1 a_2 c-\left(a_1^2+1\right) d}{c \sqrt{a_1^2+a_2^2+1}}\right)+4 \pi T\left(\frac{\sqrt{2}
\sqrt{a_1^2+1} d}{\sqrt{a_1^2 a_2^2+a_1^2+a_2^2+1}},\frac{a_1 a_2 d-\left(a_2^2+1\right) c}{d \sqrt{a_1^2+a_2^2+1}}\right) +\\
&&-2 \arctan\left(\frac{a_1
a_2 c-\left(a_1^2+1\right) d}{c \sqrt{a_1^2+a_2^2+1}}\right)-2 \arctan\left(\frac{a_1 a_2 d-\left(a_2^2+1\right) c}{d \sqrt{a_1^2+a_2^2+1}}\right)+\\
&&\pi
\text{erf}\left(\frac{\sqrt{a_1^2+1} d}{\sqrt{a_1^2 a_2^2+a_1^2+a_2^2+1}}\right)+\pi
\text{erf}\left(\frac{c}{\sqrt{a_1^2+1}}\right)+\\
&&\pi +2 \arctan\left(\frac{a_1 a_2}{\sqrt{a_1^2+a_2^2+1}}\right)\\
&&\left.\right)
\end{eqnarray}
For[count = 1, count <= 200, count++,
{a1, a2, c, d} = RandomReal[{-5, 5}, 4, WorkingPrecision -> 50];
I1 = NIntegrate[
Exp[-t^2] Erfc[a1 t - c] Erfc[a2 t - d], {t, -Infinity, Infinity}];
I2 = 1/Sqrt[\[Pi]] (\[Pi] +
2 ArcTan[(a1 a2)/Sqrt[(1 + a1^2 + a2^2)]] -
2 ArcTan[ (a1 a2 c - (1 + a1^2) d)/(
Sqrt[(1 + a1^2 + a2^2)] c)] -
2 ArcTan[ (a1 a2 d - (1 + a2^2) c)/(
Sqrt[(1 + a1^2 + a2^2)] d)] + \[Pi] Erf[c/Sqrt[
1 + a1^2]] + \[Pi] Erf[(Sqrt[(1 + a1^2)] d)/ Sqrt[
1 + a1^2 + a2^2 + a1^2 a2^2]] +
4 \[Pi] OwenT[(Sqrt[2] c)/Sqrt[
1 + a1^2], (a1 a2 c - (1 + a1^2) d)/(
Sqrt[(1 + a1^2 + a2^2)] c)] +
4 \[Pi] OwenT[(Sqrt[2] Sqrt[(1 + a1^2)] d)/ Sqrt[
1 + a1^2 + a2^2 + a1^2 a2^2], (a1 a2 d - (1 + a2^2) c)/(
Sqrt[(1 + a1^2 + a2^2)] d)]);
If[Abs[I2/I1 - 1] > 10^(-3),
Print["results do not match", {a1, a2, c, d, {I1, I2}}]; Break[]];
If[Mod[count, 10] == 0, PrintTemporary[count]];
];
Here is my trial, which is partially successful but still not fully answering to your question.
Using some coordinate change, I derived that
$$ I_{n} := \int_{0}^{\infty} \mathrm{erfc}^{n}(x) \, dx = \frac{1}{\sqrt{n}} \left( \frac{2}{\sqrt{\pi}} \right)^{n-1} \int_{T^{n-1}} \int_{0}^{\infty} s^{n-1}e^{-(|x|^{2}-1)s^{2}} \mathrm{erfc}(s) \, ds d\sigma_{x}, \tag{*} $$
where $T^{n-1} = \{ x \in [0, \infty)^{n} : x_{1} + \cdots + x_{n} = \sqrt{n} \}$ is an $(n-1)$-simplex and $d\sigma_{x}$ is the surface measure on $T^{n-1}$.
Case $n = 1$. Since $T^{0} = \{1\}$ and $d\sigma_{x} = \delta(x-1)$ is a unit mass located at $x = 1$, the equation $\text{(*)}$ gives no new information.
Case $n = 2.$ It is easy to find that
$$ \int_{0}^{\infty} s e^{-(|x|^{2}-1)s^{2}} \mathrm{erfc}(s) \, ds = \frac{1}{2|x|(|x|+1)}$$
and by suitable parametrization of the line segment $T^{1}$, we get
\begin{align*}
I_{2}
&= \sqrt{\frac{2}{\pi}} \int_{T^{1}} \frac{1}{2|x|(|x|+1)} \, d\sigma_{x} \\
&= \sqrt{\frac{2}{\pi}} \int_{0}^{\sqrt{2}} \frac{\sqrt{2} \, dt}{2\sqrt{t^{2} + (\sqrt{2}-t)^{2}} ( \sqrt{t^{2} + (\sqrt{2}-t)^{2}} + 1)} \\
&= \sqrt{\frac{2}{\pi}} \int_{0}^{\frac{\pi}{4}} \frac{d\theta}{1+\cos\theta}
= \sqrt{\frac{2}{\pi}} \tan\left(\frac{\pi}{8}\right)
= \frac{2 - \sqrt{2}}{\sqrt{\pi}}.
\end{align*}
We can also write
$$ I_{2} = \boxed{ \displaystyle \frac{2}{\sqrt{\pi}} - \frac{2\sqrt{2}}{\pi^{3/2}} \arctan(\infty) }. $$
Case $n = 3$. We have
$$ \int_{0}^{\infty} s^{2}e^{-(|x|^{2}-1)s^{2}} \mathrm{erfc}(s) \, ds = \frac{1}{2\sqrt{\pi}(|x|^{2}-1)} \left( \frac{\arctan\sqrt{|x|^{2}-1}}{\sqrt{|x|^{2}-1}} - \frac{1}{|x|^{2}} \right). $$
Now we know that $T^{2}$ is a regular triangle with side length $\sqrt{6}$. Thus by introducing polar coordinate change, we get
\begin{align*}
I_{3}
&= \frac{2}{\pi^{3/2}\sqrt{3}} \int_{T^{2}} \frac{1}{|x|^{2}-1} \left( \frac{\arctan\sqrt{|x|^{2}-1}}{\sqrt{|x|^{2}-1}} - \frac{1}{|x|^{2}} \right) \, d\sigma_{x} \\
&= \frac{4\sqrt{3}}{\pi^{3/2}} \int_{0}^{\frac{\pi}{3}} \int_{0}^{\frac{\sec\theta}{\sqrt{2}}} \frac{1}{r} \left( \frac{\arctan r}{r} - \frac{1}{r^{2} + 1} \right) \, drd\theta \\
&= \frac{4\sqrt{3}}{\pi^{3/2}} \int_{0}^{\frac{\pi}{3}} \left( 1- \frac{\arctan\left(\sec\theta / \sqrt{2}\right)}{\sec\theta / \sqrt{2}} \right) \, d\theta \\
&= \frac{4\sqrt{3}}{\pi^{3/2}} \int_{\frac{1}{\sqrt{2}}}^{\sqrt{2}} \left( 1 - \frac{\arctan u}{u} \right) \frac{du}{u\sqrt{2u^{2} - 1}}.
\end{align*}
I haven't tried the last integral, but Mathematica suggests that
$$ \int_{\frac{1}{\sqrt{2}}}^{\sqrt{2}} \left( 1 - \frac{\arctan u}{u} \right) \frac{du}{u\sqrt{2u^{2} - 1}} = \frac{\sqrt{3}\pi}{4} -\sqrt{\frac{3}{2}} \arctan\left( \sqrt{2} \right).$$
This gives
$$ I_{3} = \boxed{ \displaystyle \frac{3}{\sqrt{\pi}} - \frac{6\sqrt{2}}{\pi^{3/2}} \arctan\left( \sqrt{2} \right) }, $$
which can be numerically checked.
I am trying to generalize this calculation for higher dimensions, but it seems somewhat daunting. For example, if $n = 4$ we have to evaluate
$$ \int_{0}^{\infty} \mathrm{erfc}^{4}(x) \, dx = \frac{1}{\pi^{3/2}} \int_{T^{3}} \frac{2|x|+1}{(|x|+1)^{2}|x|^{3}} \, d\sigma_{x}. $$
Nevertheless, using the formula
$$ \int_{T^{n-1}} f(|x|) \, d\sigma_{x} = n! \int_{0}^{\sqrt{\frac{1}{n-1}}} \int_{0}^{\sqrt{\frac{n}{n-2}}s_{1}} \cdots \int_{0}^{\sqrt{\frac{3}{1}}s_{n-2}} f\left( \sqrt{1+s_{1}^{2} + \cdots + s_{n-1}^{2}} \right) \, ds_{n-1} \cdots ds_{1} $$
together with aid of Mathematica, I was able to evaluate $I_{4}$ and it was
$$ I_{4} = \boxed{ \displaystyle \frac{4}{\sqrt{\pi}} - \frac{24\sqrt{2}}{\pi^{3/2}} \arctan \left( \frac{1}{\sqrt{8}} \right) }. $$
These series of observations lead us to believe that $I_{n}$ is of the form
$$ I_{n} = \frac{n}{\sqrt{\pi}} - \frac{n!\sqrt{2}}{\pi^{3/2}} \arctan \alpha_{n} $$
for some reasonably nice $\alpha_{n}$ (with $\alpha_{1} = 0$, $\alpha_{2} = \infty$, $\alpha_{3} = \sqrt{2}$ and $\alpha_{4} = \frac{1}{\sqrt{8}}$), but inverse symbolic calculations show that this seems no longer the case for $n \geq 5$.
Indeed, for $n = 5$ we have
$$ I_{5} = \frac{5}{\sqrt{\pi}} - \frac{80\sqrt{2}}{\pi^{3/2}} \arctan \sqrt{\frac{2}{3}} + \frac{240\sqrt{2}}{\pi^{5/2}} A\left( \sqrt{\frac{5}{3}}, \sqrt{\frac{3}{2}}, \sqrt{\frac{4}{5}} \right), $$
where $A(p, q, r)$ is the generalized Ahmed's integral. I have no idea whether it will reduce to a familiar closed form or not.
Best Answer
Since $\int \exp(b x^2+c x) dx= \frac{\sqrt{\pi } e^{-\frac{c^2}{4 b}} \text{erfi}\left(\frac{2 b x+c}{2 \sqrt{b}}\right)}{2 \sqrt{b}}$ we integrate by parts and we have: \begin{eqnarray} I&=& -\frac{\sqrt{\pi } e^{-\frac{c^2}{4 b}} \text{erfi}\left(\frac{c}{2 \sqrt{b}}\right)}{2 \sqrt{b}}- \frac{2 a}{\sqrt{\pi}} \int\limits_0^\infty \frac{\sqrt{\pi } e^{-\frac{c^2}{4 b}} \text{erfi}\left(\frac{2 b x+c}{2 \sqrt{b}}\right)}{2 \sqrt{b}} \exp(-a^2 x^2) dx\\ &=& -\frac{\sqrt{\pi } e^{-\frac{c^2}{4 b}} \text{erfi}\left(\frac{c}{2 \sqrt{b}}\right)}{2 \sqrt{b}}+\frac{e^{-\frac{c^2}{4 b}}\sqrt{\pi}}{\imath \sqrt{b}}2 T(\epsilon, \imath \frac{\sqrt{b}}{a}, \imath \frac{c \sqrt{2}}{2 \sqrt{b}} ) \\ &=& -\frac{\sqrt{\pi } e^{-\frac{c^2}{4 b}} \left(4 i T\left(\frac{i c}{\sqrt{2} \sqrt{b} \sqrt{1-\frac{b}{a^2}}},\frac{i \sqrt{b}}{a}\right)-\text{erfi}\left(\frac{c}{2 \sqrt{b} \sqrt{1-\frac{b}{a^2}}}\right)+\text{erfi}\left(\frac{c}{2 \sqrt{b}}\right)\right)}{2 \sqrt{b}} \end{eqnarray} where in the second line we took a small number $0 < \epsilon << 1$ and we used the definition of the generalized Owen's T function Generalized Owen's T function and in the last line we simplified the result. Here $T(\cdot,\cdot)$ is the Owen's T function. The result is valid for $0 < b < a^2$.