Hydrogeology Journal DOI 10.1007/s10040-013-1069-x

Characterization of land subsidence induced by groundwater withdrawals in the plain of Beijing city, China Youquan Zhang & Huili Gong & Zhaoqin Gu & Rong Wang & Xiaojuan Li & Wenji Zhao Abstract The plain of Beijing city in China suffers severe land subsidence owing to groundwater overdraft. The maximum subsidence rate could reach 6cm/year through the 2000s. An integrated subsidence-monitoring program was designed, including levelling survey, borehole extensometers and multilayer monitoring of groundwater level, with the aim to understand both hydrological and mechanical processes and to characterize the land subsidence. From multilayer compaction monitoring, the major compression layers were identified. The major strata contributing to compression deformation are the second (64.5–82.3m) and third (102–117m) aquitards, which contributed around 39% of the total subsidence. Meanwhile, irrecoverable deformations were also observed in the second (82.3–102m) and third (117–148m) confined aquifers; they exhibit elasto-plastic mechanical behavior, which is attributed to the thin beds of silt or silty clay. Stress–strain analysis and oedometer tests were conducted to study the aquifer-system response to pumping and to estimate the specific storage of the major hydrogeologic units. The results reveal the creep behavior and elastoplastic, visco-elasto-plastic mechanical behavior of the aquitards at different depths. The compressibility of the aquitards in the inelastic range is about one order of magnitude larger than for the elastic range. Keywords China . Subsidence . Groundwater management . Aquitard . Mechanical behavior

Received: 26 November 2012 / Accepted: 13 October 2013 * Springer-Verlag Berlin Heidelberg 2013 Y. Zhang ()) : H. Gong : Z. Gu : X. Li : W. Zhao College of Resource Environment and Tourism, Capital Normal University, Beijing 100048, China e-mail: [email protected] R. Wang Beijing Institute of Hydrogeology and Engineering Geology, Beijing 100195, China

Introduction Land subsidence caused by groundwater extraction has been a worldwide problem (Anderssohn et al. 2008; Galloway et al. 1999). Many areas that suffered subsidence have been identified and monitored, and qualitative and quantitative characterizations have been done. Previous researchers proposed that the hydromechanical behavior of an aquitard is one of the key points to understand the aquifer-system response to pumping and to mitigate subsidence (Galloway and Burbey 2011). The mechanical behavior of aquitards has been investigated and determined through laboratory testing, modeling and stress–strain analysis in situ (Burbey 2001; Gambolati et al. 2000; Helm 1975, 1976; Lofgen 1969; Men 1999a, b; Poland and Davis 1969; Terzaghi 1925). Riley (1970) identified the elastic expansion and plastic compaction from the stress–strain diagrams obtained in situ. Two orders of magnitude difference between the elastic compressibility and plastic compressibility have been found. Stress–strain curves also reflected the lag response between applied stress change and deformation within the aquitard (Riley 1984). Men (1999b) found the creep behavior of shallow, soft clay in Shanghai from laboratory tests. They proposed the creep behavior should be considered in calculating land subsidence. Based on simulation, Burbey (2001) evaluates what properties of an aquifer system may influence the stress–strain feature and determine which factors can influence the shape of hysteresis loops. Helm introduced poroviscosity theory to study the time-dependent mechanical behavior of land subsidence (Helm 1998). The compressible material was taken as a viscous fluid (Jackson et al. 2004) and the corresponding creep behavior can be better characterized based on this theory. Several analytic and numerical models of poroviscosity proposed by Li and Helm (Helm 1998; Li and Helm 1995, 1997) can be used to calculate and evaluate the 3D deformation due to groundwater extraction. The aquitard is the main compaction layer in most cases and the compressibility is two or more orders of magnitude greater than the aquifer. Thus, the compaction of the aquifer was usually ignored or underestimated in the study of aquifer-system deformation. However, some irrecoverable deformations have been confirmed within the sandy layers of aquifer systems, and in some

cases have notably contributed to total compaction both in laboratory testing and in situ (Liu et al. 2004; Shi et al. 2007; Zhang et al. 2007). The mechanical response of the aquitard and aquifer can help improve the assessment and prediction. The mechanical features of different strata are complex but important in assessing and formulating appropriate models of aquifer-system deformation. Beijing, located in the northwest of the North China Plain, is one of the cities that have undergone typical subsidence caused by groundwater extraction, especially in the plain area of the city. Annual groundwater withdrawals in the Beijing plain have consistently exceeded natural recharge since the 1970s (Jia et al. 2004). This over-extraction of groundwater has resulted in long-term depressurization of the aquifer system, regional decline of water levels, and local subsidence bowls. In order to meet the demands of urban development, some well fields and non-public water wells have been developed in the downtown areas, and this has resulted in the partial transfer of land subsidence to these new areas. The deformation in these subsidence bowls shows a different magnitude. The maximum subsidence reached up to 1.1 m during the period 1955–2007. A project was designed to monitor land subsidence, explain the mechanism, analyze impact on the high-speed railway, mitigate further subsidence and forestall fissure hazards in Beijing’s plain area. In order to meet the monitoring needs, borehole extensometers and multilayer monitoring wells have been used in the study area. The monitoring data reveal that the deformation is much more complex than previously realized. The purpose of this article is to identify the contribution of deformation of different soil layers and to analyze how the mechanical behavior of the aquifer system responded to the groundwater pumpage. This will be helpful for formulating appropriate models of aquifer-system deformation, and for supporting hazard mitigation measures in the future.

