2.1 Variably saturated flow

ParFlow can operate in variably saturated mode using the well-known mixed form of the Richards equation (Celia et al., 1990). The mixed form of the Richards equation implemented in ParFlow is

where S_s is the specific storage coefficient [L⁻¹], S_w is the relative saturation [-] as a function of pressure head p of the fluid or water [L], t is time [T], ϕ is the porosity of the medium [-], q is the specific volumetric (Darcy) flux [LT⁻¹], k_s is the saturated hydraulic conductivity tensor [LT⁻¹], k_r is the relative permeability [-], which is a function of pressure head, q_s is the general source or sink term [T⁻¹] (includes wells and surface fluxes, e.g., evaporation and transpiration), and z is depth below the surface [L]. The Richards equation assumes that the air phase is infinitely mobile (Richards, 1931). ParFlow has been used to numerically simulate river–aquifer exchange (free-surface flow and subsurface flow; Frei et al., 2009) and highly heterogenous problems under variably saturated flow conditions (Woodward, 1998; Jones and Woodward, 2001; Kollet et al., 2010). Under saturated conditions, e.g., simulating linear groundwater movement under assumed predevelopment conditions, the steady-state saturated mode can be used.

2.2 Steady-state saturated flow

The most basic operational mode is the solution of the steady-state, fully saturated groundwater flow equation:

where q_s represents a general source or sink term, e.g., wells [T⁻¹], q is the Darcy flux [LT⁻¹], which is usually written as

where k_s is the saturated hydraulic conductivity [LT⁻¹], and P represents the 3-D hydraulic head-potential [L]. ParFlow does include a direct solution option for the steady-state saturated flow that is distinct from the transient solver. For example, ParFlow uses the solver “impes” under the single-phase, fully saturated, steady-state condition relative to the variably saturated, transient mode wherein the Richards equation solver is used (Maxwell et al., 2016). When studying sophisticated or complex phenomena, e.g., simulating a fully coupled system (i.e., surface and subsurface flow), an overland flow boundary condition is employed.

2.3 Overland flow

Surface water systems are connected to the subsurface, and these interactions are particularly important for rivers. However, these connections have been historically difficult to explicitly represent in numerical simulations. A common approach has been to use river-routing codes, like Hydrologic Engineering Center (HEC) codes, as well as MODFLOW and its River Package to determine head in the river, which is then used as a boundary condition for the subsurface model. This approach prevents feedbacks between the two models, and a better representation of the physical processes in these kinds of problems is one of the motivations for IHMs. Overland flow is implemented in ParFlow as a two-dimensional kinematic wave equation approximation of the shallow water equations. The continuity equation for two-dimensional shallow overland flow is given as

where υ is the depth-averaged velocity vector [LT⁻¹], ψ_s is the surface ponding depth [L], t is time [T], and q_s is a general source or sink (e.g., precipitation rate) [T⁻¹]. Ignoring the dynamic and diffusion terms results in the momentum equation,

which is known as the kinematic wave approximation. The S_f,i and S_o,i represent the friction [-] and bed slopes (gravity forcing term) [-], respectively, where i indicates the x and y directions (also shown in Eqs. 7 and 8) (Maxwell et al., 2015). Manning's equation is used to generate a flow depth–discharge relationship:

where n is the Manning roughness coefficient $\left[{\mathrm{TL}}^{-\mathrm{1}/\mathrm{3}}\right]$. The flow of water out of an overland flow simulation domain only occurs horizontally at an outlet controlled by specifying a type of boundary condition at the edge of the simulation domain. In a natural system, the outlet is usually taken as the region where a river enters another water body such as a stream or lake. ParFlow determines the overland flow direction through the D4 flow-routing approach. In a simulation domain, the D4 flow-routing approach allows for flow to be assigned from a focal cell to only one neighboring cell accessed via the steepest or most vertical slope. The shallow overland flow formulation (Eq. 9) assumes that the flow depth is averaged vertically and neglects a vertical change in momentum in the column of surface water. To account for vertical flow (from the surface to the subsurface or subsurface to the surface), a formulation that couples the system of equations through a boundary condition at the land surface becomes necessary. Equation (5) can be modified to include an exchange rate with the subsurface, q_e, as

which is common in other IHMs. In ParFlow, the overland flow equations are coupled directly to the Richards equation at the top boundary cell under saturated conditions. Conditions of continuity of pressure (i.e., the pressures of the subsurface and surface domains are equal right at the ground surface) and flux at the top cell of the boundary between the subsurface and surface systems are assigned. Fig. 1 demonstrates the continuity of pressure at the ground surface for flow from the surface into the subsurface. This assignment is done by setting pressure head in Eq. (1) equal to the vertically averaged surface pressure, ψ_s,

and the flux, q_e, equal to the specified boundary conditions (e.g., Neumann or Dirichlet type). For example, if Neumann-type boundary conditions are specified, which are given as

and one solves for the flux term in Eq. (10), the result is

where the ∥ψ,0∥ operator is defined as the greater of the quantities, ψ and 0. Substituting Eq. (12) for the boundary condition in Eq. (11), requiring the aforementioned flux continuity q_BC=q_e, leads to

