Nonlinear Flow Process: A New Package to Compute Nonlinear Flow in MODFLOW

A new MODFLOW package (Nonlinear Flow Process; NLFP) simulating nonlinear flow following the Forchheimer equation was developed and implemented in MODLFOW‐2005. The method is based on an iterative modification of the conductance calculated and used by MODFLOW to obtain an effective Forchheimer conductance. The package is compatible with the different layer types, boundary conditions, and solvers as well as the wetting capability of MODFLOW. The correct implementation is demonstrated using four different benchmark scenarios for which analytical solutions are available. A scenario considering transient flow in a more realistic setting and a larger model domain with a higher number of cells demonstrates that NLFP performs well under more complex conditions, although it converges moderately slower than the standard MODFLOW depending on the nonlinearity of flow. Thus, this new tool opens a field of opportunities to groundwater flow simulation with MODFLOW, especially for core sample simulation or vuggy karstified aquifers as well as for nonlinear flow in the vicinity of pumping wells.


Introduction
Groundwater flow in porous aquifers is mostly considered to be laminar and thus described by Darcy's law (Darcy 1856), that is, by a linear relationship between specific discharge and hydraulic gradient. However, in some cases where the observed flow behavior is nonlinear, Darcy's law is found to be inadequate (Scanlon et al. 2003;Kuniansky et al. 2008). As widely used groundwater modeling software such as MODFLOW-2005(Harbaugh 2005) is generally based on Darcy's law, there is a lack of tools accounting for such nonlinear flow behavior through porous media or karst aquifers. The Conduit Flow Process (CFP) released by the USGS ) was the first MODFLOW package that allowed the use of a nonlinear (turbulent) flow equation in MODFLOW. This package offers two different modes that can also be combined: the first mode is based on a hybrid approach coupling to MODFLOW a discrete pipe network where flow can be laminar or turbulent (Liedl et al. 2003); the second mode is intended for porous or vuggy karstified aquifers and allows turbulent flow simulations in preferential flow layers where a nonlinear flow law is used if the specific discharge exceeds a given threshold . Reimann et al. (2011Reimann et al. ( , 2012 modified and extended CFP mode 2 to enable the use of different nonlinear flow laws. The aforementioned packages are based on an approach where a switch from the linear Darcy law to the nonlinear equation occurs when the Reynolds number (calculated from the specific discharge and a NGWA.org Vol. 53, No. 4-Groundwater-July-August 2015 (pages 645-650) user-specified diameter of the flow pathways) exceeds a threshold defined by the user. But in contrast to the flow in a single conduit, the wide spectrum of conduit sizes present in real aquifers can be generally expected to result in a more gradual transition from laminar to turbulent flow. Parameter studies based on laboratory experiments (Chin et al. 2009) and lattice Boltzmann simulations (Sukop et al. 2013) using core samples of a vuggy Floridian Aquifer suggest that an approach where Darcy's law is entirely replaced by the nonlinear Forchheimer equation might be more appropriate than a sharp transition from laminar to turbulent flow. Beyond this, a gradual transition to turbulent flow is numerically more stable than a sharp transition. The purpose of this Methods Note is to present such an approach and its implementation into a new package that allows the simulation of nonlinear Forchheimer flow in MODFLOW. This package is not only intended for karst aquifers but can also be employed in porous aquifers where nonlinear flow occurs close to pumping wells (e.g., Forchheimer 1901;Halford 2000;Wen et al. 2011). Nonlinear Flow Process (NLFP) is an open access package and its source code can be downloaded at: http://iewarchiv.uni-graz.at/zusatz/nlfp/. The background theory and the limitations of the package are presented in the first part of the paper. Then, four benchmark tests showing the validity of NLFP are described in the second part. Finally, one additional example comparing NLFP to CFP mode 2 in a more realistic context is shown in the third part.

