Groundwater (3D) Modeling
Contaminant Transport, Uncertainty Quantification, & Vulnerabilities
Executive Overview
Groundwater flow and contaminant transport models are often utilized to support plume delineation, remedial design, risk assessment, and long-term performance predictions. However, predictive reliability is highly sensitive to structural assumptions, parameter uncertainty, and conceptual site model limitations. Groundwater models are often considered reliable when calibration residuals are within acceptable ranges. However, calibration fit does not ensure predictive reliability. Structural non-uniqueness, parameter compensation, and sensitivity to boundary conditions often remain unexamined.
Independent groundwater modeling review evaluates whether governing equations, boundary conditions, calibration procedures, and uncertainty quantification methods are scientifically defensible and appropriately constrained by available data.
Governing Flow Equation
Most numerical groundwater models are based on the groundwater flow equation derived from Darcy’s Law and mass conservation:
Darcy Flux
$$q = -K \nabla h$$
Seepage velocity
$$v = \frac{\theta}{q}$$
- = Darcy flux
- = hydraulic conductivity tensor
- = hydraulic head
- = effective porosity
$$\frac{\partial}{\partial x}\!\left(K_x \frac{\partial h}{\partial x}\right)
+\frac{\partial}{\partial y}\!\left(K_y \frac{\partial h}{\partial y}\right)
+\frac{\partial}{\partial z}\!\left(K_z \frac{\partial h}{\partial z}\right)
+ W
=
S_s \frac{\partial h}{\partial t}$$
For unconfined aquifers:
$$S_y \frac{\partial h}{\partial t}$$
- ∂/∂x, ∂/∂y, ∂/∂z = partial derivative with respect to the x, y, z-directions
- ∂h/∂x, ∂h/∂y, ∂h/∂z = hydraulic gradient in the x, y, z-directions
- W (Source/Sink Term) = volumetric flux per unit volume
- Sₛ = specific storage
- Sy = specific yield
- ∂h/∂t = rate of change of hydraulic head over time
Common Vulnerabilities:
- Over-simplified hydraulic conductivity zoning.
- Insufficient vertical discretization.
- Unrealistic recharge assumptions.
- Boundary conditions not physically supported.
Contaminant Transport Equation
Contaminant migration is generally described using the advection–dispersion equation (ADE):
$$\frac{\partial (R C)}{\partial t}
=
\nabla \cdot (D \nabla C)
–
\nabla \cdot (\mathbf{v} C)
+
S
–
\lambda C$$
=
\theta \left(
D_{xx} \frac{\partial^2 C}{\partial x^2}
+
D_{yy} \frac{\partial^2 C}{\partial y^2}
+
D_{zz} \frac{\partial^2 C}{\partial z^2}
\right)
–
\left(
q_x \frac{\partial C}{\partial x}
+
q_y \frac{\partial C}{\partial y}
+
q_z \frac{\partial C}{\partial z}
\right)
+
\theta R$$
- = dissolved concentration
- = retardation factor
- = dispersion tensor
- = seepage velocity
- = source/sink term
- = decay constant
Retardation & Sorption
$$R = 1 + \frac{\theta \rho_b}{K_d}$$
- = bulk density
- = distribution coefficient
- = porosity
Literature-derived distribution coefficient (Kd) values may not reflect site-specific geochemistry, leading to overconfident travel time projections.
Sensitivity Analysis & Uncertainty
Key uncertain parameters:
- Sorption coefficient (Kd) variability.
- Interfacial adsorption coefficients.
- Hydraulic conductivity.
- van Genuchten parameters.
- Recharge rates.
- Soil layering.
Monte Carlo approach:
$$Var(Y) \approx \frac{1}{N-1} \sum (Y_i – \bar{Y})^2$$
- Var(Y) = Sample variance of Y
- N = Number of samples
- Yᵢ = Individual realization of Y
- Ȳ = Sample mean
- Yᵢ − Ȳ = Deviation from the mean
Structural vulnerability develops when a single deterministic prediction is presented without accompanying uncertainty bounds.
Source Term Representation
For many sites, uncertainty is primarily influenced by source characterization.
$$S = \frac{M_0}{V} e^{-k t}$$
- = initial mass
- = source volume
- = mass depletion rate constant
Oversimplified assumptions about source depletion may lead to an underestimation of long-term plume persistence.
Hydrodynamic Dispersion
Dispersion tensor:
$$D = \alpha_L v + \alpha_T v + D_m$$
- = longitudinal dispersivity
- = transverse dispersivity
- = molecular diffusion
Literature dispersivity values without scale justification can artificially smooth plume fronts.
Calibration & Parameter Estimation
Model calibration seeks to minimize residuals between observed and simulated values.
Objective function:
$$\Phi = \sum_{i=1}^{n} w_i \left( h_i^{\mathrm{obs}} – h_i^{\mathrm{sim}} \right)^2$$
- = weighting factor
- = observed head
- = simulated head
A low objective function does not guarantee predictive reliability. Calibration may compensate for structural error through parameter adjustment (non-uniqueness).
Uncertainty Quantification
Predictive reliability requires quantifying parameter and structural uncertainty. Variance-based uncertainty propagation:
$$\sigma_C^2
=
\sum_i \left( \frac{\partial C}{\partial p_i} \right)^2 \sigma_{p_i}^2
+
2 \sum_i \sum_j
\frac{\partial C}{\partial p_i}
\frac{\partial C}{\partial p_j}
\operatorname{Cov}(p_i, p_j)$$
- C = Model output concentration
- σ_C = Standard deviation of prediction C
- ∂C/∂pᵢ = Sensitivity coefficient
- σ_pᵢ² = Variance of parameter pᵢ
- Cov(pᵢ , pⱼ) = Covariance between parameters pᵢ and pⱼ
Monte Carlo simulation:
$$C_j = f(p_{1j}, p_{2j}, \ldots, p_{nj})$$
- = uncertain parameters
- = parameter variance
Inverse modeling uncertainty quantification is presented in the appendix.
Sensitivity Coefficient
$$S_{p_i} = \frac{\partial C}{\partial p_i} \cdot \frac{p_i}{C}$$
High sensitivity parameters necessitate more stringent site-specific constraints.
Predictive Uncertainty Bounds
Confidence intervals may be expressed as:
$$C_{\text{upper/lower}} = C_{\text{mean}} \pm z \sigma_C$$
- = confidence multiplier (e.g., 1.96 for 95%)
Models presented without predictive uncertainty bounds may overstate confidence.
Mass Balance Evaluation
A defensible model must satisfy mass conservation:
$$\text{Total Inflow} – \text{Total Outflow} = \Delta \text{Storage}$$
Structural Uncertainty
Parameter uncertainty accounts for only a portion of the risk. Structural uncertainty includes:
- Missing heterogeneity.
- Oversimplified conceptual site model.
- Incorrect boundary conditions.
- Ignored preferential pathways.
Conclusion
Groundwater modeling is fundamentally influenced by structural assumptions and the uncertainty of parameters. Numerical calibration may yield visually satisfactory fits; however, the reliability of predictions is contingent upon thorough uncertainty quantification, robust conceptual models, and clear communication of limitations.
Mathematical Appendix: Parameter Estimation & Predictive Uncertainty
Groundwater and contaminant transport predictions inherently contain uncertainty. Modern calibration and inversion approaches quantify this uncertainty rather than masking it behind a single “best-fit” result. The following equations summarize the conceptual framework used in advanced model evaluation and defensibility review.
Weighted Least Squares Objective Function
$$
\Phi(\mathbf{p})
=
\sum_{i=1}^{n}
w_i
\left(
y_i^{obs}
–
y_i^{sim}(\mathbf{p})
\right)^2
$$
- Φ(p) = objective function
- p = parameter vector
- wᵢ = weighting factor
- yᵢ(obs) = observed data
- yᵢ(sim)(p) = model output
Model calibration process involves minimizing the objective function (Φ). A low residual error does not remove structural or predictive uncertainty. Multiple parameter combinations may produce similar fits values (non-uniqueness).
Likelihood Function (Gaussian Residual Assumption)
$$
L(\mathbf{p})
=
\prod_{i=1}^{n}
\frac{1}{\sqrt{2\pi\sigma_i^2}}
\exp
\left(
–
\frac{
\left(
y_i^{obs}
–
y_i^{sim}(\mathbf{p})
\right)^2
}{
2\sigma_i^2
}
\right)
$$
- σᵢ² = measurement error variance
- exp = exponential function
This defines the likelihood of observing measured data given the parameter set P. Larger residuals reduce likelihood. Underestimating measurement uncertainty can falsely increase model confidence.
Posterior Distribution (Bayesian Inference)
$$
P(\mathbf{p} \mid \mathbf{y})
\propto
L(\mathbf{p})
\cdot
P(\mathbf{p})
$$
- P(p | y) = posterior distribution
- L(p) = likelihood
- P(p) = prior parameter distribution
Posterior parameter distributions integrate calibration likelihood with prior parameter uncertainty. Parameter estimates remain uncertain even after calibration. This equation expresses how data reduce — but do not eliminate — uncertainty.
Multi-Gaussian Parameter (K) Field
$$
\ln K(\mathbf{x})
\sim
\mathcal{N}
\left(
\mu,
C(h)
\right)
$$
Exponential covariance model:
$$
C(h)
=
\sigma^2
\exp
\left(
-\frac{h}{\lambda}
\right)
$$
- ln K(x) = log hydraulic conductivity at location x
- μ = mean
- C(h) = covariance function
- h = separation distance
Hydraulic conductivity variability is often modeled as a spatially correlated random field.
Linearized Predictive Covariance
$$
\Sigma_{pred}
=
J
\Sigma_p
J^T
+
\Sigma_{\epsilon}
$$
- Σ_pred = predictive covariance matrix
- J = Jacobian (sensitivity matrix)
- Σ_p = parameter covariance matrix
- Σ_ε = observation error covariance
Future plume predictions become more uncertain when the model is highly sensitive to uncertain parameters. Predictive confidence must account for both parameter variability and structural sensitivity.
Monte Carlo Predictive Ensemble
$$
C_j = f(\mathbf{p}_j)
$$
Predictive variance estimate:
$$
\sigma_C^2
=
\frac{1}{N-1}
\sum_{j=1}^{N}
\left(
C_j – \bar{C}
\right)^2
$$
$$C_{\text{upper}} = \bar{C} + 1.96\,\sigma_C$$
$$C_{\text{lower}} = \bar{C} – 1.96\,\sigma_C$$
- pⱼ = sampled parameter realization
- Cⱼ = model prediction for realization j
Defensible modeling presents ranges of outcomes rather than singular deterministic values. Decision-making should consider upper-bound scenarios, especially when remedial performance is at stake.
References
Anderson, M. P., Woessner, W. W., & Hunt, R. J. (2015). Applied groundwater modeling: Simulation of flow and advective transport (2nd ed.). Academic Press.
Carrera, J., & Neuman, S. P. (1986). Estimation of aquifer parameters under transient and steady-state conditions: 1. Maximum likelihood method incorporating prior information. Water Resources Research, 22(2), 199–210.
Doherty, J., 2025. Calibration and Uncertainty Analysis for Complex Environmental Models. Second Edition. Watermark Numerical Computing, Brisbane, Australia. ISBN:978-0-9943786-1-3.
Doherty, J., & Hunt, R. J. (2010). Approaches to highly parameterized inversion: A guide to using PEST for groundwater-model calibration. U.S. Geological Survey Scientific Investigations Report 2010–5169.
Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., & Rubin, D. B. (2013). Bayesian data analysis (3rd ed.). Chapman & Hall/CRC.
Rubin, Y. (2003). Applied stochastic hydrogeology. Oxford University Press.
Zheng, C., & Bennett, G. D. (2002). Applied contaminant transport modeling (2nd ed.). Wiley-Interscience.
Groundwater Modeling FAQs
Common groundwater modeling vulnerabilities include oversimplified hydraulic conductivity zoning, unsupported source assumptions, inadequate vertical discretization, and failure to quantify uncertainty. Independent groundwater modeling review evaluates these structural risks before predictive decisions are made.
Uncertainty quantification determines the sensitivity of model predictions to parameter variability. Without uncertainty analysis, groundwater modeling projections may produce overly confident remediation time frames or plume forecasts.
Bayesian groundwater modeling constitutes a probabilistic framework that integrates prior parameter uncertainty with observed field data to generate posterior parameter distributions and predictive uncertainty bounds. Instead of producing a singular “best-fit” model, Bayesian approaches quantify the spectrum of plausible subsurface conditions and the corresponding plume dynamics.
Groundwater flow modeling predicts hydraulic head and velocity fields, whereas contaminant transport modeling predicts concentration migration using advection-dispersion equations. Both must be integrated for sound environmental decision-making.
Parameter uncertainty reflects limited knowledge of hydraulic or transport values, while structural uncertainty arises from simplifications in conceptual model design. Structural uncertainty often dominates predictive risk.
Equifinality refers to the condition where multiple parameter combinations produce similar calibration results but different predictive outcomes. Several models may fit historical data equally well while predicting varying future plume behavior.
Multi-Gaussian geostatistical modeling represents hydraulic conductivity as a spatially correlated random field rather than uniform zones. It incorporates variance and correlation length to simulate realistic heterogeneity across the subsurface.
Groundwater models often support remedy selection, plume containment design, and long-term monitoring plans. Structural vulnerabilities or unquantified uncertainty may affect confidence and long-term liability exposure.
Calibration adjusts model parameters to match observed data, while validation evaluates predictive performance under independent conditions. A well-calibrated model does not inherently produce reliable long-term predictions.
Groundwater models often exhibit non-uniqueness, indicating that various parameter combinations may yield comparable calibration results. Variations in conceptual model structure can therefore lead to divergent long-term plume predictions.