Study area Geographical location Beijing is located in the northwest of the North China Plain, which is one of the three highest-risk deltas with respect to subsidence in China (Syvitski et al. 2009). Figure 1 shows the geographical location of the study area: the plain area of Beijing city. This area is bounded by the Taihang Mountains in the west and by Jundu Mountains in the north and northeast. The south of the study area is adjacent to Hebei province and Tianjin city. Beijing’s plain area is about 6,390 km2 and it has more than 20 million inhabitants, living in the metropolitan areas. The study area is characterized by a continental monsoon climate, with an average annual precipitation of 588 mm/year (1956∼2002). Eighty percent of the precipitation is concentrated in the period mid-June to September. The average annual temperature is around 11.7 °C and Hydrogeology Journal

the maximum value may reach up to 42.6 °C in summer. The surface topography of the study area is relatively flat ground surface, with elevations ranging from around 60 m above sea level in the north to around 20 m in the southeast.

Hydrogeological setting The geology of the Beijing plain area is composed of water-bearing alluvial-pluvial and river channel deposits overlying bedrock. Tertiary and older sedimentary and volcanic rock units underlie Quaternary sediments and form the lateral and basal boundaries of the aquifer system. Relative to the alluvial-pluvial fill, basement rock is nearly impermeable. Locally Cambrian-Ordovician rocks underlie the unconsolidated sediments and have a hydraulic connection with the Quaternary aquifer systems. Five rivers—Chaobai River, Yongding River, Wenyu River,Ju River, Juma River—flow through the study area and have deposited a thick mass of Quaternary sediment of Pliocene through Holocene in age. Along the plain, several faults form the boundary between the graben and horst that are responsible for the structural framework of the basement (Jia et al. 2004). The geometry of this buried faulted basement structure provides a control on the total thickness of the unconsolidated deposits. The thickness increases from tens of meters in the mountain-front zone to several hundred meters in the central or southeast plain area. The sediment, forming the principal water-bearing units, consists of gravel, fine-to-coarse sand, silt, and clays (Jia et al. 2004). In general, the granularity of the aquifers tends to decrease from the distal fan to the proximal fan. The regional hydrostratigraphy, illustrated through representative boreholes, is shown in Fig. 2. Borehole log and cross sections show that the aquifers are composed of fineto-coarse sand and gravel in the proximal fan and fine sand, silt or clay in the mid and distal-fan areas (see Fig. 2). Figure 2 presents four Quaternary chronostratigraphic subdivisions: early Pleistocene (Q1); middle Pleistocene (Q2); late Pleistocene (Q3) and Holocene (Q4). The multilayered aquifer system may be separated into three aquifer groups, within depth intervals of 0–100, 100–200 and 200–300 m, corresponding to Q3-Q4, Q2 and Q1 geological eras, respectively. The lines presented in Fig. 2 were used to time stratigraphic contact and to help identify aquifer groups. Aquifer groups Q3-Q4 and Q2 are the major water-supplying units of the aquifer system and supply irrigation and industry-domestic water respectively (Xie et al. 2003). Aquifer group Q1 is locally distributed in the deepest basin deposit and scarcely exploited during the past decades. The three aquifer groups were further divided into five aquifers, which are separated by confining units of relatively lower permeability (see Fig. 2). This subdivision is based on lithology, consolidated condition, head and hydraulic characteristics from published descriptions (Xie et al. 2003) and from data collected during mapping as part of the study. The primary factor, hydraulic separation, DOI 10.1007/s10040-013-1069-x

Fig. 1 Location of study area: the plain of Beijing city. Dashed lines represent the boundaries of the alluvial fan (proximal-fan, mid-fan and distal-fan areas). Cross-section A–A′ is also shown in Fig. 2

is indicated by (1) lack of drawdown response in adjacent water-bearing zones during aquifer tests, and (2) differences in hydraulic head, and water-level fluctuations between individual aquifers. From the Fig. 2, the top-most confined aquifers (the first and second confined aquifers) are locally connected in the proximal fan. In general, the material of the first confined aquifer is less-compressible, coarse sand and fined-grained sand. The second confined aquifer between depths of 82.3 and 102 m is composed of fine-grain sand and little silty clay. The third confined aquifer, distributed over extensive areas in the middle to proximal fan, is the principal pumping layer for the groundwater demand of industrial and domestic use. This aquifer portion mainly consists of medium-tocoarse sand, silt and localized sand-silt mixtures. The fourth and fifth aquifers are locally distributed in the center of the thicker Quaternary deposit and have not been fully exploited in recent years. The hydrogeological cross section also illustrates the horizontal and vertical extent of the aquitards. The top-most aquitard (the first aquitard) is widely distributed over the plain area in the depth range of 35 to 49 m. This aquitard is composed mostly of clay and minor amounts of fine-tomedium sand. The second aquitard is distributed over extensive areas in the study region. This confining unit, Hydrogeology Journal

which consisted of silty clay, clay and minor amounts of silty fine sand, is located at about depth of 64.5 to 82.3 m. The third aquitard is chiefly composed of silty clay and clayey fine sand. It occurs at the depth interval 102–117 m and has a wide variety of thickness throughout the study area. The fourth aquitard is only locally distributed at the Quaternary deposit zone of north plain. This aquitard, located between 148 and 218 m deep, is chiefly formed by silt, silty clay and minor amounts of fine-to-medium sand interbeds.

Methodology Monitoring and measurement of land subsidence Land subsidence in the Beijing plain area has previously been monitored using a variety of techniques including leveling, extensometers and global positioning systems (GPS; Jia et al. 2004). The data used to map and analyze subsidence in this study were based on the conventional first-order leveling. The leveling network is now composed of 612 benchmarks, 23 major leveling lines and 40 closed loops (Jia et al. 2004). All the available benchmarks were scattered in the study area and formed an overall length of about 3,665 km. The first survey of the network was carried out in 1955. The network was resurveyed at first-order or second-order accuracy for a different DOI 10.1007/s10040-013-1069-x

