Importance of Aquitard Response Time for Groundwater Management in Multi-Aquifer Systems Subject to Land Subsidence

Groundwater management models have been widely applied to obtain optimal pumping strategies for land subsidence control, but most of them do not explicitly incorporate land subsidence variables (such as cumulative settlement and land subsidence rates) within the model constraints and neglect the transient effect due to aquitard storage. Here, three operating scenarios of a hypothetical multi-aquifer system, which include a highly compressible aquitard, were implemented with the aim of evaluating land subsidence and identifying management schemes with the support of an optimization model for groundwater management. In a 50-year management period, maximizing pumping while restricting drawdown to 10 m after year 25 stabilizes groundwater levels within the aquifer, but land subsidence continues to reach 4.8 m at year 50. The effect of reducing pumping rates and how early in the management period this is implemented is also analyzed. Restricting the pumping rate as early as year 6 leads to reduced land subsidence at year 50 by 17%. If pumping reduction is delayed, larger land subsidence rates occurred in the system (7.9, 8.3 and 9.6 cm/year in the tested cases); however, if the total settlement is evaluated as a proportion of the thickness of the aquitard, values of the order of 10% are presented. Our results highlight the importance of timely decisions for groundwater management based on the response time of the aquitards in multi-aquifer systems.


Introduction
In multi-aquifer hydrogeologic systems, when compressible aquitards are present or are interbedded within the aquifers, intensive groundwater pumping causes land subsidence, generating multiple adverse effects from extensive damage to infrastructure to a gradual depletion of the storage capacity of the hydrogeological system (Galloway and Burbey 2011; Bagheri et al. 2019).Each of these processes occurs gradually depending on the characteristics of the hydrogeological system; the time scales involved in the system's response to anthropogenic disturbances can be very large and diverse (Rousseau-Gueutin et al. 2013;Currell et al. 2016) and the time for an aquitard to reach equilibrium in the face of a disturbance can take hundreds of years (Rudolph and Frind 1991;Zapata-Norberto 2019;Zapata-Norberto et al. 2018).These diverse time scales must be considered in groundwater management.
Groundwater management models have been widely applied to evaluate measures to counteract the effects of saltwater intrusion, to evaluate social and hydro-environmental policies, and evaluate new sources of water supply (Abd-Elaty et al. 2021;Khatiri et al. 2020;Zhu 2013).Pumping strategies for subsidence control have also been investigated without explicitly incorporating land subsidence variables within the model constraints but focusing on limiting drawdown (Qin et al. 2018;Chang et al. 2007Chang et al. , 2011;;Psilovikos 2006;Steinbrügge et al. 2005;Maddock 1972).Chang et al. (2007) developed an optimal stochastic groundwater management model explicitly considering land subsidence.By neglecting transient effects within the aquitard (assuming that the time required for the aquitard to reach steady state given a perturbation in hydraulic head is smaller than the management time period), they employ an uncoupled representation of groundwater flow and consolidation and the response matrix technique for the deterministic management groundwater model.They consider random hydrogeologic parameters and, by using a first-order variance-estimation method, transfer the original deterministic management model into its stochastic form.Chang et al. (2011) extended the deterministic model to simultaneously consider inelastic and elastic compaction.In these two works the transient effect in the aquitard was neglected, i.e., they assumed that the time required for a perturbation (drawdown) to be transmitted to the full thickness of the aquitard is much less than the time step used in the optimization model.However, this is not applicable in hydrogeological systems that include highly compressible aquitards such as those found in Mexico City, where the aquitard response time can be hundreds of years (Rudolph and Frind 1991;Ortega-Guerrero et al. 1999;Zapata-Norberto et al. 2018;Zapata-Norberto 2019).
Díaz-Nigenda (2022) proposed a methodological approach to identify the transient effect in an aquitard and include it in a groundwater management model for a multiaquifer system, considering land subsidence induced by pumping; their methodology is based on the theory of Bredehoeft and Pinder (1970) and, in order to apply it to highly compressible aquitards, the propagation of the drawdown through the aquitard is computed using Duhamel's theorem.The model integrated three modules that analyze the optimization of the system, the transient propagation of drawdown through the aquitard and the calculation of land subsidence.
Here, the shortcomings of a partial vision in groundwater management are analyzed through an approach that only considers maximizing pumping and restricting drawdown in the aquifer.This approach may lead to recovery of piezometric levels (drawdown is reduced) in some areas, while increasing drawdown in other areas, causing that a greater area may be affected by subsidence.
On the other hand, the importance of considering the response time of the system is also analyzed.In multi-aquifer systems, aquifer units have a characteristic response time which typically can be orders of magnitude smaller than the characteristic response time of aquitard units, and the magnitude of the latter controls land subsidence.In addition, the importance of timely action when dealing with land subsidence and other adverse effects of groundwater overexploitation is emphasized; early detection of adverse effects is important as well as implementing management measures to mitigate the problem.

Materials and Methods
Since hydraulic diffusivity of aquitards is typically several orders of magnitude smaller than that of aquifers, for simplicity we assume that drawdown in the aquifer can be uncoupled from the drawdown in the aquitard during a time step (Neuman et al. 1982).

