Revisiting MODFLOW’s Capability to Model Flow Through Sedimentary Structures

,


Introduction
Sedimentary basins contain vast amounts of freshwater and need to be carefully managed to secure their longevity. Basins are comprised of discrete sedimentary structures (from herein "structures") of unique shape and orientation following millions of years of deposition, erosion and tectonics. Although basins hosting groundwater resources are predominantly horizontal and extensive, folding and faulting can induce tilting of structures, and erosion can cause units to have irregular lateral boundaries (Catuneanu 2006;Fossen 2016). Furthermore, deposition and consolidation cause sedimentary units to exhibit anisotropic hydraulic conductivity up to 10,000:1 at a regional scale (Michael and Voss 2009;Yager et al. 2009), with maximum conductivity oriented parallel to geologic deposition. Structures of primary interest for water supply are the connected highly permeable sand bodies related to extensive shoreline or dune systems, as well as the connected low permeability structures deposited in low energy settings which can cause compartmentalization of flow systems (Koltermann and Gorelick 1996;Oriani and Renard 2014). Geological models developed for flow modeling in sedimentary basins, where these structures are ubiquitous, must preserve their physical shape and hence connectivity. The flow model must equally ensure that connected and true flow paths are recreated (Renard and Allard 2013). Groundwater models must therefore represent the geometry and principal direction of conductivity associated with sedimentary structures, and solve for the sharp changes in flow direction that these structures generate.
Among numerical groundwater models, MODFLOW is one of the most widely used codes (Harbaugh 2005;Kumar 2019). MODFLOW is open source with a variety of GUIs and an extensive range of plugins to suit various hydrological scenarios. However, MODFLOW has historically been limited to Cartesian discretization and a two-point finite-difference scheme, which is not suitable for geometrically complex hydrogeological units (Hoaglund and Pollard 2003;Li et al. 2010;Gao et al. 2021). The two-point flux approximation (TPFA) in MODFLOW has an inherent orthogonality assumption whereby the conductivity tensor must align with cell boundaries (Harbaugh 2005). Problems arise when structures and their conductivity tensors are nonuniformly tilted in the vertical direction, and when their irregular shape becomes significantly rasterized when mapped onto a structured grid. Accordingly, the modeling of fluid flow in complex sedimentary reservoirs has been undertaken using finite-element (FE) (Yager et al. 2009;Fries et -al. 2014;Borghi et al. 2015) and finite volume (FV) (Rühaak et al. 2008;Panday et al. 2013;Stefansson et al. 2018) methods. The FE method has elaborate meshing capabilities and is mass conservative when utilizing the control volume finite element (CVFE) method (Forsyth 1990;Fung et al. 1992). Similarly, control volume finite-difference (CVFD), when the FV approach is applied to FD, can be used with flexible discretization (Ponting 1989;Heinemann et al. 1991) when combined with a multipoint flux approximation (MPFA) scheme (Edwards and Rogers 1998;Aavatsmark 2002). MPFA utilizes multiple neighboring cells to construct gradient terms, which allows more rigorous treatment of irregular cell geometries and directional anisotropy.
MODFLOW has evolved over many decades in an attempt to better represent the geometry and tensor orientation of complex flow fields. MODFLOW-2000 hydrogeologic-unit flow (HUF) and model-layer variabledirection horizontal anisotropy (LDVA) (Anderman et al. 2002) packages allow upscaling of heterogeneous layers and can solve for anisotropy in the horizontal plane. The DISU and GNC packages of MODFLOW-USG (Panday et al. 2013) allow for flexible and unstructured gridding, and corrections for discretization that does not comply with CVFD requirements. Most recently, MODFLOW 6 ) combined the flexible gridding of MODFLOW-USG and incorporated the XT3D capability which solves the groundwater flow equation using CVFD with a three-dimensional (3D) MPFA. The advances of MODFLOW 6 should therefore overcome many of the perceived limitations of MOD-FLOW codes when simulating the geometries and flow orientation induced by sedimentary structures (Hoaglund and Pollard 2003;Li et al. 2010;Gao et al. 2021). Yet, given its recent release in 2017, MODFLOW 6 and its XT3D capability has been applied in few research studies (Bennett et al. 2019;Goode and Senior 2020).
This study highlights the shortcomings of MOD-FLOW versions prior to MODFLOW 6 in the context of modeling flow through sedimentary structures and evaluates MODFLOW 6's improved capacity to address them. We highlight the limitations of older MODFLOW versions by turning our attention to structures which are geometrically complex or which have directional anisotropy, both of which are almost always present and control flow in sedimentary basins (Yager et al. 2009). We initially undertake benchmark simulations to assess traditional and contemporary discretization methods and solvers in MODFLOW, before applying our findings to a case study which compares modeled hydraulic heads and groundwater fluxes between MODFLOW 6 and prior versions. Finally, we discuss the significance of discretization methods and the use of MPFA for real world applications.