Equation (13) shows that the surface water equations are represented as a boundary condition to the Richards equation. That is, the boundary condition links flow processes in the subsurface with those at the land surface. This boundary condition eliminates the exchange flux and accounts for the movement of the free surface of ponded water at the land surface (Kollet and Maxwell, 2006; Williams and Maxwell, 2011).

Figure 1Coupled surface and subsurface flow systems. The physical system is represented in (a), and a schematic of the overland flow boundary condition (continuity of pressure and flux at the ground surface) is in (b). The equation $p={\mathit{\psi }}_{\mathrm{s}}=\mathit{\psi }$ in Fig. 1 signifies that at the ground surface, the vertically averaged surface pressure and subsurface pressure head are equal, which is the unique overland flow boundary used by ParFlow.

Download

Many IHMs couple subsurface and surface flows by making use of the exchange flux, q_e, model. The exchange flux between the domains (the surface and the subsurface) depends on hydraulic conductivity and the gradient across some interface where indirect coupling is used (VanderKwaak, 1999; Panday and Huyakorn, 2004). The exchange flux concept gives a general formulation of a single set of coupled surface–subsurface equations. The exchange flux term, q_e, may be included in the shallow overland flow continuity equation as the exchange rate term with the subsurface (Eq. 9) in a coupled system (Kollet and Maxwell, 2006).

2.4 Multiphase flow and transport equations

Most applications of the code have reflected ParFlow's core functionality as a single-phase flow solver, but there are also embedded capabilities for the multiphase flow of immiscible fluids and solute transport. Multiphase systems are distinguished from single-phase systems by the presence of one or more interfaces separating the phases, with moving boundaries between phases. The flow equations that are solved in multiphase systems in a porous medium comprise a set of mass balance and momentum equations. The equations are given by

where $i=\mathrm{1},\mathrm{\dots },n$ denotes a given phase (such as air or water). In these equations, ϕ is the porosity of the medium [-], which explains the fluid capacity of the porous medium, and for each phase, i S_i(xt) is the relative saturation [-], which indicates the content of phase i in the porous medium, υ_i(xt) represents the Darcy velocity vector [LT⁻¹], Q_i(xt) stands for the source or sink term [T⁻¹], p_i(xt) is the average pressure $\left[{\mathrm{ML}}^{-\mathrm{1}}{\mathrm{T}}^{-\mathrm{2}}\right]$, ρ_i(xt) is the mass density [ML⁻³], λ_i is the mobility [L³TM⁻¹], g is the gravity vector [LT⁻²], and x and t represent the space vector and time, respectively. ParFlow solves for the pressures on a discrete mesh and uses a time-stepping algorithm based on a mass conservative backward Euler scheme and spatial discretization (a finite-volume method). ParFlow's multiphase flow capability has not been applied in major studies; however, this capability is also available for testing (Ashby et al., 1993; Tompson et al., 1994; Falgout et al., 1999; Maxwell et al., 2016).

The transport equations included in the ParFlow package describe mass conservation in a convective flow (no diffusion) with degradation effects and adsorption included along with extraction and injection wells (Beisman et al., 2015; Maxwell et al., 2016). The transport equation is defined as follows:

where c_i,j(xt) represents the concentration fraction of contaminant [-], λ_i is degradation rate [T⁻¹], F_i(xt) is the mass concentration [L³M⁻¹], ρ_s(x) is the density of the solid mass [ML⁻³], n_I is injection wells [-], ${\mathit{\gamma }}_k^{\mathrm{I};i}\left(t\right)$ is the injection rate [T⁻¹], ${\mathrm{\Omega }}_k^{\mathrm{I}}\left(\mathit{x}\right)$ represents the area of the injection well [-], $c_{i,j}^{-k}\left(\mathit{x}t\right)$ is the injected concentration fraction [-], n_E is the extraction wells [-], ${\mathit{\gamma }}_k^{\mathrm{E};i}\left(t\right)$ is the extraction rate [T⁻¹], ${\mathrm{\Omega }}_k^{\mathrm{E}}\left(\mathit{x}\right)$ is an extraction well area [-], $i=\mathrm{0},\mathrm{\dots },n_{\mathrm{p}-\mathrm{1}}\left(n_{\mathrm{p}}\mathit{ϵ}\left\{\mathrm{1},\mathrm{2},\mathrm{3}\right\}\right)$ is the number of phases, $j=\mathrm{0},\mathrm{\dots },n_{\mathrm{c}}-\mathrm{1}$ represents the number of contaminants, c_i,j is the concentration of contaminant j in phase i, k is hydraulic conductivity [LT⁻¹], $\mathit{\chi }{\mathrm{\Omega }}_k^{\mathrm{I}}$ is the characteristic function of an injection well region, and $\mathit{\chi }{\mathrm{\Omega }}_k^{\mathrm{E}}$ is the characteristic function of an extraction well region. The mass concentration term, F_i,j, is taken to be instantaneous in time and a linear function of contaminant concentration:

where K_d;j is the distribution coefficient of the component [L³M⁻¹]. The transport–advection equation or convective flow calculation performed by ParFlow offers a choice of a first-order explicit upwind scheme or a second-order explicit Godunov scheme. The advection calculations are discretized as boundary value problems for each primary dimension over each compute cell. The discretization is a fully explicit, forward Euler, first-order accurate in time approach. The implementation of a second-order explicit Godunov scheme (second-order advection scheme) minimizes numerical dispersion and presents a more accurate computational process at these timescales than either an implicit or lower-order explicit scheme. The stability issue here is that the simulation time step is restricted via the Courant–Friedrichs–Lewy (CFL) condition, which demands that time steps are chosen as small enough to ensure that mass is not transported more than one grid cell in a single time step in order to maintain stability (Beisman, 2007).

2.5 Computational grids

An accurate numerical approximation of a set of partial differential equations is strongly dependent on the simulation grid. Integrated hydrologic models can use unstructured or structured meshes for the discretization of the governing equations. The choice of grid type to adopt is problem-specific and often a subjective choice since the same domain can be represented in many ways, but there are some clear trade-offs. For example, structured grid models, such as ParFlow, may be preferred to unstructured grid models because structured grids provide significant advantages in computational simplicity and speed, and they are amenable to efficient parallelization (Durbin, 2002; Kumar et al., 2009; Osei-Kuffuor et al., 2014). ParFlow adopts a regular, structured grid specifically for its parallel performance. There are currently two regular grid formulations included in ParFlow, an orthogonal grid and a terrain-following formulation (TFG); both allow for variable vertical discretization (thickness over an entire layer) over the domain.

2.5.1 Orthogonal grid

Orthogonal grids have many advantages, and many approaches are available to transform an irregular grid into an orthogonal grid such as conformal mapping. This mapping defines a transformed set of partial differential equations using an elliptical system with “control functions” determined in such a way that the generated grid would be either orthogonal or nearly orthogonal. However, conformal mapping may not allow flexibility in the control of the grid node distribution, which diminishes its usefulness with complex geometries (Mobley and Stewart, 1980; Haussling and Coleman, 1981; Visbal and Knight, 1982; Ryskin and Leal, 1983; Allievi and Calisal, 1992; Eca, 1996).

A Cartesian, regular, orthogonal grid formulation is implemented by default in ParFlow, though some adaptive meshing capabilities are still included in the source code. For example, layers within a simulation domain can be made to have varying thickness. Figure 2a shows the standard way that topography or any other non-rectangular domain boundaries are represented in ParFlow. The domain limits, and any other internal boundaries, can be defined using grid-independent triangulated irregular network (TIN) files that define a geometry, or a gridded indicator file can be used to define geometric elements. ParFlow uses an octree space-partitioning algorithm (a grid-based algorithm or mesh generators filled with structured grids) (Maxwell, 2013) to depict complex structure and land surface representations (e.g., topography, watershed boundaries, and different hydrologic facies) in three-dimensional space (Kollet et al., 2010). These land surface features are mapped onto the orthogonal grid, and looping structures that encompass these irregular shapes are constructed (Ashby et al., 1997). The grid cells above the ground surface are inactive (shown in Fig. 2a) and are stored in the solution vector but not included in the solution.

Figure 2Representation of orthogonal (a) and terrain-following (b) grid formulations and schematics of the associated finite-difference dependences (right). The i, j, and k are the x, y, and z cell indices.

Download

3.1 Newton–Krylov solver for variably saturated flow

The cell-centered fully implicit discretization scheme applied to the Richards equation leads to a set of coupled discrete nonlinear equations that need to be solved at each time step, and, for variably saturated subsurface flow, ParFlow does this with the inexact Newton–Krylov method implemented in the KINSOL package (Hindmarsh et al., 2005; Collier et al., 2015). Newton–Krylov methods were initially utilized in the context of partial differential equations by Brown and Saad (1990). In the approach, a coupled nonlinear system as a result of discretization of the partial differential equation is solved iteratively. Within each iteration, the nonlinear system is linearized via a Taylor expansion. After linearization, an iterative Krylov method is used to solve the resulting linear Jacobian system (Woodward, 1998; Osei-Kuffuor et al., 2014). For variably saturated subsurface flow, ParFlow uses the GMRES Krylov method (Saad and Schultz, 1986). Figure 3 is a flowchart of the solution technique ParFlow uses to provide approximate solutions to systems of nonlinear equations.

Figure 3Working flowchart of ParFlow's solver for linear and nonlinear system solutions.

Download

The benefit of this Newton–Krylov method is that the Krylov linear solver requires only matrix–vector products. Because the system matrix is the Jacobian of the nonlinear function, these matrix–vector products may be approximated by taking directional derivatives of the nonlinear function in the direction of the vector to be multiplied. This approximation is the main advantage of the Newton–Krylov approach as it removes the requirement for matrix entries in the linear solver. An inexact Newton method is derived from a Newton method by using an approximate linear solver at each nonlinear iteration, as is done in the Newton–Krylov method (Dembo and Eisenstat, 1982; Dennis Jr. and Schabel, 1996). This approach takes advantage of the fact that when the nonlinear system is far from converged, the linear model used to update the solution is a poor approximation. Thus, the convergence criteria for an early linear system solver are relaxed. The tolerance required for the solution of the linear system is decreased as the nonlinear function residuals approach zero. The convergence rate of the resulting nonlinear solver can be linear or quadratic, depending on the algorithm used. Through the KINSOL package, ParFlow can either use a constant tolerance factor or ones from Eisenstat and Walker (1996). Krylov methods can be very robust, but they can be slow to converge. As a result, it is often necessary to implement a preconditioner, or accelerator, for these solvers.

