On the possible contribution of clayey inter-layers to delayed land subsidence above producing aquifers

In recent years, measurements of land subsidence above pumped aquifers by permanent GPS and InSAR have exhibited some delay relative to drawdown ranging from months to years. The current modeling approaches accounting for water fluid dynamics and porous medium geomechanics may fail to predict such a delay and may underestimate the land settlement after the well shutdown. In the present communication, an investigation is made on the residual compaction of the intervening clayey formations as a possible contribution to retarded land subsidence. The pore pressure variation within the aquifer and its propagation in the clay are simulated by a finite element flow model, with the resulting pore pressure decline used as input data in a hypoplastic geomechanical model. A proper sensitivity analysis on (i) aquifer depth, (ii) ratio between the sandy and the clayey layers thickness and hydraulic conductivity, (iii) oedometric compressibility in first and second loading cycles, is performed for a typical geology of a Quaternary sedimentary basin. The results show that a certain fraction, up to 20 % of the overall land subsidence, can take place after the shutdown of the producing wells depending on actual basin, litho-stratigraphy and parameter values.


Introduction
A major consequence of groundwater pumping is anthropogenic land subsidence.This can be a matter of concern if the affected areas are highly urbanized, with the loss in ground elevation generating possible structural problems to the buildings.If the aquifers underlie the coastland, the environmental problems are the exposure to flooding during high tides and severe sea storms.Hence, land subsidence must be investigated, monitored and reliably predicted to implement adequate remedial measures.
Extensometric records of soil deformation combined with land subsidence measurements by permanent GPS stations and InSAR techniques have provided evidence of some delay between the well shutdown or a significant reduction of the pumping rates and the resulting land settlement (Hettema et al., 2002;Teatini et al., 2006;Wu et al., 2010;Galloway and Sneed, 2013;Chang et al., 2014).While pore pressure in the aquifers progressively recovers, the subsidence may still continue, with a delay ranging from a few months to a few years.This suggests that the clayey layers confining the aquifer may keep on depressurizing and depleting after the shutdown of the producing wells.In the present study an investigation is made on the possible correlation between the delayed subsidence and the depletion of the clayey interlayers of a multi-aquifer system.
The current standard modelling approach accounting for water fluid-dynamics and porous medium geomechanics usually does not accurately address these inter-layers (or aquitards) and tend to underestimate the land settlement after field abandonment.Before using a more complex modelling approach, such as visco-elastic or visco-plastic (Wu et al., 2010;Chang et al., 2014), a sensitivity analysis is performed with the aid of a 1-way coupled hysteretic hypoplastic model (Gambolati et al., 2001) in which the pore pressure variation propagates in the clayey inter-layer.The sensitivity analysis takes into account the main parameters controlling land subsidence: aquifer depth, ratio between the sandy and the clayey layers thickness and hydraulic conductivity, and oedometric compressibility in the first loading and second unloading-loading cycle.
Published by Copernicus Publications on behalf of the International Association of Hydrological Sciences.

Mathematical model
Aquifer dynamics and land subsidence caused by water withdrawal is theoretically described by the coupled process involving mechanics and flow in porous media (Biot, 1941).Uncoupling the flow field from the stress field is a usual assumption in groundwater hydrology.Gambolati et al. (2000) have shown that uncoupling has generally no measurable influence over any timescale of practical interest.Hence, in the present analysis a 1-way coupled approach is followed where the pore pressure variation is first simulated by a 3-D groundwater flow model, and then the land settlement is predicted by a 3-D geomechanical model with the pore pressure change used as an external distributed source of strength.
The aquifer hydrodynamics relies on the classical groundwater flow equation: where k is the hydraulic permeability tensor, h is the hydraulic head change, t is time, q the source/sink and S s the specific elastic storage defined as S s = γ (c M + ϕβ) with γ the specific weight of water, c M the oedometric compressibility of the porous medium, ϕ the porosity and β the volumetric fluid compressibility.A nonlinearity is present in Eq. (1) since S s is a function of the vertical effective stress σ z , hence h, via c M .Equation ( 1) is solved in space by the finite element method using linear tetrahedral elements and in time by an implicit finite difference scheme, with the nonlinearity addressed by a Newton-like iterative method.
The incremental pore pressure variation p = h • γ induced by water pumping is implemented into the geomechanical model to compute the medium displacements.The equilibrium equations along the coordinate directions x, y and z in an isotropic medium can be written in terms of incremental soil displacement u (Verruijt, 1969): Table 1.The values of the four parameters used in the simulations: depth H a and permeability of the producing aquifer, thickness of the sandy ( S s ) and clayey ( S c ) layer, ratio between the oedometric compressibility in first and second loading cycle and ratio between the sandy and the clayey hydraulic conductivity.
where u x , u y and u z are the components of u along the coordinate directions, λ and G are the Lamé constant and shear modulus of porous medium, respectively, and ε is the volume strain.A nonlinearity is introduced in Eq. ( 2) since λ and G are functions of σ z , and hence u, via c M as with ν the Poisson ratio.Using data based on field and lab measurements the compressibility c M is related to σ z , with a different relationship in first and second loading cycle accounting for geomechanical hysteresis (Baù et al., 2002).Equation ( 2) is solved by the finite element method using the same tetrahedral elements as the groundwater flow model.The nonlinearity is solved by a Newton-like iterative method.

