Changes in permafrost soils by freezing or melting processes pose a challenge in predicting groundwater flow for extended time periods. A further important aspect is the density-driven flow arising for example due to the dilution of the saline groundwater by the melt water from the permafrost. This work extends an existing freeze-thaw model in porous media to include saline soils and implements the fully coupled system of partial differential equations using efficient numerical methods. These include a linearly implicit time stepping scheme and an efficient linear solver based on the geometric multigrid preconditioner. For numerical experiments, we take an extension of the InterFrost benchmark. Our results highlight the advantage of adaptive time stepping and show that dependence of the freezing point on salt mass fraction leads to a faster melting process. They also demonstrate how meltwater can impact formation of salt fingers by the density-driven flow.

Salinization of coastal aquifers is a current problem in many regions worldwide. Simulation of this process provides an important tool for the forecast of drinking water resources. The consideration of the unsaturated phreatic zone has an essential influence on the accuracy of the prediction. However, for the large temporal and spatial scales of the aquifers, the presence of the unsaturated subdomains creates challenging difficulties for the numerical methods. In this work, we investigate two approaches for simulating haline density-driven flow in partially saturated aquifers. The first approach, based on the Richards equation, uses a representation of the saturation field and employs an adaptive linearly implicit time discretization scheme. The second approach explicitly represents the water table using a level-set method and relies on weak coupling combined with a standard implicit Euler scheme. Both approaches employ a finite-volume discretization for spatial representation of fluid flow and salt transport. The implementation is based on the UG4 toolkit together with the parallel groundwater flow simulator d3f++. The resulting linear systems are solved using a geometric multigrid method. We compare the aforementioned approaches with respect to theoretical and numerical aspects. The performance of both methods is evaluated in three numerical experiments. In particular we demonstrate robustness of the linearly implicit scheme and introduce a stabilization for the level-set method. The high-performance computing (HPC) potential of the approaches is assessed and demonstrated as well.

We consider a Henry-like intrusion problem, where fluid flow is driven by variations in fluid density. The Multi-Level Monte Carlo (MLMC) method is employed to estimate the mean value of a quantity of interest (QoI). The QoI is defined as the earliest time at which the mass fraction of salt exceeds a given threshold. In our setting, porosity, permeability, recharge, and fracture thickness are treated as uncertain parameters and modeled as random variables. For each realization of these parameters, the evolution of the salt mass fraction is governed by a system of nonlinear, time-dependent partial differential equations (PDEs). We demonstrate that the MLMC method can be effectively applied to this problem, significantly reducing computational costs compared to classical Monte Carlo methods. The findings of this study have the potential to enhance and accelerate the monitoring of drinking water resources and pollution dynamics.

Performance assessment and field-development planning for carbon geo-sequestration relies on reservoir simulation. The standard approach is to use regular (corner-point) grids in space and a fully implicit discretisation in time, solving for all primary unknowns at once using Newton’s method. While the former limits the physical realism of the simulation model built from the geomodel, the latter leads to a large ill-conditioned system of equations that is inefficient to solve. We overcome the regularisation issue via a hierarchy of fully unstructured, geobody-conforming, finite element meshes which can be refined until mesh convergence is achieved. To address the time discretisation issue we employ, for the first time, the linearly implicit extrapolation scheme (LIMEX) to solve the highly non-linear coupled two-phase flow and reactive transport equations. To solve the arising large sparse systems of linear equations, we apply the geometric multigrid (GMG) method that demonstrates optimal, linear complexity and allows an efficient parallelization on supercomputers. Another novel feature of our formulation is the consideration of the kinetics of mass transfer between the carbonic and aqueous phases. This approach removes the first-order dependence on mesh refinement of CO_2 that is dissolved, a characteristic feature of standard equilibrium models. We demonstrate the parallel scalability of our simulation framework with an implementation based on the UG4 platform. Proof-of-concept results accurately capture key features of CO_2 migration including filtration by capillary barriers and convective dissolution of CO_2 at the base of the plume.