3.2 Multigrid solver

Multigrid (MG) methods constitute a class of techniques or algorithms for solving differential equations (system of equations) using a hierarchy of discretization (Briggs et al., 2000). Multigrid algorithms are applied primarily to solve linear and nonlinear boundary value problems and can be used as either preconditioners or solvers. The most efficient method for preconditioning linear systems in ParFlow is the ParFlow Multigrid (PFMG) algorithm (Ashby and Falgout, 1996; Jones and Woodward, 2001). Multigrid algorithms arise from the discretization of elliptic partial differential equations (Briggs et al., 2000) and, in ideal cases, have convergence rates that do not depend on the problem size. In these cases, the number of iterations remains constant even as problem sizes grow large. Thus, the algorithm is algorithmically scalable. However, it may take longer to evaluate each iteration as problem sizes increase. As a result, ParFlow utilizes the highly efficient implementation of the PFMG in the hypre library (Falgout and Yang, 2002).

For variably saturated subsurface flow, ParFlow uses the Newton–Krylov method coupled with a multigrid preconditioner to accurately solve for the water pressure (hydraulic head) in the subsurface and diagnoses the saturation field (which is used in determining the water table) (Woodward, 1998; Jones and Woodward, 2000, 2001; Kollet et al., 2010). The water table is calculated for computational cells having hydraulic heads above the bottom of the cells. Generally, a cell is saturated if the hydraulic head in the cell is above the node elevation (cell center), or the cell is unsaturated if the hydraulic head in the cell is below the node elevation. For saturated flow, ParFlow uses the conjugate gradient method also coupled with a multigrid method. It is important to note that subsurface flow systems are usually much larger radially than they are thick, so it is common for computational grids to have highly anisotropic cell aspect ratios to balance the lateral and vertical discretization. Combined with anisotropy in the permeability field, these high aspect ratios produce numerical anisotropy in the problem, which can cause the multigrid algorithms to converge slowly (Jones and Woodward, 2001). To correct this problem, a semi-coarsening strategy or algorithm is employed, whereby the grid is coarsened in one direction at a time. The direction chosen is the one with the smallest grid spacing, i.e., the tightest coupling. In an instance in which more than one direction has the same minimum spacing, then the algorithm chooses the direction in the order of x, followed by y, and then in z. To decide on how and when to terminate the coarsening algorithm, Ashby and Falgout (1996) determined that a semi-coarsening down to a $(\mathrm{1}\times \mathrm{1}\times \mathrm{1})$ grid is ideal for groundwater problems.

3.3 Multigrid-preconditioned conjugate gradient (MGCG)

ParFlow uses the multigrid-preconditioned conjugate gradient (CG) solver to solve the groundwater equations under steady-state and fully saturated flow conditions (Ashby and Falgout, 1996). These problems are symmetric and positive definite, two properties the CG method was designed to target. While CG lends itself to efficient implementations, the number of iterations required to solve a system that results from the discretization of the saturated flow equation increases as the problem size grows. The PFMG algorithm is used as a preconditioner to combat this growth and results in an algorithm for which the number of iterations required to solve the system grows only minimally. See Ashby and Falgout (1996) for a detailed description of these solvers and the parallel implementation of the multigrid-preconditioned CG method in ParFlow (Gasper et al., 2014; Osei-Kuffuor et al., 2014).

3.4 Preconditioned Newton–Krylov for coupled subsurface–surface flows

As discussed above, coupling between subsurface and surface or overland flow in ParFlow is activated by specifying an overland boundary condition at the top surface of the computational domain, but this mode of coupling allows for activation and deactivation of the overland boundary condition during simulations in which ponding or drying occurs. Thus, surface–subsurface coupling can occur anywhere in the domain during a simulation, and it can change dynamically during the simulation. Overland flow may occur by the Dunne or Horton mechanism depending on local dynamics. Overland flow routing is enabled when the subsurface cells are fully saturated. In ParFlow the coupling between the subsurface and surface flows is handled implicitly. ParFlow solves this implicit system with the inexact Newton–Krylov method described above. However, in this case, the preconditioning matrix is adjusted to include terms from the surface coupling. In the standard saturated or variably saturated case, the multigrid method is given the linear system matrix, or a symmetric version, resulting from the discretization of the subsurface model. Because ParFlow uses a structured mesh, these matrices have a defined structure, making their evaluation and the application of a multigrid straightforward. Due to varying topographic height of the surface boundary, where the surface coupling is enforced, the surface effects add nonstructured entries in the linear system matrices. These entries increase the complexity of the matrix entry evaluations and reduce the effectiveness of the multigrid preconditioner. In this case, the matrix–vector products are most effectively performed through the computation of the linear system entries rather than the finite-difference approximation to the directional derivative. For the preconditioning, surface couplings are only included if they model flow between cells at the same vertical height, i.e., in situations in which overland flow boundary conditions are imposed or activated. This restriction maintains the structured property of the preconditioning matrix while still including much of the surface coupling in the preconditioner. Both these adjustments led to considerable speedup in coupled simulations (Osei-Kuffuor et al., 2014).

