Solution Tutorial — 2D Equilibrium Sorption with 1st Order Decay

This tutorial presents the two analytical solutions implemented in the 2D Equilibrium Sorption with 1st Order Decay applet: the exact Cleary-Ungs solution (evaluated by Gauss-Legendre quadrature) and the approximate Domenico solution (closed-form). Both describe 2D transient solute transport in a semi-infinite, homogeneous aquifer with a finite-width continuous source, linear equilibrium sorption, and optional first-order decay.

1Introduction

The physical domain is an infinite 2D horizontal aquifer. The x-axis is aligned with the uniform pore-water velocity $v$. A continuous line source of concentration $C_0$ is located at $x = 0$, extending from $-Y/2$ to $+Y/2$ in the transverse ($y$) direction. The source is on at all times ($t > 0$).

Figure 1. 2D aquifer domain: finite-width continuous source at x = 0, uniform flow in +x direction.

The model assumptions are:

Uniform, steady pore-water velocity $v$ in the $+x$ direction
Linear equilibrium sorption: $S = K_d\,c$, giving retardation factor $R = 1 + \rho K_d/\theta$
First-order aqueous-phase decay at rate $\lambda$
Molecular diffusion negligible: $D_x = \alpha_x v$, $D_y = \alpha_y v$
Homogeneous, isotropic aquifer properties

2Governing Equation

The retarded advection-dispersion equation with first-order decay in two dimensions is:

$$R\,\frac{\partial c}{\partial t} = D_x\frac{\partial^2 c}{\partial x^2} + D_y\frac{\partial^2 c}{\partial y^2} - v\frac{\partial c}{\partial x} - \lambda c$$

(1)

where $R = 1 + \rho K_d/\theta$ is the retardation factor, $\theta$ is porosity, $\rho$ is soil bulk density, $K_d$ is the linear sorption partition coefficient, $\lambda$ is the first-order decay rate [1/T], and $D_x = \alpha_x v$, $D_y = \alpha_y v$ are the longitudinal and transverse dispersion coefficients, with $\alpha_x$ and $\alpha_y$ the corresponding dispersivities.

3Initial & Boundary Conditions

The aquifer is initially solute-free, and the following boundary conditions apply:

$$c = C_0 \qquad x = 0 \text{ and } -\tfrac{Y}{2} < y < \tfrac{Y}{2}$$ $$c = 0 \qquad x = 0 \text{ and } \left(y < -\tfrac{Y}{2} \text{ or } y > \tfrac{Y}{2}\right)$$ $$c,\;\frac{\partial c}{\partial x} = 0 \qquad x \to \infty \qquad\qquad c,\;\frac{\partial c}{\partial y} = 0 \qquad y \to \pm\infty$$

(2)

A first-type (Dirichlet) boundary condition is applied at $x = 0$: the concentration is fixed at $C_0$ within the source zone and at zero outside. This differs from a flux (third-type) condition; at points very close to the source, the two conditions give different results, but they converge far from the source.

4Cleary-Ungs Exact Solution

The exact solution to Eq. (1) subject to conditions (2) is given by Cleary & Ungs (1978) and documented in Wexler (1992):

$$c(x,y,t) = \frac{C_0\,x}{4\sqrt{\pi\,\alpha_x\,v}}\,\exp\!\left(\frac{x}{2\alpha_x}\right) \int_0^{t/R} \tau^{-3/2} \exp\!\left[-\left(\frac{v}{4\alpha_x} + \lambda\right)\tau - \frac{x^2}{4\alpha_x v\,\tau}\right]$$ $$\times\left\{\mathrm{erf}\!\left[\frac{Y/2 - y}{2\sqrt{\alpha_y v\,\tau}}\right] + \mathrm{erf}\!\left[\frac{Y/2 + y}{2\sqrt{\alpha_y v\,\tau}}\right]\right\}d\tau$$

(3)

where $\mathrm{erf}(\cdot)$ is the error function. The integral must be evaluated numerically. The applet uses Gauss-Legendre quadrature as described by Wexler (1992), replacing the integral with a weighted sum over a fixed number of quadrature points (20, 60, or 96).

Accuracy note: For large Péclet numbers ($x/\alpha_x \gg 1$), the integrand becomes sharply peaked and the numerical quadrature may produce oscillations with fewer points. If oscillations appear, switch to 96 integration points.

5Domenico Approximate Solution

Domenico and co-workers developed a closed-form approximate solution based on simplifying the time-dependence of the transverse spreading. The 2D form is obtained as the limit of the 3D solution for a large vertical source (Domenico 1987; Wiedemeier et al. 1999):

$$c(x,y,t) = \frac{C_0}{4}\exp\!\left[\frac{x}{2\alpha_x}\!\left(1 - \sqrt{1 + \frac{4\lambda\alpha_x}{v}}\right)\right]$$ $$\times\;\mathrm{erfc}\!\left[\frac{x - \dfrac{v\,t}{R}\sqrt{1 + \dfrac{4\lambda\alpha_x}{v}}}{2\sqrt{\alpha_x v\,t/R}}\right]$$ $$\times\left[\mathrm{erf}\!\left(\frac{y + Y/2}{2\sqrt{\alpha_y x}}\right) - \mathrm{erf}\!\left(\frac{y - Y/2}{2\sqrt{\alpha_y x}}\right)\right]$$

(4)

where $\mathrm{erfc}(\cdot) = 1 - \mathrm{erf}(\cdot)$ is the complementary error function.

Validity Restrictions

Due to the nature of the approximations, the Domenico solution is not valid at short times or close to the source. The following conditions should be satisfied for results to be within ~5% of the true value:

Condition	Description
$t \gtrsim 5\,\alpha_x / v$	Near-steady-state condition: source zone concentration is within 5% of $C_0$
$x / \alpha_x \gtrsim 30$	Far-field condition: concentration along the x-axis at the advective front ($x = vt$) is within 5% of its true value

Warning: At short times or near the source (small $x/\alpha_x$), the Domenico solution may overestimate or underestimate concentrations significantly. Use the Cleary-Ungs solution or check the absolute difference plot to assess the error.

6Comparing the Two Solutions

The applet provides an |Dom − Cleary-Ungs| display mode that maps the absolute pointwise difference between the two solutions across the entire domain. This is the most direct way to assess where and by how much the Domenico approximation deviates from the exact result.

The difference is typically largest near the source ($x/\alpha_x \lesssim 30$) and at early times ($t \lesssim 5\,\alpha_x/v$)
In the far field at quasi-steady state, the two solutions often agree to within a few percent
The difference increases for high decay rates ($\lambda$) and strong retardation ($R$), because the Domenico approximation handles these terms differently

Tip: Run both solutions with identical parameters, then switch to the difference view to identify the spatial regions where the Domenico approximation is unreliable for your specific problem.

7Symbol Definitions

Symbol	Name	Units
$c$	Aqueous-phase solute concentration	M/L³ fluid
$C_0$	Source zone concentration at x = 0	M/L³
$v$	Pore-water velocity in the +x direction	L/T
$\alpha_x$	Longitudinal dispersivity	L
$\alpha_y$	Transverse dispersivity	L
$D_x$	Longitudinal dispersion coefficient: $D_x = \alpha_x v$	L²/T
$D_y$	Transverse dispersion coefficient: $D_y = \alpha_y v$	L²/T
$\theta$	Soil porosity	—
$\rho$	Soil bulk density	M soil / L³
$K_d$	Linear sorption partition coefficient	L³ fluid / M soil
$R$	Retardation factor: $R = 1 + \rho K_d/\theta$	—
$\lambda$	First-order aqueous-phase decay rate	1/T
$Y$	Source width in the transverse y-direction (source spans −Y/2 to +Y/2)	L
$x$	Longitudinal coordinate (positive downstream)	L
$y$	Transverse coordinate	L
$t$	Time	T
$\mathrm{erf}(\cdot)$	Error function	—
$\mathrm{erfc}(\cdot)$	Complementary error function: $\mathrm{erfc}(z) = 1 - \mathrm{erf}(z)$	—

8References

Cleary, R. W. & Ungs, M. J. (1978). Analytical models for groundwater pollution and hydrology. Princeton University Water Resources Program Report 78-WR-15. Princeton University, Princeton, NJ.
Domenico, P. A. (1987). An analytical model for multidimensional transport of a decaying contaminant species. Journal of Hydrology, 91, 49–58.
Domenico, P. A. & Robbins, G. A. (1985). A new method of contaminant plume analysis. Ground Water, 23(4), 476–485.
Wexler, E. J. (1992). Analytical solutions for one-, two-, and three-dimensional solute transport in ground-water flow systems with uniform flow. Techniques of Water-Resources Investigations of the United States Geological Survey, Chapter B7, Book 3.
Wiedemeier, T. H., Rifai, H. S., Newell, C. J., & Wilson, J. T. (1999). Natural Attenuation of Fuels and Chlorinated Solvents in the Subsurface. John Wiley & Sons, New York.

► Launch the Model ← Interface Tutorial