Delayed Drainage of a Largely Deformed Aquitard due to Abrupt Water Head Decline in Adjacent Aquifer

A governing equation of drawdown was put forward to describe the one-dimensional large-strain consolidation behavior of an aquitard without consideration of the creeping effect. An analytical solution was derived to characterize the drawdown variation in the aquitard subjected to sudden hydraulic head decline in adjacent confined aquifer. The characteristics of the groundwater dynamics and water balance in the aquitard have been analyzed based on the analytical solution. A comparison analysis of results has been made between the large-strain theory and the classical small-strain theory. The type-curve fitting method was used to determine the hydrogeological parameters, on the basis of the observed variations of aquitard deformation with time. The analytical solution was thus validated by a comparison with the observed experimental results. It is found that the water drainage of aquitard is obviously delayed in response to the water head decline in the adjacent aquifer. All delayed water release from the aquitard terminates when the consolidation time reaches the value of l/cv0. The aquitard deformation predicted by the large-strain theory is less than that given by the small-strain theory, and the prediction discrepancy of these two theories increases with the increasing soil compressibility.


Introduction
A multilayered aquifer system usually consists of aquifers of relatively high permeability alternating with aquitards of relatively low permeability in between [1]. Aquitards may extend over large land areas throughout the earth, covering entire basins or alluvial plains, and sometimes they even extend into the sea, such as the Dakota aquifer system in the United States and the Yangtze Delta in China [2][3][4]. Abundant groundwater stored in the aquitards can be used for water supply when the aquifers are developed [5]. The role of aquitards in the subsurface groundwater system has been increasingly recognized in many respects. The groundwater dynamics in the aquitard is closely related to the groundwater exploitation, groundwater contamination, underground storage utilization, and land subsidence [6,7].
Aquitards consist of fine sediments such as clay and silt in low-energy depositional environments [8,9]. Aquitards, whose compressibility is often much larger than aquifers, are capable of releasing substantial quantities of water from the original storage in the process of dissipation of excess pore water pressure [10]. The laboratory measurement and regional groundwater flow modeling have demonstrated that the permeability of aquitard is not zero [8,11]. The hydraulic conductivity in an aquitard is usually less than 10 −8 m/s, and the permeability varies over a large range from 10 −15 to 10 −23 m 2 [12].
The understanding of the groundwater dynamics and water balance is of great importance to the prediction, evaluation, and control of land subsidence, and it is also meaningful to the groundwater resource evaluation and exploitation, as well as the analysis of contaminant and heat transport in the aquifer system [13][14][15]. Many studies have recently been carried out in terms of the groundwater flow and depletion in the aquitard. A deconvolution method was proposed [16], by which the drawdown in the confining layer can be calculated on the basis of the water head in neighboring aquifers. In particular, this water head may fluctuate in an arbitrary manner, which is not necessarily caused by a prescribed pumping regime [16]. As the water head varies in the confining layer, the water flow through it can be described by Darcy's law, which gives a direct estimation of water depletion. Konikow and Neuzil [3] presented a simplified method to estimate the groundwater depletion from the confining layers, and the computational errors were generally small. An analytical approach was proposed by Li et al. [5] to evaluate the depletion of groundwater in the confining layers throughout the entire exploitation history and also within limited periods, especially when there are few data on the drawdown history of the pumped aquifer. Alternatively, numerical models can be used to analyze the groundwater depletion in the confining layers by providing abundant information of the aquifer-confining layer system [17,18]. However, it should be pointed out that the aforementioned studies were mostly done based on Terzaghi's theory of one-dimensional consolidation for saturated clays.
Since the hydraulic properties of the confining layers of low permeability are usually not well understood, it is difficult to manage effectively the groundwater resources and have an accurate description of contaminant transport in the confined aquifer systems [19]. In general, the hydraulic conductivity and specific storage of an aquitard can be determined by the laboratory and in situ experiments [20][21][22]. For example, Burbey [23] used a graphical technique to determine the hydrogeological parameters by taking advantages of the time-subsidence data in a pumping test. In terms of this graphic technique, the straight-line portion of semi-log plot of time-compaction data represents the aquitard compaction. A type-curve fitting method was proposed by Zhou et al. [4] to calculate the hydraulic conductivity and specific storativity with the observed flux in response to the sudden decline of water head in the adjacent confined aquifer. Zhuang et al. [19] proposed a type-curve fitting method to estimate the hydraulic conductivity and specific storativity of an aquitard on the basis of the aquitard compaction variation and a piecewise linear function describing the drawdown history of the underlying aquifer. Konikow and Neuzil [3] presented a general relationship between the porosity n and the hydraulic conductivity by considering the clay content of the aquitard and a relationship between n and specific storativity by considering the consolidation degree. The parameters of aquitard were assumed to be constant during the consolidation process in the aforementioned studies.
Researchers have also investigated the groundwater dynamics and water balance in the small-strain consolidated aquitard [4,6]. For the large-strain consolidated aquitard (e.g., spatially the newly formed sedimentary soft soil), the variation of hydraulic properties is nonlinear during the consolidation of aquitard [24][25][26]. The small-strain consolidation theory is not suitable when the deformation of aquitard is larger than 10% of its thickness [27][28][29]. In this paper, an analytical solution was proposed to account for the drawdown variation in an aquitard subjected to largestrain consolidation, given that the water head in the adjacent confined aquifer suffers from a sudden decline. The characteristics of the groundwater dynamics, water balance, and aquitard deformation were accordingly analyzed. Additionally, a type-curve fitting method was presented for the determination of hydrogeological parameters, and the obtained parameters were validated by the experiments.