5.1 ParFlow–Common Land Model (PF.CLM)

The Common Land Model (CLM) is a land surface model designed to complete land–water–energy balance at the land surface (Dai et al., 2003). CLM parameterizes the moisture, energy, and momentum balances at the land surface and includes a variety of customizable land surface characteristics and modules, including land surface type (land cover type, soil texture, and soil color), vegetation and soil properties (e.g., canopy roughness, zero-plane displacement, leaf dimension, rooting depths, specific heat capacity of dry soil, thermal conductivity of dry soil, porosity), optical properties (e.g., albedos of thick canopy), and physiological properties related to the functioning of the photosynthesis–conductance model (e.g., green leaf area, dead leaf, and stem area indices). A combination of numerical schemes is employed to solve the governing equations. CLM uses a time integration scheme that proceeds through a split-hybrid approach, in which the solution procedure is split into “energy balance” and “water balance” phases in a very modularized structure (Mikkelson et al., 2013; Steiner et al., 2005, 2009). The CLM described here and as incorporated in ParFlow is a modified version of the original CLM introduced by Dai et al. (2003), though the original version was coupled to ParFlow in previous model applications (e.g., Maxwell and Miller, 2005). The current coupled model, PF.CLM, consists of ParFlow incorporated with a land surface model (Jefferson et al., 2015, 2017; Jefferson and Maxwell, 2015). The modified CLM is composed of a series of land surface modules that are called as a subroutine within ParFlow to compute energy and water fluxes (e.g., evaporation and transpiration) to and out of the soil. For example, the modified CLM computes the bare ground surface evaporative flux, E_gr, as

where β (dimensionless) denotes the soil resistance factor, ρ_a represents air density [ML⁻³], u_∗ represents friction velocity [LT⁻¹], and q_∗ (dimensionless) stands for the humidity scaling parameter (Jefferson and Maxwell, 2015). Evapotranspiration for vegetated land surface, E_veg, is computed as

where r_b is the air density boundary resistance factor [LT⁻¹], q_sat (dimensionless) is saturated humidity at the land surface, and q_af (dimensionless) is the canopy humidity. The combination of q_sat and q_af forms the potential evapotranspiration. The potential evapotranspiration is divided into transpiration R_pp,dry (dimensionless), which depends on the dry fraction of the canopy, and evaporation from foliage covered by water L_w (dimensionless). L_SAI (dimensionless) is the summation of the leaf and stem area indices that estimates the total surface from which evaporation can occur. A detailed description of the equations PF.CLM uses can be found in Jefferson et al. (2015, 2017) and Jefferson and Maxwell (2015).

PF.CLM simulates variably saturated subsurface flow, surface or overland flow, and aboveground processes. PF.CLM was developed prior to the current Community Land Model (see Sect. 5.2), and the module structure of the current and early versions are different. PF.CLM has been updated over the years to improve its capabilities. PF.CLM was first done in the early 2000s as an undiversified, a column proof-of-concept model, whereby data or a message was transmitted between the coupled models via input/output files (Maxwell and Miller, 2005). Later, PF.CLM was presented in a distributed or diversified approach with a parallel input/output file structure wherein CLM is called as a set sequence of steps within ParFlow (Kollet and Maxwell, 2008a). These modifications, for example, were done to incorporate subsurface pressure values from ParFlow into chosen computations (Jefferson and Maxwell, 2015). These, to some extent, differentiate the modified version (PF.CLM) from the original CLM by Dai et al. (2003). Within the coupled PF.CLM, ParFlow solves the governing equations for overland and subsurface flow systems, and the CLM modules add the energy balance and mass fluxes from the soil, canopy, and root zone that can occur (i.e., interception, evapotranspiration, etc.) (Jefferson and Maxwell, 2015).

At the coupling interface where the models overlap and undergo online communication (Fig. 4b), ParFlow calculates and passes soil moisture and pressure heads of the subsurface to CLM, and CLM calculates and transmits transpiration from plants, canopy and ground surface evaporation, snow accumulation and melt, and infiltration from precipitation to ParFlow (Ferguson et al., 2016). In short, CLM does all canopy water balances and snow, but once the water through falls to the ground or snow melts, ParFlow takes over and estimates the water balances via the nonlinear Richards equation. The coupled model, PF.CLM, has been shown to more accurately predict root-depth soil moisture compared to the uncoupled model, i.e., the stand-alone land surface model (CLM), with capability to compute near-surface soil moisture. This increased accuracy results from the coupling of soil saturations determined by ParFlow and their impacts on other processes including runoff and infiltration (Kollet, 2009; Shrestha et al., 2014; Gebler et al., 2015; Gilbert and Maxwell, 2017). For example, Maxwell and Miller (2005) found that simulations of deeper soil saturation (more than 40 cm) vary between PF.CLM and uncoupled models, with PF.CLM simulations closely matching the observed data. Table 2 contains summaries of studies conducted with ParFlow coupled to either the original version of CLM by Dai et al. (2003) or the modified CLM (ParFlow with a land surface model).

