A new type curve method for estimating aquitard hydraulic parameters in a multi-layered aquifer system

Journal of Hydrology 527 (2015) 212–220

Contents lists available at ScienceDirect

Journal of Hydrology journal homepage: www.elsevier.com/locate/jhydrol

A new type curve method for estimating aquitard hydraulic parameters in a multi-layered aquifer system Chao Zhuang a, Zhifang Zhou a,⇑, Hongbin Zhan b,c, Guangya Wang d a

School of Earth Science and Engineering, Hohai University, Nanjing 210098, PR China Department of Geology and Geophysics, Texas A&M University, College Station, TX 77843-3115, United States c School of Environmental Studies, China University of Geosciences, Wuhan, Hubei 430074, PR China d Geological Survey of Jiangsu Province, Nanjing 210018, PR China b

a r t i c l e

i n f o

Article history: Received 4 February 2015 Received in revised form 17 April 2015 Accepted 27 April 2015 Available online 5 May 2015 This manuscript was handled by Corrado Corradini, Editor-in-Chief, with the assistance of Adrian Deane Werner, Associate Editor Keywords: Type curve method Rate of aquitard compaction Drawdown rate of change Hydraulic parameters

s u m m a r y We propose a new type curve method for estimating the hydraulic conductivity (K) and the specific storage (Ss) of an aquitard in a multi-layered aquifer system. Theoretical type curves of rate of aquitard compaction versus time are obtained through drawdown-time series in the aquitard, under a condition that the drawdown history of the underlying aquifer can be approximated as piecewise linear segments (PLS), in each of which the drawdown shows a linear increase trend with time. Comparing with previous analytical solutions, which are limited by the assumption of a constant head in the overlying unconfined aquifer, the superposition principle can be called when water levels fluctuate with time in both upper and lower aquifers adjacent to the aquitard. When applying the new PLS method to estimate K and Ss values of an aquitard, the starting point of drawdown history at adjacent aquifers should be determined at first. After that the drawdown history is approximated as PLS to obtain drawdown rate of change for each linear stage and turning points of drawdown history. The new PLS method is demonstrated at a field site in Changzhou City, Jiangsu Province of China. Estimated results show that K and Ss are, respectively, 9.80  1010 m/s and 4.46  104 m1 for a thick aquitard, and 2.34–3.50  1010 m/s and 2.28  104 m1 for a thin aquitard. Estimated field K values, as well as Ss values, through the new PLS method, are generally in accordance with the results obtained from the geologic method by Konikow and Neuzil (2007). They are also in agreements with the numerically calibrated results of Shi et al. (2008a) at the regional scale and the laboratory experimental values of in-situ soil samples collected by Geological Survey of Jiangsu Province of China. Ó 2015 Elsevier B.V. All rights reserved.

1. Introduction A multi-layered aquifer system usually consists of multiple aquifers with alternating aquitards in between (Cihan et al., 2014; Zhan and Bian, 2006; Zhou et al., 2013). It is a rather common feature that can be often found beneath alluvial plains and sedimentary basins around the world. The permeability contrasts of alternating aquifers and aquitards are often over several orders of magnitude with the permeability values of most aquitards being less than 108 m/s (Neuzil, 1986). Consequently, those aquitards are often treated either as no-flow barriers or as leaking layers for cross-formation flow to occur, under rather idealized simplification.

⇑ Corresponding author. Tel./fax: +86 25 83787140. E-mail address: [email protected] (Z. Zhou). http://dx.doi.org/10.1016/j.jhydrol.2015.04.062 0022-1694/Ó 2015 Elsevier B.V. All rights reserved.

In recent decades, the importance of aquitards has been increasingly recognized in many aspects (Hart et al., 2006; Hendry and Wassenaar, 1999; Keller et al., 1989), including their impact on recharge and contaminant transport from overlying aquifers into deeper formations and in waste disposal (Neuzil, 1994) and in land subsidence research (Poland and Davis, 1969; Poland et al., 1975). Aquitards, whose compressibility values are often much larger than aquifers, are capable of releasing substantial quantities of water from storage in the process of dissipating excess pore water pressure (Leake, 1990). However, the accuracy of estimating quantities of water released from storage relies in a large degree on the accuracy of knowing the in-situ aquitard hydraulic parameters such as the hydraulic conductivity (K) and the specific storage (Ss), which are often not easy to obtain in-situ (Konikow and Neuzil, 2007). In an aquifer system, hydrodynamic response to groundwater pumping is highly dependent on the aquitard’s, as well as the aquifer’s hydraulic properties

C. Zhuang et al. / Journal of Hydrology 527 (2015) 212–220

(Keller et al., 1989; Rudolph and Frind, 1991; Smith et al., 2013). For better simulation of underground flow, we have to obtain reliable K and Ss values of aquitards at first. K and Ss of an aquitard can be determined by laboratory and in-situ experiments (van der Kamp, 2001; Yu et al., 2013). Conventional constant head and falling head permeameters are widely used in the laboratory for determining the aquitard K value, but the operation time is usually very long, obviously because of its very low permeability value (Remy, 1973). Consolidation experiments have also been used to determine K (Keller et al., 1989) and Ss (Shaver, 1998) values of aquitards. Unfortunately, laboratory experimental results are often unrepresentative of in-situ conditions (Smith et al., 2013; van der Kamp, 2001) due to the limitation of sample volumes and ubiquitous geological heterogeneity (Eaton, 2006; Marsily et al., 2005). The latter includes fractures (Keller et al., 1986; Rudolph and Frind, 1991; van der Kamp, 2001), dikes, sand streaks, interbeds, and so on (Meriano and Eyles, 2009). In addition, aquitard samples are inevitably subject to stress perturbation during the processes of collection, transportation, and laboratory installation (Clark, 1998). In-situ testing results are representative of what happens in actual field conditions. Field pumping tests can be used, such as Wolff (1970), Neuman and Witherspoon (1972), and Burbey (2003). The drawdown-time series in a pumping test may be interpreted using conventional well hydraulic theories, such as Hantush and Jacob (1955), Hantush (1960), Papadopulos and Cooper (1967), and Neuman and Witherspoon (1969a). Slug test is another commonly used field method for obtaining in-situ K values (Butler et al., 1996), but its results only represent the localized K values that are not suitable for large scale studies (Keller et al., 1989). If possible, a good practice recommended is to check both laboratory-based and field-based K and Ss values of aquitards to identify any scale dependent (or independent) trends, which have been a focus of hydrogeological investigations for many decades. If the aquifer water level changes with time in an arbitrary manner, one may need to use the deconvolution method (Neuman and Gardner, 1989) or the harmonic analysis method (Boldt-Leppin and Jim Hendry, 2003) based on Fourier transform to deal with time series of water level changes. Jiang et al. (2013) applied the harmonic analysis method using water-level signals in adjacent aquifers. In reality, both the deconvolution method and the harmonic analysis method reflect the aquitard drawdown responses to adjacent aquifer drawdown changes, and they can only determine the hydraulic diffusivity value (K/Ss) of an aquitard, which is similar to traditional field pumping tests. And laboratory consolidation experiments are usually operated to firstly determine the Ss value and then to calculate the K value based on the obtained hydraulic diffusivity value. Several studies introduced other methods to determine K and Ss of an aquitard. For example, Burbey (2003) used a graphical technique by taking advantages of time-subsidence data in a pumping test, in which the straight-line portion of the semi-log plot of time-compaction data represented aquitard compaction. Konikow and Neuzil (2007) demonstrated a general relationship between porosity (n) and K considering the clay content of aquitard and the relation between n and Ss considering the consolidation degree, which is called the geologic method here. Zhou et al. (2013) recently calculated K and Ss values using a type curve method, in which the aquifer drawdown was approximated by a step function in respect to time. Nearly all the analytical methods assumed Ss of the aquitard being a constant value, however, in reality, Ss is non-linear and should be evaluated by analyzing the past maximum effective stress of the unit. In a multi-layered aquifer system experiencing cyclical pumping, the drawdown may fluctuate from season to season and even decrease leading to rebound at the land surface. However, in many