Mathematical Model and Analytical Solution.
A multilayered aquifer-confining layer system was considered in this study. The system consists of an aquitard of infinite lateral extent, which are bounded by two aquifers lying, respectively, above and below the aquitard. The initial thickness of aquitard is l. The aquitard is assumed to be homogeneous and saturated and to have a unified thickness (see Figure 1).
In the multilayered aquifer-confining layer system, the scalar equivalent hydraulic conductivity of each aquifer exceeds that of the aquitard by at least two orders of magnitude. The seepage in the aquitard follows the path of the least resistance; that is, water flows downward in the vertical direction and enters the aquifer with a sharp refraction at the interface between the aquitard and the aquifer. Therefore, it is reasonable to consider the water flow in the confining layer to be essentially vertical [16,30]. The Lagrangian coordinate a was used, which was assumed to be positive in vertically downward direction, with the coordinate origin located at the top surface of the aquitard (see Figure 1).
With the soil particles and pore water assumed to be incompressible, the equations for one-dimensional largestrain consolidation of saturated clay in the Lagrangian coordinate system are given as follows [31]: where σ is the total vertical stress; G s is the specific gravity of solid particles, and γ w is the unit weight of water; e = e a, t is the void ratio at time t; e 0 = e a, 0 is the initial void ratio.  Figure 1: Conceptual model of the multilayered aquifer-confining layer system.

