Efficient Simulations for Contamination of Groundwater Aquifers under Uncertainties

Accurate modeling of contamination in subsurface flow and water aquifers is crucial for agriculture and environmental protection. Here, we demonstrate a parallel method to quantify the propagation of the uncertainty in the dispersal of pollution in density‐driven flow. We solve an Elder‐like problem, where we use random fields to model the limited knowledge on the porosity and permeability. The uncertain solution, mass fraction, is approximated via low‐cost generalized polynomial chaos expansion (gPCE). Parallelization is done in both the physical and parametric spaces.

1 Density-driven groundwater flow problem with uncertain porosity and permeability Numerical simulation of the density-driven flow plays an important role in modeling how pollution can migrate through an aquifer [1][2][3][4]. Some pollutants, like salt, are soluble in fresh water and can leak into groundwater systems. Such pollutants can change the density of a fluid and induce density-driven flows within the aquifer. This results in a faster propagation of the contamination due to convection.
We applied an established gPCE technique [5][6][7][8] to estimated the propagation of input uncertainties associated with the porosity and permeability into the QoI, which could be a) the mean and the variance of the mass fraction in a sub-domain at each time step, b) the exceedance probability or quantiles in selected points, or c) probability density function in a point. These statistics could be further used for more efficient Bayesian inference, data assimilation, optimal design of the experiment and optimal control. Our overall goal is to model the dispersion of water pollution and monitor the movement of subsurface water flow. In addition, we used the novel technique of parallelization in spatial and stochastic spaces. We solved each deterministic problem and all stochastic realizations in parallel.
We consider a domain D ⊂ R d , d = 2, filled with two phases: a solid phase and a liquid phase. The solid phase is an immobile solid porous matrix with a mobile liquid in its pores. The matrix is characterized by its porosity φ : D → R and its permeability K : D → R d×d [9]. The liquid phase is a solution of salt (sodium chloride) in water. We denote the mass fraction of the brine (a saturated saltwater solution) in the liquid phase by c(t, x) : [0, +∞) × D → [0, 1]. The liquid phase has density ρ = ρ(c) and viscosity µ = µ(c).
The motion of the liquid phase through the solid matrix is modeled by the Darcy velocity q(t, x) : [0, +∞) × D → R d . In these terms, the conservation laws for the liquid phase and the dissolved salt can be written as where the tensor field D represents the molecular diffusion and the mechanical dispersion of the salt. In addition, some momentum equation must be considered for q. We assume the Darcy's law denotes the hydrostatic pressure and g = (0, . . . , −g) T ∈ R d is the gravity vector.
For simplicity, we assume an isotropic porous medium characterized by a scalar permeability We assume the linear dependence for the density: ρ(c) = ρ 0 + (ρ 1 − ρ 0 )c, where ρ 0 and ρ 1 denote the densities of pure water and the brine, respectively. Thus, c ∈ [0, 1], where c = 0 corresponds to the pure water and c = 1 to the saturated solution.
We assume also that D = φD m I, where D m is the coefficient of the molecular diffusion. The viscosity is considered to be constant. We model porosity φ as a random field, namely, φ = φ(x, θ), x ∈ D, θ = (θ 1 , ..., θ M ), and θ i are random variables. We assume K to be isotropic and dependent on φ: K = KI, K = K(φ) ∈ R.
We use a Kozeny-Carman-like dependence: K(φ) = κ KC · φ 3 1 − φ 2 , where κ KC is a scaling factor. To solve the deterministic problem (stochastic porosity is fixed) we used plugin d 3 f of the simulation framework ug4, [10,11]. This framework presents a flexible tool for the numerical solution of non-stationary and nonlinear systems of PDEs, and is parallelized and highly scalable. For the solution, Equations 1-2 are discretized by a vertex-centered finite-volume method on unstructured grids in the geometric space, as presented in [12]. The implicit Euler scheme is used for the time discretization. The obtained large sparse system of nonlinear algebraic equations in every time step is solved by the Newton method. The sparse linear system, which appears in the nonlinear iterations, is solved by the BiCGStab method with the multigrid preconditioning (V-cycle) and the ILU β -smoothers.
Numerics. We consider a 2D reservoir geometry D = (0, 600) × (0, 150) m 2 . We set the Dirichlet boundary conditions for c on the top and bottom boundaries. On the left and right sides, we impose the no-flux boundary conditions for equation (2). For equation (1), we impose the no-flux boundary conditions on the whole ∂D. We set p to 0 at the two upper corners of the domain. We assume the porosity is defined as follow φ(x, θ) = 0.09 + 0.005 · θ · (cos(πx/300) + sin(πy/150)), 1]. We computed the mean values and variances by the quasi Monte Carlo (qMC) and gPCE based surrogate methods [7,8,13,14]. The number of qMC simulations (Halton sequence) is 200, and for the gPCE of order p = 4 we used 9 Gauss-Legendre quadrature points. In Fig. 1(a) more advanced sparse grid rule may be needed to approximate the gPCE coefficients of the solution. Also, the standard gPCE may fail for more complicated porosity fields.