Table 2Selected coupling studies involving the application of ParFlow and atmospheric, land surface, and subsurface models.

“CLM” indicates that coupling with ParFlow was by the original Common Land Model or Community Land Model. “CLM (modified)” indicates that the modified version of
the Common Land Model by Dai et al. (2003) was a module for ParFlow.

Download XLSX

5.2 ParFlow in the Terrestrial Systems Modeling Platform, TerrSysMP

ParFlow is part of the Terrestrial System Modeling Platform TerrSysMP, which comprises the nonhydrostatic fully compressible limited-area atmospheric prediction model, COSMO, designed for both operational numerical weather prediction and various scientific applications on the meso-β (horizontal scales of 20–200 km) and meso-γ (horizontal scales of 2–20 km) (Duniec and Mazur, 2011; Levis and Jaeger, 2011; Bettems et al., 2015), and the Community Land Model version 3.5 (CLM3.5). Currently, it is used in direct simulations of severe weather events triggered by deep moist convection, including intense mesoscale convective complexes, prefrontal squall-line storms, supercell thunderstorms, and heavy snowfall from wintertime mesocyclones. COSMO solves nonhydrostatic, fully compressible hydro-thermodynamical equations in advection form using the traditional finite-difference method (Vogel et al., 2009; Mironov et al., 2010; Baldauf et al., 2011; Wagner et al., 2016).

An online coupling between ParFlow and the COSMO model is performed via CLM3.5 (Gasper et al., 2014; Shrestha et al., 2014; Keune et al., 2016). Similar to the Common Land Model (by Dai et al., 2003), the CLM3.5 module accounts for surface moisture, carbon, and energy fluxes between the shallow or near-surface soil (discretized or specified top soil layer), snow, and the atmosphere (Oleson et al., 2008). The model components of a fully coupled system consisting of COSMO, CLM3.5, and ParFlow are assembled by making use of the multiple–executable approach (e.g., with the OASIS3–MCT model coupler). The OASIS3–MCT coupler employs communication strategies based on the message passing interface standards MPI1/MPI2 and the Project for Integrated Earth System Modeling, PRISM, Model Interface Library (PSMILe) for parallel communication of two-dimensional arrays between the OASIS3–MCT coupler and the coupling models (Valcke et al., 2012; Valcke, 2013). The OASIS3–MCT specifies the series of coupling, frequency of the couplings, the coupling fields, the spatial grid of the coupling fields, the transformation type of the (two-dimensional) coupled fields, and simulation time management and integration.

At the coupling interface, the OASIS3–MCT interface interchanges the atmospheric forcing terms and the surface fluxes in serial mode. The lowest level and current time step of the atmospheric state of COSMO is used as the forcing term for CLM3.5. CLM3.5 then computes and returns the surface energy and momentum fluxes, outgoing longwave radiation, and albedo to COSMO (Baldauf et al., 2011). The air temperature, wind speed, specific humidity, convective and grid-scale precipitation, pressure, incoming shortwave (direct and diffuse) and longwave radiation, and measurement height are sent from COSMO to CLM3.5. In CLM3.5, a mosaic tiling approach may be used to represent the subgrid-scale variability of land surface characteristics, which considers a certain number of patches or tiles within a grid cell. The surface fluxes and surface state variables are first calculated for each tile and then spatially averaged over the whole grid cell (Shrestha et al., 2014). As with PF.CLM3.5, the one-dimensional soil column moisture predicted by CLM3.5 gets replaced by ParFlow's variably saturated flow solver, so ParFlow is responsible for all calculations relating to soil moisture redistribution and groundwater flow. Within the OASIS3–MCT ParFlow sends the calculated pressure and relative saturation for the coupled region soil layers to CLM3.5. CLM3.5 also transmits depth-differentiated source and sink terms for soil moisture including soil moisture flux, e.g., precipitation, and soil evapotranspiration for the coupled region soil layers to ParFlow. Applications of TerrSysMP in fully coupled mode from saturated subsurface across the ground surface into the atmosphere include a study on the impact of groundwater on the European heat wave of 2003 and the influence of anthropogenic water use on the robustness of the continental sink for atmospheric moisture content (Keune et al., 2016).

5.3 ParFlow–Weather Research and Forecasting models (PF.WRF)

The Weather Research and Forecasting (WRF) model is a mesoscale numerical weather prediction system designed to be flexible and efficient in a massively parallel computing architecture. WRF is a widely used model that provides a common framework for idealized dynamical studies, full-physics numerical weather prediction, air-quality simulations, and regional climate simulations (Michalakes et al., 1999, 2001; Skamarock et al., 2005). The model contains numerous mesoscale physics options such as microphysics parameterizations (including explicitly resolved water vapor, cloud, and precipitation processes), surface layer physics, shortwave radiation, longwave radiation, land surface, planetary boundary layer, data assimilation, and other physics and dynamics alternatives suitable for both large-eddy and global-scale simulations. Similar to COSMO, the WRF model is a fully compressible, conservative-form, nonhydrostatic atmospheric model that uses time-splitting integration techniques (discussed below) to efficiently integrate the Euler equations (Skamarock and Klemp, 2007).

