Analytical Solutions of One-Dimensional Contaminant Transport in Soils with Source Production-Decay

Abstract

An analytical solution in closed form of the advection-dispersion equation in one-dimensional contaminated soils is proposed in this paper. This is valid for non-conservative solutes with first order reaction, linear equilibrium sorption, and a time-dependent Robin boundary condition. The Robin boundary condition is expressed as a combined production-decay function representing a realistic description of the source release phenomena in time. The proposed model is particularly useful to describe sources as the contaminant release due to the failure in underground tanks or pipelines, Non Aqueous Phase Liquid pools, or radioactive decay series. The developed analytical model tends towards the known analytical solutions for particular values of the rate constants.

1. Introduction

The study of contaminant transport in porous media is of great importance for scientists working in environmental engineering, hydrology, geology, soil physics, chemical, and petroleum engineering. The concentration distribution in space and time, of a reactant pollutant in groundwater, can be described by the advection-dispersion equation (ADE). This equation contains terms taking into account the physical process governing the advective transport, molecular diffusion, and hydrodynamic dispersion, chemical reaction in the liquid phase, contaminant decay or production, and adsorption/desorption processes with the solid phase.

The literature contains many analytical solutions of the ADE both in open and closed form. These analytical solutions are useful for providing pollution scenarios in risk analysis, to investigate the effects of chemical-physical parameters on pollutant transport, and also to validate the numerical models. The most adopted methods to find the analytical solutions are Fourier analysis, Laplace transform, and expansion in Bessel or Hankel functions [1,2,3,4,5,6].

Exact analytical solutions of the advection-dispersion equation are seldom possible if the considered problem to solve takes into consideration some complex chemical processes such as non-linear sorption and multi-component reaction chains considered together. Some, solutions taking into account non-linear adsorption/desorption, have been developed for one-dimensional problems [7,8]. However, multi-dimensional exact analytical solutions in semi-finite and finite domain under first- or third-type boundary conditions are usually expressed in open form and include integrals to be numerically evaluated, or contain an infinite series in their expressions [9,10,11,12,13,14,15,16].

Leij et al. and Wang et al. [17,18,19] provided a series of multi-dimensional analytical solutions by using the Green Function Method (GFM) [20] for various plane and linear sources. With this method, under some particular initial and boundary conditions, a more complex solution can be constructed by combining the single directional solutions. Multi-dimensional closed form approximated analytical solutions can also be constructed by applying GFM under some simplifying hypothesis [21,22,23]. For this reason, the study and the development of one-dimensional analytical solutions in closed form is of high interest. Furthermore, a one-dimensional solution can be suitable to model contaminant in laboratory columns, pollution in aquifers where the contamination source extends through the saturated zone in either transverse y, and z directions or cases where dispersion in y and z directions can be neglected. Moreover, a one-dimensional solution as it is can be applied in calculating vertical flow of contaminants in saturated soil.

Many one-dimensional exact analytical solutions in closed form have been proposed, and scientists’ attention has been dedicated to solutions having time-dependent sources [24]. A library of one-dimensional analytical models that encloses some solutions with source decay was proposed by Van Genuchten et al. [25]. Guerrero et al. [26] proposed a Duhamel theorem based approach in order to compute one-dimensional solutions with time dependent boundary conditions and gave some solutions in closed form for exponential source decay.

We propose here a closed form analytical solution of the one-dimensional ADE in semi-finite domain for a reacting solute under first order decay and linear sorption described by a retardation factor. The source of contamination is described by a Robin time-dependent boundary condition, representing a combined production-decay release mechanism.

In [27,28,29] a technique for solving the ADE for Dirichlet, the Neumann and Cauchy problem is proposed, and the technique is valid for multi coupled multi-species processes, considers exponential dependence in time of the initial source, and utilizes Laplace transforms. Even if a scientific intersection between this approach and the proposed one is present, in this work, by considering only one component advection-dispersion-reaction process, we are able to consider more general initial data condition, like the box case, and we are able to manage the problem in a simple way, by reducing the whole problem to the calculation of an unique integral.

The solution so proposed is of simple use, is written in an analytical form that includes the Van Genuchten solutions [25] for specific boundary conditions, and is easily adoptable for developing multi-dimensional analytical solutions by using the Green Function Method.

2. Materials and Methods

2.1. Governing Equation

The one-dimensional advection-dispersion equation governing contaminant transport subject to first order decay and linear equilibrium sorption in a semi-finite or a finite domain of length L₀ can be written as:

$$R\frac{\partial C(x,t)}{\partial t}=D\frac{{\partial }^{2}C(x,t)}{\partial x^2}-v\frac{\partial C(x,t)}{\partial x}-\lambda C(x,t)$$

(1)

Subject to the following initial and boundary conditions:

$$C(x,0)=0$$

(2)

$$C(0,t)=h(t)$$

(3a)

or

$$-D\frac{\partial C(0,t)}{\partial x}+vC(0,t)=vg(t)$$

(3b)

Moreover, we can assume one of the following homogeneous Neumann conditions, whose physical meaning is to exclude material flux at the exit (finite or infinite domain):

$$\frac{\partial C(L_0,t)}{\partial x}=0$$

(4a)

or

$$\frac{\partial C(\infty ,t)}{\partial x}=0$$

(4b)