Benchmarks Purpose
Groundwater flow is directed by sedimentary structures because of their boundaries of contrasting conductivity, as well as directional anisotropy. Therefore, groundwater models must preserve the geometric shape of structures, as well as the directional alignment of principal conductivity. The purpose of this section is to benchmark MODFLOW's past and current ability to estimate flow using different types of grids (1), along structure boundaries (2) and with a dipping conductivity tensor (3). These requirements may be distilled by benchmarking MOD-FLOW's ability to estimate flow through a dipping highly conductive anisotropic layer in a low conductivity domain.

Methodology
The individual impact of grid design, geological boundaries and anisotropic conductivity can all be tested using a simple model consisting of an idealized tilted channel ( Figure 1). The idealized channel eliminates irregularities from real-life structures and allows comparison of numerically simulated hydraulic head and flow to an analytical solution. The benchmark model is in two-dimensional and can be considered as either plan or in transect, representing a sand channel or dipping sand layer, respectively. Both scenarios exhibit an identical analytical solution.
A head boundary function was imposed to induce a hydraulic gradient in the same orientation as the channel: where h(x, y) is the head at cell center coordinates at x and y (or z in transect), i u is the hydraulic gradient parallel to the direction of the channel, and θ is the orientation of the channel from the x -axis. Parameters used for the benchmark are presented in Table 1. Hydraulic conductivity of the structure is treated as isotropic for two analyses (Sections 2.3.1 and 2.3.2) and anisotropic for the last analysis (Section 2.3.3). A very low permeability domain (K domain ) surrounds the sand facies so flow is constrained to within the sand to allow an analytical solution for both flux and volumetric flow. This novel setup is particularly useful for the purpose of computational benchmarking, as it allows comparison with the analytical solution of the flow. The analytical flux (q u ) within the channel is equal to K max × i u , that is, 1 m/d in the direction of θ given our set of parameters.  Despite the hydraulic similarity of the plan and transect scenarios, the discretization options of MOD-FLOW are specific to each orientation (presented with results in Figure 3). Plan scenarios consider both rectilinear discretization and flexible triangular discretization, implemented so that they contain a similar number of cells. Flexible discretization allows cell edges to align with the channel boundary, generating realistic geometries with fewer cells, while hydraulic properties are assigned to rectilinear grid cells according to the proportion of cell area in the channel boundary. Flexible discretization was achieved using a triangle mesh generator (Shewchuk 1996) as inputs to the DISV package available in MODFLOW 6. Transect options are limited by MODFLOW's flat cell top and bottom requirement, meaning flexible grids in transect is not possible. In transect, we consider rectilinear as well as vertically offset discretization. Vertically offset represents hydraulic units as continuous layers where the depth and thickness of model layers are set to represent the geological layer. Rectilinear requires more cells than vertically offset, as well as preprocessing to assign hydraulic conductivities to cells which do not necessarily follow hydrogeologic layers (Anderson et al. 2015). Jupyter Notebooks for the benchmarks are included as Supporting Information.

