I think this is a possible solution.
For the case where we want to integrate over all of $\mathbb{R}$ (and probably for other cases, but this is simplest), it looks like we can solve the integral by going back to the initial representation of the convolution of a uniform and a Gaussian, taking the product there, and exchanging orders of integration.
That is, solve
$$
\frac{1}{(b-a)(d-c)2\pi\sigma_1\sigma_2} \int_{c}^{d}\int_{a}^{b}\int_{-\infty}^{\infty} \exp \left\{-\frac{(\tau -t)^2}{2\sigma_1^2} \right\}\exp \left\{-\frac{(\rho -t)^2}{2\sigma_2^2} \right\} dt d\tau d\rho
$$
The innermost integral is just a cross-correlation (or convolution with one of the means negated) of two Gaussians, with solution
$$
\begin{align*}
&\frac{1}{2\pi \sigma_1\sigma_2}\int_{-\infty}^{\infty} \exp \left\{-\frac{(\tau -t)^2}{2\sigma_1^2} \right\}\exp \left\{-\frac{(\rho -t)^2}{2\sigma_2^2} \right\} dt\\
&=\frac{1}{\sqrt{2\pi(\sigma_1^2+\sigma_2^2)}}\exp\left(-\frac{(\tau-\rho)^2}{2(\sigma_1^2+\sigma_2^2)} \right)
\end{align*}
$$
So now we should be able to integrate this Gaussian function twice to get the desired quantity. The antiderivative of $\text{erf}$ is $z\text{erf}(z)+\exp(-z^2)/\sqrt{\pi} + C$, which will be needed on the final integration.
However, this looks like a pain, but Mathematica gives the solution to
$$
\int _a^b\int _c^d\exp \left(-\frac{(\tau-\rho)^2}{2 \left(\sigma_1^2+\sigma_2^2\right)}\right) d\rho d\tau
$$
as
$$
\sqrt{\frac{\pi }{2}} \sqrt{\sigma_1^2+\sigma_2^2} \left(\sqrt{\frac{2}{\pi }} \sqrt{\sigma_1^2+\sigma_2^2} \left(e^{-\frac{(b-c)^2}{2 \left(\sigma_1^2+\sigma_2^2\right)}}-e^{-\frac{(a-c)^2}{2 \left(\sigma_1^2+\sigma_2^2\right)}}\right)+\sqrt{\frac{2}{\pi }} \sqrt{\sigma_1^2+\sigma_2^2} \left(e^{-\frac{(a-d)^2}{2 \left(\sigma_1^2+\sigma_2^2\right)}}-e^{-\frac{(b-d)^2}{2 \left(\sigma_1^2+\sigma_2^2\right)}}\right)+(a-c) \text{erf}\left(\frac{c-a}{\sqrt{2} \sqrt{\sigma_1^2+\sigma_2^2}}\right)+(d-a) \text{erf}\left(\frac{d-a}{\sqrt{2} \sqrt{\sigma_1^2+\sigma_2^2}}\right)-(b-c) \text{erf}\left(\frac{c-b}{\sqrt{2} \sqrt{\sigma_1^2+\sigma_2^2}}\right)+(b-d) \text{erf}\left(\frac{d-b}{\sqrt{2} \sqrt{\sigma_1^2+\sigma_2^2}}\right)\right)
$$
so I believe the final answer should be:
$$
\frac{1}{2(b-a)(d-c)}\left(\sqrt{\frac{2}{\pi }} \sqrt{\sigma_1^2+\sigma_2^2} \left(e^{-\frac{(b-c)^2}{2 \left(\sigma_1^2+\sigma_2^2\right)}}-e^{-\frac{(a-c)^2}{2 \left(\sigma_1^2+\sigma_2^2\right)}}\right)+\sqrt{\frac{2}{\pi }} \sqrt{\sigma_1^2+\sigma_2^2} \left(e^{-\frac{(a-d)^2}{2 \left(\sigma_1^2+\sigma_2^2\right)}}-e^{-\frac{(b-d)^2}{2 \left(\sigma_1^2+\sigma_2^2\right)}}\right)+(a-c) \text{erf}\left(\frac{c-a}{\sqrt{2} \sqrt{\sigma_1^2+\sigma_2^2}}\right)+(d-a) \text{erf}\left(\frac{d-a}{\sqrt{2} \sqrt{\sigma_1^2+\sigma_2^2}}\right)-(b-c) \text{erf}\left(\frac{c-b}{\sqrt{2} \sqrt{\sigma_1^2+\sigma_2^2}}\right)+(b-d) \text{erf}\left(\frac{d-b}{\sqrt{2} \sqrt{\sigma_1^2+\sigma_2^2}}\right)\right)
$$
You can obtain a result for the second integral in terms of Owen's T-function which is defined as
$$
T(x, p) = \frac{1}{2\pi}\int_{0}^{p} \frac{e^{-\frac{1}{2} x^2 (1+t^2)}}{1+t^2} \, dt \quad \left(-\infty < x, p < +\infty\right).
$$
From this definition, a simple calculation using Feynman's trick tells us that
$\require{cancel}$
\begin{align*}
\frac{d}{d\beta} 2\sqrt{\frac{\pi}{\alpha}}T\left(\sqrt{\frac{2\alpha}{\alpha+1}}\beta, \frac{1}{\sqrt{\alpha}} \right)& = \cancel{2}\sqrt{\frac{\pi}{\alpha}} \frac{1}{\cancel{2}\pi}\int_{0}^{ \frac{1}{\sqrt{\alpha}}}\frac{\partial}{\partial \beta} \frac{e^{-\frac{1}{\cancel{2}}\frac{\cancel{2}\alpha}{\alpha+1} \beta^2 (1+t^2)}}{1+t^2} dt\\
& = \frac{1}{\sqrt{\pi \alpha}}\int_{0}^{ \frac{1}{\sqrt{\alpha}}} \frac{-2\frac{\alpha}{\alpha+1}\cancel{(1+t^2)}\beta e^{-\frac{\alpha}{\alpha+1}\beta^2 (1+t^2)}}{\cancel{1+t^2}} dt\\
& = -\frac{2\beta}{\sqrt{\pi }} \frac{\sqrt{\alpha}}{\alpha+1}e^{-\frac{\alpha}{\alpha+1}\beta^2}\int_{0}^{ \frac{1}{\sqrt{\alpha}}} e^{-\left(\sqrt{\frac{\alpha}{\alpha+1}}\beta t\right)^2} \, dt\\
& \overset{\color{blue}{u = \sqrt{\frac{\alpha}{\alpha+1}}\beta t}}{=} -\frac{\color{green}{2}\cancel{\beta}}{\color{green}{\sqrt{\pi }}} \frac{\cancel{\sqrt{\alpha}}}{\alpha+1}e^{-\frac{\alpha}{\alpha+1}\beta^2}\frac{\color{blue}{\sqrt{\alpha+1}}}{\color{blue}{\cancel{\beta}}\cancel{\color{blue}{\sqrt{\alpha}}}} \color{green}{\int_{0}^{\frac{\beta}{\sqrt{\alpha+1}}} e^{-u^2} \, du}\\
& = - \frac{1}{\sqrt{\alpha+1}}e^{-\left(1-\frac{1}{\alpha+1}\right)\beta^2}\color{green}{\text{erf}\left( \frac{\beta}{\sqrt{\alpha+1}}\right)} \tag{1}
\end{align*}
Now let's define a function $I(\beta)$ in terms of the integral we want
$$
I(\beta) := \int_{0}^{\infty} e^{-\alpha x^2} \text{erf}(x+\beta) \, dx
$$
by again using Feynman's trick we get
\begin{align*}
\frac{d}{d\beta}I(\beta) = \int_{0}^{\infty}e^{-\alpha x^2} \frac{\partial}{\partial \beta }\text{erf}(x+\beta) \, dx = \int_{0}^{\infty}e^{-\alpha x^2} \frac{2}{\sqrt{\pi}}e^{-(x+\beta)^2} \, dx = \frac{2}{\sqrt{\pi}}\int_{0}^{\infty}e^{-\left(\left(\color{blue}{\sqrt{\alpha+1}}\right)^2x^2 +2\color{purple}{\beta}x +\color{green}{\beta^2}\right)} \, dx
\end{align*}
where we used the fundamental theorem of calculus to differentiate the error function, along with its definition $\text{erf}(x) = \frac{2}{\sqrt{\pi}}\int_{0}^{x} e^{-t^2} \, dt$. Now, from equation $4$ in section $3.2$ of the tables you linked we know that
$$
\int_{0}^{\infty} e^{-(\color{blue}{a}^2 t^2 + 2\color{purple}{z} t + \color{green}{c})} \, dt = \frac{\sqrt{\pi}}{2a}\text{erfc}\left( \frac{z}{a}\right)e^{-\left(c-\frac{z^2}{a^2}\right)} = \frac{\sqrt{\pi}}{2a}e^{-\left(c-\frac{z^2}{a^2}\right)}-\frac{\sqrt{\pi}}{2a}\text{erf}\left( \frac{z}{a}\right)e^{-\left(c-\frac{z^2}{a^2}\right)}
$$
using that $\text{erfc}(x) = 1 - \text{erf}(x)$ by definition. We thus get that
$$
I'(\beta) = \frac{1}{\sqrt{\alpha+1}}e^{-\left(1-\frac{1}{\alpha+1}\right)\beta^2}- \frac{1}{\sqrt{\alpha+1}}\text{erf}\left( \frac{\beta}{\sqrt{\alpha+1}}\right)e^{-\left(1-\frac{1}{\alpha+1}\right)\beta^2} \tag{2}
$$
Using the definition of the error function, notice that
$$
\frac{d}{d\beta}\left[ \frac{1}{2} \sqrt{\frac{\pi}{\alpha}}\text{erf}\left(\sqrt{\frac{\alpha}{\alpha+1}} \beta\right) \right]= \frac{1}{\sqrt{\alpha+1}}e^{-\left(1-\frac{1}{\alpha+1}\right)\beta^2}
$$
And thus, we can combine the previous equation with equations $(1)$ and $(2)$ to get that
$$
\frac{d}{d\beta}I(\beta) = \frac{d}{d\beta}\left[\frac{1}{2} \sqrt{\frac{\pi}{\alpha}}\text{erf}\left(\sqrt{\frac{\alpha}{\alpha+1}} \beta\right) + 2\sqrt{\frac{\pi}{\alpha}}T\left(\sqrt{\frac{2\alpha}{\alpha+1}}\beta, \frac{1}{\sqrt{\alpha}} \right)\right]
$$
But since we have two functions whose derivatives are equal, the functions being differentiated must be equal up to some constant $\color{blue}{C}$:
$$
I(\beta) = \sqrt{\frac{\pi}{\alpha}}\left(\frac{1}{2} \text{erf}\left(\sqrt{\frac{\alpha}{\alpha+1}} \beta\right) + 2T\left(\sqrt{\frac{2\alpha}{\alpha+1}}\beta, \frac{1}{\sqrt{\alpha}} \right)\right) + \color{blue}{C}
$$
Evaluating at $\beta =0$ we get
$$
I(0) = \sqrt{\frac{\pi}{\alpha}}\left(\frac{1}{2} \text{erf}\left(0\right) + 2T\left(0, \frac{1}{\sqrt{\alpha}} \right)\right) +C = \sqrt{\frac{\pi}{\alpha}}\left( \frac{\text{arctan}\left(\frac{1}{\sqrt{\alpha}} \right)}{\pi}\right) +C \tag{3}
$$
where we used the fact that $T(0, x)= \frac{1}{2\pi}\int_{0}^{x} \frac{1}{1+t^2} \, dt= \frac{\arctan(x)}{2\pi}$. If we now evaluate $I(0)$ using our original integral definition we can obtain the value for $C$ and complete the problem. We see that
$$
I(0) = \int_{0}^{\infty} e^{-\sqrt{\alpha}^2 x^2} \text{erf}(x) \, dx
$$
and using equation $2$ from section $4.3$ of the tables you linked we get
\begin{align}
I(0) = \frac{\sqrt{\pi}}{2\sqrt{\alpha}} - \frac{1}{\sqrt{\pi\alpha}}\arctan(\sqrt{\alpha}) = \cancel{\frac{\sqrt{\pi}}{2\sqrt{\alpha}}}-\frac{1}{\sqrt{\pi\alpha}}\left(\cancel{\frac{\pi}{2}} - \arctan\left(\frac{1}{\sqrt{\alpha}}\right) \right) =\frac{\arctan\left(\frac{1}{\sqrt{\alpha}}\right)}{\sqrt{\pi \alpha}} \tag{4}
\end{align}
where we used that $\arctan(x) + \arctan\left( \frac{1}{x}\right)= \frac{\pi}{2}$ since $ \alpha > 0$ in order for the integral to converge. Combining equations $(3)$ and $(4)$ we get that $C=0$, so we can finally conclude that for $a>0$ and $b \in \mathbb{R}$:
$$
\boxed{\int_{0}^{\infty} e^{-\alpha x^2} \text{erf}(x+\beta) \, dx = \sqrt{\frac{\pi}{\alpha}}\left(\frac{1}{2} \text{erf}\left(\sqrt{\frac{\alpha}{\alpha+1}} \beta\right) + 2T\left(\sqrt{\frac{2\alpha}{\alpha+1}}\beta, \frac{1}{\sqrt{\alpha}} \right)\right)}
$$
The integral you give in your question can be derived from the general case above. First, notice that
$$
\frac{1}{\pi}\int_{0}^{1} \frac{e^{-\frac{x^2}{2} (1+t^2)}}{1+t^2} \, dt=2T\left(x,1 \right) = \frac{1}{4} \text{erfc}\left( \frac{x}{\sqrt{2}}\right) \left(2- \text{erfc}\left( \frac{x}{\sqrt{2}}\right) \right)
$$
which can be proven by differentiating both sides of the equation by $x$ and using Feynman's trick. We thus get that
\begin{align}
\int_{0}^{\infty} e^{- x^2} \text{erf}(x+\beta) \, dx &= \sqrt{\frac{\pi}{1}}\left(\frac{1}{2} \text{erf}\left(\sqrt{\frac{1}{1+1}} \beta\right) + 2T\left(\sqrt{\frac{\cancel{2}}{\cancel{1+1}}}\beta, \frac{1}{\sqrt{1}} \right)\right)\\
&= \sqrt{\pi}\left(\frac{1}{2} \text{erf}\left( \frac{\beta}{\sqrt{2}}\right) + 2T\left(\beta, 1\right)\right)\\
&= \sqrt{\pi}\left(\frac{1}{2} \left(1-\text{erfc}\left( \frac{\beta}{\sqrt{2}}\right)\right) + \frac{1}{4} \text{erfc}\left( \frac{\beta}{\sqrt{2}}\right) \left(2- \text{erfc}\left( \frac{\beta}{\sqrt{2}}\right) \right)\right)\\
&= \sqrt{\pi}\left(\frac{1}{2}-\cancel{\frac{1}{2}\text{erfc}\left( \frac{\beta}{\sqrt{2}}\right)} + \cancel{\frac{1}{2} \text{erfc}\left( \frac{\beta}{\sqrt{2}}\right)}- \frac{1}{4} \text{erfc}\left( \frac{\beta}{\sqrt{2}}\right)\text{erfc}\left( \frac{\beta}{\sqrt{2}}\right) \right)\\
&= \frac{\sqrt{\pi}}{2}\left(1-\frac{1}{2}\text{erfc}^2\left(\frac{\beta}{\sqrt{2}}\right)\right)
\end{align}
As a final note, a nice complimentary identity I found whilst solving this problem is the integral
$$
\int_{-\infty}^{\infty} e^{-\alpha x^2} \text{erf}(\gamma x + \beta) \, dx = \sqrt{\frac{\pi}{\alpha}} \text{erf}\left( \frac{\beta \sqrt{\alpha}}{\sqrt{\gamma^2 +\alpha}}\right)
$$
which might also be useful to you.
Best Answer
I have asked this question because my intention was to calculate the Fourier transform: \begin{equation} \int_{-\infty}^{+\infty} e^{-\alpha^2t^2+\beta t} \text{erf}\left(i\alpha t + \mu \right) e^{-i\omega t}dt = \mathcal{F}\left\lbrace e^{-\alpha^2t^2+\beta t} \big(\text{erf}\left(i\alpha t + \mu\right) \right\rbrace \triangleq I(\omega) \end{equation} Maybe I have found something that can be of use to someone.
Consider the Faddeeva function (also known as Kramp function, complex error function, plasma dispersion function): \begin{equation} w(z) \triangleq e^{-z^2}\left(1+\frac{2i}{\sqrt\pi}\int_{0}^{z}e^{u^2}du\right) \equiv e^{-z^2}\left(1+\text{erf}(iz)\right) \\ \end{equation}
Its integral representation is equal to the convolution of a Gaussian with a simple pole, i.e. to the Hilbert transform of a Gaussian: \begin{equation} w(z) = \frac{i}{\pi}\int_{-\infty}^{+\infty}\frac{e^{-u^2}}{z-u}du = e^{-z^2} * \frac{i}{\pi z} = i\mathcal{H}\left\lbrace e^{-\alpha^2z^2}\right\rbrace\quad \quad \Im z >0 \\ \end{equation} This observation simplifies the calculation of its Fourier transform: \begin{equation} \begin{aligned} \mathcal{F}\left\lbrace e^{-\alpha^2t^2}\text{erf}(i\alpha t) \right\rbrace &= \mathcal{F}\left\lbrace w(\alpha t)-e^{-\alpha^2t^2}\right\rbrace \\ &=\mathcal{F}\left\lbrace e^{-\alpha^2t^2} * \frac{i}{\pi t}-e^{-\alpha^2t^2}\right\rbrace \\ &=\mathcal{F}\left\lbrace e^{-\alpha^2t^2}\right\rbrace \mathcal{F}\left\lbrace \frac{i}{\pi t}\right\rbrace-\mathcal{F}\left\lbrace e^{-\alpha^2t^2}\right\rbrace \\ &=\frac{\sqrt\pi}{\alpha}\exp\left(\frac{i\omega}{2\alpha}\right)^2\big(\text{sgn}(\omega)-1\big) %&=\frac{\sqrt\pi}{\alpha}e^{\left(\frac{i\omega}{2\alpha}\right)^2}\left(\text{sgn}(\omega)-1\right) \end{aligned} \end{equation} The Fourier transform $\mathcal{F}\left\lbrace\frac{i}{\pi t}\right\rbrace = \text{sgn}(\omega)$ is obtained by taking the Cauchy principal value of the integral.
Now we can rearrange the integrand of $I(\omega)$ in order to isolate a $e^{-\alpha^2t^2}\text{erf}(i\alpha t)$ term. By defining the variable change $t' \triangleq t + \frac{\mu}{i\alpha}$ we obtain: \begin{equation} \begin{aligned} I(\omega) &= \int_{-\infty}^{+\infty} e^{-\alpha^2t^2+\beta t} \text{erf}\left(i\alpha t + \mu \right) e^{-i\omega t}dt \\ &=\int_{-\infty}^{+\infty} e^{-\alpha^2\left(t'-\frac{\mu}{i\alpha}\right)^2+(\beta-i\omega)\left(t'-\frac{\mu}{i\alpha}\right)} \text{erf}\left(i\alpha t'\right) dt'\\ %&=\int_{-\infty}^{+\infty} e^{-\alpha^2 t'^2-2i\alpha\mu_qt'+ \mu^2+a_0t'-\beta\frac{\mu_q}{i\alpha}-i\omega t'+i\omega \frac{\mu_q}{i\alpha}+\mu^2} \text{erf}\left(i\alpha t'\right) dt'\\ %&=\int_{-\infty}^{+\infty} e^{-\alpha^2 t'^2+s_qt'+(i\omega-a_0) \frac{\mu}{i\alpha}+\mu_q^2} \text{erf}\left(i\alpha t'\right) e^{-i\omega t'} dt'\\ &= e^{\mu^2+\frac{\mu}{i\alpha}(i\omega-\beta)} \mathcal{F}\left\lbrace e^{(\beta - 2i\alpha\mu)t'}e^{-\alpha^2t'^2}\text{erf}(i\alpha t')\right\rbrace\\ &= \frac{\sqrt\pi}{\alpha}e^{\left(\frac{i\omega-\beta}{2\alpha}\right)^2} \big(\text{sgn}(i\omega-\beta+2i\alpha\mu)-1\big) \end{aligned} \end{equation}