The online ParFlow WRF coupling (PF.WRF) extends the WRF platform down to the bedrock by including highly resolved three-dimensional groundwater and variably saturated shallow or deep vadose zone flows, as well as a fully integrated lateral flow above the ground surface (Molders and Ruhaak, 2002; Seuffert et al., 2002; Anyah et al., 2008; Maxwell et al., 2011). The land surface model portion that links ParFlow to WRF is supplied by WRF through its land surface component, the Noah Land Surface Model (Ek et al., 2003); the stand-alone version of WRF has no explicit model of subsurface flow. Energy and moisture fluxes from the land surface are transmitted between the two models via the Noah LSM that accounts for the coupling interface and is conceptually identical to the coupling in PF–COSMO. The three-dimensional variably saturated subsurface and two-dimensional overland flow equations, and the three-dimensional atmospheric equations given by ParFlow and WRF, are simultaneously solved by the individual model solvers. Land surface processes, such as evapotranspiration, are determined in the Noah LSM as a function of potential evaporation and vegetation fraction. This effect is calculated with the formulation

where E(x) stands for the rate of soil evapotranspiration (length per unit time), fx represents empirical coefficient, f_avg denotes vegetation fraction, and E_pot is potential evaporation, determined depending on atmospheric conditions from the WRF boundary layer parameterization (Ek et al., 2003). The vegetation fraction is zero over bare soils (i.e., only soil evaporation), so Eq. (22) becomes

The quantity F is parameterized as follows:

where ϕ is the porosity of the medium, and S_w and S_res are relative saturation and residual saturation, respectively, from the van Genuchten relationships (van Genuchten, 1980; Williams and Maxwell, 2011). Basically, F refers to the parameterization of the interrelationship between evaporation and near-ground soil water content and provides one of the connections between Noah LSM, ParFlow, and thus WRF.

In the presence of a vegetation layer, plant transpiration (length per unit time) is determined as follows:

where C_plant(−) represents a constant coefficient between 0 and 1, which depends on vegetation species, and the G(z) function represents soil moisture, which provides other connections between the coupled models (i.e., ParFlow, Noah, and WRF). The solution procedure of PF.WRF uses an operator-splitting approach in which both model components use the same time step. WRF soil moisture information, including runoff, surface ponding effects, and unsaturated and saturated flow, which includes an explicitly resolved water table, is calculated and sent directly to the Noah LSM within WRF by ParFlow and utilized by the Noah LSM in the next time step. WRF supplies ParFlow with evapotranspiration rates and precipitation via the Noah LSM (Jiang et al., 2009). The interdependence between the energy and land balance of the subsurface, ground surface, and lower atmosphere can fully be studied with this coupling approach. The coupled PF.WRF via the Noah LSM has been used to simulate explicit water storage and precipitation within basins, surface runoff, and the land–atmosphere feedbacks and wind patterns as a result of subsurface heterogeneity (Maxwell et al., 2011; Williams and Maxwell, 2011). Studies with the coupled model PF.WRF are highlighted in Table 2.

5.4 ParFlow–Advanced Regional Prediction System (PF.ARPS)

The Advanced Regional Prediction System (ARPS) is composed of a parallel mesoscale atmospheric model created to explicitly predict convective storms and weather systems. The ARPS platform aids in effectively investigating the changes and predictability of storm-scale weather in both idealized and more realistic settings. The model deals with the three-dimensional, fully compressible, nonhydrostatic, spatially filtered Navier–Stokes equations (Rihani et al., 2015). The governing equations include the conservation of momentum, mass, water, heat or thermodynamics, turbulent kinetic energy, and the equation of state of moist air by making use of a terrain-following curvilinear coordinate system (Xue et al., 2000). The governing equations presented in a coordinate system with z as the vertical coordinate are given as follows.

Equations (26) to (29) are momentum, continuity, thermodynamics, and the equation of state, respectively. The material (total) derivative d∕dt is defined as

The variables v, ρ, T, P, g, F, and Q in Eqs. (26) to (29) represent velocity [LT⁻¹], density [ML⁻³], temperature [K], pressure $\left[{\mathrm{ML}}^{-\mathrm{1}}{\mathrm{T}}^{-\mathrm{2}}\right]$, gravity [LT⁻²], frictional force [MLT⁻²], and the diabatic heat source $\left[{\mathrm{ML}}^{-\mathrm{2}}{\mathrm{T}}^{-\mathrm{2}}\right]$, respectively (Xu et al., 1991). The ARPS model employs a high-order monotonic advection technique for scalar transport and fourth-order advection for other variables, e.g., mass density and mass mixing ratio. A split-explicit time advancement scheme is utilized with leapfrog on the large time steps, and an explicit and implicit scheme for the smaller time steps is used to inculcate the acoustic terms in the equations (Rihani et al., 2015).

