3-D land subsidence simulation using the NDIS package for MODFLOW

The standard subsidence package for MODFLOW, MODFLOW-SUB simulates aquifer-system compaction and subsidence assuming that only 1-D-vertical displacement of the aquifer system occurs in response to applied stresses such as drawdowns accompanying groundwater extraction. In the present paper, 3-D movement of an aquifer system in responses to one or more pumping wells is considered using the new aquifer-system deformation package for MODFLOW, NDIS. The simulation of aquifersystem 3-D movement using NDIS was conducted with a stress or hydraulic head dependent specific storage coefficient to simulate nonlinear deformation behavior of aquifer-system sedimentary materials. NDIS’s numerical simulation for aquifer horizontal movement is consistent with an analytic solution for horizontal motion in response to pumping from a leaky confined aquifer (Li, 2007). For purposes of comparison, vertical subsidence of the aquifer system in response to groundwater pumping is simulated by the both the NDIS and MODFLOW-SUB models. The results of the simulations show that land subsidence simulated by MODFLOW-SUB is significantly larger and less sensitive to pumping rate and time than that simulated by NDIS. The NDIS simulations also suggest that if the total pumpage is the same, pumping from a single well may induce more land subsidence than pumping from multiple wells.


Introduction
Land subsidence has been attributed to the compaction (compression) of aquifer systems caused by ground-water withdrawal.Widespread land subsidence caused by the overexploitation of groundwater is a global problem.Over the past decades, many experts have attempted to develop models to evaluate transient aquifer-system compaction and subsidence accurately and efficiently using MODFLOW (Leake and Prudic, 1991;Hoffmann et al., 2003;Leake and Galloway, 2007).Most of the developed models only account for the 1-D-vertical displacement that is known to be associated with changes in aquifer-system storage, attributed principally to aquitards-fine-grained, low permeability interbeds and confining units.However, the occurrence of horizontal movement caused by ground water withdrawal indicates that aquifer-system deformation and subsidence result from 3-D movement of the aquifer system.The NDIS numerical model presented here simulates 3-D displacement of the solid matrix caused by groundwater extraction, and is based on the principle of effective stress for the case of a constant geostatic load.
The purpose of this paper is to present the NDIS numerical model for 3-D aquifer-system movement due to groundwater withdrawal, and to evaluate NDIS with respect to the simulated aquifer-system vertical movement (subsidence) by comparing simulation results from NDIS with those from MODFLOW-SUB.One of the advantages of the approach presented in this paper is that one can simplify the processes of solving for 3-D aquifer-system deformation with stressdependent aquifer-system specific storage coefficients that change due to compaction driven by hydraulic force loading on the aquifer-system skeleton.

