Environmental Management (2013) 51:1109–1125 DOI 10.1007/s00267-013-0037-5

RESEARCH

Groundwater Exploitation Management Under Land Subsidence Constraint: Empirical Evidence from the Hangzhou–Jiaxing–Huzhou Plain, China Guoliang Cao • Dongmei Han • Jessa Moser

Received: 4 April 2012 / Accepted: 25 February 2013 / Published online: 20 April 2013 Ó Springer Science+Business Media New York 2013

Abstract Land subsidence caused by extensive groundwater pumping has become a factor which cannot be ignored in the sustainable exploitation of groundwater resources. The Hangzhou–Jiaxing–Huzhou Plain is one of the locations with China’s most severe land subsidence problems; the region has experienced dramatic land subsidence since the 1960s. Historical records of groundwater extraction, hydraulic head, and land subsidence show the latter to be the result of continual and excessive extraction of groundwater from deep confined aquifers. This study reconstructs land subsidence using an integrated regional groundwater flow and land subsidence model. The model is calibrated using land subsidence and groundwater level measurements from 1996 to 2007. Simulation results reproduce the cones of depression for groundwater heads and nadirs of land subsidence reasonably well. The calibrated model is used to evaluate the efficacy of land subsidence prevention plans from 2008 to 2010, and to predict future land subsidence over the next decade considering several groundwater exploitation scenarios. The results show the main cause of land subsidence to be inelastic compaction of the aquifer system resulting from G. Cao Center for Water Research, College of Engineering, Peking University, Beijing 100871, China G. Cao  J. Moser Department of Geological Sciences, University of Alabama, Tuscaloosa, AL 35487, USA D. Han (&) Key Laboratory of Water Cycle & Related Land Surface Processes, Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Sciences, Beijing 100101, China e-mail: [email protected]

continuously declining water levels. The model reveals that while the area of land subsidence will continue to extend, the rate of this extension may be significantly decreased by reduction of groundwater extraction. If the current land subsidence prevention and reclamation plans are continued and surface water diversion projects implemented, though land subsidence cannot be halted, the rate at which it is occurring can be effectively reduced. Keywords Regional land subsidence  Groundwater over-exploitation  Inelastic compaction  Numerical model  Groundwater withdrawal reduction

Introduction Growing populations, urbanization, and rising demands from industrial and agricultural sectors have rendered land subsidence caused by excessive groundwater exploitation a concern of global significance (e.g., Poland and Davis 1969; Poland 1984; Yamamoto 1995; Galloway and others 1999; Hu and others 2004; Teatini and others 2005). As such, land subsidence mitigation has become an important consideration in the management of groundwater resources. Typical methods of subsidence prevention include control of groundwater withdrawal and groundwater artificial recharge (Galloway and Burbey 2011). Numerical models have proven beneficial in evaluating the evolution of land subsidence due to groundwater pumping, and are at present the most powerful predictive tool in assessing potential future land subsidence developments (Calderhead and others 2011). Regulatory efforts to manage groundwater extraction have been shown to effectively recover lowered groundwater head and retard the development of land subsidence

123

1110

(Daito and others 1991; Gallardo and others 2009). However, feasible groundwater exploitation strategies on a basin scale are complicated by local differences in groundwater demands from domestic, industrial, and agricultural sectors across regions. The issue of land subsidence in susceptible groundwater basins where groundwater demand is high will inevitably persist unless the optimal groundwater scenarios are designed to adequately represent social, economic, and hydrogeologic conditions. To evaluate the ability of potential groundwater management alternatives to achieve the goal of preventing further land subsidence, groundwater flow models coupled with aquifer-system compaction models have been proven useful (e.g., Onta and Gupta 1995; Xu and van der Gun 1995; Larson and others 2001; Teatini and others 2006; Ortiz-Zamora and Ortega-Guerrero 2010). Although optimization algorithms have been proposed for management of groundwater pumping explicitly considering land subsidence in hypothetical groundwater pumping problems (Chang and others 2007), these are not adequate in practice as they fail to address issues such as the complexity of the aquifer system, the variation of spatial distributions in groundwater pumping, and the impact of potential water management strategies. The Hangzhou–Jiaxing–Huzhou (HJH) Plain, as a part of the Southern region of the Yangtze River Delta, refers to the plain around the cities of Hangzhou, Jiaxing, and Huzhou. This area has suffered widespread surface water contamination since the 1980s. The primary cause of contamination is uncontrolled industrial and domestic sewage disposal, which has rendered approximately 90 % of the surface water in the HJH Plain seriously contaminated and undrinkable (Li and others 2009). This surface water contamination situation along with the fact that surface water supply networks have been limited to cities and suburban areas (Zhang 2005) has led to increased use of groundwater as a secondary water supply source to sustain the rapid economic development in this region over the last three decades. Consequently, the HJH Plain has become one of the three most subsiding regions in China (Xu and others 2008). Land subsidence caused by extensive groundwater withdrawal in the HJH Plain has been investigated in previous studies (Zhao and others 2004, 2006; Li and others 2006). Negative effects include flooding associated with lowering land surfaces, complications in the discharging of pollutants in urban areas, river navigation problems, and waterlogging of soil (Zhao and others 2003, 2004). Historical records of groundwater exploitation, hydraulic head, and land subsidence in the HJH Plain show that the primary cause of subsidence is continuous overexploitation of groundwater from deep, confined aquifers.

123

Environmental Management (2013) 51:1109–1125

This necessitates that land subsidence be a critical consideration in determining regulatory constraints for groundwater pumping. It also underscores the importance of the development of a regional groundwater flow model coupled with a land subsidence model, which will serve as a much needed aid in the management of water resources. However, no previous modeling work has been done to reconstruct or predict the occurrence of land subsidence in the HJH Plain. The main purposes of this study are: (1) to develop a coupled regional groundwater flow and land subsidence model that can be used to reconstruct groundwater level changes in a confined aquifer system, as well as the subsequent compaction of the aquifer and aquitard matrices; (2) to evaluate the efficacy of current groundwater regulations in preventing land subsidence in the HJH Plain; (3) to assess future environmental impacts of current groundwater regulations, along with several alternative groundwater management scenarios, over the next decade.