Capability to Model Flow Through Dipping Isotropic Layer
The benchmark compares numerical head and flux using three different grid construction methods: rectilinear (plan and transect), flexible triangular (plan) and vertically offset (transect). Here, we verify whether flow in dipping layers is robustly solved for all grid construction methods by measuring the flux value at the center of the channel far from the edges, where boundary effects are minimal, and comparing this to the analytical solution. Those results are visualized in Figure 2.
We can observe that the analytical flux, 1 m/d at 30 • , is reproduced with the rectilinear grid ( Figure 2a). It should be noted that MPFA essentially reduces to TPFA ( Figure 2b) for isotropic rectilinear grids as this is already accurate . Without MPFA, the flexible grid produces perturbed heads and fluxes ( Figure 2c) caused by the distorted grids and cell center connections which are not perpendicular to cell faces. However, this issue is resolved with MPFA ( Figure 2d).
Vertically offset discretization produces a +15% error in flux magnitude, and more significantly a nearly horizontal head gradient and flux direction, rather than at 30 • and aligned with the channel (Figure 2e). The flux direction error is not a problem of head boundary values given that the homogeneous case, that is when K domain is equal to K channel , produces the expected diagonal flow (Appendix S1). Preliminary analysis (Provost, A. and Langevin, C, personal communication, August 2022) suggests that the primary source of error is insufficient connectivity, and therefore insufficient hydraulic communication, between cells in adjacent grid layers in the channel. An interpolation scheme such as XT3D is unable to compensate for the error (Figure 2f) because the issue is largely one of connectivity, and not simply a matter of irregular geometry. The error in the flux magnitude and direction is particularly evident in the heterogeneous case of our benchmark due to the relatively large vertical offsets between cells and the hydraulic isolation of the channel from the rest of the model domain, which serves to suppress the vertical flows between grid layers that are necessary for an accurate solution within the channel. Further investigation is needed to verify this hypothesis and develop a practical approach to incorporating the appropriate grid connectivity into MODFLOW models with large vertical offsets between cells.
The outcome of this benchmark is that diagonal head gradients are well represented for rectilinear and flexible grids, but vertically offset grids are not yet suitable for steeply dipping heterogeneous layers where a diagonal head gradient is imposed.

Capability to Model Flow Along Dipping Layer Boundary
Rectilinear and vertically offset discretizations, which are the only options for MODFLOW versions previous to MODFLOW 6, result in the rasterization of structure boundaries (Figure 3c). We investigate errors due to this approximation by drawing our attention to fluxes at the channel boundary, expecting flow to be parallel with the boundary. Therefore, we examine fluxes at the channel boundary for all discretization scenarios (Figure 3a through 3c), as well as translating this to implications for total flow by integrating fluxes (Figure 3d). We also investigate the effect of resolution by increasing the number of cells used in the simulation.
As expected from its original purpose, flexible discretization has no boundary effect, with fluxes at the boundary staying consistent (Figure 3a), and the volumetric flow conforming to the analytical solution even for coarse resolutions (green line in Figure 3d). Despite its visually rasterized representation, the vertically offset grid does not suffer from boundary effects (Figure 3b), an artifact of inappropriate hydraulic connectivity. However, the overestimated flux as previously identified causes the volumetric flow through the channel to be equally overestimated and does not improve with resolution (blue line in Figure 3d). The rectilinear grid, however, produces irregular flux magnitudes and directions for cells along the channel boundary (Figure 3c), which we coin as the "staircase effect" given the way that flow must make its way along the boundary in a convoluted way. The retardation of fluxes at the boundary subsequently reduces the volumetric flow through the channel, which improves with resolution due to a smaller percentage of channel cells being located along the boundary (orange line in Figure 3d). A grid convergence study was done for the rectilinear scenario showing that convergence of volumetric flow to the analytical follows a power law of exponent −1, consistently found with the staircase effect (Kereyu and Gofe 2016). Therefore, to restrict the error of the flow through a 30 • tilting channel to within 5% and 10%, 20 and 10 cells per channel width, respectively, are necessary. While the rectilinear staircase error is resolution-dependent, the vertically offset error is not. Therefore, if rectilinear discretization is not able to adequately refine dipping structures due to computational limitations, then vertically offset would be the less erroneous option. The transition point of the required rectilinear resolution to better vertically offset for the benchmark is 5 cells per channel width, the point where the orange and blue line crosses in Figure 3d. The transition point would depend on the angle of the structure.
This benchmark highlights one of the benefits of flexible gridding, which represents accurate fluxes and flows even at a coarse scale, whereas rectilinear grids require more refinement due to rasterization of boundaries.