where C(x,t) is the contaminant concentration [ML⁻³], D is the dispersion coefficient [L²T⁻¹], v is the average space velocity [LT⁻¹], λ is the first order decay constant [T⁻¹], R is the retardation factor [-], h(t) and g(t) are general time-dependent functions [ML⁻³]. Equations (1)–(4) are valid for a saturated homogeneous porous medium, constant and uniform fluid flow, and dispersion coefficients where the source of contamination is written as a Dirichlet or Robin boundary conditions.

2.2. The Proposed Analytical Solution

Here we consider R = 1 (thus, no adsorption phenomena are considered), even though the reader can generalize the proposed solution to arbitrary values of R by substituting to the velocity ν the retarded velocity ν_r = ν/R, to the longitudinal dispersion coefficient D the new coefficient D_r defined as D_r = D/R, and finally to the decay coefficient λ the new one λ_r = λ/R [30].

The analytical solution is derived for a general inlet distribution subject to Robin boundary conditions and semi-finite domain.

Let us consider Equation (1) subject to the following initial and boundary conditions:

$$C(x,0)=0$$

(5)

$$-D\frac{\partial C(0,t)}{\partial x}+vC(0,t)=vg(t)$$

(6)

And to the following constraint:

$$\frac{\partial C(\infty ,t)}{\partial x}=0$$

(7)

Let us introduce the new variable c(x,t):

$$C(x,t)=\exp \left(\frac{vx}{D}\right)\exp \left(-\lambda t\right)c(x,t)$$

(8)

In this way, it is possible to shift the problem subject to Robin condition to a problem subject to Neumann condition.

By deriving Equation (8) with respect to x, we obtain:

$$\frac{\partial C(x,t)}{\partial x}=\frac{v\exp \left(\frac{vx}{D}\right)\exp \left(-\lambda t\right)c(x,t)}{D}+\exp \left(\frac{vx}{D}\right)\exp \left(-\lambda t\right)\left(\frac{\partial c(x,t)}{\partial x}\right)$$

(9)

By making the substitution of Equation (8) in the first term of the right-hand side of Equation (9) this equation can be rewritten as:

$$\frac{\partial C(x,t)}{\partial x}=\frac{vC(x,t)}{D}+\exp \left(\frac{vx}{D}\right)\exp \left(-\lambda t\right)\left(\frac{\partial c(x,t)}{\partial x}\right)$$

(10)

By evaluating Equation (10) in x = 0:

$$\frac{\partial C(0,t)}{\partial x}-\frac{vC(0,t)}{D}=\exp \left(-\lambda t\right)\left(\frac{\partial c(0,t)}{\partial x}\right)$$

(11)

Knowing that:

$$-D\frac{\partial C(0,t)}{\partial x}+vC(0,t)=vg(t)$$

(12)

Doing some passages, Equation (12) becomes:

$$\frac{\partial C(0,t)}{\partial x}-\frac{vC(0,t)}{D}=-\frac{vg(t)}{D}$$

(13)

By setting $\frac{vg(t)}{D}=f(t)$ we obtain:

$$\frac{\partial C(0,t)}{\partial x}-\frac{vC(0,t)}{D}=-f(t)$$

(14)

By comparing Equations (11) and (14) we obtain that:

$$-f(t)=\exp \left(-\lambda t\right)\left(\frac{\partial c(0,t)}{\partial x}\right)$$

(15)

and finally:

$$\left(\frac{\partial c(0,t)}{\partial x}\right)=-\frac{f(t)}{\exp \left(-\lambda t\right)}$$

(16)

So, the initial general problem described by Equation (1) can be transformed into the new general problem described by:

$$\frac{\partial c(x,t)}{\partial t}=D\frac{{\partial }^{2}c(x,t)}{\partial x^2}-v\frac{\partial c(x,t)}{\partial x}$$

(17)

with the new boundary conditions given by Equation (16), that is a one-dimensional problem subject to second-type boundary conditions. Solving with Laplace transform in time, the general solution is the following:

$$c(x,t)=-\frac{v{n}_1\left(x,t\right)}{2}+\frac{\sqrt{D}\exp \left(-\frac{vx}{2D}\right)\exp \left(\frac{-v^{2}t}{4D}\right)n_2\left(x,t\right)}{\sqrt{\pi }}$$

(18)

where:

$$n_1\left(x,t\right)={\displaystyle \underset{0}{\overset{t}{\int }}f(q)\exp \left(\lambda q\right)erfc\left(\frac{v\sqrt{t-q}+\frac{x}{\sqrt{t-q}}}{2\sqrt{D}}\right)}dq$$

(19)

$$n_2\left(x,t\right)={\displaystyle \underset{0}{\overset{t}{\int }}\frac{f(q)\exp \left(\lambda q\right)\exp \left(\frac{q{v}^2}{4D}\right)\exp \left(-\frac{x^2}{4D(t-q)}\right)}{\sqrt{t-q}}}dq$$

(20)

Equations (18) and (19) are related by the following equality:

$$\frac{\partial n_1\left(x,t\right)}{\partial x}=-\frac{\exp \left(-\frac{t{v}^2}{4D}\right)\exp \left(-\frac{vx}{2D}\right)n_2\left(x,t\right)}{\sqrt{\pi D}}$$

(21)

The integral general solution finally is:

$$c(x,t)=-\frac{v\hspace{0.17em}{n}_1\left(x,t\right)}{2}-D\frac{\partial n_1\left(x,t\right)}{\partial x}$$

(22)