Basic equation
A governing equation for groundwater flow through a saturated anisotropic and hetergenous porous material is expressed by (Bear, 1972): where h is hydraulic head; K is hydraulic conductivity, a diagonal matrix of rank two; W denotes a volumetric flux per unit volume and represents sources and/or sinks of water mass within an aquifer system; and S s (= S sk + S sw ) is the specific storage coefficient, where S sk and S sw are the specific storage coefficients of the soil skeleton and pore water, respectively, and S sk can be either a constant for linear elastic materials or a stress-dependent variable for inelastic materials.The individual solid particles are assumed to be incompressible.The term W denotes a volumetric flux per unit volume and represents sources and/or sinks of water mass within an aquifer system.For an aquifer system that has steady flow without the term W , Eq. ( 1) can be further simplified to the Laplace equation.The US Geological Survey (USGS) developed the 3-D program MODFLOW that numerically simulates steady and transient groundwater flow.To simulate land subsidence in response to groundwater withdrawal, USGS further developed the MODFLOW-SUB package simulating 1-D aquifer-system compaction and subsidence based on vertical stress and strain of the aquifer system.In contrast, a new 3-D model, NDIS developed from DIS (Zhang, 2009) similarly uses the parent program MODFLOW-96 (McDonald and Harbaugh, 1988) and the bulk flux relation through permeable, deformable and saturated porous materials.The bulk flux relation is defined by Helm (1979Helm ( , 1984)): where q b is the bulk flux; v s is the velocity of solids; q is the transient specific discharge that is defined by q = n (v s − v w ); and v w is the fluid (water) velocity.Note that a bolded variable with an arrow represents a vector throughout this paper.For steady flow, Eq. ( 2) can be expressed by the following form: where the subscripts st denote the steady flow.If no aquifer movement is assumed once an equilibrium condition (or steady flow) is reached, i.e., v s−st (x, y, z) = 0, and the bulk flux is assumed to be constant, i.e., q b (x, y, z, t) = q b−st (x, y, z), combining Eqs.
(2) and (3) results in: The expression (4) indicates that aquifer velocity can be written in terms of the specific discharge for transient and steady flow.Moreover, on the right side of Eq. ( 4), the terms q(x, y, z, t) and q st (x, y, z) can be found by running MOD-FLOW for transient flow and steady flow, respectively.It should be pointed out that the assumption of v s−st (x, y, z) = 0 suggests that deformation due to secondary consolidation is negligible and the assumption of the constant bulk flux implies a constant pumping rate.The former means that the aquifer system has no further displacement or deformation once the unsteady flow becomes steady one that is under the equilibrium condition, and the latter can be proved to be true for a constant pumpage (Li and Helm, 2010).Therefore, the velocity field of solids v s for transient flow can be computed as a function of space and time from Eq. ( 4).Using Eq. ( 4), the displacement field of solids, u s can be further computed through integrating the velocity field of solids over the time as follows: where subscript 0 stands for the initial value of displacement and time.The velocity and displacement fields of the solid skeleton vary with time.

Stress-dependent specific storage coefficient
Based on an exponential relation between the effective volume stress and volume strain of the soil skeleton, the specific storage coefficient can be found to be a function of volume stress written in the following form (Li and Ding, 2013): where α is the compressibility of the soil skeleton; S sk0 (= γ w α 0 ) is the initial specific storage coefficient, where the symbol stands for "by definition", γ w is the water unit weight and α 0 is the initial compressibility of the aquifer system; σ v0 is the initial effective volume stress; σ v is the current effective volume stress and the subscript 0 denotes the initial value of a variable throughout this paper.Because the effective volume stress is equal or larger than its initial value (S sk ≤ S sk0 at t ≥ 0), suggests that in response to changes in hydraulic head due to groundwater extraction, an aquifer system with elastic materials may deform larger than that with inelastic materials because of the strain-hardening effect.Equation ( 6) can be written in terms of transient hydraulic head, h through Terzaghi's principle (1925) of effective stress: Hubbert's potential (1940): h = z + p/γ w in terms of hydraulic head in the following form: where σ v is the total volume stress and is normally assumed to be constant; z denotes the elevation from a chosen datum; p is the pore water pressure.Note that S sk (h) in Eq. ( 7) can be used in Eq. ( 1) because MODFLOW-96 for NDIS has been modified to solve the nonlinear partial differential equation (PED) for hydraulic head.For simplicity, the constant S ske0 is applied to Eq. ( 1) for elasticity when σ v ≤ σ v−pre and the head-dependent S skv (h) = S skv0 σ v0 /σ v for inelasticity when σ v >σ v−pre that is preconsolidation stress.