Geological and Hydrological Setting Study Site The HJH Plain is one of the most developed economic regions in China (Jean and Tang 2008), as well as an important center for rice production (half of the farmland in the area is cultivated with rice). Geographically, it is located in the southern part of the Yangtze River Delta, including the northern half of the Zhejiang province in between the Yangtze and Qiantang rivers, south of the Shanghai city and Jiangsu province (Fig. 1). The total area of the HJH Plain is approximately 6,500 km2. The topography is relatively flat, gently sloping from 3–5 m above mean sea level (a.s.l.) in the south, to 1–2 m (a.s.l.) in the north, where it abuts the Taihu Lake. The average elevation is 2–4 m (a.s.l.). The climate in the HJH Plain is subtropical monsoon with a mean annual temperature range of 15–17 °C and precipitation of *1,200 mm, with significant seasonal variations (higher in summer and autumn, and lower in winter and spring) (Jiang and others 1998). The averaged annual pan evaporation is *900 mm (Jiang and others 1998). The HJH Plain is characterized by a dense natural surface water network, which accounts for *10 % of the total area of the plain. Surface water from these rivers serves as the primary source of water in the HJH Plain (with groundwater as the second most important water source). However, due to rapid urbanization and industrial development, these rivers have been severely contaminated by direct sewage disposal. In 1990, approximately 50 % of the total surface water available in the HJH Plain met China’s

Environmental Management (2013) 51:1109–1125

1111

Fig. 1 Location of the Hang-Jia-Hu Plain. Cross-section line (A–A0 ) shown in Fig. 2

national water quality standard for drinking water. However, currently, there is no river or lake in this region can serve as water supply for domestic or agricultural use because of failure to meet national water quality standards. Water shortage due to reduced water quality has caused the HJH Plain to rely more heavily on groundwater to meet its rising water demands. Hydrogeological Settings The HJH Plain is predominantly comprised of Quaternary sediments ranging in thickness from approximately 50 m in the southwest to more than 300 m in the northeast (Jiang and others 1998). Geological and hydrogeologic surveys revealed a multi-layered aquifer system beneath the plain (Jiang and others 1998; Liu and others 2006); the Quaternary deposits in the plain may be divided into a shallow unconfined aquifer underlain by three confined aquifers (hereafter referred to as I, II, and III in descending order) (Fig. 2). The corresponding geologic units of the aquifers are a Holocene series (Q4), late Pleistocene series (Q3), middle Pleistocene series (Q2) and an early Pleistocene series (Q1). Interbedded aquitards between the confined aquifers are called as aquitard 1 and aquitard 2. The geometry of the uppermost late Pleistocene confined aquifer (confined aquifer I) is controlled by three paleochannels spreading toward the northeast. It is composed of medium to fine-grained sand and sand with

gravel in the upper stream, and of medium to fine-grained sand and fine sand in the lower stream. This aquifer reaches a thickness of 7–18 m, with a depth of aquifer top ranging from 23 to 60 m. Flow in the middle Pleistocence aquifer (confined aquifer II) is also heavily influenced by paleochannels. The confined aquifer is primarily composed of sand and gravel, and is divided into several subaquifers in the Jiaxing area. The depth of top of this confined aquifer varies from 60 to 120 m. Confined aquifer II is much thicker than confined aquifer I, ranging from 10 to 50 m. This is the most exploited aquifer unit in this area. The early Pleistocene confined aquifer III is mainly composed of medium to coarse-grained sand and sand with gravel, with low permeability clays in the lower part of the aquifer. The depth of top of this confined aquifer ranges from 140 to 170 m. The thickness of this aquifer increases from 8 m in the southwest to as much as 117 m in the northeast. The aquitards in this region consist primarily of clay and silty clay (Zhang and others 2007; Shi and others 2008), and are distributed widely in the HJH Plain. It is noted that the west boundary of the second aquitard move northeast as the aquifer depth increases (Fig. 2). The Geological Services of the Zhejiang Province has developed a number of contour maps of depth of top and thickness of the aquifers by homogenizing lithologic information from bore holes. These contour maps are used to determine the thickness of the interbedded aquitards.

123

1112

Environmental Management (2013) 51:1109–1125

Fig. 2 Hydrogeological cross-section along A–A0 line (Fig. 1) showing the main formations of multiaquifer system underlying the Hang-Jia-Hu Plain Table 1 Description of hydraulic properties for the formations in the HJH after model calibration Layer number (from top to bottom)

Description

Average thickness (m)

Hydraulic conductivity (m/day)

Specific storage (1/m)

Layer 1

Aquifer I

11

7.4

Layer 2

Aquitard 1

18

1.6 9 10

-6

Elastic skeletal-specific storage (1/m)

Inelastic skeletal-specific storage (1/m)

1.4 9 10-6

4.2 9 10-5

5.0 9 10-4

-6

-6

3.0 9 10-4

1.4 9 10

2.0 9 10

Layer 3 Layer 4

Aquifer II Aquitard 2

25 27

19.8 1.0 9 10-6

1.4 9 10 1.4 9 10-6

2.4 9 10 3.0 9 10-6

4.5 9 10-4 2.0 9 10-4

Layer 5

Aquifer III

31

9.3

1.4 9 10-6

1.4 9 10-5

2.1 9 10-4

Hydraulic Properties Several pumping tests were performed in approximately 45 wells across the HJH Plain, the results were generally analyzed by the Theis–Jacob method (ZIGS 2004 report). The estimated transmissivity T ranges from 500 to 2,000 m2/day in the paleochannels, and between 100 and 500 m2/day in the alluvial aquifers. The horizontal hydraulic conductivity Kh is derived from the ratio of transmissivity and aquifer thickness. The estimated hydraulic conductivities vary over 1–30, 5–35, and 1–50 m/day for the aquifer I, II and III, respectively. The conductivity of clay layers is difficult to measure by in situ test. Laboratory permeability test of 46 samples from three boreholes indicate that the conductivity ranges between 10-7 and 10-5 m/day (Table 1). The elastic skeletal-specific storage and inelastic skeletalspecific storage can be obtained from consolidation experiment using both the elastic (Cs) and the inelastic (Cc) compression indices (Calderhead and others 2011). The estimated elastic compression indices vary over 0.01 to 0.08 and 0.02 to 0.09 for the aquifers and aquitards, respectively. The estimated inelastic compression indices vary over 0.1 to 0.3 and 0.1 to 0.5 for the aquifers and aquitards, respectively

123

-6

-5

(Shi and others 2008). These estimated compression indices and porosity estimates then can be used to constrain the specific elastic and inelastic storage given known effective stress. The elastic and inelastic storage coefficients can be calculated through stress–strain analysis of the observed groundwater level and compaction variations in extensometers (Riley 1969). That the HJH Plain has only a single group of extensometers (with only a short time series of available data) makes it difficult to calculate the elastic and inelastic storage coefficients of the aquifer system using Riley’s stress–strain analysis method (Riley 1969). The estimated specific elastic and inelastic storage through this method is on magnitude of 10-6 to 10-5 and 10-5 to 10-4 1/m, respectively (Zhang and others 2007).

Groundwater Exploitation and Land Subsidence Pumping and Land Subsidence Rates up to Present Time Groundwater pumping from the deep aquifer of the HJH Plain can be traced back to 1914 when a well of *170 m