3. Results and Discussion

3.1. Linear Combination of Exponential Inlet Distribution as a Robin (Third-Type) Boundary Conditions in Semi-Finite Domain

By setting:

$$g(t)=C_{0}y(1-e^{-{\lambda }_{p}t})+C_0{e}^{-{\lambda }_{s}t}$$

(23)

where $C_0$ is the initial concentration [ML⁻³] in the source and $y=\raisebox{1ex}{$C_0{}_{\infty }$}\!\left/ \!\raisebox{-1ex}{$C_0$}\right.$ is the ratio between the residual source concentration at the steady state and the initial concentration; ${\lambda }_p$ is the production constant [T⁻¹] and ${\lambda }_s$ is the decay constant [T⁻¹].

Equations (19) and (20) become:

$$n_1\left(x,t\right)={\displaystyle \underset{0}{\overset{t}{\int }}\frac{v}{D}{C}_{0}y(1-e^{-{\lambda }_{p}q})+C_0{e}^{-{\lambda }_{s}q}\exp \left(\lambda q\right)erfc\left(\frac{v\sqrt{t-q}+\frac{x}{\sqrt{t-q}}}{2\sqrt{D}}\right)}dq$$

(24)

$$n_2\left(x,t\right)={\displaystyle \underset{0}{\overset{t}{\int }}\frac{\frac{v}{D}{C}_{0}y(1-e^{-{\lambda }_{p}q})+C_0{e}^{-{\lambda }_{s}q}\exp \left(\lambda q\right)\exp \left(\frac{q{v}^2}{4D}\right)\exp \left(-\frac{x^2}{4D(t-q)}\right)}{\sqrt{t-q}}}dq$$

(25)

by solving Equations (22)–(25) with Equation (8) substitution, we get the final solution:

$$C\left(x,t\right)=C_{0}y\left(A\left(x,t\right)+B\left(x,t\right)\right)+C_{0}E\left(x,t\right)$$

(26)

$$\begin{array}{ll}A(x,t)=& \frac{v}{v+u}\exp \left(\frac{(v-u)x}{2D_x}\right)erfc\left(\frac{x-tu}{2\sqrt{D_{x}t}}\right)+\frac{v}{v-u}\exp \left(\frac{(v+u)x}{2D_x}\right)erfc\left(\frac{x+tu}{2\sqrt{D_{x}t}}\right)+\\ & +\frac{v^2}{2\lambda D_x}\exp \left(\frac{vx}{D_x}\right)\exp \left(-\lambda t\right)erfc\left(\frac{x+tv}{2\sqrt{D_{x}t}}\right)\end{array}$$

(27)

$$\begin{array}{ll}B(x,t)=& -\exp \left(-{\lambda }_{p}t\right)\left(\frac{v}{v+w}\exp \left(\frac{(v-w)x}{2D_x}\right)erfc\left(\frac{x-tw}{2\sqrt{D_{x}t}}\right)+\frac{v}{v-w}\exp \left(\frac{(v+w)x}{2D_x}\right)erfc\left(\frac{x+tw}{2\sqrt{D_{x}t}}\right)\right)+\\ & -\frac{v^2}{2\left(\lambda -{\lambda }_p\right)D_x}\exp \left(\frac{vx}{D_x}\right)\exp \left(-\lambda t\right)erfc\left(\frac{x+tv}{2\sqrt{D_{x}t}}\right)\hfill \end{array}$$

(28)

$$\begin{array}{ll}E(x,t)=& \exp \left(-{\lambda }_{s}t\right)\left(\frac{v}{v+W}\exp \left(\frac{(v-W)x}{2D_x}\right)erfc\left(\frac{x-tW}{2\sqrt{D_{x}t}}\right)+\frac{v}{v-W}\exp \left(\frac{(v+W)x}{2D_x}\right)erfc\left(\frac{x+tW}{2\sqrt{D_{x}t}}\right)\right)+\\ & +\frac{v^2}{2\left(\lambda -{\lambda }_s\right)D_x}\exp \left(\frac{vx}{D_x}\right)\exp \left(-\lambda t\right)erfc\left(\frac{x+tv}{2\sqrt{D_{x}t}}\right)\hfill \end{array}$$

(29)

$$u=\sqrt{4D_x\lambda +v^2},\hspace{0.17em}w=\sqrt{4D_x(\lambda -{\lambda }_p)+v^2},\hspace{0.17em}W=\sqrt{4D_x(\lambda -{\lambda }_s)+v^2}$$

That is the one-dimensional analytical solution here proposed in semi-finite domain. This solution can be derived, by doing numerous passages, from the general equation given by Srinivasan and Clement [28,29]. By setting zero initial condition and three species transport problem, the general equation given by Srinivasan and Clement [29] becomes:

$$c_i\left(x,t\right)={\displaystyle \sum _{i1=1}^i\left[\left({\displaystyle \prod _{i2=i1+1}^i{y}_{i2}{k}_{i2-1}}\right){\displaystyle \sum _{i2=i1}^i{\displaystyle \sum _{i3=1}^{i1}\left\{G_2^1\right\}}}\right]}\hspace{1em}i=1,2,3$$

with boundary conditions in line with [29].

It is possible to see that the proposed analytical solution reduces to known solutions for particular values of the production/decay constants.

Case 1: For ${\lambda }_p=0$ and ${\lambda }_s=0$ we have:

$$g(t)=C_0$$