Capability to Model Flow Through Dipping Anisotropic Layer
Simulation of flow through dipping anisotropic layers requires every component of the hydraulic conductivity tensor (Aavatsmark et al. 1998), which MODFLOW codes prior to MODFLOW 6 were not designed to handle. Instead, MODFLOW 6's MPFA scheme, activated using the XT3D option, uses a larger stencil and weighted interpolation to allow for the consideration of the full tensor. Here we benchmark MODFLOW 6 and prior MODFLOW versions capability for modeling flow in a dipping anisotropic channel. Given that anisotropy is typically considered between the horizontal and vertical direction, and that intra-cell fluxes computed on vertically offset grids are still subject to significant loss of accuracy in steeply dipping channels, we focus on the transect scenario, with the rectilinear grid. The previous simulations have used an isotropic conductivity tensor, and now we replace this with a varying anisotropic tensor in the direction of the channel (30 • ).
We find that the flux magnitude when using TPFA is extremely diminished, up to 0.28% of the analytical flux for 1000:1, while MPFA matched the analytical flux to 99.97% (Table 2).
We extend this particular benchmark to include multiple dips and anisotropy ratios for TPFA and MPFA, and then recover the effective conductivity ellipses based on outputted fluxes. Please refer to Appendix S2 for methodology and Jupyter Notebook S3. The effective conductivity ellipses (Figure 4a) show that MPFA matches the input conductivity, whereas the TPFA ellipses are diminished, at its worst at 45 • . The TPFA method essentially truncates the inputted conductivity tensor by setting the off-diagonal components of the K tensor to zero, and uses the x and y components of the ellipse radius at angle theta for the diagonal components ). As we increase the anisotropy ratio, we increase the magnitude of the error (Figure 4b). Effective conductivity drops steeply to only a fraction of its input  Table 2 Computed Fluxes, q mag (m/d), Flux Direction, q theta ( • ) and Model Run Time (s) for the Benchmark in Figure 1 with an Anisotropic Conductivity Tensor  value for high anisotropy, for example a 3 • dip with 1000:1 anisotropy produces an equivalent conductivity of only 27% of its actual value. This benchmark highlights that MODFLOW versions prior to MODFLOW 6 are far from robust for modeling flow through highly anisotropic and dipping layers, and that MODFLOW 6's XT3D capability must be used. However, it should be noted that at high anisotropy ratios, MPFA can be susceptible to spatial oscillations in the solution, and users of the option should critically evaluate results to check there are no sources or sinks present that are unrealistic ). An example that focuses on correcting this issue, in this case for the dispersion tensor rather than the conductivity tensor, is for MT3DMS, the transport package used with MODFLOW (Yan and Valocchi 2020).

Case Study
Here, we present a case study where we apply our recommendations for modeling structure with MODFLOW on a real basin with an existing groundwater model. The site exhibits multiple superimposing aspects associated with sedimentary structures, including irregular shaped structures in plan and dipping anisotropic layers, as well as sharp changes in flow direction. The purpose is to showcase the resulting differences in modeled flow patterns between using traditional and contemporary MODFLOW discretization methods, and once again assert the significance of MODFLOW 6's new packages.