Environmental Management (2013) 51:1109–1125

1113

Fig. 3 Historical groundwater withdrawals from the confined aquifers of the HJH Plain for a different locales and b the distinct aquifers

depth was drilled for drinking water. Extensive groundwater exploitation from the aquifers of the HJH Plain began in the middle of the 1960s, as groundwater development in the area dramatically increased from 1964 to 2007 (Fig. 3). Over that time period, approximately 3.4 billion m3 of groundwater was extracted from the confined aquifer system for industrial, municipal and agricultural purposes. The HJH Plain has been affected by land subsidence since 1964 (Hu and others 2004). By the end of 2007, the maximum cumulative subsidence had reached approximately 1,100 mm, and the total affected area with cumulative subsidence over 50 mm was approximately 4,200 km2 (Wu and others 2008). Over time the evolution of subsidence in the HJH Plain showed a nonlinear relationship with groundwater pumping (Fig. 4a). The overall evolution of land subsidence can be divided into four obvious periods prior to 2001: 1964–1972, 1973–1988, 1989–1996, and 1997–2000. These periods are outlined below. After 2001, the relationship between groundwater withdrawals and land subsidence became more complex. From 1964 to 1972, groundwater pumping across the HJH increased slowly and was concentrated in urban areas where average pumping rates were *20 million m3/year, which was primarily dedicated to municipal use. The measured subsidence rate was 3–15 mm/year from 1964 to 1972, with an average of 10 mm/year. Increases in land subsidence during this period were significant compared to the slow increase in total pumping rate relative to the subsequent period, 1973–1988. During this period, the

Fig. 4 Relationship between a land subsidence rate at the subsidence center in the metropolitan area of Jiaxing and local groundwater development, and b cumulative land subsidence in the land subsidence center in the metropolitan area of Jiaxing and average water levels in aquifer II and III

extent of land subsidence was restricted to the metropolitan area of Jiaxing. Groundwater pumping increased sharply beginning in 1973, and annually increasing rates averaging 13 % (largely due to municipal and agricultural use) resulted in the development of a cone of depression in the potentiometric

123

1114

surfaces of the confined aquifers in Jiaxing. The land subsidence rate in Jiaxing increased greatly, from 16 mm/ year in 1973 to *50 mm/year in 1988. During this period, the extent of subsidence had extended beyond the metropolitan area of the Jiaxing to suburban areas. By 1989, the environmental issues ensuing groundwater exploitation had attracted local government attention, and measures were taken to begin constraining pumping. However, groundwater use in rural areas still continued to increase, exceeding that of urban areas, where withdrawal had began to decrease. In the metropolitan area of Jiaxing, groundwater withdrawal decreased from *15 to *10 million m3. From 1989 to 1996, total groundwater pumping increased from *120 to 150 million m3. Concurrent land subsidence rates in Jiaxing show a decreasing trend, from 46 to 27 mm/year. During this period, land subsidence began to occur in other urban areas besides the metropolitan area of Jiaxing. It was not until 1997, when government regulations tightened, leading to reduced groundwater extraction rates and closure of a number of illegally constructed wells, that groundwater pumping began to show decreasing trends in both urban and rural areas. In 2000, the total pumping across the HJH Plain was reduced to *110 million m3, and *3 million m3 in the metropolitan area of Jiaxing. The decreasing trend in land subsidence rates in Jiaxing shows a more gradual slope than that in the previous period (Fig. 4a). From 1997 to 2000, land subsidence spread widely across the entire plain region, and the land subsidence center shifted from Jiaxing to Haiyan and Pinghu. From 2001 to 2005, groundwater exploitation began to increase again in response to the rapid deterioration of surface water quality, particularly in 2002 and 2003. The land subsidence rate increased dramatically in 2003. These increases were curtailed by severe governmental groundwater exploitation prohibition measures enacted in 2006 and 2007. In 2006 and 2007, the land subsidence rate began decreasing again, accompanying drastic reductions of groundwater withdrawal across the entire plain. Monitoring Groundwater Level and Land Subsidence Monitoring of local land subsidence in the HJH Plain has been carried out since the 1950s. In 1953, a monitoring network with *20 survey benchmarks was constructed in the metropolitan area of Jiaxing, with leveling done at second-order accuracy. Second-order leveling surveys were repeated in 1973 and 1981. The coverage of the leveling network was extended by adding to *40 leveling benchmarks in 1983. Extensive land leveling was conducted in the metropolitan area of Jiaxing from 1988 to 1994 by the Geo-Environmental Monitoring Station of the Zhejiang Province (GEMSZR), covering an area of 200 km2 (Jiang

123

Environmental Management (2013) 51:1109–1125

2008). Land subsidence monitoring was halted in 1995 and resumed in 1998, at which time the monitoring network was extended to a coverage area of 300 km2 (Jiang 2008). Wider coverage of the leveling survey was first carried out by the Department of Water Resources of the Zhejiang Province (DWRZJ). The monitoring network constructed by the DWRZJ covered more critically known land subsidence areas, and the second and third-order leveling was done in 1974 and only repeated in 1994. In 1999, the coverage of the GEMSZR monitoring network had extended to 3,500 km2 across the HJH Plain, and measurement of land subsidence thereafter was repeated every year (ZIGS 2004 report). During the last decade, traditional leveling was integrated with borehole extensor meters and Global Positioning System (GPS) technology. In April of 2005, the HJH Plain’s first vertical extensometer was installed in the metropolitan area of Jiaxing, beginning the detection of compression within the various layers of the aquifer system. The earliest steps toward the development of a regional groundwater level monitoring network were taken at the beginning of the 1980s. However, historical groundwater level measurements are rather scarce and scattered. Presently, there are a total of 86 wells available for monitoring groundwater level, 2 of which are located in the first aquifer, 53 in the second, 18 in the third, and the remaining 13 functioning as multi-aquifer monitoring wells. The groundwater monitoring network is maintained and measured monthly by GEMSZR. The evolution of the hydraulic heads in the aquifers has been divided into four main periods (ZIGS 2004 report). The first period corresponds to the pre-development conditions prior to the beginning of major groundwater extraction in the early 1950s. In this period static water levels in all aquifers were close to the ground. In the second period, from the 1950s to the early 1980s, groundwater levels decreased significantly. Lowest groundwater water levels in the center of groundwater depression cones were 22, 38 and 20 m below ground surface in 1980 in the aquifer I, II and III, respectively. The groundwater depression cones extended rapidly into regional scale in the period from the 1980s to 2000. Groundwater levels in the aquifer I were observed to be 10–40 m, with 35–40 m in the depression center below groundwater surface in 2000. An annual decline in the aquifer piezometric surface of about 1.2 m is estimated for the period of 1980–2000 in the aquifer II. Groundwater levels in aquifer II was averaged 31 m across the aquifer, with 43 m in depression cone center below groundwater surface in 2000. In the aquifer III, the annual decline in the piezometric surface is estimated 1.2 m in the 1980s and 0.8 m in the 1990s. In 2000, the average groundwater level in the aquifer III was 33 m, with 43 m in the depression center below ground surface.