213

long-term intensively exploited deep aquifers that were weakly affected by seasonal precipitation, the drawdown may be continuously increasing with time, such as the Su-Xi-Chang (SXC) area of China (Zhang et al., 2010). This special situation is what we exactly consider here. The objective of this study is to determine K and Ss of an aquitard, considering such a situation. The new method proposed in this study is demonstrated by a field application to the aquifer system beneath Qingliang Primary School at Changzhou City, Jiangsu Province of China. The findings of the method are compared with the geological method of Konikow and Neuzil (2007), and they are also compared with the numerically calibrated results of K by Shi et al. (2008a) at regional scale and with laboratory experimental K and Ss values of in-situ soil samples collected by Geological Survey of Jiangsu Province of China. 2. Theory and methods 2.1. Theory When multiple pumping wells are operating, drawdown at a fixed location is the superposition of drawdowns from all those pumping wells, each of which may have a different pumping intensity and history from others. This situation obviously cannot be handled by the conventional analytical solutions developed for a single pumping well, for instance by Hantush (1960) for confined systems or Neuman and Witherspoon (1969b) for unconfined systems. Sometimes, the drawdown may show a linear increase with time within a period of time. For instance, the water-table declined about 1 m/yr in the area of Mexico City, which partially contributed to the total subsidence of 7.5 m in 100 years at the city center (Custodio, 2002). The groundwater level had been decreasing at an average rate of 0.83 m/yr from 1978 to 1996 in Luancheng County, piedmont of Mt. Taihang in the North China Plain (Zhang et al., 2003), and the James River Lowland region of eastern South Dakota experienced head declines averaging about 7 m/yr between 1909 and 1915 (Konikow and Neuzil, 2007). We propose an aquifer system, as shown schematically in Fig. 1, with following assumptions: (1) The aquitard, with thickness, l, is of horizontally infinite extent and bounded by overlying and underlying aquifers. (2) The aquitard is homogeneous and horizontally isotropic. (3) The initial hydraulic head in the aquitard is uniform and the water level in the overlying undisturbed aquifer is constant. (4) The entire system remains saturated and the flow in the aquitard follows Darcy’s law. (5) The hydraulic conductivity of either aquifer is at least two orders of magnitude larger than that of the aquitard, for ensuring vertical flow direction in the aquitard. Neuman and Witherspoon (1969b) pointed out that the errors induced by assumption (5) were usually less than 5%. Sepúlveda (2008) found that there was almost no difference between above assumption (5) and actual three-dimensional flow in the aquitard.

Fig. 1. Schematic diagram of the studied multi-layered aquifer system.

214

C. Zhuang et al. / Journal of Hydrology 527 (2015) 212–220

Let the z-axis be vertical, positive upward, with its origin at the interface of the underlying aquifer and the aquitard (see Fig. 1). Drawdown history of the underlying aquifer is approximated as piecewise linear segments (PLS) composed of N stages, in each of which the drawdown shows a linear increase with time and the drawdown rate of change is defined as ui (>0) (1 6 i 6 N). Based on above assumptions, when time t is within the ith range (ti1  ti), the dimensionless flow governing equation of the aquitard can be given as: 2

@sD @ sD ¼ 2 ; @tD @zD

ð1Þ

subject to the following dimensionless initial and boundary conditions:

sD ðzD ; 0Þ ¼ 0ð0 6 zD 6 1Þ;

ð2Þ

sD ð1; tD Þ ¼ 0;

ð3Þ "

sD0 ¼ sD ð0; t D Þ ¼ bi tD  tD;1

# i1 X bi ai1  bj ðaj  aj1 Þ ;

ð4Þ

j¼1

where ai = tD,i/tD,1 or ai = ti/t1 (a0 = 0, a1 = 1), and the dimensionless terms are defined as: 2

sD ¼ Ks=uSs l ;

zD ¼ z=l;

2

tD ¼ Kt=Ss l

2

t D;i ¼ Kt i =Ss l ;

ð5Þ

where s, z, t and ti are respectively the aquitard drawdown, the vertical coordinate, the time, and the turning point at the junction between the ith and the (i + 1)th linear stage of the aquifer drawdown history, and the subscript ‘‘D’’ is the symbol of dimensionless variable. bi = ui/u is a positive constant, and u is a reference drawdown rate of change that numerically equals 1, thus bi is essentially the dimensionless or normalized ui and numerically equals the latter. The illustration of approximating the drawdown history as PLS in dimensionless form of Eq. (4) is shown in Fig. 2. If bi = bi+1 for all i values, the N stages of drawdown change reduce to one stage under a constant drawdown rate of change. When we approximate the drawdown history of underlying aquifer as PLS, aquifer drawdown sD is a continuous function in respect to tD, but its derivative with respect to tD is not continuous at the turning points tD,i. To overcome this problem, we have to derive the drawdown distribution in the aquitard sequentially, namely making the final drawdown distribution in the (i  1)th stage as the initial condition for the ith stage. Using the variable separation approach the analytical solution for the initial-boundary value problem of Eqs. (1)–(4) is derived as