The PF.ARPS forms a fully coupled model that simulates spatial variations in aboveground processes and feedbacks, forced by physical processes in the atmosphere and below the ground surface. In the online coupling process, the ARPS land surface model forms the interface between ParFlow and ARPS to transmit information (i.e., surface moisture fluxes) between the coupled models. ParFlow, as a component of the coupled model, replaces the subsurface hydrology in the ARPS land surface model. Thus, ARPS is integrated into ParFlow as a subroutine to create a numerical overlay at the coupling interphase (specified layers of soil within the land surface model in ARPS) with the same number of soil layers at the ground surface within ParFlow. The solution approach employed is an operator splitting that allows ParFlow to match the ARPS internal time steps. ParFlow calculates the subsurface moisture field at each time step of a simulation and passes the information to the ARPS land surface model, which is used in each subsequent time step. At the beginning of each time step, the surface fluxes from ARPS that are important to ParFlow include the evapotranspiration rate and spatially variable precipitation (Maxwell et al., 2007). PF.ARPS has been applied to investigate the effects of soil moisture heterogeneity on atmospheric boundary layer processes. PF.ARPS keeps a realistic soil moisture that is topographically driven and shows a spatiotemporal relationship between water depth, land surface, and lower atmospheric variables (Maxwell et al., 2007; Rihani et al., 2015). A summary of current studies involving PF.ARPS is included in Table 2.

5.5 ParFlow–CrunchFlow (ParCrunchFlow)

CrunchFlow is a software package developed to simulate multicomponent multidimensional reactive flow and transport in porous and/or fluid media (Steefel, 2009). Systems of chemical reactions that can be solved by the code include kinetically controlled homogenous and heterogeneous mineral dissolution reactions, equilibrium-controlled homogeneous reactions, thermodynamically controlled reactions, and biologically mediated reactions (Steefel and Lasaga, 1994; Steefel and Yabusaki, 1996). In CrunchFlow, discretization of the governing coupled partial differential equations that connect subsurface kinetic reactions and multicomponent equilibrium, flow, and solute transport is based on finite volume (Li et al., 2007, 2010). The coupling of reactions and transport in CrunchFlow that are available at runtimes is performed using two approaches. These are briefly discussed below.

The first is a global implicit or one-step method approach based on a backwards Euler time discretization, with a global solution of the coupled reactive transport equations using Newton's method. This global implicit scheme solves the transport and reaction terms simultaneously (up to two-dimensional) (Kirkner and Reeves, 1988; Steefel, 2009). The second is a time or operator splitting of the reaction and transport terms, which is based on an explicit forward Euler method: the sequential noniterative approach, SNIA (in which the transport and reaction terms are solved) (Steefel and Van Cappellen, 1990; Navarre-Sitchler et al., 2011). The stability criterion associated with the explicit approach is that the simulation time step is restricted via the Courant–Friedrichs–Lewy (CFL) condition, under the circumstance that the transportation of mass does not occur over multiple grid cells, but a single grid cell in a time step. Thus, a small time step must be used to ensure this condition holds. This small step size may lead to simulations that will demand much more time to solve Beisman (2007), so more processors are used in order to decrease the processor workload and decrease the solution time of the simulation. The coupling of fully saturated flow to the reactive transport calculations and coupling between a partially saturated flow and transport (flow and diffusion) can be done successively. However, these simulations require calculations of the flow and liquid saturation fields with a different model.

ParCrunchFlow is a parallel reactive transport model developed by combining ParFlow with CrunchFlow. ParCrunchFlow was designed to only be applicable for subsurface simulation. The coupled model relies on ParFlow's robustness ability to efficiently represent heterogeneous domains and simulate complex flow to provide a more realistic representation of the interactions between biogeochemical processes and nonuniform flow fields in the subsurface than the uncoupled model. ParFlow provides a solution of the Richards equation to ParCrunchFlow, which is not present in the biogeochemical code CrunchFlow. ParCrunchFlow employs an operator-splitting method for reactive transport, in which the transport and reaction terms are decoupled and calculated independently. Online coupling between the models is achieved through a sequential noniterative approach, whereby the reaction terms in CrunchFlow's operator-splitting solver get connected to ParFlow's advection terms. ParCrunchFlow takes advantage of the multidimensional advection capability of ParFlow instead of CrunchFlow's advective–dispersive transport capabilities (up to two-dimensional). A steady-state governing differential equation for reaction and advection (with no dispersion and diffusion terms) in a single-phase system is given by

where C_i is the concentration of species i, v represents the velocity of flow, R_i indicates the total reaction rate of species i, and N_tot represents the total species number. In the coupling process, the advection terms are calculated by ParFlow's transport solver through a first-order explicit upwind scheme or a second-order explicit Godunov scheme. Low-order upwind weighting schemes can introduce numerical dispersion, which can impact the simulated reactions, and a comparison of several upwinding schemes can be found in Benson et al. (2017). CrunchFlow calculates the reaction terms using the Newton–Raphson method. For example, in the coupled model ParCrunchFlow, the ParFlow code assigns all hydrological parameters, undertakes the functions relating to parallelization including domain decomposition and message transmission, and solves for pressure and flow fields. The CrunchFlow module is then used to evaluate all reaction terms and conversions between mobile and immobile concentrations. A sequence of simulations of a floodplain aquifer, comprising a biologically mediated reduction of nitrate, has been performed with ParCrunchFlow. The simulations demonstrate that ParCrunchFlow more realistically represents the changes in chemical concentrations seen in most field-scale systems than CrunchFlow alone (summarized in Table 2) (Beisman, 2007; Beisman et al., 2015).