Geofluids
where u is the total pore water pressure; u hs is the hydrostatic pore water pressure; u h is the long-term departure of pore water pressure from the hydrostatic pore water pressure due to a hydraulic gradient across the soil layer, and p is the excess pore water pressure. The equilibrium of fluid phase requires And the effective stress principle is given as where σ ′ is the effective vertical stress. Darcy's law is given as where v w and v s are the velocities of the pore water and soil particle relative to the datum plane; k v is the hydraulic conductivity of clay in the vertical direction. Equation (5) shows that the water flow velocity can be decomposed into two components, including one part arising from a long-term hydraulic gradient across the consolidating soil layer and the other part relating to the excess pore water pressure gradients. The governing equation of continuity of pore water flow is The governing equation of consolidation is just the combination of (5) and (6).
The nonlinear variations of the compressibility and permeability of soil in large-strain consolidation process are assumed to, respectively, abide by (8) [32] and (9) [26], and soil creeping here is not taken into account.
where c v0 = k v0 / m v1 γ w is the consolidation coefficient of aquitard, k v0 is the initial hydraulic conductivity of the aquitard, and m v and m v1 are the coefficients of compressibility of the aquitard at small-and large-strain states, respectively. It should be noted that the small-strain consolidation theory does not consider the variation of hydraulic properties of consolidating aquitard, and it is not suitable when the deformation of aquitard is larger than 10% of its thickness, while the large-strain consolidation theory considers the nonlinear variations of the compressibility and permeability of the aquitard. The relationship between void ratio e and total pore water pressure u is given as where σ 0 ′ is the initial effective vertical stress. The decrease of excess pore water pressure equals the increase of drawdown when the total stress is constant. Therefore, the relation of void ratio e with the drawdown s (i.e., the hydraulic head decline) of the aquitard is obtained as where S s = m v1 γ w is the specific storativity of the aquitard.
Substituting (8), (9), and (11) into (7) gives The water head is assumed not to vary in the aquifer above the aquitard, and it is assumed to decline by φ in the confining aquifer underlying the aquitard. Equation (12) is solved at the following initial and boundary conditions.
The analytical solution of drawdown variation with the large-strain consolidation theory is facilitated by the use of the variable transformation and separation (Appendix).  (17) gives The cumulative fluxes at the surface and bottom of the aquitard are, respectively, determined to be The settlement at the top surface of the aquitard is then given as Furthermore, the final settlement of the aquitard is identified by letting t ⟶ ∞ The cumulative fluxes at the top and bottom of the aquitard are given in a dimensionless form, respectively, as where t denotes the dimensionless time. According to the water balance principle, the water released from the entire thickness range of the aquitard is equal to the difference of the cumulative fluxes per unit horizontal area at the top and bottom of the aquitard.
The water balance results of the aquitard are shown in Figure 2. The curves of a = 0 and a = 1 represent the cumulative fluxes at the top and bottom surfaces of the aquitard, respectively. The curve Q t represents the cumulative water release of the aquitard. Three characteristic evolution stages are observed in Figure 2. At the initial stage (i.e., 0 < t < 0 1), the cumulative flux is zero at the cross section a = 0. The cumulative flux increases significantly at the cross section a = 1, owing to the cumulative water release of the aquitard. At the middle stage (i.e., 0 1 < t < 0 7), the cumulative water release of the aquitard increases slowly. The cumulative flux increases slightly, and the increase rate tends to be stable at the cross section a = 0. Also, the increase rate of cumulative flux tends to be steady at the cross section a = 1, indicating that the underlying aquifer is recharged not only by the released water of the aquitard but also by the leakage from the upper aquifer. At the last stage (i.e., t < 0 7), the cumulative flux at the cross section a = 0 is larger than the cumulative water release of the aquitard. The main supply source at the cross section a = 1 is the cumulative flux at the cross section a = 0. The curves of a = 0 and a = 1 are parallel to each other, indicating that the water release tends to be completed.

Identification of the Parameters of the Aquitard.
To determine the hydrogeological parameters of the aquitard while the water head declines by φ in the adjacent confined aquifer, we express (23) in a dimensionless form for simplicity It is assumed that the soil particles and pore water are incompressible. Therefore, the settlement at the top 4 Geofluids surface of the aquitard equals the cumulative water release of the aquitard. The logarithmic forms of (23) and (27) are, respectively, expressed as The standard curve of S 0, t is presented in Figure 3. The second terms at the right side of (30) and (31) are constant. Therefore, in the logarithmic plot, the test curve of the flux S 0, t is analogous to the type curve of the dimensionless flux S at a = 0 (see Figure 4 for reference). The curve of S 0, t is superposed onto the standard curve of S, and meanwhile, the axes of the two graphs are kept to be parallel to each other. There is a horizontal difference of c v /l 2 and a vertical difference of 1/S ∞ between these two curves. The matched point is chosen to be the intersection point of the standard curve with the test curve of S 0, t . The coordinates of matched point S, S, t, and t are then substituted into (32), (33), and (34) to determine the hydrogeological parameters.