( sD ðzD ; t D Þ ¼  p23

"

bi tD  t D;1

i 1 X X 1 bj n3

(

n¼1

j¼1

i1 X bi ai1  bj ðaj  aj1 Þ j¼1

#) ð1  zD Þ

;  ) exp n2 p2 /ði  jÞðtD  aj t D;1 Þ  2 2  sin npzD  exp n p ðtD  aj1 t D;1 Þ

ð6Þ

where the function /(x) is defined as

 /ðxÞ ¼

0; x ¼ 0 1; x – 0

ð7Þ

:

The details of derivation of Eq. (6) are listed in Appendix A. Eq. (6) implies that sD contains an infinite series of exponentially decayed terms with respect to tD. Such series precisely represent the excess pore water pressure in the aquitard, and they suggest that aquitard drawdown lags behind drawdown variation in the underlying aquifer, showing a delayed drainage in aquitard, because the hydraulic diffusivity (which equals K/Ss) of aquitard is much less than that of aquifer (Liu and Helm, 2008). When tD approaches infinity, each item of the series in Eq. (6) approaches zero and sD reaches its steady-state solution. Our exercise of Eq. (6) indicates that as long as the dimensional duration of the ith stage exceeds the hydraulic response time of the aquitard, which can be approximately obtained as t  Ssl2/K (Alley, 2002), the delayed drainage in aquitard at this particularly stage is almost non-existent. Since we assume ui being positive, the continuous increase in aquifer drawdown leads to the consolidation head in the aquitard being always less than previous levels due to time lag, and the compaction being continuous and inelastic throughout its history. Dissipating excess pore water is usually accompanied by aquitard compaction. And cumulative compaction is equal to its depletion from storage (Konikow and Neuzil, 2007). Usually, the measurement of aquitard compaction is not conducted from beginning, thus one is unable to know the cumulative compaction at a later time. However, determining the compaction rate by fitting the time-compaction data is straightforward, thus is illustrated in the field application part of this study. Based on the water balance equation of Zhou et al. (2013), the cumulative compaction is the integration of drawdown multiplying Ss over the entire aquitard thickness. The compaction rate, f, is the derivative of the cumulative compaction with respect to time. So the dimensionless compaction rate, fD, is f D ðt D Þ ¼ bi 

8 9 < exp½ð2n  1Þ2 p2 ðt D  aj1 t D;1 Þ = h i ; 2 2 2: ; p /ði  jÞðt  a t Þ /ði  jÞexp ð2n  1Þ ð2n  1Þ D j D;1 n¼1

i 1 X 8X bj 2

p

j¼1

1

ð8Þ

where fD is defined as:

f D ¼ 2f =uSs l:

ð9Þ

The detailed derivation of Eq. (8) is also included in Appendix A. Eq. (8) implies that fD also consists of an infinite series of exponentially decayed terms with respect to tD. It is easy to obtain the steady-state dimensionless compaction rate of the ith stage as f D jtD !1 ¼ bi , and its corresponding dimensional form as

Fig. 2. Illustration of approximating the drawdown history as piecewise linear segments in dimensionless form.

f jt!1 ¼ ui Ss l=2. Up to now, we only consider water level fluctuations in the underlying aquifer. In fact, water levels in overlying aquifer may also fluctuate with time. If this is true, one can call the principle of superposition and simply add the aquitard response to either aquifer drawdown fluctuation together. So the constant ui(=biu) used in above analysis essentially represents the sum of drawdown rate of change of the overlying and underlying aquifers. Eq. (8) shows that the relationship of fD versus tD depends on a series of bi or ai values. For convenience, we only take the situation of three linear stages for discussion (Fig. 3). It is easy to see the delayed aquitard compaction in Fig. 3, in which a2 = 2, a3 = 5,

215

C. Zhuang et al. / Journal of Hydrology 527 (2015) 212–220

101

101

0.1

100

0.15

α2=2, α3=5

0.05

0.2

fD

fD

α2=2, α3=5

0.05

100

tD,1=0.02

0.15

0.2

tD,1=0.02 βi=2 βi=1

βi=1

βi=0.5

(a)

βi=0.5 10-1 -2 10

0.1

10-1

100

10-1 -2 10

tD

(b) 10-1

100

tD

Fig. 3. Type curves of dimensionless compaction rate, fD, versus dimensionless time, tD, with a series of bi and ai values. (a) Dimensionless drawdown rate of change for each stage with b1 = 0.5, b2 = 1, b3 = 0.5, (b) Dimensionless drawdown rate of change for each stage with b1 = 0.5, b2 = 1, b3 = 2.

b1 = 0.5, and b2 = 1 for both Fig. 3a and 3b, and b3 = 0.5 for Fig. 3a and b3 = 2 for Fig. 3b. If bi < bi1, taking b2 = 1 and b3 = 0.5 for example, fD decreases rapidly with tD when it turns into the ith stage (see Fig. 3a); if bi > bi1, taking b2 = 1 and b3 = 2 for example, fD increases rapidly with tD when it turns into the ith stage (see Fig. 3b). If the duration of the ith stage is long enough, fD approaches bi eventually. The logarithmic forms of fD and tD derived from Eq. (5) and Eq. (9) are

log f D ðt D Þ ¼ log f ðtÞ þ log

log t D ¼ log t þ log

K 2

Ss l

:

2

uSs l

;

ð10Þ

ð11Þ

The second term on the right side of Eqs. (10) or (11) is independent of time. If plotted in a double-logarithmic plot, the measured time-dependent compaction rate curve of aquitard is similar in shape to one of the type curves shown in Fig. 3. And there is a shift of 2/uSsl along the horizontal axis and a shift of K/Ssl2 along the vertical axis. Therefore, the type curve method can be used to determine aquitard K and Ss values. On the measured curve and the matched type curve, coordinates f, fD, t, and tD of an chosen match-point are substituted into Eqs. (5) and (9) to estimate the hydraulic conductivity

K¼

2ltD f ; utfD

ð12Þ

and the specific storage

Ss ¼

2f

ulfD

:

ð13Þ

If the measured time-dependent compaction rate has reached the constant value in the ith stage, Eqs. (12) and (13) can be respectively simplified to

K ¼ 2lf =ui t;

ð14Þ

Ss ¼ 2f =ui l:

ð15Þ

When applying the type curve method to estimate K and Ss values of an aquitard, the starting point of aquifer drawdown history should be determined at first, then approximate the aquifer drawdown history as PLS to obtain the bi values and the turning points (ti) and corresponding ai values. Since the degree of time lag of drawdown

in the aquitard is strictly dependent on the drawdown history of the system, the starting point of the drawdown history must be reasonably determined or evaluated. If the early monitoring data of water level is not available, one should refer to the groundwater exploitation history to estimate the starting point. It is noteworthy that natural geologic formations are heterogeneous with complex patterns of spatial variability (Copty et al., 2008). In addition, the aquitard compaction process may result in substantial volumetric changes (Rudolph and Frind, 1991) and subsequent alternation in K and Ss values (Bredehoeft and Hanshaw, 1968; Rudolph and Frind, 1991). Therefore, the estimated results based on this new type curve method should be regarded as bulk average values over the entire aquitard thickness and through the entire compaction history. Fig. 1 implies that only one aquitard can be evaluated, however, it is obvious that the new PLS method can be used to evaluate multiple aquitards in a multi-layered aquifer system. 3. Field application 3.1. Background The proposed new method for estimating aquitard K and Ss values will be tested using a field site near Changzhou City, Jiangsu Province of China. The geographic location of the field site is shown in Fig. 4. This field site experienced substantial ground subsidence due to extensive groundwater pumping. Changzhou City is part of Yangtze Delta, and is located in southern Jiangsu Province with its north adjacent to Yangtze River and its southeast adjacent to Tai Lake. Regional ground elevation of the city is 3–6 m above mean sea level (m.s.l.). There is only one extensometer (F 0 0) in Suzhou, Wuxi and Changzhou cities in Jiangsu Province (SXC area) that has long-term observation data of ground subsidence. F00 was constructed at the center of regional groundwater depression cone in 1983 and it started service since 1984. Its ground elevation is about 4.4 m above m.s.l. 3.2. Hydrogeological units Changzhou City is underlain by a multi-layered aquifer system composed of Quaternary sediments. There are four main aquifers in this region, including an unconfined aquifer and three confined aquifers. Fig. 5 shows the stratigraphic distribution at the F00 site (without the unconfined aquifer). The first confined aquifer (CA1)

216

C. Zhuang et al. / Journal of Hydrology 527 (2015) 212–220

Fig. 4. Location of F00 site in Changzhou City, Jiangsu Province, China.

Depth interval (m)

Cross section

0.46 0.46

Specific Hydraulic storagae conductivity Hydrogeological Stratigraphy (Ss,m-1) (K,m/s) _ _ CA1 3.89E-09 Aquitard 1 (U) 1.40E-03 Aquitard 1 (U) 1.14E-03 4.46E-09

0.40

1.20E-03

3.21E-09

Aquitard 1 (M)

0.42

6.94E-03

1.00E-09

0.40

9.64E-03

1.00E-09

Aquitard 1 (M) Aquitard 1 (M)

0.42

6.40E-04

1.00E-09

Aquitard 1 (M)

0.42

1.85E-3

8.20E-10

Aquitard 1 (U)

Silty sand and fine sand

_

_

_

Clay, silty clay and silt

0.39

7.93E-04

1.33E-09

_

_

_

0.37

6.96E-04

1.00E-09

Soil types

18.55~28.55

Silty sand Silty clay Silty clay

28.55~35.70 35.70~43.95 43.95~50.45 50.45~61.20

Silty clay with interbedded siltstone Clay Silty clay and clay

61.20~71.32

Clay with interbeded silty clay

71.32~92.80

Soft plastic silty clay

5.48~18.55

Porosity (n) _

CA2

92.80~107.80 107.80~117.70 117.70~143.60

Silty-fine sand with interbedded clay

117.70~143.60

Hard plastic silty clay

Aquitard 2

CA3 Aquitard 3 (U)

Fig. 5. Stratigraphic description at F00 site. Data synthesized in this figure were collected by Geological Survey of Jiangsu Province of China. U denotes upper, M denotes middle, L denotes lower and read 1.40E03 as1.40  103.

is composed of silty sand and buried from 5.48 m to 18.55 m below ground surface. The second confined aquifer (CA2) is mainly composed of silty and fine sand and its burial depth is 92.80–107.80 m below ground surface. Due to its moderate burial depth, high water storage capacity and good water quality, CA2 is the major groundwater supply layer in the SXC area (Hu et al., 2009). The relatively high permeability of this aquifer, whose hydraulic conductivity is

5.8  105 m/s (Shi et al., 2007; Zhang et al., 2007), is another reason for choosing this one for water supply. The third confined aquifer (CA3) is buried from 117.65 m to 143.60 m below ground surface, and is mainly composed of silty fine sand with interbedded clay. Between the confined aquifers are (upper) aquitard 1 and (lower) aquitard 2, whose thicknesses are 74.25 m and 9.85 m, respectively. It is noteworthy that aquitard 1 is essentially a

217

C. Zhuang et al. / Journal of Hydrology 527 (2015) 212–220

10 0 -10 -20

Groundwater level, m

spatially heterogeneous layer at the F00 site. However, the drawdown variation in aquitard 1 in response to adjacent aquifer drawdown changes is unknown due to lack of piezometers in this aquitard, so this aquitard is simplified as homogeneous in this study for the sake of analytical treatment to determine the bulk average hydraulic parameter values over the entire thickness. We extend the porosity range of aquitard 1 in Fig. 5 to 0.35–0.5 for covering all potential porosity values of sub-layers in aquitard 1, considering its high heterogeneity. Primary sedimentary material of aquitard 1 is silty clay. Sample granularity analysis for aquitard 1 shows that grain sizes of 0.05–0.1 mm and less than 0.01 mm constitute 39.5% and 50.0% by weight, respectively (Shi et al., 2007). Aquitard 2 is composed of clay and silty clay with an average porosity of 0.39.

CA1 β11=0.7

-30 -40

β12=2.64

β23=2.62

-50

CA3

CA2 -60 -70

DRC of CA1 DRC of CA2 DRC of CA3

-80

β13=1

-90 1950 1955 1960 1965 1970 1975 1980 1985 1990 1995 2000

3.3. Groundwater exploitation and water level fluctuations in aquifers Groundwater exploitation in Changzhou City can be traced back to 1919, when the earliest deep well was developed. The original water levels for all aquifers were about 0.6 m below ground surface (Li et al., 2006). The unconfined aquifer is not developed yet and its water level remains almost unchanged over years. CA1 has a similar situation with the unconfined aquifer. Observation data show that water level of CA1 is greatly affected by precipitation and is in close connection with the unconfined aquifer. Being cut across by Yangtze River, CA1 could be recharged directly by Yangtze River. CA1 has not been exploited in SXC area due to its poor water quality, and its water level remains almost unchanged over years. Groundwater overexploitation of CA2 in Changzhou City began in late 1970s and continued into late 1990s. Before 1970, the average extraction rate was less than 105 m3/d. It was, respectively, 2  105 m3/d in 1976 and 2–2.36  105 m3/d during 1982–1989. From 1990 to 1997, it dropped to 1.50  105 m3/d. According to the observation data, water level variation in CA2 was not directly related to precipitation, but was closely correlated with groundwater extraction rates. The water level of CA2 has declined continuously for several decades due to overexploitation. In early 1950s, it was close to ground surface and the extraction rate was about 3.65  106 m3/yr. Since 1957, the number of deep pumping wells substantially increased, and the extraction rate reached 4  107 m3/yr in middle 1960s. Water levels at the center of the depression cone in CA2 were 20 m and 64 m below ground surface in 1967 and 1981, respectively (Wang et al., 2009), and they dropped over 62.3 m from 1960s to 1994 (Wu et al., 2009). The water level dropped to 83.3 m below ground surface in 1994 (Shi et al., 2007; Zhang et al., 2007). Groundwater exploitation in the SXC area was strictly regulated since 1995 and prohibited after 2000. At the F00 site, the water level of CA2 has already increased since 1994 (see Fig. 6), and the drawdown history of CA2 can be approximated into three continuous stages, namely 1950–1965, 1965–1983 and 1983–1994. The average drawdown rate of change in each stage is 0.7 m/yr, 2.64 m/yr and 1.0 m/yr, respectively (see Fig. 6). CA3 is deeply buried and distributed mostly in the northern Changzhou City. Groundwater recharge was quite difficult for CA3 which mainly received limited lateral recharge from other areas and a small amount of recharge through leakage from fractured deep bedrocks. Water level of CA3 was originally close to ground surface as well. During the groundwater exploitation periods of CA2, the water level of CA3 was actually higher than that of CA2. This was mainly because of the lower exploitation intensity of CA3. Although groundwater exploitation was limited since 1995 and prohibited after 2000, water level of CA3 did not rise until to 2000 (see Fig. 6).

Time, year Fig. 6. Groundwater level fluctuations in confined aquifers at F00 site. DRC denotes drawdown rate of change.

3.4. Aquitard compaction histories According to Terzaghi’s consolidation theory, when groundwater level declines, the drop in pore water pressure increases the effective stress on soil grains, leading to deformation of underground units. Land subsidence in the SXC area was closely related to groundwater exploitation. Land subsidence in Changzhou City was first observed in 1960s (Shi et al., 2008a; Wang et al., 2009), and developed rapidly since 1970s. At the F00 site, the total land subsidence during 1984–1994 was 464.55 mm and the cumulative compaction of aquitard 1 was 270.34 mm (see Fig. 7a), which is about 58.19% of the total subsidence; during 1984–2001, aquitard 1 contributed 57.15% of the total land subsidence (Shi et al., 2007; Wu et al., 2009, 2008; Zhang et al., 2007). The sub-units buried from 35.40 m to 92.80 m were the main consolidation units (Shi et al., 2007, 2008b). During this 11-year period, the compaction rate of aquitard 1 gradually declined from 30.97 mm/yr to 21.36 mm/yr (see Fig. 7b), reflecting the delay phenomenon of aquitard compaction. The cumulative compaction of aquitard 2 during the same 11-year period was 29.16 mm (see Fig. 7a). Aquitard 2 had an average compaction rate of 2.94 mm/yr during this period (see Fig. 7b). The compaction of aquitard 2 stopped after 1997, and it is speculated that the time for aquitard 2 to reach constant compaction rate was 2–3 yr. 3.5. Parameter estimation and verification For CA2, we only consider the drawdown history before 1994. We choose 1950 as the starting point, and 1965 and 1983 as two turning points of such a drawdown history. Since the water level of CA1 remains unchanged over years, the compaction of aquitard 1 is entirely caused by water level fluctuations in CA2. For such a case, based on the principle of the new PLS method, we have t1 = 15 yr, t2 = 33 yr, t3 = 44 yr, and further obtain a1 = 1, a2 = 2.2, a3 = 2.93. The average drawdown rates of change for three stages of the approximated PLS are 0.7 m/yr, 2.64 m/yr and 1.0 m/yr, respectively, leading to b11 ¼ 0:7, b12 ¼ 2:64 and b13 ¼ 1, where the superscript 1 represents aquitard 1. For CA3, the water level data before 1989 are scarce. However, the sum of drawdown rate of change of CA2 and CA3 for the third stage, interconnecting with aquitard 2, is 2.62 m/yr, so we have b23 ¼ 2:62, where the superscript 2 represents aquitard 2. From Fig. 8 we find that the measured compaction rate curve of aquitard 1 matches the type cure of tD,1 = 0.18. On the matched curve, we focus on the intermediate part of the measured data

218

C. Zhuang et al. / Journal of Hydrology 527 (2015) 212–220

400

35

Aquitard 1 Aquitard 2

30

300 250 200

R2 = 0.9980

150

a

R2 = 0.9998

100

Compaction rate, mm/yr

Cumulative compaction, mm

350

25 20

Aquitard 1 Aquitard 2

15

b

10 5

50 0 1984 1986 1988 1990 1992 1994 1996 1998 2000

0

34

35

36

37

38

39

40

41

42

43

44

Time, year

Time, year

Fig. 7. Time series of aquitard compaction at the F00 site. (a) Cumulative compaction. (b) Compaction rate.

1.0

1

10

2

0.1

0.22 0.18 0.14

Porosity, n (dimensionless)

10

0.26 0.3

0.05 10

f,mm/yr

fD

tD,1=0.02

β2=2.64

0

10

1

0.6

0.4

0.2

sil

s

d an

t

ys

ilt

β3=1 β1=0.7

10

0.8

Aquitard 1 (Laboratoy result) Aquitard 1 (This study) Aquitard 2 (Laboratory result) Aquitard 2 (This study) t Shi et al. (2008a) ten t Range of variation t on c en ten o n cont lay c c h lay clay lay hig c hc e y g r i rat me h ve e de so o on m tst sil

0.0 -16 10

α2=2.2, α3=2.93

-14

-12

10

-10

10

-8

10

-6

10

10

-4

10

Hydraulic conductivity, K (m/s)

-1

10 -2 10 0 10 0

10 -1

tD

10 1

10 0 10 2

t,year

Fig. 9. Comparison of estimated hydraulic conductivity values in this study with experimental results obtained from permeameter test and that of Shi et al. (2008a) in the geologic method calibrated from Konikow and Neuzil (2007). Laboratory K values are obtained through falling head permeameter tests.

Fig. 8. Observed data of aquitard 1 matched with the type curves. The solid lines represent type curves and the symbols represent the observed data.

1.0

0.8

Aquitard 1 (Laboratory result) Aquitard 1 (This study) Aquitard 2 (Laboratory result) Aquitard 2 (This study)

Range of variation 0.6

0.4

Water compre ssibility

Porosity, n (dimensionless)

for better match. A match point, the solid red dot in Fig. 8, with coordinates of tD = 0.49, t = 39 yr, fD = 1.6, f = 26.5 mm/yr, is selected. Substituting these values into Eqs. (12) and (13) yields the estimated hydraulic conductivity, K1 = 9.80  1010 m/s, and the specific storage, Ss1 = 4.46  104 m1, of aquitard 1. The estimated hydraulic conductivity of aquitard 2 is K2 = 2.34– 3.50  1010 m/s according to Eq. (14), and its estimated specific storage is Ss2 = 2.28  104 m1 according to Eq. (15) on the basis that the time for aquitard 2 to reach its constant compaction rate is 2–3 yr. Based on soil types listed in Fig. 5 and the granular constitution analysis mentioned above, the hydraulic conductivity (K) values of aquitard 1 and 2 are supposed to vary within the parallelogram outlined with dot-dash lines as shown in Fig. 9. Irreversible deformation (see Fig. 7a) suggests that these two aquitards are subject to normal consolidation and plastic deformation. Their specific storage (Ss) values are supposed to vary within irregularly closed dot-dash curves as shown in Fig. 10. After the model calibration by Shi et al. (2008a), the range of K values for aquitards in the SXC area and Shanghai City is 1095  1010 m/s. Fig. 9 shows that the parallelogram completely contains the range of Shi et al. (2008a) and basically contains the experimental results with exception at one point. The estimated

ted

da

li so

on

rc ve

O

d ate lly olid a rm cons No

0.2

0.0 -7 10

-6

10

-5

10

-4

-3

10

10

-2

10

Specific storage, Ss (m ) -1

Fig. 10. Comparison of estimated specific storage values with laboratory results and that of Shi et al. (2008a) in the geologic method calibrated from Konikow and Neuzil (2007). Ignoring the compressibility of water, experimental specific storage values are calculated as Ss = qgav/(1 + e0), where q is the density of water, g is the acceleration of gravity, e0 is the initial void ratio, and av is the compressibility of aquitard obtained from laboratory consolidation tests.

219

C. Zhuang et al. / Journal of Hydrology 527 (2015) 212–220

K values are five to six orders of magnitude smaller than those of confined aquifers, indicating that the flow in aquitards is essentially vertical. Fig. 10 shows that the estimated Ss results are within the closed dot-dash curves, and they are generally smaller than the experimental results, which is because the aquitards have experienced consolidation for a very long time and the samples are inevitably subject to stress perturbation due to many factors such as sample collection, transportation and laboratory installation. Above all, the estimated parameters of aquitard 1 are the bulk values of the entire thickness and through the entire deformation history. In order to obtain a better definition of a layered aquitard, one must observe drawdown fluctuations not only in adjacent aquifers but also in the aquitard itself, combining with soil compaction measurement of each sub-unit.

41372253) and the Program for Non-profit Industry Financial Program for MWR (Grant No. 201301083). The authors like to thank Zhaofeng Li, Zhi Dou, Zhou Chen, Qiaona Guo, Yong Huang and Hu Zheng for their suggestions on the manuscript. Appendix A Make the final drawdown distribution in the (i  1)th stage, which is defined as sD(zD, tD,i1), as the initial condition for the ith stage. It is noteworthy that at tD,0 = 0, sD(zD, 0) = 0, which is Eq. (2). For the ith stage, the governing equation, the upper and lower boundary conditions of flow in aquitard are Eqs. (1), (3) and (4), respectively, with tD,i1 6 tD 6 tD,i. For such a case, the analytical solution for drawdown in aquitard is easily derived using the variable separation approach as