Background and Theory
The Forchheimer equation (Forchheimer 1901;Chin et al. 2009) extends the linear Darcy law by a quadratic term: where i is the hydraulic gradient [-], q is the specific discharge [L T −1 ], a [T L −1 ] and b [T 2 L −2 ] are the Forchheimer parameters. It reproduces Darcy's law in the limit of small discharges q where a is the inverse of the hydraulic conductivity K 0 [L T −1 ]: In return, it turns into a quadratic flow law in the limit of large discharge. Equation 1 can be transformed to Darcy's law with a discharge-dependent hydraulic conductivity: where is the effective Forchheimer hydraulic conductivity [L T −1 ]. As MODFLOW uses conductance instead of hydraulic conductivity to compute flow, Equation 4 has to be written in terms of conductance too. The linear conductance [L 2 T −1 ] is calculated by MODFLOW using the distance l [L] and the area of the flow cross-section A [L 2 ] between two nodes: In analogy, the Forchheimer conductance [L 2 T −1 ] can be defined as: Since the Forchheimer conductance is dependent on the specific discharge it must be calculated iteratively by a fixed point scheme. Analogous to the treatment of unconfined flow in MODFLOW, the conductance is updated in each iteration step based on reconstructing the specific discharge from the Forchheimer conductance and the hydraulic gradient of the previous iteration step: Here C new F and C old F denote the Forchheimer conductances [L 2 T −1 ] at the new and the previous iteration step, respectively, and h old [-] the head difference between two nodes from the previous iteration step.
Thus, Equation 7 is used to modify the conductance in MODFLOW at each iteration. The computation of C new F occurs only after the read and prepare procedure, where the input data are read and boundary conditions, layer type and no-flow cells are identified. As the inverse of the hydraulic conductivity K 0 is already available, NLFP requires only the ratio b/a as input parameter from the user; typical values of b/a for various materials are provided, for example, by Chin et al. (2009), Sidiropoulou et al. (2007, or Zeng and Grigg (2006). Setting b/a zero corresponds to using the standard Darcy equation. It is noteworthy that in contrast to the CFP package, where the parameters controlling the nonlinear flow are constant within each layer, different values of b/a can be specified for each cell in NLFP to account for inhomogeneity within the layers. In addition to this aquifer parameter, a criterion for convergence of the Forchheimer conductance can be specified by the user, as a convergence of the hydraulic heads alone does not guarantee that fluxes and Forchheimer conductances have also converged. This phenomenon was observed in the second benchmark test and is in general likely to occur under fixed-head boundary conditions.
At present, the NLFP package can be used only with the Block Centered Flow package (BCF). An adaptation to the Layer Flow Package (LPF) will be subject of future work. NLFP was successfully tested with the three solvers DE45, SIP, and PCG2; the latter showed the fastest convergence for the test runs. Just as CFP mode 2, the current version of NLFP computes nonlinear flow within preferential flow layers, but leaves the vertical conductance unchanged. The implementation of Equation 7 for the vertical flow component will be also subject of future code development.
As can be seen from Equation 7 the Forchheimer conductance depends on the head difference between the respective cell and its neighbor. Since head values are not available for inactive or dry cells, NLFP checks and accounts for the status of the neighboring cells before trying to compute the Forchheimer conductance. The package is also compatible with each of the four existing MODFLOW layer types (strictly confined, strictly unconfined, unconfined with constant or variable transmissivity). The NLFP package was compiled with the latest available version of CFP (which can be found on the USGS website). As the source code was written in a single file, a separate compilation of NLFP with the latest version of MODFLOW-2005 is also possible.

Benchmark Tests
In order to verify the correct implementation of the Forchheimer equation in the NLFP source code, four benchmark scenarios are considered.

First Test: Specified Flow in a Confined Aquifer
The first model is composed of 1 row and 11 columns. The size of the cells is constant 10 m × 10 m; the thickness of the single confined layer is also 10 m. The hydraulic conductivity is set to K 0 = 1/a = 1 m/s and b/a = 0.25 s/m. On the left-hand boundary, a constanthead is set to 11 m, while on the right-hand boundary a constant discharge is specified using the WELL package of MODFLOW. The discharge is varied from 10 −6 m/s to 10 3 m/s in steady-state simulations. This scenario allows a direct comparison of the hydraulic gradient calculated by NLFP with that obtained from the Forchheimer equation (Equation 1). As can be seen from Figure 1 NLFP agrees well with the Forchheimer equation over the entire range of hydraulic gradients investigated. For low hydraulic gradients the Forchheimer equation approaches the linear Darcy law, whereas at high hydraulic gradients the nonlinear term of the Forchheimer equation dominates (i = bq 2 ) and thus the graph approaches a power law with an exponent of 2. The latter corresponds to well-known turbulent flow laws such as the Darcy-Weisbach equation (e.g., Reimann et al. 2011).

Second Test: Fixed Gradient in a Confined Aquifer
The second test example comprises a single row of 11 confined cells with constant-head values defined at both boundaries. The size of the cells is 1 m width × 10 m length; the thickness of the single confined layer is 1 m, yielding a model cross-section of 1 m 2 . The other hydraulic parameters are identical to the previous scenario. A sequence of steady-state simulations is conducted varying the hydraulic gradient by several orders of magnitude. It is easily recognized that the exact solution where the head increases linearly from one boundary to the other can be achieved for any value of the Forchheimer conductance, provided that it is the same everywhere. Thus, the correct head values are obtained, while Forchheimer conductance and flux still change from step to step. Therefore, the additional convergence criterion imposed for the Forchheimer conductance is necessary, in particular, for constant-head boundary conditions. Figure 1 shows the results computed by NLFP compared with those obtained from the Forchheimer equation

Third Test: Specified Flow in an Unconfined Aquifer
To verify that NLFP also computes nonlinear flow for unconfined conditions correctly, a corresponding analytical solution is derived. The problem is simplified by considering only one-dimensional (1D) steady-state flow. The model setting, hydraulics parameters and boundary conditions correspond to those of the first test case, that is, Combining Equations 1 and 9 we obtain: Using the polynomial division Equation 10 can be rewritten as: Integrating and substituting the constant-head boundary condition x = 0, h = h 0 finally give the analytical solution: Although this equation cannot explicitly be solved for h, it is sufficient for validating the numerically calculated hydraulic head distribution.
The analytical solution provided by Equation 12 and the respective numerical results are shown in Figure 2 for three different values of b/a (but with a kept constant as 1 s/m) and a specified discharge of 10 m 2 /s. As can be seen NLFP with b/a = 0 s/m (represented by squares in Figure 2) is able to reproduce the water table obtained with the Darcy equation (solid curve in Figure 2). The NLFP package also shows good agreement with the analytical solution for a value of b/a = 0.25 s/m. Finally, the analytical solution of the Forchheimer equation with b/a = 20.25 s/m (represented by circles) is shown. The strong deviation between analytical and numerical solution in this highly nonlinear regime points out the need for an appropriate discretization. Large spatial variations in the water level require an opposite variation in the specific discharge and thus strong variations in the Forchheimer conductance. Therefore, the discretization must be fine enough to capture this variation. The numerical results obtained with the cell size reduced from 10 m to 5 m, 2 m, and 1 m confirm that a finer discretization indeed brings the numerical results close to the analytical solution ( Figure 2). This proves that the numerical implementation is also correct when nonlinear Forchheimer flow is calculated under unconfined conditions.

Fourth Test: Pumping Well in a 2D Confined Aquifer
In order to investigate the validity of the NLFP package in 2D applications, a radial symmetric pumping well benchmark is implemented. The model grid is taken from Kinzelbach and Rausch (1995) and represents only a quarter of the considered pumped field. The cell discretization is fine in the vicinity of the pumping well (see Figure 3) and coarse at the other model boundaries. The model dimensions are 100 m × 100 m and the thickness of the confined aquifer is 1 m. The hydraulic conductivity is set to 10 −3 m/s and b/a to 1000 s/m. The pumping rate is 0.0025 m 3 /s and represents only a quarter of the total pumped discharge. Similar to the previous benchmark, only steady-state flow is considered. Results are presented for MODFLOW and NLFP in Figure 3 and are compared with the drawdown obtained from the corresponding analytical solution of the Forchheimer equation (Bear 1979): where h(r) and H are the hydraulic heads at the distance r [L] from the pumping well and at the radius of well influence R [L], respectively, Q w the pumping rate of the well [L 3 T −1 ], T the transmissivity [L 2 T −1 ], and B the aquifer thickness [L]. As can be seen from Figure 3, NLFP matches well the analytical solution, the differences between the diagonal and the longitudinal head profile being due to the radial flow toward the well. Thus, NLFP computes correctly nonlinear flow close to pumping wells. In summary, the four benchmark tests show that the NLFP package is able to simulate Forchheimer flow under both confined and unconfined conditions as well as for 2D problems. However, the scenarios considered above are highly simplified. The performance of NLFP in a more complex case thus is considered in the next part.

Application Example
The computational performance of the NLFP package is examined using a more realistic scenario taken from  Mayaud et al. (2013). The selected catchment is 3 km long and 1.5 km wide. The model domain is discretized with cells of constant size (10 m × 10 m) and represents a simplified confined karst system that is drained by a single conduit (length 2300 m) created by defining high permeability contrasts with its surrounding matrix. An artificial event with allogenic and autogenic recharge is applied (Figure 4). A first steady-state stress period is used to define the initial conditions for the subsequent transient simulation. The model scenario is implemented to investigate how far a conduit constriction located in the middle of the conduit (represented by reducing the hydraulic conductivity of individual model cells within a sequence of highly conductive cells) can influence the discharge behavior of a karst spring (for more details see Mayaud et al. 2013). In subsequent model runs, the value of b/a is varied from 0 s/m (standard MODFLOW) up to 500.25 s/m.  Figure 4 shows the spring hydrographs obtained with b/a = 500.25 s/m together with those obtained using the standard (laminar) MODFLOW and MODFLOW-CFP mode 2 . The model always converged and the computation time varied from 55 s for linear flow simulated with the standard MODFLOW-2005 to 79 s for NLFP with b/a = 500.25 s/m. The spring discharge simulated by NLFP showed an increasingly damped response of peak flow and baseflow, which is similar to the behavior of CFP mode 2 described by Mayaud et al. (2013). was tested in four benchmark scenarios and shown to work in both confined and unconfined conditions. The package was also successfully applied to a more realistic setting with a larger model domain and a higher number of cells. Depending on the choice of the Forchheimer parameter b/a, which determines the nonlinearity of the equation, NLFP was found to converge only moderately slower than the standard MODFLOW. In general, the smooth transition from linear Darcy flow to nonlinear flow, which results from the use of the Forchheimer equation, is expected to be numerically more robust than the abrupt transition implemented in other nonlinear flow packages for MODFLOW. Compared to CFP, where one value for an entire layer must be specified for both the critical Reynolds number and a representative conduit diameter, the NLFP package requires only one additional aquifer parameter which can be chosen spatially variable. Thus, this new tool opens a new field of modeling opportunities related to nonlinear flow, for example, in (vuggy) karstified aquifers or in the vicinity of pumping wells.