New Semi-Analytical Solutions for Advection–Dispersion Equations in Multilayer Porous Media

Abstract

A new semi-analytical solution to the advection–dispersion–reaction equation for modelling solute transport in layered porous media is derived using the Laplace transform. Our solution approach involves introducing unknown functions representing the dispersive flux at the interfaces between adjacent layers, allowing the multilayer problem to be solved separately on each layer in the Laplace domain before being numerically inverted back to the time domain. The derived solution is applicable to the most general form of linear advection–dispersion–reaction equation, a finite medium comprising an arbitrary number of layers, continuity of concentration and dispersive flux at the interfaces between adjacent layers and transient boundary conditions of arbitrary type at the inlet and outlet. The derived semi-analytical solution extends and addresses deficiencies of existing analytical solutions in a layered medium, which consider analogous processes such as diffusion or reaction–diffusion only and/or require the solution of complicated nonlinear transcendental equations to evaluate the solution expressions. Code implementing our semi-analytical solution is supplied and applied to a selection of test cases, with the reported results in excellent agreement with a standard numerical solution and other analytical results available in the literature.

Appendix: Numerical Solution

In this appendix, we briefly outline the finite volume scheme used to obtain a numerical solution to the multilayer transport model (1)–(6), as mentioned in Sect. 4.

The interval [0, L] is discretized using a uniform grid consisting of n nodes with the kth node located at $x = (k-1)h =: x_{k}$, where $k = 1,\ldots ,n$ and $h = L/(n-1)$. The number of nodes n is chosen to ensure that a node coincides with each interface ($x = \ell _{i}$, $i=1,\ldots ,m-1$). Let $\overline{c}_{k}(t)$ denote the numerical approximation to c(x, t) at $x = x_{k}$ and

$$\begin{aligned} J_{i,k}&= D_{i}\frac{\overline{c}_{k}-\overline{c}_{k-1}}{h} - v_{i}\frac{\overline{c}_{k-1}+\overline{c}_{k}}{2},\\ S_{i,k}&= -\mu _{i}\overline{c}_{k} + \gamma _{i}. \end{aligned}$$

The discrete system takes the form:

$$\begin{aligned} \mathbf {M}\frac{\text {d}\mathbf {c}}{\text {d}t} = \mathbf {F}(\mathbf {c}),\quad \mathbf {c}(0) = \mathbf {c}_{0}, \end{aligned}$$

(38)

where $\mathbf {M}\in \mathbb {R}^{n\times n}$, $\mathbf {c} = (\overline{c}_{1},\ldots ,\overline{c}_{n})^{T}\in \mathbb {R}^{n}$ and $\mathbf {F} = (F_{1},\ldots ,F_{n})^{T}\in \mathbb {R}^{n}$. The initial solution vector $\mathbf {c}_{0}\in \mathbb {R}^{n}$ gets its entries from the initial condition (2) with first entry $f_{1}$, last entry $f_{m}$ and kth entry ($k = 2,\ldots ,n-1$) equal to $f_{i}$ if $x_{k}\in (\ell _{i-1},\ell _{i})$ or $(f_{i}+f_{i+1})/2$ if $x_{k} = \ell _{i}$, where $i = 1,\ldots ,m-1$ is the interface index. The form of $\mathbf {M}$ depends on the choice of boundary conditions at the inlet and outlet:

$$\begin{aligned} \mathbf {M} = \left\{ \begin{array}{ll} \mathbf {I}, &\quad{} \text {if }\;b_{0}\ne 0\text { and }b_{L}\ne 0,\\ \mathbf {I} - \mathbf {e}_{1}\mathbf {e}_{1}^{T}, &\quad{} \text {if }\;b_{0} = 0\text { and }b_{L}\ne 0,\\ \mathbf {I} - \mathbf {e}_{n}\mathbf {e}_{n}^{T}, &\quad{} \text {if }\;b_{0} \ne 0\;\text { and }\;b_{L} = 0,\\ \mathbf {I} - \mathbf {e}_{1}\mathbf {e}_{1}^{T} - \mathbf {e}_{n}\mathbf {e}_{n}^{T}, &\quad{} \text {if }\;b_{0} = 0\;\text { and }\;b_{L} = 0,\\ \end{array}\right. \end{aligned}$$

where $\mathbf {I}$ is the $n\times n$ identity matrix and $\mathbf {e}_{k}$ is the kth column of $\mathbf {I}$. The components of $\mathbf {F}$ are defined as follows:

$$\begin{aligned} F_{1} = a_{0}\overline{c}_{1} - g_{0}(t) \end{aligned}$$

if $b_{0} = 0$,

$$\begin{aligned} F_{1} = \frac{J_{1,2} + \frac{D_{1}}{b_{0}}g_{0}(t) + (v_{1} - D_{1}\frac{a_{0}}{b_{0}})\overline{c}_{1} + \frac{h}{2}S_{1,1}}{\frac{h}{2}R_{1}} \end{aligned}$$

if $b_{0}\ne 0$,

$$\begin{aligned} F_{k} = \frac{J_{i,k+1} - J_{i,k} + hS_{i,k}}{hR_{i}} \end{aligned}$$

if $x_{k}\in (\ell _{i-1},\ell _{i})$,

$$\begin{aligned} F_{k} = \frac{\theta _{i+1}J_{i+1,k+1} - \theta _{i}J_{i,k} + \frac{h}{2}\left( \theta _{i}S_{i,k} + \theta _{i+1}S_{i+1,k}\right) }{\frac{h}{2}\left( \theta _{i}R_{i} + \theta _{i+1}R_{i+1}\right) } \end{aligned}$$

if $x_{k} = \ell _{i}$ and $i = 2,\ldots ,m-1$,

$$\begin{aligned} F_{n} = a_{L}\overline{c}_{n} - g_{L}(t) \end{aligned}$$

if $b_{L} = 0$, and

$$\begin{aligned} F_{n} = \frac{\frac{D_{m}}{b_{L}}g_{L}(t) - (v_{m} + D_{m}\frac{a_{L}}{b_{L}})\overline{c}_{n} - J_{m,n} + \frac{h}{2}S_{m,n}}{\frac{h}{2}R_{m}} \end{aligned}$$

if $b_{L}\ne 0$. The initial value problem (38) is solved using MATLAB’s in-built ode15s solver with the default tolerances and options of $\texttt {Mass} = \mathbf {M}$ and $\texttt {MassSingular} = \texttt {true}$ if $b_{0}$ and/or $b_{L}$ is zero. The interval of integration (tspan) is chosen to return the solution at the appropriate times shown in Figs. 2 and 3.