4. Conclusions A new method for the drawdown variation in an aquitard was derived, under the condition that the drawdown history of adjacent aquifers could be approximated as PLS. Type curves of dimensionless rate of aquitard compaction versus dimensionless time were derived from developed analytical solutions and a new type curve method was proposed for estimating the vertical hydraulic conductivity and the specific storage of an aquitard. This new PLS method was demonstrated by a field application in a multi-layered aquifer system beneath Qingliang Primary School at Changzhou City, Jiangsu Province of China. The main conclusions are as follows:

" sD ðzD ; t D Þ ¼ bi tD  t D;1

# i1 X bi ai1  bj ðaj  aj1 Þ ð1  zD Þ j¼1

1 1 X X þ Gn;i ðt D Þ sin npzD þ Hn;i ðtD Þ sin npzD ; n¼1

where

Gn;i ðt D Þ ¼ 

 2bi  1  exp½n2 p2 ðt D  t D;i1 Þ ; 3 3 n p

 Z Hn;i ðtD Þ ¼ 2

1

sD ðzD ; t D;i1 Þ sinnpzD dzD 

0

(1) The aquitard drawdown response accompanied by the aquitard compaction lags behind water level variations in adjacent aquifers. When the duration of one of the linear stages exceeds the hydraulic response time of the aquitard, the delay phenomenon in this stage is virtually non-existent. (2) When applying the new PLS method to estimate K and Ss values of an aquitard, the starting point of the drawdown history of adjacent aquifers should be determined at first. After that, the principle of superposition is called to approximate the aquifer drawdown history as PLS to obtain the drawdown rate of each stage and the turning points. (3) The estimated K values of aquitards 1 and 2 at the selected field site are in general agreement with the results of the geologic method developed by Konikow and Neuzil (2007) and are in accordance with experimental results and within the numerically calibrated range of Shi et al. (2008a). The estimated Ss values of aquitards 1 and 2 generally agree with the results of the geologic method as well. But they are usually smaller than the experimental results, because of the long-term consolidation and stress perturbation during the processes of sample collection, transportation, and laboratory installation. (4) The estimated aquitard K and Ss results based on the new type curve method developed in this study should be regarded as the bulk average values over the entire aquitard thickness and through the entire compaction history. To obtain a better definition of a layered aquitard, long-term water level data should be observed not only in the aquifers but also in the aquitards, and should be integrated with compaction measurement of each sub-unit.