Methodology for Land Subsidence Analysis
Following Díaz-Nigenda (2022), we employ the response matrix method to compute drawdown in the aquifer at N c selected control points, due to pumping at N well wells.We build a model using Modflow (Harbaugh 2005) and its graphic interface ModelMuse (Winston 2019) to represent a hydrogeologic system composed of an aquifer of thickness b, overlain by an aquitard of thickness b'.From the Modflow model we compute response coefficients m,j,k [L −2 T], which represent the drawdown at the m th control point due to a unit pumping rate at the j th well during the k th time step.
The drawdown within the aquifer, s (m) at the m th control point (observation well) and at the k th time step is computed by where, s (m,k) is the drawdown at the m th control point, [L]; Q (j,k) is the pumping rate of the j th well [L 3 T −1 ] and ΔQ j,p+1 is the change in the pumping rate of the j th well observed in the To compute the propagation of drawdown through the aquitard, Díaz-Nigenda (2022) employ Duhamel's Theorem: where s ′ m z, t k , is the drawdown in the aquitard at control point m and at elevation z.Equation (2) allows extending the methodology proposed by Chang et al. (2007) to highly , where storage effects are important.s ′ m is the drawdown in a onedimensional aquitard (in the vertical direction) due to an instantaneous unit drawdown, s 0 = 1 , at the interface with the aquifer, while the other end remains at a constant head, given by the analytic solution (Hanshaw and Bredehoeft 1968) where, b ′ is the aquitard thickness; K ′ is the hydraulic conductivity of the aquitard; and Ss ′ is the specific storage of the aquitard; t * is dimensionless time (Bredehoeft and Pinder 1970).Equation ( 2) is discretized as Based on (3), Bredehoeft and Pinder (1970) propose that for t * < 0.2 the system is in a transient state, while for larger values of t * the flow within the aquitard reaches a steady state and storage effects can be neglected.Díaz-Nigenda (2022) suggested that transient and storage effects are significant even for values of t * close to 0.4 and that only for t * ≥ 0.5 flow in the aquitard reaches steady state.Therefore, for t * < 0.5 transient and storage effects need to be considered through our methodology.
The thickness of the aquitard b' is divided into N b' segments of constant length Δb � ; then total settlement at control point m and time t k , m t k is computed by where, z p is the elevation of the centroid of the p th segment; w is density of water, g is the acceleration of gravity, and Gand are Lame constants.For simplicity we only consider irreversible consolidation, that is, if drawdown is positive m t k is computed by ( 7), other- wise it is zero.

Numerical Example
A numerical experiment was designed considering a hypothetical multi-aquifer system (Fig. 1) composed by an aquitard that overlies a leaky aquifer.Left and right boundaries are set to constant head (Dirichlet) while upper and bottom boundaries are set to no flow.The aquifer has three zones, each one homogeneous and isotropic but with different hydraulic conductivity, K, and specific storage, S s (Table 1); three extraction wells (A, B, C) and 8 observation wells/control points (a, b, c, d, e, f, g and h) were considered.

Definition of Scenarios
We consider a hypothetical scenario of an aquifer being developed.Pumping rates start from zero and are increased rapidly during the first 5 years to reach 85% of a total annual (a) Controlling drawdown.A maximum allowed drawdown s max is defined, which is enforced at all control points; thus, pumping rates are reduced to ensure compliance in some wells but in others the pumping rates are increased.This strategy is applied from year 26 onwards.Pumping rates that maximize the pumped volume in the system Z subject to drawdown restrictions are computed by solving where N y is the number of years.(b) Intervention time.Pumping rates are reduced or increased by the same procedure, but this reduction or increase is enforced at earlier times.Different intervention times are employed to illustrate that early action is needed to deal with the "memory" effect of land subsidence.
Table 2 shows some of the characteristics of the analysis scenarios.Figure 2 shows the pumping rates considered in the study cases.

Results and Discussion
The hydrogeologic system in Fig. 1 has two relevant and contrasting time scales.The first one corresponds to the time required for a perturbation in the piezometric level due to pumping in each well to reach the nearest lateral Dirichlet boundary.Once simulation time is larger than this time scale, drawdown and hydraulic head in the aquifer tend to stabilize.We can approximate this time scale by t c = L 2 S s K (Currell et al. 2016), where L is the distance from the pumping well to the nearest Dirichlet boundary.In our example, t c ranges from 0.5 to 3 days for the three pumping wells.The second relevant time scale corresponds to the time for a step perturbation in hydraulic head / drawdown to completely propagate through the aquitard and reach steady state.In this case t c = b � 2 S � se K � , which leads to t c = 101.5 years; however, based on the analytic solution given by (3), one can assume that the aquitard is close to steady state for half of that time period, that is, t * = 0.5 , which leads to 50.75 years, hence the 50-year planning horizon adopted in the test cases.