Results and Discussion
The proposed analytical solutions in this paper for calculating the drawdown and depletion of the aquitard are evaluated by solving a large-strain consolidation problem. The drawdown of underlying confined aquifer is 10 m, and the water head of ground surface water (or unconfined aquifer) is assumed to be unchanged. A comparison analysis is also performed between the results obtained by the large-strain and small-strain consolidation theories. The large-strain theory herein refers to the analytical solution developed in this paper, and the small-strain consolidation theory is the Terzaghi's one-dimensional consolidation theory. The initial thickness of the aquitard is 10 m. The aquitard parameters used in the large-strain consolidation analysis are k v0 = 10 −9 m/s, S s = 0 01, 0.02, 0.03, 0.04 m −1 . The corresponding compressibility coefficients are 1, 2, 3, and 4 MPa −1 . In the small-strain analysis, the values of specific storativity and hydraulic conductivity are the same as those used in the large-strain consolidation analysis.

Geofluids
The spatial distributions of drawdown in the aquitard with different specific storativity values (S s = 0 01, 0.02, 0.03, 0.04 m −1 ) at t = 4 years were calculated using the largeand small-strain theories. The results in Figure 5 indicate that the drawdown in the aquitard increases with the decreasing specific storativity. The excess pore water pressure dissipates more quickly at a relatively small specific storativity. It is also observed that the drawdown obtained by the large-strain theory is less than that obtained by the small-strain theory. The dissipation of excess pore water pressure in the large-strain consolidation is slower than that in the small-strain consolidation. The difference of predicted drawdown between the large-and small-strain theories increases with the increasing specific storativity.
The drawdown at different positions (a = 1, 3, 5, 7, 9 m) in the aquitard is plotted in Figure 6. The drawdown increases initially and arrives at a stable value eventually. The variation rate of drawdown is large at the beginning and then gradually decreases to zero. It is also interesting to note that the drawdown increases with the increasing a value at a given time. The pore water releases more quickly at relatively high a value and specific storativity. The drawdown predicated by the large-strain consolidation theory is smaller than that obtained by the small-strain consolidation theory, and the prediction difference between the large-and smallstrain consolidation theories is insignificant. The void ratio decreases with the increasing effective stress in the aquitard. Meanwhile, the hydraulic conductivity and specific storativity of the aquitard decrease, which is ascribed to the decrease of void ratio.
The spatial distribution of water flow velocities in the aquitard is determined by (18) at different specific storativity values (S s = 0 01, 0.02, 0.03, 0.04 m −1 ) (see Figure 7). The results suggest that the water flow velocity in the upper part of the aquitard is lower than that in the lower part. The water flow velocity in the upper part is relatively small given a large specific storativity, but the case is opposite for the water flow velocity in the lower part. The water flow velocity predicated by the large-strain theory is smaller than that predicted by the small-strain theory, and the prediction difference increases with the increasing a value and specific storativity.
The variation of water flow velocity at different positions in the aquitard is described in Figure 8. The water flow velocity is zero at the beginning in the upper part of the aquitard and then increases to a steady value. For the lower part, the flow velocity increases initially and then decreases to the same steady value. Due to the decreasing hydraulic conductivity and specific storativity of the aquitard considered in the large-strain theory, the pore water flow velocity is lower than that obtained by the small-strain theory, and the predicted water flow velocities are the same when the water flow becomes stable in the aquitard.