ðA1Þ

n¼1

ðA2Þ

 2 sD ð0; tD;i1 Þ np

exp½n2 p2 ðt D  t D;i1 Þ;

ðA3Þ

In order to obtain the complete analytical solution for tD in any linear stage, sD(zD, tD,i1) should be determined sequentially. For the first stage, Hn,1(tD) = 0, and the drawdown in aquitard for this stage is

sD ðzD ; t D;1 Þ ¼ b1 t D;1 ð1  zD Þ þ

1 X Gn;1 ðtD;1 Þ sin npzD :

ðA4Þ

n¼1

Substituting Eq. (A4) into Eq. (A3) yields

Hn;2 ðt D Þ ¼ 

 2b1  exp½n2 p2 ðtD  tD;1 Þ  exp½n2 p2 t D  : n3 p3

ðA5Þ

Substituting Eq. (A5) into Eq. (A1) yields the drawdown distribution for the second linear stage, and letting tD = tD,2 yields the drawdown distribution for the second stage

sD ðzD ; t D;2 Þ ¼ b2 ðtD;2  t D;1 Þ þ b1 tD;1 ð1  zD Þ þ

1 X Gn;2 ðtD;2 Þ n¼1

1 X  sin npzD þ Hn;2 ðtD;2 Þ sin npzD :