Fig. 2 Hydrogeological cross-section A–A′ of the study area (the location is indicated in Fig. 1)

purpose from 1950s to now. The pffiffiffiffiaccuracy requirement for the leveling is around 1:8 K mm misclosure in any double run, where K is the distance between two neighboring benchmarks in km. All the results provided the basis for mapping subsidence and analyzing time-series evolution of land subsidence. The multiple extensometers and multilevel monitoring wells can provide valuable information about compaction, pore pressure and hydraulic head of the aquifer system, and make it convenient to characterize hydraulic and mechanical properties of the hydrogeologic units at different depth. A representative borehole-extensometer station, MLC3, was selected to collect data and analyze the deformation response of individual layers. Ten extensometers were respectively anchored in ten boreholes at different depths according to the stratigraphic variation. The extensometers have a minimum depth of 2.4 m and a maximum depth of about 238 m. Multi-level hydraulic head data were collected from the different monitoring wells at the same station. Long-term water levels were provided by the cooperator. Other original data were obtained from unpublished archived files and re-analyzed as part of this study.

Stress–strain analysis Riley (1970) developed a graphical method to analyze the deformation response to stress change of the aquifer Hydrogeology Journal

system by plotting applied stress versus strain (or hydraulic head versus deformation). Based on in-situ data, Riley showed that aquifer-system elastic and inelastic skeletal specific storages can be estimated from stress–strain diagrams. The skeletal specific storage (Ssk) is characterized numerically by strain per unit effective stress. It can be estimated graphically from the slope of the change in measured strain of the compacting interval (Δb/b) versus the monitored head change (Δh), and can be expressed as (Riley 1970): S sk ¼

Δb bΔh

Previous researchers have noted that head changes below the pre-consolidation stress (head) typically cause an elastic response of the compacting interval, which appears on strain–stress diagrams as hysteresis loops characterized by a similar slope (Epstein 1987; Hanson 1989; Pope and Burbey 2004; Sneed 2001; Sneed and Galloway 2000). The elastic skeletal specific storage (Sske) is the linear slope of the recovery limb of the hysteresis loops. Head changes that exceed the pre-consolidation stress (head) will cause an inelastic response of the compacting interval, which appears on strain–stress diagrams as a linear trend in multiple loops (Epstein DOI 10.1007/s10040-013-1069-x

1987; Hanson 1989; Riley 1970). The inelastic skeletal specific storage (Sskv) is the slope of the curve that connects the dominant hysteresis loops. To characterize the deformation-response pattern, stress–strain relations of individual stratigraphic units were analyzed during the study period 2004–2007.

Oedometer tests Previous researchers attributed the residual compaction to the delayed drainage and fluid-pressure equilibration of the low-permeability aquitards. In order to further analyze the mechanism of the lagging response, oedometer tests were performed to measure the deformation response to different loads. The secondary compressions for the soils were investigated by conducting a series of conventional one-dimensional (1D) consolidation tests for undisturbed soil samples taken from different locations and from a depth of 10–263 m. The specimens were tested in stainless steel rings, 5 cm diameter and 2 cm high. Each specimen was tested with 7day load increments during which the vertical deformation was measured with time. The specimen at different depths was tested to a different maximum vertical effective stress, which was determined by the pre-consolidation pressure of the specimen. Some time-compression curves obtained from oedometer tests could give a clear indication of an inflection point where the primary consolidation is assumed to end and the secondary consolidation is to start. Casagrande’s method was used to identify this inflection point on the settlement-log time plots. The idea behind the Casagrande’s logarithm of time fitting method is described in (Casagrande 1936). The test results of the second aquitard show an obvious creep behavior. Furthermore, the magnitude and characteristics of secondary compression of second aquitard were analyzed.

Results and discussion

2003). The aquifers below 200 m are not greatly exploited (Jia et al. 2004). Removal of groundwater was initiated in 1935, with the intensive withdrawal beginning in the 1950s (Bin 1964). Beginning in 1980, total withdrawal increased to around 2×109 m3/year. Groundwater overdraft has primarily occurred in the suburban area. Since 1990, relative stable pumping rates have been maintained at about 2.5× 109 m3/year (Xie et al. 2003). The net groundwater depletion, obtained by water-balance method, reached a peak in 1999 at 1.56×109 m3/year. However, since the late 1990s, the annual groundwater extraction has consistently exceeded the estimated natural recharge rate of around 1.53×109 m3/year. Currently, the groundwater volume pumped from the aquifer system is being depleted much faster than natural recharge. The reduced recharge is mainly caused by the declined in precipitation, especially after 1999 (see Fig. 3). The long-term negative effects of continued groundwater over-extraction have included the water-level decline of the aquifer system and regional subsidence. From 1970, prior to development, to 2003, the maximum water-level decline has reached up to 40 m in the north– northeast area of the Beijing plain area (Jia et al. 2004). A comparison of declining water level between 1970 and 2003, with the spatial distribution of public-pumping centers for the 1990–2000 period, was conducted (Fig. 4). The areas of maximum (30–40 m) water-level decline are concentrated in the central plain and spatially offset from the major public-pumping centers located in the north, northwest and west areas. Most of the public water wells are located in proximal-fan and mid-fan areas which consist mostly of coarse-grained deposits. The potentiometric surface slopes approximately from west to east in the western plain area and north to south in the northern plain area (Jia et al. 2004). One possible reason for the offset phenomenon is that public-pumping centers intercept groundwater flow, sustaining the water level in the down-gradient areas of the aquifer system. Another reason, as noted by Bin (1964), is that locally maximum water-level decline is related to the collective effect of extensive small non-public water wells in these areas