Modeling set-up
The geometry of the porous system used as test problem for the sensitivity analysis is shown in Fig. 1.The cylindrical domain of the geomechanical model is characterized by a radius R = 10 km and height H = 5 km, with zero displacements prescribed on the bottom and the outer boundaries.An aquifer is located at depth H a and confined by two clayey layers.An accurate discretization of 1 m thick elements is used to simulate the pore pressure propagation at the sand-clay interface and the in less permeable units along the vertical direction.A zero pore pressure is prescribed on the side surface and the natural condition is prescribed on the top and bottom boundaries.The aquifer permeability decreases with depth as summarized in Table 1 with typical values of a Quaternary sedimentary basin.
A 10-year constant water withdrawal rate Q * is prescribed from a 500 m-radius cylindrical aquifer volume to simulate the pore pressure decline experienced by several producing wells.The fluid-dynamic simulations last for 10 years after the cessation of the production.The value Q * is calibrated to have a uniform 100 bar pressure decline in the produced volume at the end of the withdrawal period.
The sensitivity analysis takes into account four parameters: -the depth H a of the sandy layer; -the ratio between the sandy and the clayey layer thickness, S s and S c respectively, with S s + 2 S c = 100 m; -the ratio C r between the oedometric compressibility in first e second loading cycle; -the ratio R k between the sandy and the clayey layers permeability.
The values of the four parameters used in the simulations are summarized in Table 1.

Numerical results
Figure 2 shows the pore pressure variation along the symmetry axis for the case with H a = 1000 m, S c = 45 m, R k = 1 × 10 −6 and C r = 8: the pore pressure recovery in the clay layers is significantly delayed relative to the aquifer.The pore pressure in the aquitards may still decline after the shutdown of the producing wells, possibly causing a delayed land subsidence.
The results of the geomechanical model are presented in terms of the maximum vertical land displacement normalized with the respect to the value at the end of the water abstraction (t = 10 years).A normalized displacement larger than 1 highlights a delay between the pore pressure in the aquifer and the land settlement.The results obtained for the investigated H a values and S s = 10 m are shown in Figs.3-6.
Note that the sudden increase of the aquifer stiffness at the well shutdown generates an abrupt increase of land subsidence although the pore pressure recovers in the aquifer.This 3-D deformation process, which is mainly accounted for the Poisson ratio and, subordinately, the depth of the depleted layer, was already pointed out by Ferronato et al. (2001).
Small values of R k and S s and large values of H a and C r enhance a delayed land settlement because the pore pressure propagates in the clay more slowly.For all the investigated H a values the delay between the pore pressure recovery and the land settlement is maximum with S s = 10 m, C r = 8 and R k = 1 × 10 −6 .With these configurations, a 9 and 20 % of the final land subsidence takes place after the cessation of withdrawal for H a = 1000 and 2000 m, respectively.With H a = 500 m subsidence is almost unvaried after the well closure.

Conclusions
A set of numerical simulations are performed using a 3-D 1way coupled groundwater flow and geomechanical model to investigate the possible contribution of clayey inter-layers to the delayed land subsidence measured above deep pumped aquifers.
The parameters addressed by the analysis include: aquifer depth H a , thickness of the sand S s and of the clay S c units, ratio between the oedometric compressibility in first e second loading cycle C r , and ratio between the aquifer and the aquitard permeability R k .The range of variability of these parameters is typical of a Quaternary sedimentary basin.
The results show that the delay between the pore pressure variation and land settlement can partly be accounted for by the delayed pressure propagation in the clayey layers confining the aquifer.The delay increases if: -the aquifer depth H a increases; -the sandy layer thickness S s decreases relative to the aquitard thickness; -the ratio between the oedometric compressibility in first and second loading cycle C r rises; -the ratio between the sandy and the clayey layers permeability R k decreases.
Within the investigated configurations the maximum subsidence occurring after the well shutdown amounts to 20 % and therefore must be taken into account for reliable subsidence evaluation and prediction.
Coupling this process with a visco-plastic mechanical behaviour of the depleted layers might cause an enhancement of the land subsidence delay.

Figure 1 .
Figure 1.Geomechanical model used in the numerical simulations.

Figure 2 .
Figure 2. Pressure change predicted by the groundwater flow model along the symmetrical axis in the case with H a = 1000 m, S c = 45 m, R k = 1 × 10 −6 and C r = 8.

Figure 3 .
Figure 3. Maximum vertical displacement U z normalized with respect to the displacement at the end of the pumping period (10 years) U z (10) in the case with H a = 500 m and S s = 10 m.In the left panel R k = 1 × 10 −6 and C r = 3, 5 and 8.In the right panel C r = 8 and R k = 1 × 10 −2 , 1 × 10 −4 and 1 × 10 −6 .

Figure 4 .
Figure 4. Maximum vertical displacement U z normalized with respect to the displacement at the end of the pumping period (10 years) U z (10) in the case with H a = 1000 m and S s = 10 m.In the left panel R k = 1 × 10 −6 and C r = 3, 5 and 8.In the right panel C r = 8 and R k = 1 × 10 −2 , 1 × 10 −4 and 1 × 10 −6 .

Figure 5 .
Figure 5. Maximum vertical displacement U z normalized with respect to the displacement at the end of the pumping period (10 years) U z (10) in the case with H a = 2000 m and S s = 10 m.In the left panel R k = 1 × 10 −6 and C r = 3, 5 and 8.In the right panel C r = 8 and R k = 1 × 10 −2 , 1 × 10 −4 and 1 × 10 −6 .

Figure 6 .
Figure 6.Maximum vertical displacement U z normalized with respect to the displacement at the end of the pumping period (10 years) U z (10) in the case with S s = 10 m C r = 8, R k = 1 × 10 −6 and H a = 500 m, 1000 m and 2000 m.