Perth Regional Aquifer Modeling System
The Perth Basin contains vast amounts of groundwater which supplies approximately 70% of water supply for Perth, Western Australia (De Silva et al. 2013). It is a sedimentary rift basin comprising around 12 km of sediments and hosting multiple major aquifers (Davidson 1995). Of note is the superficial unconfined aquifer which overlies the Leederville and Yarragadee confined aquifer systems. These three systems are separated by two major confining units, the Kardinya and South Perth Shales, respectively, which are often thick and extensive, but not conformable at all locations. The aquifers are bounded to the east at the Darling Fault by Archean crystalline rocks, and groundwater flows predominantly west toward the ocean, although there is some mounding in recharge zones which influences local flow patterns. Recharge to the superficial aquifer is primarily through rainfall, with leakage replenishing the underlying confined aquifers where aquitards are absent. The Swan River is a major drainage pathway running North-East to South-West (see Supporting Information, Figure S1 for a locality plan). The Perth Regional Aquifer Modeling System (PRAMS 3.5.2) is a coupled recharge and groundwater model, initially built in the early 2000 s for the purpose of managing the complex groundwater needs of Perth (De Silva et al. 2013;Siade et al. 2017). PRAMS covers a surface area of 9100 km 2 and utilizes a modified version of MODFLOW-2000 that incorporates a landscape model known as the Vertical Flux Manager (Silberstein et al. 2009), which employs the WAVES software for simulating flow though the unsaturated zone and ultimately groundwater recharge (Dawes et al. 2012). PRAMS 3.5.2 is discretized into 500 × 500 m cells in plan and consists of 454 rows and 214 columns with approximately 600,000 active cells and runs for 13 years with monthly time steps resulting in a CPU runtime of approximately 2 to 3 h on most common desktop computers (Siade et al. 2020). Thirteen vertically offset layers represent the main sedimentary units, which are mostly extensive with some pinch-outs. The Kings Park formation is an intrusive unit which incises multiple layers and is represented as a change in conductivity in the incised layer. The structure of the basin results in dipping of up to 15 • , such as the South Perth Shale around the Joondalup area.

Methodology
Our aim is to compare the current PRAMS model which uses MODFLOW-2000, with its MODFLOW 6 counterpart to investigate the impact of available discretization options and numerical schemes on modeled flow patterns.
We extract two submodels from PRAMS, one in plan and the other in transect. In plan, the submodel is centered on the intrusive Kings Park formation as it represents irregular geometry and an ideal candidate to apply flexible gridding. The transect submodel runs 11.5 km in a north-south direction across a major dip. The submodels are created by extracting arrays of horizontal and vertical conductivity from PRAMS, as well as head values and recharge for a typical snapshot in time. The models are run in steady-state to observe how model discretization affects long-term flow patterns. More details on methodology specific to each orientation are discussed below.

Plan Submodel
The plan submodel was created from PRAMS Layer 7 which represents the Wanneroo Member of the Leederville aquifer, which is incised from the west by the Kings Park formation. The Wanneroo Member has a relatively high horizontal conductivity (estimated at 6 m/d), with the Kings Park Formation being less conductive at 1 m/d.
Three discretization scenarios in plan were created ( Figure 5 and Table 3): (1) the original PRAMS discretization; (2) flexible triangular, as recommended from the benchmarks; and (3) fine rectilinear, which was shown from our benchmarks to accurately represent complex flow and corresponds to our reference case. The flexible gridding scenario can be locally refined and coarsened in different ways, however here we aim to have a similar number of cells as PRAMS. Heads along all four boundaries were extracted from the PRAMS submodel and marked as constant heads, and linearly interpolated for cell centers along the flexible and rectilinear model boundaries.
Fluxes for the flexible grid are streamlined along the formation boundary and better represent fine rectilinear despite using approximately the same number of cells as PRAMS ( Figure 5). Rasterization of the boundaries of the Kings Park Formation using PRAMS causes the formation to artificially narrow resulting in a misleading flowpath (point A) and induces flow that is not streamlined around the formation edge due to the staircase effect (point B). Instead, the flexible grid represents the flow regime around the formation edge and allows for more refinement flexibility, which is key for maximizing numerical accuracy in critical areas such as along model boundaries, geological boundaries and pumping zones.