Base Scenario
Figure 2A depicts drawdown in control point PC-a, as well as the pumping rate in the nearest well, A. Drawdown rates in PC-a increase as the pumping rate keeps increasing during the first 6 years; the effect of Dirichlet / constant head boundaries is evident as drawdown  2C, maximum drawdown rates in the first 6 years lead also to maximum subsidence rates (close to 59 cm/year during year 5); total land subsidence is about 5.3 m and more than 50% of this deformation occurs after year 9, continuing during the entire simulation period (rate of subsidence is 0.6 cm/year during year 50), even though drawdown within the aquifer is small after 10 years.
Similar behavior occurs at PC-b (Fig. 3A-C) and PC-c (not shown) in response to pumping wells B and C, respectively.However, in these two areas total drawdown in the aquifer is smaller than in PC-a and therefore subsidence rates and total subsidence are smaller; total settlement at PC-b is 1.8 m (Fig. 3B) and the maximum subsidence rate is 20 cm/year at year 5 (Fig. 3C).

Controlling Drawdown
This scenario is integrated with pumping rates that maximize the pumped volume in the system starting from year 26, while limiting drawdown in the aquifer to less than or equal to 10 m at the control points.Within this scenario, annual pumping rate is reduced at well A and is increased at wells B and C (Table 2).

Pumping Rate Reduction at Well A
This case illustrates the option of controlling land subsidence by reducing pumping rate.Figure 2D shows that reducing the pumping rate at year 26 leads to recovering piezometric levels in that area of the system (from 13.31 m to 10.405 m, Table 2).Although piezometric level in the aquifer recovers quickly at year 26 (Fig. 2D), recovery within the aquitard takes more time (Fig. 2E) and relatively large subsidence rates are maintained at least 6 more years after drawdown within the aquifer has stabilized (Fig. 2F).This extra time needs to be considered in planning and management.This case leads to a total settlement at year 50 of 4.8 m (Fig. 2F), 0.5 m less than the base case (5.3 m, Fig. 2C).

Intervention Time
In this section we analyze the case where the pumping rate is reduced at well A at an early time with the aim of controlling or reducing land subsidence at control point PC-a. Figure 4A depicts the evolution of cumulative settlement while Fig. 4B shows annual rate of subsidence, for four cases: the base scenario, when pumping rate is reduced at year 26 to control drawdown to 10 m, and when this reduction in pumping rate is implemented at years 10 and 6.Although maximum subsidence rates are of a similar order of magnitude (between 58 and 64 cm/year) for all four cases and occur between years 5 and 6 (Fig. 4B), total settlement at the end of the simulation period is 5.3, 4.8, 4.1 and 4.0 m for each case, respectively (Fig. 4A).Evidently, reducing the pumping rate is more effective at Fig. 4 Comparison of total settlement in the aquitard (point PC-a) with different implementation times of optimal pumping as a measure to control land subsidence diminishing total settlement if implemented early, ideally as soon as the negative impact is observed.This might be of the utmost importance in cases of coastal aquifers or near large rivers, as the risk of flooding would increase due to land subsidence.

Discussion and Conclusions
Management strategies to control land subsidence are often seen as stabilizing piezometric levels and, in intensively pumped or overexploited hydrogeologic systems, reducing the disparity in the subsurface water balance.The results obtained here highlight part of the complexity of controlling land subsidence, and the importance of considering the characteristic time scales of the process.It is shown that, in highly compressible aquitards, the characteristic time scale that controls land subsidence can be more than one order of magnitude larger than the characteristic response time of the aquifer.Selection of the management time horizon should consider both time scales.Economic evaluation of the negative effects of land subsidence must consider the fact that ground deformation will continue long after piezometric levels in the aquifer have been stabilized (by implementing managed artificial recharge, for example).Employing the methodology developed by Díaz-Nigenda (2022), which computes the delayed response of an aquitard to pumping in an aquifer by means of Duhamel's theorem, allows accounting for the disparity in the characteristic response times of the aquifer and the aquitard in a computationally efficient manner.Our work expands the applicability of computationally efficient modeling schemes such as Chang et al. (2007Chang et al. ( , 2011) ) to cases when highly compressible aquitards are present in multiaquifer hydrogeologic systems.
It is also shown that identifying the problem and acting early is important and leads to a significant reduction in the observed total land settlement.In this context, an effective monitoring system registering ground deformation (by topographic surveys or remote sensing techniques such as InSAR) and piezometric levels in the aquifer and the aquitards, is desirable.
Implementing optimal pumping in a multi-aquifer system to control land subsidence needs to consider the restriction of subsidence rates and total settlement, and economic costs of damage to infrastructure, in addition to other restrictions to diminish negative impacts due to intensive groundwater pumping.Of course, all actions should consider all groundwater stakeholders to avoid social, economic, and environmental conflicts.

Fig. 2
Fig. 2 Comparison of pumping rates at well A, and aquitard drawdown and land subsidence with the implementation of the base and control drawdown scenarios at control point PC-a

Fig. 3
Fig. 3 Comparison of pumping rates at well B, and aquitard drawdown and land subsidence with the implementation of the base and control drawdown scenarios in control point PC-b

Table 2
Annual pumping rates and drawdown of the test cases