(30)

$$W=w=u$$

(31)

$$E\left(x,t\right)=-F\left(x,t\right)$$

(32)

The solution, for this particular case, is:

$$C\left(x,t\right)=C_{0}A\left(x,t\right)$$

(33)

where:

$$\begin{array}{cc}A(x,t)=& \frac{v}{v+u}\exp \left(\frac{(v-u)x}{2D_x}\right)erfc\left(\frac{x-tu}{2\sqrt{D_{x}t}}\right)+\frac{v}{v-u}\exp \left(\frac{(v+u)x}{2D_x}\right)erfc\left(\frac{x+tu}{2\sqrt{D_{x}t}}\right)+\\ & +\frac{v^2}{2\lambda D_x}\exp \left(\frac{vx}{D_x}\right)\exp \left(-\lambda t\right)erfc\left(\frac{x+tv}{2\sqrt{D_{x}t}}\right)\hfill \end{array}$$

(34)

That is identical to case C6 described in [25], by setting $C_i=0$ and $\gamma =0$.

Case 2: For ${\lambda }_p\to \infty$ we have:

$$g(t)=C_{0}y+C_0{e}^{-{\lambda }_{s}t}$$

(35)

and the solution can be obtained by computing ${\lim }_{{\lambda }_p\to \infty }C(x,t)$.

We can easily observe that:

$$\begin{array}{ll}\underset{{\lambda }_p\to \infty }{\lim }B(x,t)=& \underset{{\lambda }_p\to \infty }{\lim }\left(-\exp \left(-{\lambda }_{p}t\right)\left(\frac{v}{v+w}\exp \left(\frac{(v-w)x}{2D_x}\right)erfc\left(\frac{x-tw}{2\sqrt{D_{x}t}}\right)+\frac{v}{v-w}\exp \left(\frac{(v+w)x}{2D_x}\right)erfc\left(\frac{x+tw}{2\sqrt{D_{x}t}}\right)\right)\right)+\\ & +\underset{{\lambda }_p\to \infty }{\lim }\left(-\frac{v^2}{2\left(\lambda -{\lambda }_p\right)D_x}\exp \left(\frac{vx}{D_x}\right)\exp \left(-\lambda t\right)erfc\left(\frac{x+tv}{2\sqrt{D_{x}t}}\right)\right)\end{array}$$

(36)

is equal to zero, hence the solution becomes:

$$C\left(x,t\right)=C_{0}yA\left(x,t\right)+C_{0}E\left(x,t\right)$$

(37)

where:

$$\begin{array}{ll}A(x,t)=& \frac{v}{v+u}\exp \left(\frac{(v-u)x}{2D_x}\right)erfc\left(\frac{x-tu}{2\sqrt{D_{x}t}}\right)+\frac{v}{v-u}\exp \left(\frac{(v+u)x}{2D_x}\right)erfc\left(\frac{x+tu}{2\sqrt{D_{x}t}}\right)+\\ & +\frac{v^2}{2\lambda D_x}\exp \left(\frac{vx}{D_x}\right)\exp \left(-\lambda t\right)erfc\left(\frac{x+tv}{2\sqrt{D_{x}t}}\right)\hfill \end{array}$$

(38)

$$\begin{array}{ll}E(x,t)=& \exp \left(-{\lambda }_{s}t\right)\left(\frac{v}{v+W}\exp \left(\frac{(v-W)x}{2D_x}\right)erfc\left(\frac{x-tW}{2\sqrt{D_{x}t}}\right)+\frac{v}{v-W}\exp \left(\frac{(v+W)x}{2D_x}\right)erfc\left(\frac{x+tW}{2\sqrt{D_{x}t}}\right)\right)+\\ & +\frac{v^2}{2\left(\lambda -{\lambda }_s\right)D_x}\exp \left(\frac{vx}{D_x}\right)\exp \left(-\lambda t\right)erfc\left(\frac{x+tv}{2\sqrt{D_{x}t}}\right)\hfill \end{array}$$

(39)

That is identical to case C14 described in [25], by setting $C_i=0$ and $\gamma =0$.

3.2. Linear Combination of Exponential Inlet Distribution as a First-Type Boundary Conditions in Semi-Finite Domain

Model presented in Section 3.1 is subject to Robin boundary conditions. It is important to remark that, even if these boundary conditions well describe mass conservation, a large number of software are implemented by using model based on first-type or Dirichlet boundary conditions in which the concentration at the source is specified. However, solutions derived on the basis of first-type boundary conditions do not conserve mass if concentrations are taken as volume-average i.e., the amount of solute per unit of volume.

One dimensional systems subject to Dirichlet boundary conditions can be mass conservative if concentrations are interpreted as flux-average i.e., the quantity of solute per fluid unit volume that passes the cross-section area in a unit of time [31]. Several papers highlight and describe the importance to do a distinction between the resident or volume-average concentration C_r and the flowing or flux-averaged concentration C_f in one-dimensional solute transport [32,33,34].

The above mentioned studies gave the following correlation between C_r and C_f:

$$C_f=C_r-\frac{D_x}{v}\frac{\partial C_r}{\partial x}$$

(40)

Equation above gives the relationship between flux-averaged concentration C_f and volume-average concentration C_r.

It was proved that this formulation leaves one-dimensional transport equation invariant, it follows that transport equation can be written both in term of C_r and C_f [35,36]. By considering Equation (1) with following boundary and initial conditions:

$$C(x,0)=0$$

(41)

$$C(0,t)=h(t)$$

(42)

$$\frac{\partial C(\infty ,t)}{\partial x}=0$$

(43)

$$h(t)=C_{0}y(1-e^{-{\lambda }_{p}t})+C_0{e}^{-{\lambda }_{s}t}$$

(44)

it is possible to obtain a solution following passages in Appendix A.

The solution is:

$$C_f\left(x,t\right)=C_{0}y\left(a(x,t)+b(x,t)\right)+C_{0}e(x,t)$$

(45)

where:

$$\begin{array}{ll}a(x,t)=& \left(\frac{1}{2}erfc\left(-\frac{1}{2}\frac{tu-x}{\sqrt{D_{x}t}}\right)+\frac{\sqrt{D_x}\exp \left(-\frac{1}{4}\frac{{\left(tu-x\right)}^2}{D_{x}t}\right)}{\sqrt{\pi t}\left(u+v\right)}\right)\exp \left(\frac{1}{2}\frac{x\left(v-u\right)}{D_x}\right)\\ & +\left(\frac{1}{2}erfc\left(\frac{1}{2}\frac{tu+x}{\sqrt{D_{x}t}}\right)+\frac{\sqrt{D_x}\exp \left(-\frac{1}{4}\frac{{\left(tu+x\right)}^2}{D_{x}t}\right)}{\sqrt{\pi t}\left(v-u\right)}\right)\exp \left(\frac{1}{2}\frac{x\left(u+v\right)}{D_x}\right)\hfill \\ & +\frac{1}{2}\frac{\exp \left(\frac{vx}{D_x}\right)\exp \left(-t\lambda \right)\exp \left(-\frac{1}{4}\frac{{\left(tv+x\right)}^2}{D_{x}t}\right)v}{\lambda \sqrt{\pi t\sqrt{D_x}}}\hfill \end{array}$$

(46)

$$\begin{array}{ll}b(x,t)=& \left(-\frac{1}{2}erfc\left(-\frac{1}{2}\frac{tw-x}{\sqrt{D_{x}t}}\right)-\frac{\sqrt{D_x}\exp \left(-\frac{1}{4}\frac{{\left(tw-x\right)}^2}{D_{x}t}\right)}{\sqrt{\pi t}\left(v+w\right)}\right)\exp \left(\frac{1}{2}\frac{x(v-w)}{D_x}\right)\exp \left(-t{\lambda }_p\right)\\ & +\left(-\frac{1}{2}erfc\left(\frac{1}{2}\frac{tw+x}{\sqrt{D_{x}t}}\right)-\frac{\sqrt{D_x}\exp \left(-\frac{1}{4}\frac{{\left(tw+x\right)}^2}{D_{x}t}\right)}{\sqrt{\pi t}\left(v-w\right)}\right)\exp \left(\frac{1}{2}\frac{x\left(v+w\right)}{D_x}\right)\exp \left(-t{\lambda }_p\right)\hfill \\ & -\frac{1}{2}\frac{\exp \left(\frac{vx}{D_x}\right)\exp \left(-t\lambda \right)\exp \left(-\frac{1}{4}\frac{{\left(tv+x\right)}^2}{D_{x}t}\right)v}{\left(\lambda -{\lambda }_p\right)\sqrt{\pi t\sqrt{D_x}}}\hfill \end{array}$$

(47)

$$\begin{array}{ll}e(x,t)=& -\left(-\frac{1}{2}erfc\left(-\frac{1}{2}\frac{tW-x}{\sqrt{D_{x}t}}\right)-\frac{\sqrt{D_x}\exp \left(-\frac{1}{4}\frac{{\left(tW-x\right)}^2}{D_{x}t}\right)}{\sqrt{\pi t}\left(v+w\right)}\right)\exp \left(\frac{1}{2}\frac{x(v-W)}{D_x}\right)\exp \left(-t{\lambda }_s\right)\\ & -\left(-\frac{1}{2}erfc\left(\frac{1}{2}\frac{tW+x}{\sqrt{D_{x}t}}\right)-\frac{\sqrt{D_x}\exp \left(-\frac{1}{4}\frac{{\left(tW+x\right)}^2}{D_{x}t}\right)}{\sqrt{\pi t}\left(v-W\right)}\right)\exp \left(\frac{1}{2}\frac{x\left(v+W\right)}{D_x}\right)\exp \left(-t{\lambda }_s\right)\hfill \\ & +\frac{1}{2}\frac{\exp \left(\frac{vx}{D_x}\right)\exp \left(-t\lambda \right)\exp \left(-\frac{1}{4}\frac{{\left(tv+x\right)}^2}{D_{x}t}\right)v}{\left(\lambda -{\lambda }_s\right)\sqrt{\pi t\sqrt{D_x}}}\hfill \end{array}$$

(48)

It is possible to see that the proposed analytical solution for the first type boundary conditions reduces to known solutions for particular values of the production/decay constants.

Case 1: For ${\lambda }_p=0$ and ${\lambda }_s=0$ we have:

$$h(t)=C_0$$

(49)

$$W=w=u$$

(50)

$$C\left(x,t\right)=C_{0}a\left(x,t\right)$$

(51)

where:

$$a(x,t)=\frac{1}{2}\left(erfc\left(-\frac{1}{2}\frac{tu-x}{\sqrt{D_{x}t}}\right)\exp \left(\frac{1}{2}\frac{x\left(v-u\right)}{D_x}\right)+erfc\left(\frac{1}{2}\frac{tu+x}{\sqrt{D_{x}t}}\right)\exp \left(\frac{1}{2}\frac{x\left(u+v\right)}{D_x}\right)\right)$$

(52)

That can be written in the following form that is equal to case C5 described in [25], with $C_i=0$ and $\gamma =0$.

$$C(x,t)=\frac{C_0}{2}\exp \left(\frac{vx}{2D_x}\right)\left\{\exp \left(-\frac{ux}{2D_x}\right)erfc\left[\frac{x-ut}{2\sqrt{D_{x}t}}\right]+\exp \left(\frac{ux}{D_x}\right)erfc\left[\frac{x+ut}{2\sqrt{D_{x}t}}\right]\right\}$$

(53)

Case 2: For ${\lambda }_p\to \infty$ we have:

$$h(t)=C_{0}y+C_0{e}^{-{\lambda }_{s}t}$$

(54)

and the solution can be obtained by computing ${\lim }_{{\lambda }_p\to \infty }C(x,t)$ that is:

$$\begin{array}{ll}{\lim }_{{\lambda }_p\to \infty }b(x,t)={\lim }_{{\lambda }_p\to \infty }& \left(-\frac{1}{2}erfc\left(-\frac{1}{2}\frac{tw-x}{\sqrt{D_{x}t}}\right)-\frac{\sqrt{D_x}\exp \left(-\frac{1}{4}\frac{{\left(tw-x\right)}^2}{D_{x}t}\right)}{\sqrt{\pi t}\left(v+w\right)}\right)\exp \left(\frac{1}{2}\frac{x(v-w)}{D_x}\right)\exp \left(-t{\lambda }_p\right)\\ & +\left(-\frac{1}{2}erfc\left(\frac{1}{2}\frac{tw+x}{\sqrt{D_{x}t}}\right)-\frac{\sqrt{D_x}\exp \left(-\frac{1}{4}\frac{{\left(tw+x\right)}^2}{D_{x}t}\right)}{\sqrt{\pi t}\left(v-w\right)}\right)\exp \left(\frac{1}{2}\frac{x\left(v+w\right)}{D_x}\right)\exp \left(-t{\lambda }_p\right)\\ & -\frac{1}{2}\frac{\exp \left(\frac{vx}{D_x}\right)\exp \left(-t\lambda \right)\exp \left(-\frac{1}{4}\frac{{\left(tv+x\right)}^2}{D_{x}t}\right)v}{\left(\lambda -{\lambda }_p\right)\sqrt{\pi t\sqrt{D_x}}}\end{array}$$

(55)

is equal to zero, hence the solution becomes:

$$C\left(x,t\right)=C_{0}ya\left(x,t\right)+C_{0}e\left(x,t\right)$$

(56)

where:

$$a(x,t)=\frac{1}{2}\left(erfc\left(-\frac{1}{2}\frac{tu-x}{\sqrt{D_{x}t}}\right)\exp \left(\frac{1}{2}\frac{x\left(v-u\right)}{D_x}\right)+erfc\left(\frac{1}{2}\frac{tu+x}{\sqrt{D_{x}t}}\right)\exp \left(\frac{1}{2}\frac{x\left(u+v\right)}{D_x}\right)\right)$$

(57)

$$e(x,t)=\frac{1}{2}erfc\left(-\frac{1}{2}\frac{tW-x}{\sqrt{D_{x}t}}\right)\exp \left(\frac{1}{2}\frac{x(v-W)}{D_x}\right)\exp \left(-t{\lambda }_s\right)+\frac{1}{2}erfc\left(\frac{1}{2}\frac{tW+x}{\sqrt{D_{x}t}}\right)\exp \left(\frac{1}{2}\frac{x\left(v+W\right)}{D_x}\right)\exp \left(-t{\lambda }_s\right)$$

(58)

That is identical to case C13 described in [25] by setting $C_i=0$ and $\gamma =0$.

3.3. Consecutive Reactions at the Source as a Robin (Third-Type) Boundary Conditions in Semi-Finite Domain

By considering consecutive reactions at the source written as:

$$A\to B\to C$$

It is possible to simulate the release of contaminant B that produces at the source from the decay of a pool of contaminant A decaying during time into product C. This solution can be particularly useful for waste radioactive decay at the source or PCE (tetrachloroethylene) to TCE (trichloroethylene) degradation in soils.

$$g(t)=C_{0}K\left(-e^{-{\lambda }_{p}t}+e^{-{\lambda }_{s}t}\right)$$

(59)

where $C_0$ is the initial concentration [ML⁻³] and $K=\frac{k_1}{k_2-k_1}$ is the global chemical kinetic rate.

$$C\left(x,t\right)=C_{0}K\left(A\left(x,t\right)+B\left(x,t\right)\right)$$

(60)

$$\begin{array}{ll}A(x,t)=& -\exp \left(-{\lambda }_{p}t\right)\left(\frac{v}{v+w}\exp \left(\frac{(v-w)x}{2D_x}\right)erfc\left(\frac{x-tw}{2\sqrt{D_{x}t}}\right)+\frac{v}{v-w}\exp \left(\frac{(v+w)x}{2D_x}\right)erfc\left(\frac{x+tw}{2\sqrt{D_{x}t}}\right)\right)+\\ & -\frac{v^2}{2\left(\lambda -B\right)D_x}\exp \left(\frac{vx}{D_x}\right)\exp \left(-\lambda t\right)erfc\left(\frac{x+tv}{2\sqrt{D_{x}t}}\right)\hfill \end{array}$$