Transect Submodel
The transect submodel investigates differences in modeled head and flux for dipping anisotropic layers. The transect submodel contains a maximum dip of 7 • and maximum anisotropy ratio of 30,000:1 based on the PRAMS calibration. Three scenarios were considered in transect: (1) original PRAMS discretization which uses a vertically offset grid with TPFA; (2) coarse rectilinear discretization with rotated conductivity tensors and MPFA, as recommended from the benchmarks; and (3) fine rectilinear discretization with rotated conductivity tensors and MPFA as our reference case as it accurately represents flow (Table 3). PRAMS contains 13 vertically offset layers, each layer representing a significant geological sequence in the Perth Basin. Hydraulic properties are overlayed onto the coarse and fine rectilinear models, with careful refinement so that the true thickness of geological layers is well represented and not over rasterized. Computation of the 9000 cells in the coarse rectilinear model is expensive when compared to the 299 cells used in PRAMS. However, this number could be reduced by refining in the Z direction as needed, instead of applying a blanket Z direction discretization step.
Conductivity tensor angles used in the rectilinear scenarios are calculated using the dip of the bottom of each geological formation, although this could have been done instead for the center of the layer, depending on the interpreted bedding angle. Generally however, the formation bottom best represents depositional angle given that formation bottoms are generally conformable. The left and right boundaries are assigned constant heads, the bottom a no-flow condition and recharge is added along the top. Boundary heads are linearly interpolated for the rectilinear scenarios. Modeled fluxes for each method are compared by integrating fluxes vertically across the transect for cells pertaining to each of the 13 original layers.
There are striking differences in flux results between PRAMS and coarse rectilinear, particularly around sloping units ( Figure 6). First, we can observe, by comparing transects showing vertical flux (Figure 6a and 6b), that the rectilinear scenario introduces a wider range in vertical fluxes within the Wanneroo and Yarragadee aquifers, and that flow direction is influenced by the shape of the aquitard. The PRAMS model significantly under estimates vertical fluxes, with almost no vertical flux modeled in the Wanneroo layer, yet greater than 0.00015 m/d for some sections in the rectilinear scenario ( Figure 6c). Second, horizontal flux in the primary aquifer of interest, the Wanneroo Member, is approximately 10 to 20% more in PRAMS than the rectilinear scenario ( Figure 6d). As in our benchmarks, these results show that using TPFA on vertically offset grids can significantly underestimate vertical fluxes. This may result in contaminants being modeled as moving primarily horizontally, which may not be accurate for sloping layers.

Conclusion
This paper revisits MODFLOW's capability to model flow through sedimentary structures, which are typically irregular in geometry with anisotropic conductivity. Benchmark simulations have highlighted the limitations of MODFLOW 2005 and prior. The rasterization of hydraulic properties onto a structured grid causes a "staircase error" whereby modeled flow is retarded along interfaces between sedimentary structures where there is a sharp change in conductivity. We showed that the flow through a 30 • tilted sand channel is underestimated by 8% due to the staircase error. More important is the gross under-estimation of flow through anisotropic dipping layers, associated with folded or tilted sedimentary structures, due to significant truncation of the conductivity tensor in the absence of an MPFA scheme. When MODFLOW-2005 is used to model flow through the same 30 • tilted sand channel, flux within the channel is reduced to only 24%, 3%, and 0.3% of its actual value for anisotropy ratios of 10:1, 100:1, and 1000:1, respectively. However, we show that MODFLOW 6 has overcome some of these issues by including flexible discretization (DISV package) and an MPFA scheme (XT3D option). This has given MODFLOW similar capability as existing CVFE codes, such as SUTRA (Provost and Voss 2019), allowing robust modeling of flow through dipping anisotropic layers.
Persisting issues however include the pervading assumption that cells have flat tops and bottoms in transect, which produces the staircase error as well as requiring a large number of cells to adequately represent thin dipping layers using a regular grid. We show that vertically offset grids, while a numerically efficient way to represent gently dipping layers, can still be subject to significant loss of accuracy in head and flow estimates when the dip is steep, especially if the dipping layer is much more conductive than the surrounding material. We hypothesize that vertically offset grids to not provide sufficient connectivity between adjacent model layers to accurately simulate the steeply dipping, highly heterogeneous case.
The main recommendations from this research are that (1) vertically offset grids are further investigated as to their source of error, and that appropriate grid connectivity be incorporated in MODFLOW models with large vertical offsets between cells, (2) modelers of significantly dipping and anisotropic layers should be prepared to employ an expensive but necessary approach of overlaying properties onto a regular grid in transect and applying the XT3D option, or risk grossly underestimating groundwater flow magnitude; (3) given the increased computation required for MPFA solutions (refer Table 2), further research should focus next on automated refining and upscaling procedures, such as that suggested by Sbai (2020), and (4) FV models with unstructured discretization in 3D, such as the use of tetrahedron cells, would be a logical direction forward to truly represent the geometry and flow pathways generated by sedimentary structures. Figure S1. Locality plan for case study.