Relationship between the land subsidence and groundwater pumping The growth of the population and accelerated urbanization in Beijing has increased its reliance on groundwater resources as its water supply over the past four decades. As the available surface-water supplies have become insufficient to satisfy the developing demand, local government turned to the subsurface for additional water resources. Most of the groundwater supplies in the Beijing plain area come from a zone of unconfined and principally confined aquifers lying at depths of 0–200 m. The largest source of groundwater, the unconfined aquifer and the first confined aquifer at the maximum depth of 80–100 m, is used for irrigation. These two aquifers account for about 68 % of the total pumpage in the aquifer system. The section located at 100–200 m depth supply about 32 % groundwater for municipal and industry needs (Xie et al. Hydrogeology Journal

Fig. 3 Volume of natural recharge, total pumpage and net groundwater depletion in Beijing from 1995 to 2005 DOI 10.1007/s10040-013-1069-x

Fig. 4 Water-level decline of the second confined aquifer in the Beijing plain, 1970–2003. ‘0–10 m’ means ‘0 to less than 10 m’

which consist of fine-grained deposits with low permeability.

Overall land subsidence and its relationship with hydraulic head Over extraction in this region has caused regional decline of piezometric levels. Decreasing pore-fluid pressure may result in effective stress increases and compaction of the strata of soil (Terzaghi 1925). Contours of historical subsidence, measured from first-order level surveying, are developed for 1955–2007 (see Fig. 5). Limited or no deformation occurred in the proximal fan due to the strata formation and the corresponding poor compressibility. However, in the mid- and distal-fan areas, the strata consists mostly of fine sand, silty clay and clay for which permeability is relatively low and compressibility is high. Therefore, the overdraft or drawdown of water in this region easily causes notable subsidence compared to the proximal-fan area. Deformation maps reveal four principal subsidence bowls SC1, SC2, SC3 and SC4 in the study area. The maximum land subsidence of 1 m in SC3 and the minimum displacement of 0.75 m in SC1 have been measured. Subsidence in these bowls is concentrated in the area deposited with thick Quaternary sediments. The monitored results also show a good agreement between the deformation of the strata and the decline of regional piezometric levels. With development, accelerated and extended periods of land subsidence are associated with the increased use of groundwater in Hydrogeology Journal

alluvial aquifer systems. The comparison of land subsidence over time at four sites is displayed in Fig. 5. Deformation occurrence in benchmark sites BM3 and BM4 began in 1986, and continued through 1999 with stable subsidence rates of 2.5 and 2.2 cm/year respectively. BM1 and BM2 are both located in the eastern area of Beijing and subsidence since 1966 has been documented. Subsidence bowl SC2 was relatively stable with a subsidence rate of 1.5 cm/year during the period 1996–1999. However, subsidence of site BM1 appeared to occur relatively fast during 1966– 1982 and gradually slowed between 1983 and 1999 (see Fig. 6). Time series deformation was analyzed from measurements by level survey and the multi-layer groundwater levels at BM1 site. The comparison of subsidence and multi-layer head changes appears in Fig. 7. The water level at the eastern suburban area of Beijing continues to decline after 1966. Many of the industry wells were replaced or closed to enhance water management around 1982; however, the trends of water-level change for the period after 1982 generally declined faster than during the period before 1982. This difference may be attributed to the decreased storativity with time. Reduced subsidence rate in response to lower groundwater extraction has also been found in the eastern suburban area. From 1966 to 1982, the elevation in the eastern suburban area was reduced by 31 mm, for an average subsidence rate of 1.9 mm/year. At the same site, the subsidence rate decreased by 73 % from 1966–1982 to 1982–1999 (see DOI 10.1007/s10040-013-1069-x

Fig. 5 Subsidence contour map for the period 1955–2007. SC1, SC2, SC3 and SC4 represent subsidence bowls

Fig. 6). The subsidence rate is not as easily related to the trends of declining water level. The reasons for reduced subsidence rate during the period of relative rapid head decline in this area remain uncertain because of a lack of detailed analytical data regarding stress history, storage capacity and the vertical and horizontal distribution of pore pressure changes in the aquitards.

Vertical variation of deformation The contribution of individual hydrogeologic units to the total deformation was estimated based on the methods and

Fig. 6 Time series land subsidence for benchmark sites BM1, BM 2, BM 3 and BM 4, indicated in Fig. 5 Hydrogeology Journal

data as explained in the ‘Methodology’ section. At the MLC3 site (see Figs. 5 and 8), the accumulated subsidence was about 13.2 cm from April 2004 to July 2007. Among the total subsidence, 78 % compaction came from Holocene and late Pleistocene strata overlying the fourth confined aquifer; 39 % of the compression occurred in the late Pleistocene strata including the second aquitard and third aquitard; and 13 % of the compression came from Holocene top soil. It should be noted that the main pumping aquifer, the third confined aquifer, contributed 8 % compaction to the total displacement, which may be attributed to the compressibility of discontinuous silt interbeds within the aquifer.

Fig. 7 Time series groundwater level depth within confined aquifers and land subsidence at benchmark site BM1 DOI 10.1007/s10040-013-1069-x

Fig. 8 Contribution of compression for the different hydrogeologic units at MLC3

extraction for irrigation water. During routine operational pumping, this aquifer expands and compacts through cycles of water-level recovery and drawdown. As shown in Fig. 9a, compression of the first confined aquifer begins immediately after the water level begins to fall, and compaction virtually ceases soon after the water level starts to rise. The measured deformations were highly variable over time, reflecting the effects of seasonal fluctuations of the effective stresses applied in this aquifer (see Fig. 9b). The stress–strain curve presents annual loops corresponding to the loading and unloading cycles. The irrecoverable compaction is small and the displacement is typically small-magnitude generally recoverable elastic deformation.

First confined aquifer The annual cycles of groundwater level rise and fall, resulting from pumping and natural recharge, are identified in Fig. 9. Groundwater level fluctuates around 36 m during the study period. Figure 9a shows the compaction history between depths of 48.5 and 64.5 m, an interval that corresponds approximately to the zone of principal

Second confined aquifer The deformation and groundwater level of this aquifer show a good correlation; however, the deformation has a delayed response to the head change (see Fig. 10a). This delay is likely attributable to the composition of this aquifer which includes silt-clay interbeds or clay lenses in the sequence of deposit. Due to the relative lower vertical permeability of the clay interbed or the lens, compared to the sandy deposits, pore-water-pressure variation delays the head change in the adjacent sandy material. The residual compaction will occur even though heads in this

Fig. 9 a Deformation and groundwater level depth in the first confined aquifer (48.5–64.5 m) during the study period. b Graph of change in effective vertical stress vs. vertical strain from April 2004 to July 2007. Elastic skeletal specific storage, Sske, is the inverse slope of line A–B

Fig. 10 a Deformation and groundwater level depth in the second confined aquifer during the study period. b Graph of change in effective vertical stress vs. vertical strain from April 2004 to July 2007. Elastic skeletal specific storage, Sske, is inverse slope of line A–B. Inelastic skeletal specific storage, Sskv, is inverse slope of line C–D

Deformation of the aquifers

Hydrogeology Journal

DOI 10.1007/s10040-013-1069-x

aquifer have recovered. Irrecoverable, permanent deformation in the annual uploading/unloading cycle can was observed in Fig. 10. All these features indicate the elastoplastic response of the aquifer to the head change accompanying pumping. This aquifer produced about 5 % of the total compaction potentially attributable to the seasonal drawdown observed in this aquifer.

Third confined aquifer Figure 11a shows the history of compaction at depths between 117 and 148 m, an interval that corresponds approximately to the zone of principal pumpage from the third confined aquifer. The compaction and groundwater level of the third confined aquifer show a good correlation. Irrecoverable compaction, which corresponds to plastic deformation, was also observed in this aquifer. The compaction shown in Fig. 11 is the net change in thickness of a sequence of deposits, including the reaction of sandy layers that respond elastically, as well as the discontinuous silt interbed or clay lens in which the response is inelastic. The stress–strain curve shows elasto-plastic behavior of the third confined aquifer. Open hysteresis loops occurred during the study period, unlike in the second confined aquifer. A lack of closure of several loops suggests that interbeds in this aquifer may not be dissipating all excess pore-water pressure for these larger drawdowns before the next seasonal fluctuation resumes.

Fig. 11 a Deformation and groundwater level depth in the third confined aquifer during the study period. b Graph of change in effective vertical stress vs. vertical strain from April 2004 to July 2007. Elastic skeletal specific storage, Sske, is inverse slope of line A–B. Inelastic skeletal specific storage, Sskv, is inverse slope of line C–D Hydrogeology Journal

Fourth and fifth confined aquifers There are no corresponding multilayer monitoring data in this section; therefore, deformation features of the fourth and fifth aquifers are not analyzed in this report.

Deformation of the aquitards First aquitard In the underlying confined aquifer, water level decreased consistently during the period 2004–2007. The extensometer measurements reveal both elastic and plastic response. No obvious residual compaction occurred in the study period. Seasonal displacements reflect elastic seasonal deformation superimposed on irrecoverable plastic compaction attributed to the high compressibility of the clay layer. It can be seen from Fig. 12b that each major episode of stress increase was accompanied by additional permanent compaction.

Second aquitard The deformation is mainly caused by the groundwater extraction in the adjacent confined aquifer. Figure 13a shows the deformation in the second aquitard and the change in water level in the underlying aquifer. The deformation curve, which contained some recovery,

Fig. 12 a Deformation of the first aquitard and groundwater level depth in the first confined aquifer during the study period. b Graph of change in effective vertical stress vs. vertical strain from April 2004 to July 2007. Elastic skeletal specific storage, Sske, is inverse slope of line A–B. Inelastic skeletal specific storage, Sskv, is inverse slope of line C–D DOI 10.1007/s10040-013-1069-x

secondary settlement may lead to some errors in values for the second aquitard; thus, the creep deformation is another reason for the lag phenomenon in addition to the hydrodynamic consolidation of the second aquitard. Generally, the stress–strain diagrams and experimental results indicated that this aquitard has a visco-elastoplastic mechanical behavior and the visco-plastic compaction contributed to most of the deformation. This aquitard is the primary compaction layer which attributed about 20 % compaction for total subsidence. From 2004 through 2007, 2.6 cm of compaction has been recorded in the 64.5–82.3 m interval. The annual specific unit compaction rate in these deposits is about 0.5 mm of compaction per meter of thickness within the same period.

Fig. 13 a Deformation of the second aquitard and groundwater level depth in the second confined aquifer. b Graph of change in effective vertical stress vs. vertical strain from April 2004 to July 2007. Elastic skeletal specific storage, Sske, is inverse slope of line A–B. Inelastic skeletal specific storage, Sskv, is inverse slope of line C–D

Third aquitard The deformation in the third aquitard is chiefly caused by the head change of the third confined aquifer below, which is the main pumping aquifer. Time-series comparisons of compaction and water level exhibit a seasonal deformation signal superimposed on a long-term compaction trend (Fig. 16a). The hydraulic head of the confined aquifer decreased continually and caused the corresponding effective stress exceeding the pre-consolidation stress. Accordingly, the compaction is primarily inelastic, not fully recovering upon decrease in effective stress. The permanent compaction, present at the end of each major episode of unloading, is irreversible plastic deformation. The stress–strain curve in Fig. 16b shows that, despite water-level change associated with the seasonal pumping schedule and recharge, residual inelastic compaction is continuing at diminishing rates in this confining unit. The residual compaction results in small or no expansion through most or perhaps all of the periods of groundwater recovery, followed by accelerated compaction during the periods of decrease in hydraulic head. This feature can be explained in two ways. One is the hydrodynamic consolidation theory of soil mechanics (Terzaghi 1925) and the aquitard drainage model (Riley 1970). The other

suggests less compaction or relative uplift from June to March; there is accelerated compaction in the period midMarch to mid-June. As stresses decreased during the period of water-level recovery, compaction ceased, and a slight expansion of this aquitard can be observed. The large-magnitude deformation observed in the second aquitard is generally caused by inelastic compaction of the highly compressible clay layer (Fig. 13b). Despite the porefluid pressure fluctuation, within an established range in the annual cycles, each episode of stress increase was accompanied by permanent compaction. In addition to plastic compaction, it can be seen that the large-magnitude residual compaction occurred in the episode of stress decrease consistent with the absence of a closed loop in the yearly stress cycles. The laboratory tests reveal that the soil retains a normal consolidated state and has a creep behavior. Figure 14 shows the creep curves of the soil sample from the second aquitard.The time for the beginning of secondary consolidation, tp, was identified in Fig. 15 using the Cassagrande method (Casagrande 1936). The ratio of the secondary compression (creep) to total deformation is around 16 %. Moreover, despite this rather average value of ratio for the compressible layer, the secondary compression represents a significant part of the total subsidence because of the relatively short duration of the primary and secondary consolidation stages. As shown in Fig. 15, the compaction of the second aquitard shows a tendency to approach Fig. 14 Creep curves of the clay sample from the second aquitard stabilization during the 1-week testing period. Clearly, ignoring under pressure of 200/600 kPa Hydrogeology Journal

DOI 10.1007/s10040-013-1069-x

elastic deformation within this aquitard is relatively small compared to the permanent and residual compaction. The inverse slope of the stress–strain curve also indicates the high compressibility of this aquitard both in elastic and inelastic range compared to other confining units.

reason is the secondary compression or creep of the aquitard. The creep feature was confirmed by the oedometer test conducted on the samples acquired from the aquitard. The Casagrande method was used to determine the creep stage. The data reveal that the aquitard is normally consolidated soil and the creep deformation of this confining unit can approach 14% of total deformation. Collectively, these observations indicate that this aquitard is chiefly responsible for the plastic and viscous behavior. Similar to the second aquitard, the

Fourth aquitard The soil layer was considered to be overconsolidated. The magnitude of secondary compression is lower in overconsolidated clays than normally consolidated clays. The behavior of the soil mass is not greatly affected by creep. Figure 17a presents the water-level change and corresponding compaction within the fourth aquitard. The long-term trend of water level and compaction are similar; whereas, the time for highest and lowest recorded compaction typically has a few days delay compared to water level (see Fig. 17a). Stress–strain analysis indicates the deformation of this aquitard is elastic–plastic behavior. The permanent compaction, corresponding to the plastic deformation, contributed most of total compression within the aquitard. Residual compaction can be found during the episode of stress unloading. The compaction continued at diminishing rates through most or perhaps all of the period of head recovery and stress relief, even though the water level recovered with a relative large magnitude. In addition to the plastic deformation, a small-magnitude elastic response also can

Fig. 16 a Deformation variation of third aquitard and groundwater level depth in the third confined aquifer. b Graph of change in effective vertical stress vs. vertical strain from April 2004 to July 2007. Elastic skeletal specific storage, Sske, is inverse slope of line A–B. Inelastic skeletal specific storage, Sskv, is inverse slope of line C–D

Fig. 17 a Deformation variation of fourth aquitard and groundwater level depth in the fourth confined aquifer. b Graph of change in effective vertical stress vs. vertical strain from April 2004 to July 2007. Elastic skeletal specific storage, Sske, is inverse slope of line A–B. Inelastic skeletal specific storage, Sskv, is inverse slope of line C–D

Fig. 15 Time-compression curve derived from consolidation test

Hydrogeology Journal

DOI 10.1007/s10040-013-1069-x

Table 1 Summary of aquifer-system storage components and referred values from San Joaquin Valley S skv  m−1

First confined aquifer Second confined aquifer Third confined aquifer First aquitard Second aquitard Third aquitard Fourth aquitard San Joaquin Valley, California

−4

1.3×10 6.0×10−5 1.2×10−4 3.6×10−4 4.7×10−4 2.8×10−5 4.6×10−4 – 2.2×10−3

be found from the comparison in Fig. 17b. Generally, the compaction contributed by the plastic process is greater than two orders of magnitude of the elastic deformation.

Estimation of aquifer-system skeletal specific storage Aquifer-system skeletal specific storage can be estimated based on graphical methods (Riley 1970) as explained in the ‘Methodology’ section. Applying Riley’s method to the data from the extensometer MLC3 resulted in the graphs in Figs. 9b, 10b, 11b, 12b, 13b, 16b and 17b. The estimated elastic skeletal specific storage (Sske) and inelastic skeletal specific storage (Sskv) corresponding to the inverse slope of line A–B and line C–D were shown in these figures. Table 1 presented all the specific storage values and a comparison with San Joaquin Valley, California, USA. Elastic skeletal specific storage for aquitards fall within the ranges of values estimated from San Joaquin Valley (Ireland et al. 1980). The estimated inelastic skeletal-specific storage for the aquifer systems spans about an order of magnitude. Inelastic skeletal-specific storage for the first, second and third aquitards have the same order as the values of San Joaquin Valley. However, the ratios of elastic skeletal specific storage to inelastic skeletal specific storage values (around 4.7 %) for the second and third aquitards are smaller than the values (around 25.3 %) estimated from the San Joaquin Valley. The relative lower ratios further indicate that these two aquitards suffered predominantly inelastic compaction. Inelastic specific storage of the fourth aquitard is about one order of magnitude smaller than estimated values for the second and third aquitards. The smaller values of elastic and inelastic skeletal specific storage may also indicate that the degree of consolidation of this aquitard is high. Inelastic skeletal specific storage of the second and third confined aquifers can reach up to the magnitude of aquitards. This may be attributed to the composition of thinly bedded silt or silty clay in these two aquifers.

Conclusion Long-term extraction of groundwater has caused severe land subsidence over the extensive plain area in Beijing. Integrated subsidence-monitoring datasets, including levelling survey, borehole extensometers and multilayer monitoring of groundwater level, were collected and reanalyzed to understand aquifer-system response. Using stress–strain analyses and Hydrogeology Journal

S ske  m−1

1.7×10−5 7.4×10−6 5.6×10−6 1.7×10−5 1.7×10−5 2.2×10−5 5.1×10−6 5.9×10−6–2.2×10−5

oedometer tests, hydraulic and mechanical properties of the hydrogeologic units at different depths were characterized, which can help improve the assessment and prediction of land subsidence in the future.

1. The level survey revealed four primary subsidence bowls with different displacement rates from 1955 to 2007. The maximum subsidence rate approached up to 6 cm/year during the study period. Both the deformation and water-level change of the aquifer system show a good correlation. 2. From multilayer compaction monitoring, the major compression layers were identified. The major strata contributing to compression deformation are the second and third aquitards, which contributed around 39 % of total subsidence. Meanwhile, irrecoverable deformations were also observed in the second and third confined aquifers. It was revealed that the third confined aquifer contributes a notable compaction, 8 %, to the total subsidence. 3. Stress–strain analyses show that mechanical features of different strata are complex. Compared to the aquitards, the aquifers present poor compressibility and mostly respond to head change accompanying pumping. In contrast to previous findings, plastic effects were also found in second and third confined aquifers, which may be attributed to the thin interbed or clay lens within the water-bearing sand layer. Even though the mechanical behavior is the same, i.e. visco-elastic–plastic, the dominant effect may vary among the different aquitards depending on the strata formation and the stress history. The results of the oedometer tests show that the creep deformation occurring mainly in the normally consolidated soil and the visco behavior is obvious in the second aquitard. The ratio of the secondary compression (creep) to total deformation is around 16 %. Moreover, despite this rather average value of ratio for the compressible layer, the secondary compression represents a significant part of the total compaction because of the relatively short duration of the primary and secondary consolidation stages. 4. Meanwhile, aquifer-system specific storage was also estimated based on a graphical method (Riley 1970). The ratio of elastic skeletal specific storage to inelastic skeletal specific storage values of the second and third aquitards were around 4.7 %. The relative lower ratio and higher inelastic specific storage further indicate that the second and third aquitards suffered predominantly inelastic compaction. Inelastic DOI 10.1007/s10040-013-1069-x

skeletal specific storage of the fourth aquitard is about one order of magnitude smaller than estimated values for the second and third aquitards. The smaller values of elastic and inelastic skeletal specific storage may also indicate that the degree of consolidation of this aquitard is high. Acknowledgements The authors appreciate helpful review comments and suggestions on this article provided by Thomas Burbey, Devin Galloway and an anonymous reviewer, which greatly improved the quality of manuscript. Draft versions of the manuscript were reviewed by John Hess and Meijing Zhang. This research was partially supported by the Beijing Institute of Hydrogeology and Engineering Geology. The authors appreciate the assistance of Ye Chao, Liu Jiurong, Jia Sanman and Shen Yuanyuan who provided the water-level data and borehole extensometer data. This report is funded by the National Nature Science Foundation of China grant Nos. 41201376, 41130744, 41171335, Beijing Natural Science Foundation (No. 8133050, KZ201010028030), National Basic Research Program (973) of China (2006CB504400) and program (BJYRS-ZT-01-01).

References Anderssohn J, Wetzel H-U, Walter TR et al (2008) Land subsidence pattern controlled by old alpine basement faults in the Kashmar Valley, northeast Iran: results from InSAR and levelling. Geophys J Int 174:287–294. doi:10.1111/j.1365-246X.2008.03805.x Bin D (1964) Beijing hydrogeology survey report. First Hydrology Institute 0–96, Chinese Geological Survey, Beijing Burbey TJ (2001) Stress–strain analyses for aquifer-system characterization. Ground Water 39:128–136. doi:10.1111/ j.1745-6584.2001.tb00358.x Casagrande A (1936) The determination of the pre-consolidation load and its practical significance. Proceedings of the International Conference on Soil Mech and Found Eng (ICSMFE) Cambridge, MA, 22–26 June 1936, vol 3. Harvard University, Cambridge, MA, pp 60–64 Epstein VJ (1987) Hydrologic and geologic factors affecting land subsidence near Eloy, Arizona. US Geol Surv Water Resour Invest Rep 87–4143. http://pubs.usgs.gov/wri/1987/4143/ report.pdf. Accessed 7 July 2013 Galloway DL, Burbey TJ (2011) Review: regional land subsidence accompanying groundwater extraction. Hydrogeol J 19:1459– 1486. doi:10.1007/s10040-011-0775-5 Galloway D, Jones DR, Ingebritsen SE (1999) Land subsidence in the United States. US Geol Surv Circ 1182. http://pubs.usgs.gov/circ/ circ1182/pdf/circ1182_intro.pdf. Accessed 7 July 2013 Gambolati G, Teatini P, Baú D et al (2000) Importance of poroelastic coupling in dynamically active aquifers of the Po River Basin, Italy. Water Resour Res 36:2443–2459. doi:10.1029/2000wr900127 Hanson RT (1989) Aquifer-system compaction, Tucson Basin and Avra Valley, Arizona. US Geol Surv Water Resour Invest Rep 88–4172. http://pubs.usgs.gov/wri/1988/4172/report.pdf. Accessed 7 July 2013 Helm DC (1975) One-dimensional simulation of aquifer system compaction near Pixley, California: 1. constant parameters. Water Resour Res 11:465–478. doi:10.1029/WR011i003p00465 Helm DC (1976) One-dimensional simulation of aquifer system compaction near Pixley, California: 2. stress-dependent parameters. Water Resour Res 12:375–391. doi:10.1029/ WR012i003p00375 Helm D (1998) Poroviscosity in land subsidence case histories and current research: Proceedings of the Dr. J.F. Poland Symposium on Land Subsidence: Special Publication No. 8, Assoc. of Eng. Geologists, Washington, DC Ireland RL, Poland JF, Riley FS (1980) Land subsidence in the San Joaquin Valley, California, as of 1980. US Geol Surv Prof Pap Hydrogeology Journal

437-I. http://www.calwater.ca.gov/Admin_Record/C040893.pdf. Accessed 7 July 2013 Jackson J, Helm D, Brumley J (2004) The role of poroviscosity in evaluating land subsidence due to groundwater extraction from sedimentary basin sequences. Geofis Int 43(4):689–695 Jia SM, Guo M, Ye C et al (2004) Investigation of land subsidence in Beijing. Beijing Institute of Hydrogeology and Engineering Geology, Beijing, pp 0–212 Li J, Helm D (1995) A general formulation for saturated aquifer deformation under dynamic and viscous conditions. Proceedings of the Fifth International Symposium on Land Subsidence, Hague, Netherlands, October 1995. IAHS Publ. 234, pp 323– 332. http://ks360352.kimsufi.com/redbooks/a234/ iahs_234_0323.pdf. Accessed 7 July 2013 Li J, Helm D (1997) Numerical formulation of dynamic behavior within saturated soil characterized by elasto-viscous behavior with an application to Las Vegas Valley. In: Yuan J (ed) Computer method and advances in geomechanics. Proceedings of the 9th International Conference of the International Association for Computer Method and Advances in Geomechanics, vol 2. Wuhan, China, Balkema, Amsterdam, pp 911–916 Liu C-H, Pan Y-W, Liao J-J et al (2004) Characterization of land subsidence in the Choshui River alluvial fan, Taiwan. Environ Geol 45:1154–1166. doi:10.1007/s00254-004-0983-6 Lofgen BE (1969) Field measurement of aquifer-system compaction, San Joaquin Valley, California, USA. In: Tison LJ (ed) Land subsidence. Proceedings of the Tokyo Symposium, Sept 1969, IAHS Publ. 88, pp 272–284 http://iahs.info/redbooks/ a088/088031.pdf. Accessed 7 July 2013 Men F (1999a) Preliminary investigations on rheological properties of clay and ground settlement in Shanghai City (I) (in Chinese). J Nat Disasters 8:117–126 Men F (1999b) Preliminary investigations on rheological properties of clay and ground settlement in Shanghai City (II) (in Chinese). J Nat Disasters 8:123–132 Poland J, Davis G (1969) Land subsidence due to withdrawal of fluids. Rev Eng Geol 2:187–269 Pope JP, Burbey TJ (2004) Multiple-aquifer characterization from single borehole extensometer records. Ground Water 42:45–58 Riley FS (1970) Analysis of borehole extensometer data from central California. In: Tison LJ (ed) Land subsidence. Proceedings of the Tokyo Symposium, Sept 1969, IAHS Publ. 88, pp 423–431. http:// aquadoc.typepad.com/files/088031.pdf. Accessed 7 July 2013 Riley FS (1984) Developments in borehole extensometry. In: Johnson AI, Carbognin L, Ubertini L (eds) Land subsidence. Proceedings of the Third International Symposium on Land Subsidence, Venice, Italy, March 1984, IAHS Publ. 151, pp 169–186. http://iahs.info/redbooks/ a151/iahs_151_0169.pdf. Accessed 7 July 2013 Shi X, Xue Y, Wu J et al (2007) Characterization of regional land subsidence in Yangtze Delta, China: the example of Su-XiChang area and the city of Shanghai. Hydrogeol J 16:593–607. doi:10.1007/s10040-007-0237-2 Sneed M (2001) Hydraulic and mechanical properties affecting groundwater flow and aquifer-system compaction, San Joaquin Valley, California. US Geol Surv Open-File Rep 01–35. http:// pubs.usgs.gov/of/2001/ofr01-35/pdf/p1-4.pdf. Accessed 7 July 2013 Sneed M, Galloway DL (2000) Aquifer-system compaction: analyses and simulations—the Holly site. Edwards Air Force Base, Antelope Valley, California. US Geol Surv Water Resour Invest Rep 00–4015. http://pubs.usgs.gov/wri/2000/wri004015/ pdf/wrir004015.pdf. Accessed 7 July 2013 Syvitski JPM, Kettner AJ, Overeem I et al (2009) Sinking deltas due to human activities. Nat Geosci 2:681–686. doi:10.1038/ngeo629 Terzaghi K (1925) Settlement and consolidation of clay. Engineering News-Record, 26 Nov 1925, pp 874–878 Xie Z, Wang Y, Zhang S (2003) Investigation and evaluation of groundwater resource and water environment around the capital area. Beijing Institute of Geological Survey, Beijing, 244 pp Zhang Y, Xue Y-Q, Wu J-C et al (2007) Stress–strain measurements of deforming aquifer systems that underlie Shanghai, China. Environ Eng Geosci 8:217–228 DOI 10.1007/s10040-013-1069-x