(61)

$$\begin{array}{ll}B(x,t)=& \exp \left(-{\lambda }_{s}t\right)\left(\frac{v}{v+W}\exp \left(\frac{(v-W)x}{2D_x}\right)erfc\left(\frac{x-tW}{2\sqrt{D_{x}t}}\right)+\frac{v}{v-W}\exp \left(\frac{(v+W)x}{2D_x}\right)erfc\left(\frac{x+tW}{2\sqrt{D_{x}t}}\right)\right)+\\ & +\frac{v^2}{2\left(\lambda -S\right)D_x}\exp \left(\frac{vx}{D_x}\right)\exp \left(-\lambda t\right)erfc\left(\frac{x+tv}{2\sqrt{D_{x}t}}\right)\hfill \end{array}$$

(62)

$$w=\sqrt{4D_x(\lambda -B)+v^2},\hspace{0.17em}W=\sqrt{4D_x(\lambda -S)+v^2}$$

3.4. Consecutive Reactions at the Source as First-Type Boundary Conditions in Semi-Finite Domain

By following passages in Appendix A it is possible to derive, also in this case, the solution for consecutive reactions in semi-finite domain subject to first-type boundary conditions.

$$C(x,0)=0$$

(63)

$$C(0,t)=h(t)$$

(64)

$$\frac{\partial C(\infty ,t)}{\partial x}=0$$

(65)

$$h(t)=C_{0}K\left(-e^{-{\lambda }_{p}t}+e^{-{\lambda }_{s}t}\right)$$

(66)

$$C_f\left(x,t\right)=C_{0}K\left(a(x,t)+b(x,t)\right)$$

(67)

$$\begin{array}{ll}a(x,t)=& \left(-\frac{1}{2}erfc\left(-\frac{1}{2}\frac{tw-x}{\sqrt{D_{x}t}}\right)-\frac{\sqrt{D_x}\exp \left(-\frac{1}{4}\frac{{\left(tw-x\right)}^2}{D_{x}t}\right)}{\sqrt{\pi t}\left(w+v\right)}\right)\exp \left(\frac{1}{2}\frac{x\left(v-w\right)}{D_x}\right)\exp (-t{\lambda }_p)\\ & +\left(-\frac{1}{2}erfc\left(\frac{1}{2}\frac{tw+x}{\sqrt{D_{x}t}}\right)+\frac{\sqrt{D_x}\exp \left(-\frac{1}{4}\frac{{\left(tw+x\right)}^2}{D_{x}t}\right)}{\sqrt{\pi t}\left(v-w\right)}\right)\exp \left(\frac{1}{2}\frac{x\left(w+v\right)}{D_x}\right)\exp (-t{\lambda }_p)\hfill \\ & -\frac{1}{2}\frac{\exp \left(\frac{vx}{D_x}\right)\exp \left(-t\lambda \right)\exp \left(-\frac{1}{4}\frac{{\left(tv+x\right)}^2}{D_{x}t}\right)v}{(\lambda -{\lambda }_p)\sqrt{\pi t\sqrt{D_x}}}\hfill \end{array}$$

(68)

$$\begin{array}{ll}b(x,t)=& -\left(-\frac{1}{2}erfc\left(-\frac{1}{2}\frac{tW-x}{\sqrt{D_{x}t}}\right)-\frac{\sqrt{D_x}\exp \left(-\frac{1}{4}\frac{{\left(tW-x\right)}^2}{D_{x}t}\right)}{\sqrt{\pi t}\left(v+w\right)}\right)\exp \left(\frac{1}{2}\frac{x(v-W)}{D_x}\right)\exp \left(-t{\lambda }_s\right)\\ & -\left(-\frac{1}{2}erfc\left(\frac{1}{2}\frac{tW+x}{\sqrt{D_{x}t}}\right)-\frac{\sqrt{D_x}\exp \left(-\frac{1}{4}\frac{{\left(tW+x\right)}^2}{D_{x}t}\right)}{\sqrt{\pi t}\left(v-W\right)}\right)\exp \left(\frac{1}{2}\frac{x\left(v+W\right)}{D_x}\right)\exp \left(-t{\lambda }_s\right)\hfill \\ & +\frac{1}{2}\frac{\exp \left(\frac{vx}{D_x}\right)\exp \left(-t\lambda \right)\exp \left(-\frac{1}{4}\frac{{\left(tv+x\right)}^2}{D_{x}t}\right)v}{\left(\lambda -{\lambda }_s\right)\sqrt{\pi t\sqrt{D_x}}}\hfill \end{array}$$

(69)

3.5. Finite Release at the Source as a Robin (Third-Type) Boundary Conditions in Semi-Finite Domain

By setting:

$$g(t)=C_{\infty }{e}^{-\lambda t}H(\Delta -t)$$

(70)

where $H(\Delta -t)$ is the Heaviside function, and $C_{\infty }$ is the residual source concentration, ${\lambda }^{}$ is the first order reaction constant.

The solution is:

$$C(x,t)=e^{\frac{vx}{D_x}}{e}^{-\lambda t}(A(x,t)H(\Delta -t)+B(x,t)H(t-\Delta ))$$

(71)

$$\begin{array}{ll}A(x,t)=& \left(-\frac{D_x{C}_{\infty }}{v+u}erfc\left(\frac{tu-x}{2\sqrt{D_{x}t}}\right)+2\frac{D_x{C}_{\infty }}{v+u}\right)\exp \left(t\lambda -\frac{(v+u)x}{2D_x}\right)\\ & -\frac{C_{\infty }(u+v)}{4\lambda }\exp \left(t\lambda +\frac{(u-v)x}{2D_x}\right)erfc\left(\frac{tu+x}{2\sqrt{D_{x}t}}\right)\hfill \\ & +\frac{v{C}_{\infty }}{2\lambda }erfc\left(\frac{tu+x}{2\sqrt{D_{x}t}}\right)\hfill \end{array}$$

(72)

$$\begin{array}{ll}B(x,t)=& \left(\frac{v{C}_{\infty }}{2\lambda }erfc\left(-\frac{(t-\Delta )v+x}{2\sqrt{D_x(t-\Delta )}}\right)-\frac{v{C}_{\infty }}{\lambda }\right)\exp (\Delta \lambda )+\frac{D_x{C}_{\infty }}{v+u}\exp \left(t\lambda -\frac{(v+u)x}{2D_x}\right)\\ & \left(-erfc\left(\frac{x+(\Delta -t)u}{2\sqrt{D_x(t-\Delta )}}\right)-erfc\left(-\frac{x-tu}{2\sqrt{D_{x}t}}\right)+2\right)+\frac{D_x{C}_{\infty }}{u-v}\exp \left(t\lambda +\frac{(u-v)x}{2D_x}\right)\hfill \\ & \left(-erfc\left(-\frac{x+(t-\Delta )u}{2\sqrt{D_x(t-\Delta )}}\right)-erfc\left(\frac{x+tu}{2\sqrt{D_{x}t}}\right)+2\right)+\frac{v{C}_{\infty }}{2\lambda }erfc\left(\frac{x+tv}{2\sqrt{D_{x}t}}\right)\hfill \end{array}$$

(73)

3.6. Example of Application

The sequential decay of multi-species contaminants as nitrogen, chlorinated solvent, and radionuclide is of large importance in soil contamination. A lot of researchers have developed a large number of models involving sets of advective-dispersive transport equations coupled by first-order decay [37,38] but analytical solutions in closed form that considers both production and decay at the source are scarcely available [22,23]. In order to illustrate the model and the effect of the boundary conditions, the following example of application has been considered: one-dimensional transport of an intermediate radionuclide, subject to decay chain during its movement into groundwater, where it is possible to consider its production-decay also at the source.

$${}^{238}Pu\to {}^{234}U\to {}^{230}Th$$

Input parameters, here used, are kept from [38,39] except for input concentration value and are summarized in

Table 1. Radionuclide decay parameters.

Parameters	234U
Velocity, v [m year−1]	100
Dispersion coefficient [m2 year−1]	4000
First order decay constant λ [year−1]	0.0000028
Production constant λp [year−1]	0.079
Source decay constant λs [year−1]	0.0010028
Initial concentration of ${}^{238}Pu$ [mg m3]	5 × 104

.

It is important to notice that in the proposed model the values of parameters λ and λ_s can be completely independent.

A time of 100 years was taken into consideration [38,39,40].

Results of simulation involving Equations (60)–(62) are reported in Figure 1, Figure 2 and Figure 3. Results of simulation regarding Equations (67)–(69) are shown in Figure 4 and Figure 5.

Figure 2 and Figure 4 represent concentration profiles of U²³⁴ at x = 0 for Robin and Dirichlet boundary conditions respectively. As it is possible to see, main difference is in the initial value at t = 0. It is important to notice that the limit for t tending to zero of Equation (67) when x = 0, does not satisfy boundary conditions: it is expected that concentration of U²³⁴ at x = 0 for t equal to zero equals g(0) for Dirichlet boundary conditions, so the solution in Equation (67) is valid for all x > 0.

Figure 3 and Figure 5 represent concentration of U²³⁴ at t = 100 years along x for the two models, according to radionuclide decay. The profile computed under Robin boundary conditions reaches a maximum in concentration at a certain x (Robin conditions well describe mass conservation) while at the same t the solution computed under Dirichlet boundary conditions shows lower values of concentration.

4. Conclusions

The contaminant transport in porous media has an important role in environmental engineering science, hydrology, geology, and petroleum engineering. Despite the fact that the Advection-Dispersion Equations describe the physical problem, it is difficult to derive analytical solutions if chemical processes also take place. Analytical solutions are helpful for testing numerical models and for screening level environmental risk assessment.

The 1D–model proposed in this paper and its analytical solution can be suitable to define tank or pipeline failure, or pollution by DNAPLs (Dense NonAqueous Phase Liquids) or by radioactive contaminants, because the source is described as a time dependent function with a combined contaminant production-consumption, and moreover, it is given as a Robin boundary condition. The comparison of the proposed solution with the closed form analytical solution obtained for a time dependent production-decay source expressed as a Dirichlet boundary condition shows that the use of the Dirichlet boundary condition can underestimate concentration profiles.

The analytical solution is consistent with previous ones, since it tends towards known solutions for certain values of the parameters (limit cases) and it can be derived from a more general one for some particular conditions at the source.