Numerical modeling
The finite-difference numerical simulation is based on a conceptual model that consists of five thin aquitards (thickness is 1 m for each) and two aquifers (thickness is 30 m for each).The top layer (Layer 1) is an unconfined aquifer and the bottom layer (Layer 7) is a confined aquifer.Layers 2-6 are aquitards that serve as a hydraulic separator between two aquifers.The schematic of the seven layer aquifer system is shown in Fig. 1 in which Q is the pumping rate, and H is the total thickness (65 m) of the aquifer system.The option of "no delay" in MODFLOW-SUB is chosen for simulation of the aquifer-aquitard system.
The lateral boundaries of the conceptual model are specified as constant-head boundaries.The initial head is assumed to be zero.Summary properties of the aquifer system are given in Table 1.
To simulate 3-D land movement at a basin-wide scale in an efficient and accurate way, the seven-layer aquifer system is discretized horizontally, and each layer is represented by a grid of 37 rows by 37 columns.The total number of elements (model cells) equals 9583.
Two cases are simulated for purposes of comparison.For one case, a single well is simulated at row 19 and column 19 in the confined aquifer in Layer 7, and for the second case of three pumping wells are also simulated.In the three-well case, the centre well is located at the same location as for the one-well case, and the other two wells are located on a line passing through the centre well.As indicated in Eq. ( 7), S skv requires specification of the four parameters S skv0 , σ v , σ v0 and z, respectively.The parameter S skv0 is given in Table 1 and the values of the other three parameters in each layer are calculated and shown in Fig. 1, where the stress parameters are evaluated using the weight of the overburden material (geostatic load) and the initial head (h 0 ) where σ v0 = σ v − γ w (h 0 − z).The soil unit weight is 21.0 kN m −3 for the aquifer material and 18.5 kN m −3 for the aquitard material.The empirical formulae are introduced to find K 0 , the coefficient of earth pressure at rest, so that the volume stress can be evaluated from the vertical stress.In this paper, K 0 = 0.19+0.233log(PI) for clay and K 0 = 1− sinϕ for sand are used, respectively, where PI and ϕ are the plastic index and the internal friction angle, respectively.The total volume stress can be further calculated using the relation σ v = σ z (1+2K 0 )/3.For simplicity, the preconsolidation head is not introduced to each model layer although NDIS is capable to compute both elastic and inelastic deformation.
Three different values of hydraulic diffusivity of the confining layer, c v = K/S s (1000, 2000, and 4000 m 2 day −1 ) are chosen along with three applied pumping rates (100, 1000, and 10 000 m 3 day −1 ).Three pumping periods (10, 100, and 1000 days) are simulated to evaluate the displacement sensitivity to various parameters.

Resulsts and analysis
Subsidence (cumulative model layer compaction) simulated by NDIS and MODFLOW-SUB is shown in Fig. 2a and b.
A pumping rate of 1000 m 3 day −1 and a pumping period 100 days for both the single-well and three-well are simulated to evaluate the displacements.In the three-well case, the center well is located at the same location as for the onewell case, and the other two wells are located on each side of the center well along a line passing through it.The results show that the maximum subsidence for both models occurs near the well.Subsidence decreases rapidly with distance from the pumping well similar to the nature of the drawdown.The results indicate that the subsidence simulated by MODFLOW-SUB is significantly larger than that simulated by NDIS (Fig. 2).It is interesting to note that even though the total pumping rate in both the one-well and three-well cases are the same (i.e., 1000 m 3 day −1 ), land subsidence in response to three pumping wells is much less than that due to a single well.However, for the three-well case at the location between two pumping wells subsidence simulated by NDIS is slightly larger than that by MODFLOW-SUB (Fig. 2b).For horizontal movement (Figs. 3 and 5), the numerical solutions for the aquifer system with a leaky confining unit are consistent with the analytic solution developed by Li (2007).Namely, the aquifer movement has two zones.One is a compression zone and the other is an extension zone.The position of the boundary between the two zones is a function of time.This is because in the vicinity of the pumping well, the displacement of the solid matrix approaches zero as the well screen is assumed to be permeable only to water.The inward moving grains tend to accumulate in the surrounding area of the well and form a compression zone around the well that extends outward from the well as pumping continues.
The aquifers (Layers 1 and 7) show significant horizontal movement compared with the aquitards where there is negligible horizontal motion.Large radial strain occurs under compression near the well and decreases to zero at a transient point of the maximum radial displacement (Fig. 3).Similar to the relation between vertical displacement and distance from the pumping well, for horizontal displacement beyond the point of maximum horizontal motion the horizontal displacement gradually decreases along the radial direction.The displacement curves for all layers eventually approach the constant head boundary where strain becomes zero.The peak horizontal displacements in the confined (Layer 7) and unconfined aquifer (Layer 1) occur at about 300 and 700 m from the pumping well, respectively (Fig. 3).
Sensitivity analyses were conducted for different pumping rates and discharge periods to evaluate movement sensitivity using the NDIS and MODFLOW-SUB models.The results show that increasing pumping rate and time cause increases in displacements.Figure 4a and b show that subsidence simulated by NDIS is less sensitive than that by MODFLOW-SUB to changes in the pumping rate and time.
Horizontal land movement can result in the formation of earth fissures.Infrastructure associated with human development located in area affected by land-surface movements can result in a variety of potential problems.fined aquifer (Layer 7).The peak value of the horizontal land subsidence in Layer 7 is located 300 m from the center of the well.This location is also a transition point between extensional and compressional strain.
Figure 6 shows the impact of hydraulic diffusivity on the vertical and horizontal movement (Layer 7) using the NDIS model.It is interesting to note that the vertical displacement at any location is decreasing along with increasing hydraulic diffusivity.However, this is true for horizontal displacement only when horizontal distance is less than 1000 m.When the horizontal distance larger than 1000 m, the horizontal displacement is increasing with increasing hydraulic diffusivity.
The plots in Fig. 7 compare the sensitivity of vertical displacement to hydraulic diffusivity for the NDIS and MODLFOW-SUB models.In general, subsidence decreases with increasing hydraulic diffusivity for the both models.However, the results show subsidence simulated by NDIS is less sensitive to hydraulics diffusivity than that simulated by MODFLOW-SUB.