We investigate the applicability of the well-known multilevel Monte Carlo (MLMC) method to the class of density-driven flow problems, in particular the problem of salinisation of coastal aquifers. As a test case, we solve the uncertain Henry saltwater intrusion problem. Unknown porosity, permeability and recharge parameters are modelled by using random fields. The classical deterministic Henry problem is non-linear and time-dependent, and can easily take several hours of computing time. Uncertain settings require the solution of multiple realisations of the deterministic problem, and the total computational cost increases drastically. Instead of computing of hundreds random realisations, typically the mean value and the variance are computed. The standard methods such as the Monte Carlo or surrogate-based methods are a good choice, but they compute all stochastic realisations on the same, often, very fine mesh. They also do not balance the stochastic and discretisation errors. These facts motivated us to apply the MLMC method. We demonstrate that by solving the Henry problem on multi-level spatial and temporal meshes, the MLMC method reduces the overall computational and storage costs. To reduce the computing cost further, parallelization is performed in both physical and stochastic spaces. To solve each deterministic scenario, we run the parallel multigrid solver ug4 in a black-box fashion.

We use the Multi Level Monte Carlo method to estimate uncertainties in a Henry-like salt water intrusion problem with a fracture. The flow is induced by the variation of the density of the fluid phase, which depends on the mass fraction of salt. While the fracture’s location is fixed, its aperture is uncertain. In our setting, porosity and permeability vary spatially and recharge is time-dependent. So we introduce three random variables, one controlling both the porosity and permeability fields, one for the fracture width and one for the intensity of recharge. For each realization of these uncertain parameters, the evolution of mass fraction and pressure fields is modeled using a system of non-linear, time-dependent PDEs with a solution discontinuity at the fracture. These uncertainties propagate, affecting the distribution of salt concentration, a key factor in water resource quality. We show that the MLMC method can be successfully applied to this problem. It significantly reduces the computational cost compared to classical Monte Carlo methods by effectively balancing discretisation and statistical errors, and by evaluating multiple scenarios over different spatial and temporal mesh levels. The deterministic PDE solver, using the ug4 library, runs in parallel to compute all stochastic scenarios.

We consider a class of density-driven flow problems. We are particularly interested in the problem of the salinization of coastal aquifers. We consider the Henry saltwater intrusion problem with uncertain porosity, permeability, and recharge parameters as a test case. The reason for the presence of uncertainties is the lack of knowledge, inaccurate measurements, and inability to measure parameters at each spatial or time location. This problem is nonlinear and time-dependent. The solution is the salt mass fraction, which is uncertain and changes in time. Uncertainties in porosity, permeability, recharge, and mass fraction are modeled using random fields. This work investigates the applicability of the well-known multilevel Monte Carlo (MLMC) method for such problems. The MLMC method can reduce the total computational and storage costs. Moreover, the MLMC method runs multiple scenarios on different spatial and time meshes and then estimates the mean value of the mass fraction. The parallelization is performed in both the physical space and stochastic space. To solve every deterministic scenario, we run the parallel multigrid solver ug4 in a black-box fashion. We use the solution obtained from the quasi-Monte Carlo method as a reference solution.

We are solving a problem of salinisation of coastal aquifers. As a test case example, we consider the Henry saltwater intrusion problem. Since porosity, permeability and recharge are unknown or only known at a few points, we model them using random fields and random variables. The Henry problem describes a two‐phase flow and is non‐linear and time‐dependent. The solution to be found is the expectation of the salt mass fraction, which is uncertain and time‐dependent. To estimate this expectation, we use the well‐known multilevel Monte Carlo (MLMC) method. The MLMC method takes just a few samples on computationally expensive (fine) meshes and more samples on cheap (coarse) meshes. Then, by building a telescoping sum, the MLMC method estimates the expected value at a much lower computational cost than the classical Monte Carlo method. The deterministic solver used here is the well‐known parallel and scalable ug4 solver.

We are solving a problem of salinization of coastal aquifers.As a test case example, we consider the Henry saltwater intrusion problem.Since porosity, permeability and recharge are unknown or only known at a few points, we model them using random fields.The Henry problem describes a two-phase flow and is nonlinear and time-dependent.The solution to be found is the expectation of the salt mass fraction, which is uncertain and time-dependent.To estimate this expectation we use the well known multilevel Monte Carlo (MLMC) method.The MLMC method takes just a few samples on computationaly expensive (fine) meshes and more samples on cheap (coarse) meshes.Then, by building a telescoping sum, the MLMC method estimates the expected value at a much lower cost than the classical Monte Carlo method.The deterministic solver used here is the well-known parallel and scalable ug4 solver.

Accurate modeling of contamination in subsurface flow and water aquifers is crucial for agriculture and environmental protection. Here, we demonstrate a parallel algorithm to quantify the propagation of uncertainty in the dispersal of pollution in subsurface flow. Specifically, we consider the density-driven flow and estimate how uncertainty from permeability and porosity propagates to the solution. We take a two-dimensional Elder-like problem as a numerical benchmark, and we use random fields to model our limited knowledge on the porosity and permeability. We use the well-known low-cost generalized polynomial chaos (gPC) expansion surrogate model, where the gPC coefficients are computed by projection on sparse tensor grids. The numerical solver for the deterministic problem is based on the multigrid method and is run in parallel. Computation of high-dimensional integrals over the parametric space is done in parallel too.

The pollution of groundwater, essential for supporting populations and agriculture, can have catastrophic consequences. Thus, accurate modeling of water pollution at the surface and in groundwater aquifers is vital. Here, we consider a density-driven groundwater flow problem with uncertain porosity and permeability. Addressing this problem is relevant for geothermal reservoir simulations, natural saline-disposal basins, modeling of contaminant plumes and subsurface flow predictions. This strongly nonlinear time-dependent problem describes the convection of a two-phase flow, whereby a liquid flows and propagates into groundwater reservoirs under the force of gravity to form so-called “fingers”. To achieve an accurate numerical solution, fine spatial resolution with an unstructured mesh and, therefore, high computational resources are required. Here we run a parallelized simulation toolbox ug4 with a geometric multigrid solver on a parallel cluster, and the parallelization is carried out in physical and stochastic spaces. Additionally, we demonstrate how the ug4 toolbox can be run in a black-box fashion for testing different scenarios in the density-driven flow. As a benchmark, we solve the Elder-like problem in a 3D domain. For approximations in the stochastic space, we use the generalized polynomial chaos expansion. We compute the mean, variance, and exceedance probabilities for the mass fraction. We use the solution obtained from the quasi-Monte Carlo method as a reference solution.

The generation of detailed three dimensional meshes for the simulation of groundwater flow in thin layered domains is crucial to capture important properties of the underlying domains and to reach a satisfying accuracy. At the same time, this level of detail poses high demands both on suitable hardware and numerical solver efficiency. Parallel multigrid methods have been shown to exhibit near optimal weak scalability for massively parallel computations of density driven flow. A fully automated parameterized algorithm for prism based meshing of coarse grids from height data of individual layers is presented. Special structures like pinch outs of individual layers are preserved. The resulting grid is used as a starting point for parallel mesh and hierarchy creation through interweaved projected refinement and redistribution. Efficiency and applicability of the proposed approach are demonstrated for a parallel multigrid based simulation of a realistic sample problem.

This work summarizes solution strategies for discrete systems occurring in the simulation of processes in the subsurface. The focus is on scalable solvers for large and coupled systems. The goal of this work is to enable researchers to select suitable algorithms and parameter settings to efficiently solve their problems. The work provides an overview of existing methods, highlighting their features, potential, and also frequent pitfalls. Numerical examples are provided for single phase flow, density driven flow and poroelasticity.Aspects of multiphase flow are discussed briefly; a detailed discussion of reactive transport is beyond the scope of the article. Future trends are discussed.

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.

A 3d regional density-driven flow model of a heterogeneous aquifer system at the German North Sea Coast is set up within the joint project NAWAK (“Development of sustainable adaption strategies for the water supply and distribution infrastructure on condition of climatic and demographic change”). The development of the freshwater-saltwater interface is simulated for three climate and demographic scenarios.Groundwater flow simulations are performed with the finite volume code d 3 f++ (distributed density driven flow) that has been developed with a view to the modelling of large, complex, strongly density-influenced aquifer systems over long time periods.