$\operatorname{erf}(x)$ is an odd function, therefore,
$$
\begin{align}
\int_{-\infty}^\infty(\operatorname{erf}(a+x)+\operatorname{erf}(a-x))\;\mathrm{d}x
&=\lim_{L\to\infty}\;\int_{-L}^L(\operatorname{erf}(x+a)-\operatorname{erf}(x-a))\;\mathrm{d}x\\
&=\lim_{L\to\infty}\;\int_{-L+a}^{L+a}\operatorname{erf}(x)\;\mathrm{d}x-\lim_{L\to\infty}\;\int_{-L-a}^{L-a}\operatorname{erf}(x)\;\mathrm{d}x\\
&=\lim_{L\to\infty}\;\int_{L-a}^{L+a}\operatorname{erf}(x)\;\mathrm{d}x-\lim_{L\to\infty}\;\int_{-L-a}^{-L+a}\operatorname{erf}(x)\;\mathrm{d}x\\
&=4a\tag{1}
\end{align}
$$
since $\lim\limits_{x\to\infty}\operatorname{erf}(x)=1$ and $\lim\limits_{x\to-\infty}\operatorname{erf}(x)=-1$.
Furthermore,
$$
\begin{align}
\int_{-\infty}^\infty\left(\operatorname{erf}(a+x)\operatorname{erf}(a-x)+1 \right)\;\mathrm{d}x
&=\int_{-\infty}^\infty(\operatorname{erf}(a+x)+1)(\operatorname{erf}(a-x)+1)\;\mathrm{d}x\\
&-\int_{-\infty}^\infty(\operatorname{erf}(a+x)+\operatorname{erf}(a-x))\;\mathrm{d}x\tag{2}
\end{align}
$$
To evaluate
$$
\begin{align}
\int_{-\infty}^\infty(\operatorname{erf}(a+x)+1)(\operatorname{erf}(a-x)+1)\;\mathrm{d}x
&=\frac{4}{\pi}\int_{-\infty}^\infty\int_{-\infty}^{a+x}\int_{-\infty}^{a-x}e^{-s^2-t^2}\;\mathrm{d}s\;\mathrm{d}t\;\mathrm{d}x
\end{align}
$$
note that $s\le a+x$ and $t\le a-x$; i.e. $s-a\le x\le a-t$ and $s+t\le2a$. Thus,
$$
\begin{align}
\frac{4}{\pi}\int_{-\infty}^\infty\int_{-\infty}^{a+x}\int_{-\infty}^{a-x}e^{-s^2-t^2}\;\mathrm{d}s\;\mathrm{d}t\;\mathrm{d}x
&=\frac{4}{\pi}\int\int_{s+t\le2a}\int_{s-a}^{a-t}e^{-s^2-t^2}\;\mathrm{d}x\;\mathrm{d}s\;\mathrm{d}t\\
&=\frac{4}{\pi}\int_{-\infty}^\infty\int_{-\infty}^\infty(2a-s-t)_+\;e^{-s^2-t^2}\;\mathrm{d}s\;\mathrm{d}t
\end{align}
$$
Change variables: $u=(s+t)/\sqrt{2}$ and $v=(s-t)/\sqrt{2}$ so that $s=(u+v)/\sqrt{2}$ and $t=(u-v)/\sqrt{2}$:
$$
\begin{align}
\frac{4}{\pi}\int_{-\infty}^\infty\int_{-\infty}^\infty(2a-s-t)_+\;e^{-s^2-t^2}\;\mathrm{d}s\;\mathrm{d}t
&=\frac{4}{\pi}\int_{-\infty}^\infty\int_{-\infty}^\infty(2a-\sqrt{2}u)_+\;e^{-u^2-v^2}\;\mathrm{d}u\;\mathrm{d}v\\
&=\frac{4}{\sqrt{\pi}}\int_{-\infty}^{\sqrt{2}a}(2a-\sqrt{2}u)\;e^{-u^2}\;\mathrm{d}u\\
&=4a(\operatorname{erf}(\sqrt{2}a)+1)-\frac{4}{\sqrt{\pi}}\int_{-\infty}^{\sqrt{2}a}\sqrt{2}u\;e^{-u^2}\;\mathrm{d}u\\
&=4a(\operatorname{erf}(\sqrt{2}a)+1)-\frac{2\sqrt{2}}{\sqrt{\pi}}\int_{-\infty}^{\sqrt{2}a}\;e^{-u^2}\;\mathrm{d}u^2\\
&=4a(\operatorname{erf}(\sqrt{2}a)+1)+\frac{2\sqrt{2}}{\sqrt{\pi}}\;e^{-2a^2}
\end{align}
$$
Therefore,
$$
\int_{-\infty}^\infty\left(\operatorname{erf}(a+x)+1\right)\left(\operatorname{erf}(a-x)+1\right)\;\mathrm{d}x
=4a\left(\operatorname{erf}(\sqrt{2}a)+1\right)+\frac{2\sqrt{2}}{\sqrt{\pi}}\;e^{-2a^2}\tag{3}
$$
Thus, the convolution of $\operatorname{erf}(x)+1$ with itself is $2x(\operatorname{erf}(x/\sqrt{2})+1)+\frac{2\sqrt{2}}{\sqrt{\pi}}e^{-x^2/2}$.
Subtract $4a$ from $(3)$ using $(1)$ and $(2)$ to get
$$
\int_{-\infty}^\infty\left(\operatorname{erf}(a+x)\operatorname{erf}(a-x)+1 \right)\;\mathrm{d}x
=4a\operatorname{erf}(\sqrt{2}a)+\frac{2\sqrt{2}}{\sqrt{\pi}}\;e^{-2a^2}\tag{4}
$$
My guess is you want either $(3)$ or $(4)$.
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]];
];
Best Answer
Differentiating under the integral sign using Mathematica 11.3 gives
$$\frac{\partial}{\partial a}\frac{\partial}{\partial b}\left( \frac{e^{-\frac{x^2}{2}} \text{erf}(a x) \text{erf}(b x)}{\sqrt{2 \pi }}\right)=\frac{2 \sqrt{2} x^2 e^{-a^2 x^2-b^2 x^2-\frac{x^2}{2}}}{\pi ^{3/2}}$$
Integrating the new function is then possible in Mathematica
$$\text{Integrate}\left[\frac{2 \sqrt{2} x^2 e^{-a^2 x^2-b^2 x^2-\frac{x^2}{2}}}{\pi ^{3/2}},\{x,-\infty ,\infty \},\text{Assumptions}\to a>0\land b>0\right]=\frac{4}{\pi \left(2 a^2+2 b^2+1\right)^{3/2}}$$
Finally reverse the earlier differentiations with respect to $a$ and $b$ to get
$$\int \int \frac{4}{\pi \left(2 a^2+2 b^2+1\right)^{3/2}}\,da\,db=\frac{2 \tan ^{-1}\left(\frac{2 a b}{\sqrt{2 a^2+2 b^2+1}}\right)}{\pi }$$
In Mathematica code this is: