Joint Clarification of Contaminant Plume and Hydraulic Transmissivity via a Geostatistical Approach Using Hydraulic Head and Contaminant Concentration Data

Abstract

To enable proper remediation of accidental groundwater contamination, the contaminant plume evolution needs to be accurately estimated. In the estimation, uncertainties in both the contaminant source and hydrogeological structure should be considered, especially the temporal release history and hydraulic transmissivity. Although the release history can be estimated using geostatistical approaches, previous studies use the deterministic hydraulic property field. Geostatistical approaches can also effectively estimate an unknown heterogeneous transmissivity field via the use of joint data, such as a combination of hydraulic head and tracer data. However, tracer tests implemented over a contaminated area necessarily disturb the in situ condition of the contamination. Conversely, measurements of the transient concentration data over an area are possible and can preserve the conditions. Accordingly, this study develops a geostatistical method for the joint clarification of contaminant plume and transmissivity distributions using both head and contaminant concentration data. The applicability and effectiveness of the proposed method are demonstrated through two numerical experiments assuming a two-dimensional heterogeneous confined aquifer. The use of contaminant concentration data is key to accurate estimation of the transmissivity. The accuracy of the proposed method using both head and concentration data was verified achieving a high linear correlation coefficient of 0.97 between the true and estimated concentrations for both experiments, which was 0.67 or more than the results using only the head data. Furthermore, the uncertainty of the contaminant plume evolution was successfully evaluated by considering the uncertainties of both the initial plume and the transmissivity distributions, based on their conditional realizations.

Appendix: Normalization of the Prior Model

The geostatistical estimation of the transmissivity was implemented via the normalization of the prior model using the following procedure. More details are given in Kitanidis and Lee (2014).

List of symbols
Parameter		Dimension
All
t	Time
x	Space	m × 1
mt	Number of unknowns (source intensity)
m	Number of unknowns (contaminant plume/transmissivity)
nz	Number of observations (concentration)
nφ	Number of observations (head)
δ(x)	Dirac delta function
Estimation of initial contaminant plume distribution
s	Discretized unknown source intensity	mt × 1
z0	Discretized unknown initial contaminant plume	m × 1
${\varvec{z}}_{0}^{*}$	Observation (initial concentration)	nz × 1
Hs,${\varvec{H}}_{s}^{*}$	Jacobian matrix (for whole domain/observations)	$m \times m_{t} ,n_{z} \times m_{t}$
vz	Observation error	nz × 1
$\sigma_{{R_{z} }}$	Standard deviation of error
Rz	Error covariance matrix	nz × nz
Xs	Known drift matrix	mt × ps
βs	Unknown drift coefficients	ps × 1
Qs	Generalized prior covariance matrix of unknown function	mt × mt
θs	Structure parameter of covariance Qs
${\varvec{V}}_{{\hat{s}}}$	Covariance matrix of estimated source intensity	mt × mt
${\varvec{V}}_{{\hat{z}}}$	Covariance matrix of estimated initial contaminant plume	m × m
Estimation of hydraulic transmissivity distribution
n	Number of observations
r	Discretized unknown log-transmissivity	m × 1
z(t)	Discretized unknown contaminant plume at time t	m × 1
y	Observation	n × 1
$\user2{\varphi }$	Observation (head)	$n_{\varphi } \times 1$
${\varvec{z}}^{\user2{*}}$	Observation (transient concentration)	$n_{z} \times \left( {t_{0} , \ldots ,t_{end} } \right)$
$\overline{\user2{t}}$	Observation (mean travel time)	$n_{z} \times 1$
H	Jacobian matrix	n × m
v	Observation error	n × 1
$\sigma_{R}$	Standard deviation of error
R	Error covariance matrix	n × n
X	Known drift matrix	m × p
β	Unknown drift coefficients	p × 1
Q	Generalized prior covariance matrix of unknown function	m × m
θr	Structure parameter of covariance Q
K	Rank of approximation of Q
${\varvec{V}}_{{\hat{\user2{r}}}}$	Covariance matrix of estimated transmissivity	m × m
δr	Finite difference interval
εr	relative machine precision
$\lambda_{i}$	ith eigenvalue of Q
Vi	ith eigenvector of Q	m × 1
δls	Finite difference interval for line search
Conditional realization
$N_{{z_{0} }}$	Number of realizations (initial contaminant plume)
Nr	Number of realizations (transmissivity)
${\varvec{s}}_{{\varvec{u}}} , {\varvec{s}}_{{\varvec{c}}}$	Unconditional/conditional realization of s	mt × 1
${\varvec{z}}_{0_{{\varvec{c}}}}$	Conditional realization of z0	m × 1
${\varvec{r}}_{{\varvec{u}}} , {\varvec{r}}_{{\varvec{c}}}$	Unconditional/conditional realization of r	m × 1
Physical model
u	Groundwater velocity	2 × 1
T	Transmissivity	m × 1
Q f	Pumping rate
V	Actual groundwater velocity	2 × 1
Rf	Retardation factor
λf	Decay constant
D	Dispersion tensor	2 × 2
ε	Porosity
$\alpha_{L} , \alpha_{T}$	Dispersivity (longitudinal and transverse)
Dm	Molecular diffusion coefficient
τ	Tortuosity
Lp	Plume length

First, the drift matrix X is replaced with its normalized and isomorphic matrix U such that

$$ \begin{array}{*{20}c} {{\varvec {U}} = \left\{ {\begin{array}{*{20}c} {{\varvec{X}}/\sqrt m } & {\left( {p = 1} \right)} \\ {{\varvec{US}}_{X} {\varvec{V}}_{X}^{T} } & {\left( {p > 1} \right)} \\ \end{array} ,} \right.} \\ \end{array} $$

(28)

where ${\varvec{S}}_{X} \in {\mathbb{R}}^{p \times p}$ is a diagonal matrix of the singular values and the columns of ${\varvec{V}}_{X}^{{}} \in {\mathbb{R}}^{p \times p}$ are the orthonormal eigenvectors of ${\varvec{X}}^{T} {\varvec{X}}$. Then, U is used to compute the detrending matrix ${\varvec{P}} \in {\mathbb{R}}^{m \times m}$

$$ \begin{array}{*{20}c} {{\varvec {P}} = {\varvec {I}} - {\varvec {U}}{\varvec{U}}^{T} .} \\ \end{array} $$

(29)

The next step is to detrend the low-rank covariance Q [Eq. (20)]. First, Z_Q is replaced with ${\varvec{PZ}}_{Q}$. Then, the singular value decomposition of Z_Q is calculated such that

$$ \begin{array}{*{20}c} {{\varvec{Z}}_{Q} = {\varvec{U}}_{Z} {\varvec{S}}_{Z} {\varvec{V}}_{Z}^{T} ,} \\ \end{array} $$

(30)

where ${\varvec{U}}_{Z} \in {\mathbb{R}}^{m \times K}$ is the unitary matrix; ${\varvec{S}}_{Z} \in {\mathbb{R}}^{K \times K}$ is a diagonal matrix of the singular values; and the columns of ${\varvec{V}}_{Z}^{{}} \in {\mathbb{R}}^{K \times K}$ are the orthonormal eigenvectors of ${\varvec{Z}}_{Q}^{T} {\varvec{Z}}_{Q}$. Using ${\varvec{C}} = {\varvec{U}}_{Z}^{T} {\varvec{QU}}_{Z}$, Q is replaced with its detrended and isomorphic matrix PQP

$$ \begin{array}{*{20}c} {{\varvec {PQP}} \approx {\varvec{U}}_{Z} {\varvec {C}}{\varvec{U}}_{Z}^{T} .} \\ \end{array} $$

(31)

Then, HQ and ${\varvec{HQH}}^{T}$ can be approximated as

$$ \begin{array}{*{20}c} {{\varvec {HQ}} \approx \left( {{\varvec{HU}}_{Z} } \right){\varvec {C}}{\varvec{U}}_{Z}^{T} \equiv {\varvec {BC}}{\varvec{U}}_{Z}^{T} ,} \\ \end{array} $$

(32)

$$ \begin{array}{*{20}c} {{\varvec {HQ}}{\varvec{H}}^{T} \approx \left( {{\varvec{HU}}_{Z} } \right){\varvec {C}}\left( {{\varvec{HU}}_{Z} } \right)^{T} = {\varvec {BC}}{\varvec{B}}^{T} .} \\ \end{array} $$

(33)

Finally, $\overline{\user2{r}}$ is updated iteratively by the following linear equation system corresponding to Eq. (6)

$$ \begin{array}{*{20}c} {\left( {\begin{array}{*{20}c}{\varvec{\varSigma}}& {{\varvec{HX}}} \\ {\left( {{\varvec{HX}}} \right)^{T} } & 0 \\ \end{array} } \right)\left( {\begin{array}{*{20}c} {{\varvec{\varLambda}}^{T} {\varvec{A}}_{p} } \\ {{\varvec{MA}}_{p} } \\ \end{array} } \right) = \left( {\begin{array}{*{20}c} {{\varvec{HQA}}_{p} } \\ {{\varvec{X}}^{T} {\varvec{A}}_{p} } \\ \end{array} } \right), {\varvec{A}}_{p} = \left( {{\varvec{U}}_{Z} , {\varvec{U}}} \right),} \\ \end{array} $$

(34)

$$ \begin{array}{*{20}c} {\overline{\user2{r}} = {\varvec{A}}_{p} \left( {{\varvec{\varLambda}}^{T} {\varvec{U}}_{Z} } \right)^{T} \left( {{\varvec{y}} - h\left( {\overline{\user2{r}}} \right) + \user2{H\overline{r}}} \right).} \\ \end{array} $$

(35)

Once the optimal solution $\hat{\user2{r}}$ is obtained, the posterior covariance ${\varvec {V}}_{{\hat{r}}}$ can be approximately calculated as

$$ \begin{array}{*{20}c} {{\varvec{V}}_{{\hat{r}}} = {\varvec{U}}_{Z} {\varvec {C}}{\varvec{U}}_{Z}^{T} - {\varvec {X}}\left( {{\varvec{MA}}_{p} } \right){\varvec{A}}_{p}^{T} - {\varvec {Q}}{\varvec{H}}^{T} \left( {{\varvec{\varLambda}}^{T} {\varvec{A}}_{p} } \right){\varvec{A}}_{p}^{T} .} \\ \end{array} $$

(36)