ðA6Þ

n¼1

In the same way, we can derived the analytical solution for drawdown distribution for the ith stage, namely Eq. (6). The aquitard compaction rate, f, is the first derivative of the cumulative compaction with respect to time t, yielding

f ¼

Ss

Rl

0

sðz; tÞdz : dt

ðA7Þ

Transforming Eq. (A7) into the dimensionless form as expressed in Eq. (9) yields Acknowledgements

fD ¼ This study was financially supported by the National Natural Science Foundation of China (Grant Nos. 41172204 and

2

R1 0

sD ðzD ; t D ÞdzD : dt D

Substituting Eq. (6) into Eq. (A8) yields Eq. (8).

ðA8Þ

220

C. Zhuang et al. / Journal of Hydrology 527 (2015) 212–220

References Alley, W.M., 2002. Flow and storage in groundwater systems. Science 296 (5575), 1985–1990. Boldt-Leppin, B.E.J., Jim Hendry, M., 2003. Application of harmonic analysis of water levels to determine vertical hydraulic conductivities in clay-rich aquitards. Groundwater 41 (4), 514–522. Bredehoeft, J.D., Hanshaw, B.B., 1968. On the maintenance of anomalous fluid pressures: I. Thick sedimentary sequences. Geol. Soc. Am. Bull. 79 (9). Burbey, T.J., 2003. Use of time-subsidence data during pumping to characterize specific storage and hydraulic conductivity of semi-confining units. J. Hydrol. 281 (1–2), 3–22. Butler, J.J., McElwee, C.D., Liu, W., 1996. Improving the quality of parameter estimates obtained from slug tests. Ground Water 34 (3), 480–490. Cihan, A., Zhou, Q., Birkholzer, J.T., Kraemer, S.R., 2014. Flow in horizontally anisotropic multilayered aquifer systems with leaky wells and aquitards. Water Resour. Res. 50 (1), 741–747. Clark, J.I., 1998. The settlement and bearing capacity of very large foundations on strong soils: 1996 R.M. Hardy keynote address. Can. Geotech. J. 35 (1), 131–145. Copty, N.K., Trinchero, P., Sanchez-Vila, X., Sarioglu, M.S., Findikakis, A.N., 2008. Influence of heterogeneity on the interpretation of pumping test data in leaky aquifers. Water Resour. Res. 44 (11), W11419. Custodio, E., 2002. Aquifer overexploitation: what does it mean? Hydrogeol. J. 10 (2), 254–277. Eaton, T.T., 2006. On the importance of geological heterogeneity for flow simulation. Sed. Geol. 184 (3–4), 187–201. Hantush, M.S., 1960. Modification of the theory of leaky aquifers. J. Geophys. Res. 65 (11), 3713–3725. Hantush, M.S., Jacob, C.E., 1955. Non-steady radial flow in an infinite leaky aquifer. Eos Trans. AGU 36 (1), 95–100. Hart, D.J., Bradbury, K.R., Feinstein, D.T., 2006. The vertical hydraulic conductivity of an aquitard at two spatial scales. Groundwater 44 (2), 201–211. Hendry, M.J., Wassenaar, L.I., 1999. Implications of the distribution of dD in pore waters for groundwater flow and the timing of geologic events in a thick aquitard system. Water Resour. Res. 35 (6), 1751–1760. Hu, J., Shi, B., Inyang, H.I., Chen, J., Sui, Z., 2009. Patterns of subsidence in the lower Yangtze Delta of China: the case of the Suzhou-Wuxi-Changzhou Region. Environ. Monit. Assess. 153 (1), 61–72. Jiang, Z., Mariethoz, G., Taulis, M., Cox, M., 2013. Determination of vertical hydraulic conductivity of aquitards in a multilayered leaky system using water-level signals in adjacent aquifers. J. Hydrol. 500, 170–182. Keller, C.K., Kamp, G.V.D., Cherry, J.A., 1986. Fracture permeability and groundwater flow in clayey till near Saskatoon, Saskatchewan. Can. Geotech. J. 23 (2), 229–240. Keller, C.K., Van Der Kamp, G., Cherry, J.A., 1989. A multiscale study of the permeability of a thick clayey till. Water Resour. Res. 25 (11), 2299–2317. Konikow, L.F., Neuzil, C.E., 2007. A method to estimate groundwater depletion from confining layers. Water Resour. Res. 43 (7), W07417. Leake, S.A., 1990. Interbed storage changes and compaction in models of regional groundwater flow. Water Resour. Res. 26 (9), 1939–1950. Li, C., Tang, X., Ma, T., 2006. Land subsidence caused by groundwater exploitation in the Hangzhou–Jiaxing–Huzhou Plain, China. Hydrogeol. J. 14 (8), 1652–1665. Liu, Y., Helm, D.C., 2008. Inverse procedure for calibrating parameters that control land subsidence caused by subsurface fluid withdrawal: 2. Field application. Water Resources Res. 44 (7), W07424. Marsily, G.D., Delay, F., Gonçalvès, J., Renard, P., Teles, V., Violette, S., 2005. Dealing with spatial heterogeneity. Hydrogeol. J. 13 (1), 161–183. Meriano, M., Eyles, N., 2009. Quantitative assessment of the hydraulic role of subglaciofluvial interbeds in promoting deposition of deformation till (Northern Till, Ontario). Quatern. Sci. Rev. 28 (7–8), 608–620. Neuman, S.P., Gardner, D.A., 1989. Determination of aquitard/aquiclude hydraulic properties from arbitrary water-level fluctuations by deconvolution. Ground Water 27 (1), 66–76. Neuman, S.P., Witherspoon, P.A., 1969a. Theory of flow in a confined two aquifer system. Water Resour. Res. 5 (4), 803–816.

Neuman, S.P., Witherspoon, P.A., 1969b. Applicability of current theories of flow in leaky aquifers. Water Resour. Res. 5 (4), 817–829. Neuman, S.P., Witherspoon, P.A., 1972. Field determination of the hydraulic properties of leaky multiple aquifer systems. Water Resour. Res. 8 (5), 1284– 1298. Neuzil, C.E., 1986. Groundwater flow in low-permeability environments. Water Resour. Res. 22 (8), 1163–1195. Neuzil, C.E., 1994. How permeable are clays and shales? Water Resour. Res. 30 (2), 145–150. Papadopulos, I.S., Cooper, H.H., 1967. Drawdown in a well of large diameter. Water Resour. Res. 3 (1), 241–244. Poland, J.F., Davis, G.H., 1969. Land subsidence due to withdrawal of fluids. Rev. Eng. Geol. 2, 187–270. Poland, J.F., Lofgren, B., Ireland, R., Pugh, R., 1975. Land subsidence in the San Joaquin Valley, California, as of 1972. US Government Printing Office. Remy, J.P., 1973. The measurement of small permeabilities in the laboratory. Geotechnique 23 (3), 454–458. Rudolph, D.L., Frind, E.O., 1991. Hydraulic response of highly compressible aquitards during consolidation. Water Resour. Res. 27 (1), 17–30. Sepúlveda, N., 2008. Three-dimensional flow in the storative semiconfining layers of a leaky aquifer. Groundwater 46 (1), 144–155. Shaver, R.B., 1998. The determination of glacial till specific storage in North Dakota. Groundwater 36 (4), 552–557. Shi, X.-Q., Xue, Y.-Q., Ye, S.-J., Wu, J.-C., Zhang, Y., Yu, J., 2007. Characterization of land subsidence induced by groundwater withdrawals in Su-Xi-Chang area. China Environ. Geol. 52 (1), 27–40. Shi, X., Wu, J., Ye, S., Zhang, Y., Xue, Y., Wei, Z., Li, Q., Yu, J., 2008a. Regional land subsidence simulation in Su-Xi-Chang area and Shanghai City, China. Eng. Geol. 100 (1–2), 27–42. Shi, X., Xue, Y., Wu, J., Ye, S., Zhang, Y., Wei, Z., Yu, J., 2008b. Characterization of regional land subsidence in Yangtze Delta, China: the example of Su-Xi-Chang area and the city of Shanghai. Hydrogeol. J. 16 (3), 593–607. Smith, L.A., van der Kamp, G., Jim Hendry, M., 2013. A new technique for obtaining high-resolution pore pressure records in thick claystone aquitards and its use to determine in situ compressibility. Water Resour. Res. 49 (2), 732–743. van der Kamp, G., 2001. Methods for determining the in situ hydraulic conductivity of shallow aquitards – an overview. Hydrogeol. J. 9 (1), 5–16. Wang, G.Y., You, G., Shi, B., Yu, J., Tuck, M., 2009. Long-term land subsidence and strata compression in Changzhou, China. Eng. Geol. 104 (1–2), 109–118. Wolff, R.G., 1970. Field and laboratory determination of the hydraulic diffusivity of a confining bed. Water Resour. Res. 6 (1), 194–203. Wu, J., Shi, X., Xue, Y., Zhang, Y., Wei, Z., Yu, J., 2008. The development and control of the land subsidence in the Yangtze Delta, China. Environ. Earth Sci. 55 (8), 1725–1735. Wu, J.-C., Shi, X.-Q., Ye, S.-J., Xue, Y.-Q., Zhang, Y., Yu, J., 2009. Numerical simulation of land subsidence induced by groundwater overexploitation in Su-Xi-Chang area, China. Environ. Earth Sci. 57 (6), 1409–1421. Yu, L., Rogiers, B., Gedeon, M., Marivoet, J., Craen, M.D., Mallants, D., 2013. A critical review of laboratory and in-situ hydraulic conductivity measurements for the Boom Clay in Belgium. Appl. Clay Sci. 75–76, 1–12. Zhan, H., Bian, A., 2006. A method of calculating pumping induced leakage. J. Hydrol. 328 (3–4), 659–667. Zhang, X., Pei, D., Hu, C., 2003. Conserving groundwater for irrigation in the North China Plain. Irrig. Sci. 21 (4), 159–166. Zhang, Y., Xue, Y.-Q., Wu, J.-C., Ye, S.-J., Wei, Z.-X., Li, Q.-F., Yu, J., 2007. Characteristics of aquifer system deformation in the Southern Yangtse Delta, China. Eng. Geol. 90 (3–4), 160–173. Zhang, Y., Xue, Y.-Q., Wu, J.-C., Shi, X.-Q., Yu, J., 2010. Excessive groundwater withdrawal and resultant land subsidence in the Su-Xi-Chang area, China. Environ. Earth Sci. 61 (6), 1135–1143. Zhou, Z., Guo, Q., Dou, Z., 2013. Delayed drainage of aquitard in response to sudden change in groundwater level in adjacent confined aquifer: analytical and experimental studies. Chin. Sci. Bull. 58 (25), 3060–3069.