Geofluids
The cumulative water release within unit area of the aquitard equals the settlement at the top of the aquitard, which is calculated by (23). As indicated in Figure 9, the depletion of the aquitard increases gradually to a steady value. The cumulative water release in unit area of the aquitard, which is given by the large-strain theory, is, respectively, 0.48, 0.91, 1.30, and 1.65 m 3 , while S s is, respectively, 0.01, 0.02, 0.03, and 0.04 m −1 . The depletion of the aquitard increases with the compressibility coefficient. In addition, the predicted cumulative water release of the large-strain theory is obviously less than that of the small-strain theory, and the prediction difference of aquitard depletion increases with the increasing specific storativity value.

Laboratory Apparatus.
A physical model test has been carried out to study the groundwater dynamics, water balance, and deformation characteristics of the aquitard in response to an abrupt water head decline in the adjacent confined aquifer. A schematic description of the testing device is given in Figure 10. It is composed of the main body of consolidation container, deformation monitoring device, flow monitoring device, and water supply device.
This apparatus can measure the deformation and flux of clay layer during testing. The consolidation container is a sealed organic glass cylinder with the height, outer, and inner diameters being 50 cm, 19 cm, and 18 cm, respectively. The soil specimen consists of three parts: the lower sand layer (4.5 cm in thickness), the middle layer filled with soil for testing (17.4 cm in thickness), and the upper sand layer (5.7 cm in thickness). The middle layer representing an aquitard is silty clay, and the basic physical and mechanical properties are given in Table 1.

Test
Procedures. The water head in the clay layer is initially constant. The water head in the lower sand layer declines suddenly by φ at the beginning and then remains unchanged during testing. We measured the drainage flux at the bottom of the aquifer and the subsidence of clay layer. The testing steps are given as follows.
(1) Initial inspection of the testing device. The consolidation container must be ensured to be sealed, and the inlet and outlet pipes must not be obstructed. The testing instrument must be ensured to operate normally (2) Soil sample preparation. The inner wall of consolidation container is cleaned and smeared with a layer of grease or Vaseline; hot water is injected into the consolidation container through the overflow tank and the water table is raised slowly; no bubble is allowed in the filter layer during the water injection process; the consolidation container is successively filled with sand layer, clay layer, and sand layer, and the device to measure the settlement is placed at the top of and below the clay layer after each soil layer is consolidated under its self-weight for a duration of 24 hours (3) Sample saturation. Turn off all the valves except the exhaust valve, and pump out the air in the container with a vacuum pump; increase the negative pressure within the consolidation container gradually at a given pressure sequence (e.g., -0.02, −0.04, −0.06, −0.08, and− 0.1 MPa); the time interval for maintaining each pressure level is 12 hours (4) After the soil layers are saturated, water is added slowly to reach the target elevation; open the valve of inlet pipe, and supply water into the device; mount the water supply tank at the designed height, and fill in it with water; discharge the excess water from the overflow nozzle; adjust the position of dial gauge, fasten it, and record the initial reading; observe the   Table 2. The deformation of clay layer was also measured by the dial gage.
The data of settlement of the clay layer obtained in the test are given in Table 3.
The measured values of S 0, t are plotted against time in Figure 4. The type curve of S 0, t is then superposed onto this figure, by keeping the axes of the two graphs parallel to each other (see Figure 4). The matched point is the intersection point of the type curve S 0, t with the experimental curve of S 0, t , and its coordinates in the two systems are S = 9 83 mm, S = 0 41, t = 76 min, and t = 0 033. The parameters of the clay layer can be obtained by substituting these coordinates into (32)-(34). The hydraulic conductivity k v0 , specific storativity S s , and hydraulic diffusivity c v0 are 4.25 × 10 −4 cm/min, 3.23 × 10 −3 cm −1 , and 0.1315 cm 2 /min, respectively.
The flux q l, t within the observation time at the bottom of the clay layer is obtained by substituting the hydrogeological parameters estimated by the type-curve fitting method into (18). The main advantage of the type-curve fitting method lies in its ability to make full use of all pumping test data and to improve the accuracy of calculation. Figure 11 compares the predicted and measured fluxes. It is seen that q l, t predicted by the analytical solution in this study agrees well with the experimental results.
Based on the one-dimensional Terzaghi's consolidation theory (i.e., the small-strain consolidation theory), the hydrogeological parameters of the aquitard can also be determined using the type-curve fitting method [33]. The hydraulic conductivity, specific storativity, and hydraulic diffusivity are determined to be k v = 3 63 × 10 −4 cm/min, S s = 2 76 × 10 −3 cm −1 , and c v = 0 1315 cm 2 /min. The hydraulic diffusivity estimated by the large-strain theory is the same as that given by the small-strain theory. The hydraulic conductivity and specific storativity determined   8 Geofluids by the large-strain theory are larger than those determined by the small-strain theory. The values of hydraulic conductivity and specific storativity estimated by the largestrain theory in fact refer to the hydraulic parameters of the clay layer at the beginning and they decrease with the decreasing void ratio. The compressibility coefficient at the bottom of the clay layer in the large-strain theory is a constant (=2.8 × 10 −3 cm −1 ) as estimated by (8). The hydraulic conductivity at the bottom of the clay layer is 3.2 × 10 −4 cm/min as estimated by (9), which is close to that obtained by Darcy's law (k v = 2 9 ×10 −4 cm/min). The equivalent hydraulic conductivity given by the largestrain theory is 3.68 × 10 −4 cm/min [26], which is close to that predicted by the small-strain theory.