Summary and conclusions
The summary and conclusions of this paper are given below: -A new NDIS module for 3-D aquifer movement has been applied.The NDIS model is incorporated to the parent program MODFLOW-96 that was modified to solve a nonlinear partial differential equation for hydraulic head.
-Subsidence based on a hypothetical conceptual model is simulated using both the NDIS and MODFLOW-SUB models for purposes of comparison.The results indicate land subsidence simulated by MODFLOW-SUB is significantly larger than that simulated by NDIS.
-Two cases are considered for the conceptual model, one case has a single pumping well and the other case has three pumping wells, with two of the wells on opposite ends of a line intersecting, and at equal distances from, a center pumping well.The two cases are simulated with both the NDIS and MODLFOW-SUB models.The results show that for the same total discharge rate, the maximum land subsidence in response to a single pumping well is larger than that from the three pumping wells distributed as described above.
-Given the prescribed hydromechanical properties of the aquifers and aquitards in the conceptual model, the predominant horizontal land movement occurs in the aquifers, particularly the pumped aquifer.Horizontal motion in the aquitard's is relatively insignificant compared with the aquifer's.The characteristics of horizontal movement are consistent with the analytic solution of horizontal movement for an aquifer system with a leaky aquitard (Li, 2007).Namely, groundwater discharge radially induces two zones of horizontal deformation, compression and extension zones.The locations of the peak displacement are different in each of the two proc-iahs.net/372/437/2015/Proc.IAHS, 372, 437-442, 2015 model layers representing the pumped and non-pumped aquifers and are a function of time.
-Sensitivity analysis using NDIS indicates that vertical displacement at a location reduces with increasing diffusivity at any location but this is true for horizontal displacement only within a 1000 m radius.Farther than 1000 m from the pumping well, the horizontal displacement increases with increasing hydraulic diffusivity.
-Comparison of sensitivity analysis between the NDIS and MODLFOW-SUB models shows vertical movement decreases with increasing of pumping rates and time for both the models.However, land subsidence simulated by NDIS is less sensitive to both the pumping rate and time than that simulated by MODFLOW-SUB.

Figure 1 .
Figure 1.A conceptual model of an aquifer-aquitard system and the values of σ v and σ v0 at different depths.

Figure 2 .Figure 3 .
Figure 2. Comparison of land subsidence predicted by the NDIS and MODFLOW-SUB after 100 days pumping; (a) Land subsidence for a single well, (b) Land subsidence for three wells after 100 days pumping.

Figure 7 .
Figure 7.Comparison between NDIS and MODFLOW-SUB models for subsidence sensitivity to hydraulic diffusivity, c v is diffusivity (= K/S s , m 2 day −1 ).

Table 1 .
Properties of the conceptual model.