DORiE: A Discontinuous Galerkin Solver for Soil Water Flow and Passive Solute Transport Based on DUNE

DORiE is a software framework for modelling and simulating transient water flow and passive solute transport in unsaturated porous media with discontinuous Galerkin (DG) methods using the Distributed and Unified Numerics Environment (DUNE). It supports computations with higher-order function spaces on structured and unstructured, non-conforming grids and features advanced numerical routines like flux reconstruction and adaptive grid refinement. DORiE wraps high numerical precision and an accurate physical process representation into an application tailored to soil physicists, which is controlled entirely via configuration files and a command line interface. The application is available as Docker Image from DockerHub1, and can be compiled from source.


Summary
DORiE is a software framework for modelling and simulating transient water flow and passive solute transport in unsaturated porous media with discontinuous Galerkin (DG) methods using the Distributed and Unified Numerics Environment (DUNE). It supports computations with higher-order function spaces on structured and unstructured, non-conforming grids and features advanced numerical routines like flux reconstruction and adaptive grid refinement. DORiE wraps high numerical precision and an accurate physical process representation into an application tailored to soil physicists, which is controlled entirely via configuration files and a command line interface. The application is available as Docker Image from DockerHub 1 , and can be compiled from source.

Background
Soil water flow and its transport of solutes are key ecosystem processes governing groundwater recharge, root water uptake, and freshwater buffering. The spatial scale of soil horizons and pedons, in which the Richards equation applies, is critical for these processes (Vogel, 2019). At this scale, the intricate soil architecture and general multi-scale heterogeneity features both smooth transitions and sharp boundaries. In conjunction with the highly non-linear hydraulic properties of the soil, simulating soil water flow therefore demands robust, flexible, and performant numerical solvers.
Discontinuous Galerkin finite element methods are high-order, locally conservative discretization schemes suitable for solving a wide range of partial differential equations. Through recent advances in numerics and the availability of suitable software and hardware, they have become feasible for problems with a wide range of length scales. Their formulations can be applied to various grid element geometries and non-conforming meshes, and typically outperform finite volume and continuous Galerkin methods in terms of numerical accuracy and efficiency (Di Pietro & Ern, 2012).
DG methods for linear advection-diffusion equations are well established (Ern et al., 2009). The feasibility of locally adaptive DG methods for solving the Richards equation has also been studied extensively, e.g. by Li et al. (2007). Nevertheless, DG methods are rarely used to solve subsurface flow problems because they are difficult to formulate and their increased efficiency is counteracted by a larger computational cost for the same mesh (Miller et al., 2013). Farthing & Ogden (2017) mention several successful applications of DG methods on two-phase flow problems and the Richards equation but feature no application with DG approximation scheme in their "list of representative research and production Richards' equation codes". DuMuX, as one particular framework for solving flow and transport equations in porous media with DUNE, only offers a Richards model in quasi two-phase flow description and variations of finite volume discretization schemes (Flemisch et al., 2011).

Feature Overview
DORiE is a DG solver for the Richards equation and the convection-dispersion equation in unsaturated porous media, and a module of the DUNE framework (Blatt et al., 2016). Primarily using DUNE-PDELab (Bastian et al., 2010), it implements cell-centered finite volume (CCFV) and symmetric weighted interior penalty discontinuous Galerkin (SWIP-DG) discretization schemes for both soil water flow and solute transport with various options for controlling upwinding, DG interior penalty, and DG penalty factors (Di Pietro & Ern, 2012;Ern et al., 2009). For solving equation systems, DORiE uses the algebraic multigrid solver optimized for DG methods introduced by Bastian (2014). On unstructured grids, it supports local adaptive grid refinement (h-refinement) based on an a-posteriori diffusive flux error estimator of the water flux (Ern & Stephansen, 2008). Since solving the transport problem requires a locally conservative water flux in the convection term, DORiE applies local (div) flux reconstruction using Raviart-Thomas finite elements as described by Ern et al. (2007) to the water flow solution and optionally to the solute transport solution. Additionally, the framework encapsulates the Gaussian random field generator of the DUNE-Randomfield module for creating pseudo-random soil architectures (Klein, 2016).
We especially focused on convenient modelling of soil architectures and boundary conditions via common file formats, which is intended to make these features available to scientists with little to no programming background. The simulation and random field routines are controlled via INI configuration files which reference additional input files. DORiE uses YAML files for specifying hydraulic parameters and possibly transient boundary conditions, and HDF5 files (The HDF Group, 1997) as input files for field data like initial conditions, grid mappings, and scaling fields. While rectangular grids are generated within DORiE, unstructured grids can be built with the open-source generator GMSH (Geuzaine & Remacle, 2009) and then loaded into the application. The simulation output is provided as VTK files (Schroeder et al., 2006) which can be evaluated by third-party software or with the Python module included in DORiE.

Research and Use Cases
DORiE has been coupled to a data assimilation framework for estimating hydraulic states, parameters, and soil heterogeneity using ensemble Kalman filters based on synthetic and real data. It supports a tight integration of both code bases and allows for its models being instantiated right inside the data assimilation framework. This enables the efficient adjustment of model states and parameters in place without the need to transfer data files or call external executables. A recent study uses the resulting software for investigating the feasibility and effectiveness of representing a two-dimensional, heterogeneous flow field as one-dimensional field with reduced heterogeneity (Bauser et al., 2020). Additionally, DORiE has been used for student exercises in several iterations of the "Physics of Terrestrial Systems" lecture at the Department of Physics and Astronomy, Heidelberg University, targeted at physics students in the master of science curriculum.

Showcases
We demonstrate the capabilities of DORiE in two exemplary showcases depicted in Figure 1 and  Simulation results for solute transport in heterogeneous soil consisting of one coarse-and one fine-grained material. The material architecture is created with the included random field generator. The lower boundary condition is a water table with constant solute concentration. At the upper boundary, we model evaporation and a no-flow boundary for solute. Top: Steady-state soil water content and water flux (direction and magnitude visualized by arrows) throughout the simulation. Bottom: Total solute concentration in the soil after some 58 d of simulation time. Note that the solute accumulates at the surface.

Figure 2:
Simulation results for infiltration into sandy soil with small-scale heterogeneity created with the included random field generator. The lower boundary condition is a water table and we model constant heavy rainfall at the upper boundary. The initial condition is hydraulic equilibrium. The grid is visualized in grey. Top: Hydraulic conductivity for a completely water-saturated domain. The grid is shown in its coarsest state with which the simulation is started. Bottom: Soil water content, water flux (direction and magnitude visualized by arrows), and refined grid after 4 h of simulation time.