Conclusion
The governing equation of drawdown was proposed for the description of the one-dimensional large-strain consolidation of saturated and homogeneous clays, with the creeping effect being ignored. An analytical solution accounting for the drawdown variation in a largely deformed aquitard has been derived, which is subjected to abrupt water head decline in the adjacent confined aquifer. Meaningful results have been obtained on the basis of the assumptions in (8) and (9) relating to the compressibility and permeability of soil. The characteristics of the groundwater dynamics and water balance of the aquitard have been examined with analytical solution. A type-curve fitting method is presented for determining the hydrogeological parameters. The main conclusions are given below: (1) When the water head in the confined aquifer declines by φ, the response of groundwater dynamics in the aquitard is found to be obviously delayed. The dissipation of excess pore water pressure predicted by   9 Geofluids the large-strain theory is slower than that predicted by the small-strain theory (2) An expression for computing the water flow velocity in the aquitard has been derived. Then, an analytical solution has been proposed to account for the cumulative water release and settlement of the aquitard. A water balance equation has also been established to describe the interrelations among the recharge, discharge, and water release of the aquitard. All delayed drainage of the aquitard is expected to be completed when the consolidation time t arrives at l 2 /c v0 . The soil compressibility is one important factor affecting the discrepancy of predicted deformations between the large-and small-strain consolidation theories. This discrepancy is found to increase with the increasing soil compressibility (3) The hydrogeological parameters of the aquitard (i.e., the hydraulic conductivity, specific storativity, and hydraulic diffusivity) are estimated by fitting the curves of analytically predicted and experimentally observed settlement. The hydraulic conductivity and specific storativity determined by the largestrain theory are higher than those given by the small-strain theory. In the large-strain theory, the hydraulic conductivity and specific storativity of the aquitard are considered to decrease with the decreasing void ratio. The hydrogeological parameters estimated by the small-strain theory are constants which are assumed to be on average equivalent to their true values in the consolidation process. The flux at the bottom of the clay layer, which is obtained from the analytical solution, is in good agreement with the experimental observation, implying that the estimated parameters are able to characterize the groundwater flow in the aquitard The hydrogeological parameters of the aquitard can be estimated by the type-curve fitting method, together with the data of flux per unit horizontal area at the top and bottom of the aquitard. The type-curve fitting method can also be used to measure in situ hydrogeological parameters.