Environmental Management (2013) 51:1109–1125

Current Groundwater Exploitation Prohibition Measures In 1997, in response to growing concerns of over decreasing groundwater levels and related land subsidence problems, the Jiaxing local government drafted a Land Subsidence Prevention and Reclamation Plan (Jiaxing-LSPRP) for the period of 1997–2005. The main goal of the Jiaxing-LSPRP was to prevent over-pumpage of groundwater, thereby alleviating land subsidence, particularly in the metropolitan area of Jiaxing. In 2006, the second stage of the Jiaxing-LSPRP (2006–2010) was introduced, and the Land Subsidence Prevention and Reclamation Plan of Zhejiang province (Zhengjiang-LSPRP) was drafted to prevent land subsidence across the entire Zhejiang Province. From 2006 to 2010, a series of regulations were formulated in Zhejiang Province to promote the execution of the Jiaxing-LSPRP and Zhejiang-LSPRP. The most noteworthy measures toward advancing both the Jiaxing-LSPRP and Zhejiang-LSPRP initiatives were taken in 2010, when several severe measures were introduced to reduce groundwater pumpage and close illegal groundwater pumping operations. By the end of 2010, a total of *1,200 deep wells were closed, accounting for *85 % of planned closures. The ultimate objectives of these measures and regulations have been to constrain the extension of groundwater depression cones and reduce the subsidence rate to less than 10 mm/year in most areas of the HJH Plain. As an alternative to using contaminated local surface water for water supply, several surface water diversion projects were planned and implemented, including the Taihu Water Diversion (TWD) project (Yan and Lin 2008) and the Qiandaohu Water Diversion (QWD) project (the 12th Five-Year Water Development Plan of Hangzhou, http://www.hangzhou.gov.cn, accessed in December, 2011). The TWD project is comprised of two general networks, the eastern portion covering Jiashan and Pinghu, and the western portion covering the other administrative regions of Jiaxing. By the end of 2010, the water supply capacity of the TWD network had achieved 0.4 km3/year. The QWD project has been proposed to provide water to the area of the Hangzhou City from 2011 to 2020. The model used in this study was calibrated based on various subsidence datasets from 1996 to 2007, used to evaluate the LSPRPs (based on datasets from 2008 to 2010), and then used to predict future land subsidence during the upcoming QWD period of 2011–2020. Future Pumping Rates Scenarios In 2009, a special action taken in Jiaxing to close illegal pumping wells resulted in the closure of 396 wells. Total

1115

groundwater pumping was reduced to 60 million m3. By the end of 2010, a total of 885 wells had been closed, and groundwater pumping sharply decreased to 10 million m3. These dramatic decreases in groundwater exploitation will likely be difficult to maintain given economic demands. With the delicate balance between sustainable resource exploitation and human demand, predictive models herein considered have been designed to simulate more gradual rates of reduction in groundwater usage. Following are the three scenarios considered in this study. Scenario A This scenario assumes current groundwater pumping practices are continued. Under the current primary surface water diversion project, the water supply capacity of the TWD project is dependent upon the water quality of the Taihu Lake. However, the Taihu Lake has been suffering from water quality degradation (i.e., the blue-green algae outbreak in May, 2007) (Shi and others 2012). Therefore, considering the possibility of sudden pollution of the Taihu Lake, this scenario considers the best-fit relationship between current surface water diversion and groundwater pumping. Scenario B In 2009, as a result of well closure measures and the use of surface water from the TWD, average groundwater withdrawal decreased by 20 % in contrast to the previous year. This scenario considers the gradual increase of surface water use and decrease of groundwater pumping from 2011 to 2020. The model assumes an annual decrease in groundwater pumping of 20 % (based on the aforementioned 2009 value) across the entire HJH Plain. Scenario C This scenario assumes the planned QWD network is able to supply surface water (for domestic use) to Jiaxing from 2011 to 2015. This allows for an annual reduction in groundwater withdrawal of 25 and 20 % for confined aquifers II and III (respectively), within the urban areas of Jiaxing. After 2015, the model holds the pumping schedule for 2015 constant from 2016 to 2020.

Numerical Simulation of Groundwater Flow and Land Subsidence Numerical Simulation of One-Dimensional Compaction Land subsidence is simulated by a modular package of MODFLOW-2000 (Harbaugh and others 2000) called the Subsidence (SUB) Package (Hoffmann and others 2003). The SUB package is based on the theory of one-dimensional (vertical) consolidation, developed by Terzaghi (1925), which assumes that deformation is vertical and that the total load (total stress) does not vary. The SUB package calculates

123

1116

Environmental Management (2013) 51:1109–1125

vertical compaction based on changes in effective stress, assuming change in hydraulic head produces an equal but opposite change in effective stress. This assumption introduces error in cases of shallow, unconfined aquifers, but is reasonable when dealing with deep, confined aquifers (Leake and Prudic 1991). The compaction of the compacting layer, Db, is calculated as: Db ¼ Ssk bDh

where Sske is the elastic skeletal-specific storage; Sskv is the inelastic, or virgin, skeletal-specific storage, and r0zzðmaxÞ is the preconsolidation stress. Storage changes and the corresponding compaction are computed at every time-step. The expression that accounts for total flux q (flow per unit area) into or out of elastic and/ or inelastic skeletal storage at cell i for time-step m is: S0m S0 k ðhm  H m1 Þ þ kem ðH m1  hm1 Þ m Dt Dt

ð3Þ

where ( S0m k

¼

S0ke hm [ H m1 S0kv hm  H m1

ð4Þ

Variables hm and Hm-1 are the head at the end of timestep m and the preconsolidation head at the end of timestep m - 1 (respectively). S0k is the skeletal storage coefficient; S0ke is the elastic skeletal storage coefficient; and S0kv is the inelastic skeletal storage coefficient. The preconsolidation head is the lowest head the aquifer system has experienced. In the SUB package, the preconsolidation head is used to switch between elastic and inelastic storage properties (Leake 1990). When the water level drops below the preconsolidation head, the storage coefficient changes from elastic to inelastic values, making compaction inelastic and unrecoverable. If inelastic compaction has occurred, a new value for preconsolidation head is recorded. The package assumes elastic and inelastic storage coefficients to be constant. This assumption causes problems in shallow, unconfined aquifers, but is satisfactory in deep, confined aquifers (Leake and Prudic 1991). The SUB package simulates the delay in release of water from

123

o2 h S0s oh ¼ oz2 Kv0 ot

ð5Þ

where S0s is the specific storage of the interbed; Kv0 is the vertical hydraulic conductivity of the interbed.

ð1Þ

where Ssk is the skeletal-specific storage, b is the original thickness of the sediments, and Dh is change in hydraulic head. To account for the marked change of the skeletalspecific storage when the effective stress exceeds the preconsolidation stress, two separate values are often used: ( Sske r0zz \r0zzðmaxÞ Ssk ¼ ð2Þ Sskv r0zz  r0zzðmaxÞ

qm i ¼

compressible interbeds using the one-dimensional diffusion equation:

Model Construction Model Discretization and Boundary Conditions The finite difference grid for the model was divided into five vertical layers. Although the best way to generate a representative 3D geological model presenting the sand and clay layers, it generally requires a large number of borehole logs (Calderhead and others 2011). The lack of borehole logs makes it difficult and impractical to reliably simulate the compaction behaviors of all clay interbeds in the confined aquifers. The available top elevation and thickness contours for each aquifer are used to determine the top elevation and thickness of each aquitard. The hydrogeologic properties of the confined aquifer (e.g., vertical hydraulic conductivity and storage coefficients) are chosen such that their effect on land subsidence is equivalent to the composite effect of all interbeds. This method is most commonly implemented to simulate the compaction of multiaquifer system (e.g., Teatini and others 2006; Shi and others 2008). Each model layer consists of 230 rows and 320 columns, with a uniform horizontal spacing of 400 m. Layers 1, 3 and 5 represent confined aquifers I, II and III, respectively, while layers 2 and 4 represent the two aquitards sandwiched between the aquifers (Fig. 5). Based on available groundwater level records, land subsidence data, and the seasonal withdrawal pattern of groundwater in the HJH Plain, the simulation was carried out under transient conditions for the 12 year period of 1996–2007. The total simulation period was divided into 48 stress periods with 3 months in each stress period. The confined aquifer system in the HJH Plain is separated by a widely distributed low permeable layer of soft mud sediments from the ground surface (Fig. 2). The low permeable formations virtually obstruct the infiltration from precipitation into the confined aquifers and the connection between rivers and the confined aquifers where unconfined aquifer is absent. The interaction between the rivers (including the Great Channel) and the unconfined shallow groundwater varies seasonally. In dry season when groundwater level is generally higher than river stages, shallow groundwater discharges to rivers. Shallow groundwater can receive leakage from rivers in flooding season. Therefore, the flux into the top boundary of the

Environmental Management (2013) 51:1109–1125

1117

Fig. 5 a Plane view of and b vertical model discretization along the cross-section B–B0

model is leakage through connection with the overlying unconfined shallow aquifer, and the top boundary is prescribed as specified flux boundary (zero in areas with unconfined aquifer absent). The bottom of aquifer system contacts with Pre-Quaternary low-permeability bedrocks directly and therefore is assumed no-flow boundary. Because groundwater levels in the confined aquifers II and III are lower than the water stages in the Taihu Lake and the Hangzhou Bay, and are generally lower than the groundwater levels in adjacent aquifers in northeast, these aquifers receive lateral flux from the Taihu Lake, Hangzhou Bay, and adjacent aquifers. Little flux is also allowed through permeable boundaries contacting with mountain front sediments in the upper stream area. As the lateral flow boundaries in the study area are poorly defined (no groundwater level observations along these boundaries are available), the model uses specified flux boundaries to make a reasonable representation of field conditions. The specified flux in the lateral boundaries is modified further in the model calibration procedure with water budget constraints. The depth to bedrock in some local areas along the coast line of the HJH Plain is small, with some scattered bedrock outcrops in Haiyan and Pinghu. This implies

limited connection of groundwater in the HJH Plain to the sea in these areas. In other regions the groundwater regime extends offshore; and specified flux is assumed along these boundaries, with no density effect taken into account because seawater intrusion has never been reported in the HJH Plain. Initial Conditions The initial starting point of the transient run is set to the end of 1995 at which point the relative compaction restarts from zero. Despite the fact that the system has been under stress since the 1920s and leveling survey has begun since the 1950s, sufficient historical groundwater levels and pumping data were not available until the 1990s. The groundwater system was in a transient state during the 1990s owing to continuous decline in groundwater levels. Steady-state simulations using 1995 stresses and groundwater level constraints cannot capture the ongoing responses. Therefore, there is little choice but to begin the simulation with initial conditions derived from measured groundwater levels. Although the specified initial conditions of the groundwater system is inconsistent, to some

123

1118

degree, with the conservation equations of the simulation code, it is considered an adequate starting point. The water level distributions of all aquifers are used directly as the initial heads, with values based on potentiometric maps for 1995. Initial heads between the adjacent aquifers are linearly interpolated as inputs for the aquitards. Modifications are made in initial water level distributions to equilibrate the model before the actual transient run begins. Model preconsolidation heads for all aquifers were estimated by comparison of historical records of cumulative land subsidence and groundwater levels at the subsidence center located in Jiaxing (Fig. 4b). From 1964 to 1970, cumulative land subsidence slowly increased with concurrent slow decreases in groundwater levels. From 1971 to 1980, groundwater level decreases were followed by significant increases of cumulative compaction. After 1981, cumulative land subsidence continuously increased despite water levels remaining nearly constant. This suggests that compaction after 1981 is inelastic and may be delayed, which is the primary cause of permanent, irreversible land subsidence. In consideration of this, the initial model heads were set at preconsolidation head levels. Heads were updated in each simulation time-step, allowing successive simulations with decreasing water levels to result in permanent compaction (Wilson and Gorelick 1996). Preconsolidation heads for the aquitards are linearly interpolated between the heads in adjacent aquifers under the assumption that, at the onset of subsidence, heads at the top and bottom of the aquitards are in equilibrium with those in the adjacent aquifers.

Environmental Management (2013) 51:1109–1125

Other parameters determined during model calibration were the leakage from the overlying unconfined aquifer to the confined aquifer system, and the lateral boundary flux. These parameters and boundary conditions were manually adjusted until the model-calculated values sufficiently matched measured water levels and subsidence amounts, and the calibrated parameter zonation and parameter values were deemed reasonable. The integrated flow and land subsidence model successfully reconstructed the historical water level dynamics and land subsidence development trend in the HJH Plain. The simulated subsidence centers coincide closely with the groundwater depression cones and the areas with the greatest clay thickness. The development of land subsidence shows a strong relationship to groundwater pumping from confined aquifers both temporally and spatially. The center of subsidence has migrated from the metropolitan area of Jiaxing to surrounding urban areas. Simulation results indicate that the primary compaction occurring within the confined aquifer system is localized in the aquitards adjacent confined aquifer II. As groundwater pumpage and its distribution (temporally and spatially) are prominent factors during model calibration, considerable efforts were made to estimate the model pumping distribution based on county-level groundwater exploitation data. Estimated pumping amounts and distributions were adjusted slightly during the calibration process in some pumping cells.

Results and Discussion History Match of Water Level and Land Subsidence

Model Calibration The numerical model was calibrated against historical records of groundwater level and land subsidence with a trial-and-error procedure. Very little groundwater is extracted from aquifer I, and little data are available on water levels within this aquifer, so model calibration of water levels focused on aquifers II and III. Three major hydraulic parameters were determined during the model calibration process: hydraulic conductivity, elastic storage coefficient and inelastic storage coefficient. Based on lithologic facies and paleochannel distributions, the aquifers were divided into a number of zones with K values ranging from 5 to 30 m/day. The hydraulic conductivities for the aquitards are estimated at 10-5 to 10-7 m/day (Li and others 2006), with one uniform value assigned to each individual aquitard. The initial values of the elastic and inelastic storage coefficients are thus assigned as the values estimated by Zhang and others (2007) for the corresponding aquifer system. The hydraulic parameter zonations and values were further adjusted during model calibration.

123

The simulated and measured contour map water levels for confined aquifers II and III, as of 2007, are shown in Fig. 6. These show fairly good agreement, spatially, between observed and calculated water level features, such as the maximum drawdown points and locations of depression cones. Residuals of groundwater levels for aquifer II and III (Fig. 7) show that simulated groundwater levels higher than measured over areas around the aquifer boundaries with no observation wells located. Four wells (two from confined aquifer II and the other two from confined aquifer III) were selected that had long time series observation datasets to examine the temporal trend between observed head values and calculated heads values. Calculated heads at selected observation wells show good agreement with observed hydrographs over time (Fig. 8). Groundwater levels in the confined aquifers decrease overall (except of well OBS#12 which shows seasonal water level variations), with few locations showing recovery after declines in groundwater withdrawal, beginning in 2007 (Fig. 8).

Environmental Management (2013) 51:1109–1125

1119

Fig. 7 Groundwater level residuals (measured–simulated) in a aquifer II and b aquifer III at the end of 2007. The selected wells showing simulated and measured dynamics in Fig. 9 are labeled

Fig. 6 Contour maps of simulated and measured water levels in a aquifer II and b aquifer III at the end of 2007 (contour interval is 5 m)

Sensitivity Analysis

Due to the short time series of subsidence monitoring data from the single group of extensometers in the plain, historical contour maps of land subsidence were used to calibrate the subsidence model. There is a reasonable agreement between simulated and observed land subsidence at the end of 2007 in terms of magnitude, subsidence centers, and the regional spatial extent of subsidence (Fig. 9). These simulation results are acceptable considering the complicated characteristics of the aquifer system and limited available data. Measures of goodness of fit were made between observed and calculated water levels for all monitoring wells in aquifers II and III, and between observed and calculated cumulative land subsidence at the end of simulations at all leveling benchmarks (Fig. 10). The mean error (ME) between observed and calculated water levels was 0.35 m, and the root mean square error (RMSE) was 5.64 m for a total of 1,624 observations. The ME and RMSE between observed and calculated cumulative subsidence were 22 and 114 mm (respectively) for a total of 216 benchmarks.

The sensitivity of calibrated model response to model inputs was evaluated to help assess model reliability. The values of vertical hydraulic conductivity, elastic and inelastic storage coefficients, and initial preconsolidation heads were varied individually over ranges that reflect reasonable parameter uncertainty. As sand elastic and clay inelastic storativity could range over several orders of magnitude, those parameters were varied over large ranges than other parameters. Since the land subsidence rate is the factor of greatest concern in this study, the model sensitivity analysis was conducted on the above listed parameters using a land subsidence rate greater than 10 mm/year by the end of 2007 as criteria (Fig. 11a). The root mean square errors of cumulative land subsidence at bench marks for different parameter sets are presented to reflect the reliability of the model (Fig. 11b). The results of this analysis indicate that the model is more sensitive to changes in inelastic coefficients and preconsolidation heads than the other two parameters. Increasing vertical hydraulic conductivity or elastic storage coefficients produces a positive trend in land subsidence

123

1120

Environmental Management (2013) 51:1109–1125

Fig. 8 Observed and simulated piezometric heads at four selected observation wells (see labeled observation wells in Fig. 8)

(this same covariance is seen as values decrease). Increases in either inelastic storage coefficients or initial preconsolidation heads have a similar effect on land subsidence (with the same covariance seen as values decrease). The model is more sensitive to decreases than increases in the inelastic storage coefficient, and, conversely, more sensitive to increases than decreases in preconsolidation heads. In the early stages of increasing inelastic storage coefficients, increases in inelastic storage coefficients result in increased rates of subsidence. However, in the latter stages of increasing inelastic storage coefficients, subsidence rates begin to decrease, as more water is released from storage in the aquitards and interbeds, groundwater drawdown decreases. In the case of decreasing inelastic storage coefficients, changes in inelastic storage coefficients are more sensitive on water level than land subsidence, making the subsidence rate increase at first and then decline. The available extensometer monitoring compactions in the HJH and adjacent regions show that the second and third aquifers and adjacent aquitards are the primary consolidation layers (e.g., Shi and others 2007; Zhang and others 2007). One possible reason can explain the considerable inelastic compaction of the land layers is that sand layers in the HJH contain certain content of clay materials and clay interbeds (Xue and others 2008). This also can explain the relatively higher inelastic specific storage for the aquifers verified through model calibration (Table 1). Because the groundwater level in the HJH Plain is below mean sea level, multiplying initial preconsolidation heads

123

by values greater than one effectively lowers (decreases) the preconsolidation head and multiplying by values less than one effectively raises (increases) the preconsolidation head. Lowering the preconsolidation head, results in more elastic subsidence and decreases irreversible land subsidence rates. Decreasing preconsolidation heads effectively decreases the volume of water released resultant aquifersystem compaction, and accelerates the water level drawdown, ultimately causing increases in land subsidence rates. Increases in preconsolidation head show little effect on simulations since the model updates the preconsolidation head to the new lowest water level value for each step. Evaluation of Groundwater Pumping Management Strategies The lowest predicted water levels at the centers of water level depression cones in confined aquifers II and III from 2008 to 2020 are shown in Fig. 12. Aquifers II and III show significant rebounds in water levels in 2010 resultant reductions in groundwater pumping. In the case of Scenario A, water levels continue to increase in 2011, then show a trend of continuous decline in aquifers II and III. However, this decrease is gradual, showing overall declines of 1 m in aquifer II and 1.5 m in aquifer III over the span of 2011–2020. This can be explained as the result of considerably reduced groundwater pumpage in 2010; pumping in 2010 is approximately one tenth the average annual pumpage from 1996 to 2007.

Environmental Management (2013) 51:1109–1125

1121

Fig. 9 Cumulative land subsidence at the end of 2007 as a measured and b simulated by the land subsidence model (contour interval is 100 mm)

The water levels at the centers of depression cones are -45 and -53 m in aquifers II and III, respectively. These water levels are still higher than values before the groundwater pumping prohibition measures taken in 2010. Scenarios B and C show similar trends of general rebound in water levels beginning in 2011. In 2020, total groundwater

pumping is reduced to *2.0 million m3 in Scenario B and *5.0 million m3 in Scenario C. This results in rebound of water levels at the centers of depression cones of 7 m in aquifer II and *15 m in aquifer III. All three scenarios show essentially the same trend of decreases in area with land subsidence rates [10 mm/year

123

1122

Environmental Management (2013) 51:1109–1125

Fig. 10 a Comparison between measured and simulated heads at observation wells in aquifers II and III, and b comparison between measured and simulated cumulative land subsidence at all leveling benchmarks in 2007. Dotted lines indicate the 95 % prediction interval

Fig. 12 Lowest predicted water levels for all three pumping scenarios from 2008 to 2020 for a confined aquifer II and b aquifer III

Fig. 11 Sensitivity of calibrated model response to changes in selected model input parameters with respect to a land subsidence rate and b root mean square error of cumulative land subsidence at bench marks in 2007

across the HJH Plain (Fig. 13). The cumulative areal extent of land subsidence rates [10 mm/year across the Plain, for all scenarios, decreases to *430 km2 (*7 % of the entire plain) in 2010, and is expected to asymptotically approach zero by 2020. However, increase in the extent and

123

magnitude of regional land subsidence is still observed. In all three scenarios, the total area across the plain that is affected by cumulative land subsidence rates [100 mm is *3,700 km2 (*60 % of the entire plain) by 2020 (Fig. 13). Planned cumulative groundwater withdrawal over the prediction period 2008–2020 for Scenarios A, B and C is 352, 246, and 250 million m3, respectively. The planned pumping rates for the three scenarios from 2011 to 2020 are 18, 7, and 7.5 million m3/year, respectively. This suggests that in cases of water stress, groundwater usage

Environmental Management (2013) 51:1109–1125

1123

Another issue is the representation of the lateral boundaries and groundwater pumpage distributions. The boundary of the study area is partially defined along administrative borders, and, as a result, no observation wells are present along the boundaries. Therefore, the lateral boundaries of the regional flow and land subsidence model are defined primarily through model calibration, the impact of which should be further analyzed. Pumping data in the study area are estimated by groundwater management authorities with respect to each aquifer unit. However, the depth of pumping wells is poorly established. More detailed information about local pumping wells and certainly improve simulation of local land subsidence development, particularly near subsidence centers.

Conclusions

Fig. 13 Predicted area subject to land subsidence rates greater than a 10 mm/year and b 100 mm cumulative for different pumping scenarios

can be increased with negligible consequences to 8–18 million m3/year. Model Limitation Although great efforts toward model calibration have been taken in this study, several issues could be addressed in the future to increase model efficacy. A crucial one is ignoring the delay between decreases in groundwater level and clay compaction. The consolidation and subsidence process depends on the rate of drainage of pore water from the compressible layers and is controlled by the hydraulic diffusivity. Thus, a delay term is often required in the simulation of compaction process (Hoffmann and others 2003). However, the short period coverage of extensometer monitoring data currently available in the HJH make it difficult to estimate the delay between decrease in groundwater level and compaction. More available extensometer monitoring data in future can be used to improve the model with respect to delay process in clay layers. What is more, the numerical model would be significantly strengthened by more frequent collection of subsidence data and the installation of a larger extensometer network. The current lack of subsidence data from extensometers leaves annually collected land subsidence data from leveling investigations as the most complete dataset for use in calibrating storage coefficients. This complicates calibration, as seasonal land subsidence rebound data cannot be determined when data are only collected at yearly intervals.

A regional three-dimensional groundwater flow model coupled with a one-dimensional aquifer-system compaction model was developed to simulate the historical land subsidence generated by groundwater extraction, and to predict future land subsidence in the HJH Plain considering various management scenarios. Although the coupled simulation model is limited in its representation of the complex aquifer and groundwater flow system due to limited available data, it adequately matches the observed development of land subsidence related to the concurrent groundwater pumping in the HJH Plain. The numerical model provides a viable method for evaluating future water sustainability scenarios with the consideration of inextricably linked land subsidence issues. Several potential groundwater management alternatives were tested using the calibrated model, and future land subsidence evolution and corresponding available groundwater yield were estimated. This study shows reduction of groundwater pumping from the confined aquifers, particularly aquifer II (which is currently the most taxed aquifer unit), to be the most effective method of curtailing future development of land subsidence. Continual and effective implementation of current land subsidence prevention plans, reclamation plans, and the scheduled surface water diversion projects can reduce land subsidence due to groundwater pumping to a minimal concern across the HJH Plain over the next decade. Acknowledgments This research was supported by the strategic science and technology project of the Exploring Advanced Discipline (No. 2012QY007) in the Institute of Geographic Sciences and Natural Resources Research, China Academy of Sciences. The authors want to thank Jiankang Zhao and Mengjie Wu, from the Geo-Environment Monitoring Station of Zhejiang Province for providing the data necessary for this study. Reviewers and journal editors are thanked for their useful comments in helping to improve the original manuscript.

123

1124

References Calderhead AI, Therrien R, Rivera A, Martel R, Garfias J (2011) Simulating pumping-induced regional land subsidence with the use of InSAR and field data in the Toluca Valley, Mexico. Adv Water Resour 34:83–97 Chang YL, Tsai TL, Yang JC, Tung YK (2007) Stochastically optimal groundwater management considering land subsidence. J Water Resour Plan Manag 133:486–498 Daito K, Mizuno M, Ueshita K (1991) Control of groundwater withdrawal for preventing land subsidence in the Owari Plain, Japan. The 4th International Symposium on Land Subsidence, Houston, pp 533–542 Gallardo AH, Marui A, Takeda S, Okuda F (2009) Groundwater supply under land subsidence constrains in the Nobi Plain. Geosci J 13:151–159 Galloway DL, Burbey TJ (2011) Review: regional land subsidence accompanying groundwater extraction. Hydrogeol J 19:1459–1486 Galloway D, Jones DR, Ingebritsen SE (1999) Land subsidence in the United States: U.S. Geological Survey Circular 1182 Harbaugh AW, Banta ER, Hill MC, McDonald MG (2000) MODFLOW-2000, The U. S. Geological Survey modular groundwater model—user guide to modularization concepts and the ground-water flow process. U.S. Geological Survey Open-File Report 2000–92 Hoffmann J, Leake SA, Galloway DL, Wilson AM (2003) MODFLOW-2000 ground-water model—user guide to the subsidence and aquifer-system compaction (SUB) package. U.S. Geological Survey Open-File Report 03-233 Hu RL, Yue ZQ, Wang LC, Wang SJ (2004) Review on current status and challenging issues of land subsidence in China. Eng Geol 76:65–77 Jean DB, Tang Z (2008) A GIS based DRASTIC model for assessing groundwater vulnerability in shallow aquifers in Hangzhou– Jiaxing–Huzhou Plain, China. Res J Appl Sci 3:550–559 Jiang D (2008) Land subsidence analysis of Hangzhou–Jiaxing– Huzhou Plain in Zhejiang Province. Bull Surv Mapp 7:13–15 Jiang Y, Zhao J, Zhu C (1998) Analysis on land subsidence and reasonable exploitation of groundwater in Jiaxing City, Zhejiang Province. J Geo Hazard Control 9:36–41 Larson KJ, Bas¸ ag˘aog˘lu H, Marin˜o MA (2001) Prediction of optimal safe ground water yield and land subsidence in the Los Banos– Kettleman City area, California, using a calibrated numerical simulation model. J Hydrol 242:79–102 Leake SA (1990) Interbed storage changes and compaction in models of regional groundwater-flow. Water Resour Res 26:1939–1950 Leake SA, Prudic DE (1991) Documentation of a computer program to simulate aquifer-system compaction using the modular finitedifference ground-water flow model: U.S. Geological Survey Techniques of Water-Resources Investigations, book 6, chap A2, p 68 Li C, Tang X, Ma T (2006) Land subsidence caused by groundwater exploitation in the Hangzhou–Jiaxing–Huzhou Plain, China. Hydrogeol J 14:1652–1665 Li C, Ma T, Wang Y, Zhu X (2009) The approach and countermeasures to solve water resource problems and reconstructing a good ecological environment in Yangtze Delta Region. Int J Water 5:1–14 Liu X, Wu J, Xu J (2006) Characterizing the risk assessment of heavy metals and sampling uncertainty analysis in paddy field by geostatistics and GIS. Environ Pollut 141:257–264

123

Environmental Management (2013) 51:1109–1125 Onta PR, Gupta AD (1995) Regional management modeling of a complex groundwater system for land subsidence control. Water Resour Manag 9:1–25 Ortiz-Zamora D, Ortega-Guerrero A (2010) Evolution of long-term land subsidence near Mexico City: review, field investigations, and predictive simulations. Water Resour Res 46:W1513 Poland JF (1984) Guidebook to studies of land subsidence due to ground-water withdrawal. UNESCO, New York Poland JF, Davis GH (1969) Land subsidence due to withdrawal of fluids. In: Varnes DJ, Kiersch GA (eds) Reviews in engineering geology II. Geological Society of America, New York, pp 187–269 Riley FS (1969) Analysis of borehole extensometer data from central California. Land Sub 2:423–431 Shi X, Xue Y, Ye S, Wu J, Zhang Y, Yu J (2007) Characterization of land subsidence induced by groundwater withdrawals in Su-XiChang area, China. Environ Geol 52:27–40 Shi X, Wu J, Ye S, Zhang Y, Xue Y, Wei Z, Li Q, Yu J (2008) Regional land subsidence simulation in Su-Xi-Chang area and Shanghai City, China. Eng Geol 100:27–42 Shi X, Fang R, Wu J, Xu H, Sun Y, Yu J (2012) Sustainable development and utilization of groundwater resources considering land subsidence in Suzhou, China. Eng Geol 124:77–89 Teatini P, Ferronato M, Gambolati G, Bertoni W, Gonella M (2005) A century of land subsidence in Ravenna, Italy. Environ Geol 47:831–846 Teatini P, Ferronato M, Gambolati G, Gonella M (2006) Groundwater pumping and land subsidence in the Emilia-Romagna coastland, Italy: modeling the past occurrence and the future trend. Water Resour Res 42:W1406 Terzaghi K (1925) Earthworks mechanics based on soil physics (Erdbaumechanik auf bodenphysikalischer Grundlage). Franz Deuticke, Vienna Wilson AM, Gorelick S (1996) The effects of pulsed pumping on land subsidence in the Santa Clara Valley, California. J Hydrol 174:375–396 Wu M, Zhao J, Liu S, Shen H, Xue D (2008) Prevention and control of land subsidence in Zhejiang. Zhejiang Land Res 3:54–57 Xu T, Van der Gun J (1995) Predicting land subsidence with a constant-parameter coupled model for groundwater flow and aquitard compaction: The Markerwaard case. IAHS Publ 234:369–379 Xu YS, Shen SL, Cai ZY, Zhou GY (2008) The state of land subsidence and prediction approaches due to groundwater withdrawal in China. Nat Hazards 45:123–135 Xue Y, Wu J, Zhang Y, Ye S, Shi X, Wei Z, Li Q, Yu J (2008) Simulation of regional land subsidence in the southern Yangtze Delta. Sci China Ser D: Earth Sci 51:808–825 Yamamoto S (1995) Recent trend of land subsidence in Japan. IAHS Publ 234:487–492 Yan B, Lin H (2008) Research on diversion water scheme of improving environmental water quality in the Hangjiahu area. China Rural Water Hydropower 9:33–35 (in Chinese with English abstract) Zhang GF (2005) Situation of prohibition and limitation extraction groundwater aspect in Zhejiang and its countermeasure. Zhejiang Hydrotech 5:1–3 (in Chinese with English abstract) Zhang Y, Xue YQ, Wu JC, Ye SJ, Wei ZX, Li QF, Yu J (2007) Characteristics of aquifer system deformation in the Southern Yangtse Delta, China. Eng Geol 90:160–173 Zhao J, Wu M, Liu S (2003) Status and the prevention and treatment countermeasures of land subsiding in the coastal plain. J Geo

Environmental Management (2013) 51:1109–1125 Hazards Environ Preserv 14:16–20 (in Chinese with English abstract) Zhao J, Wu M, Liu S (2004) Study on land subsidence and monitoring network in the Hangjiahu plain of Zhejiang Province. J Geo Hazards Environ Preserv 1:16–20 (in Chinese with English abstract) Zhao J, Wu M, Liu S, Shen H (2006) Groundwater exploitation and land subsidence in the coastal plain, Zhejiang Province. Geol J China Univ 12:185–194 (in Chinese with English abstract)

1125 Zhejiang Institute of Geological Survey (ZIGS) (2004) Groundwater sustainable exploitation and geological disasters survey and assessment for Hang-Jia-Hu region in Yangtze River Delta (unpublished report)

123