Three dimensional numerical modeling of land subsidence in Shanghai

Shanghai city has been suffering land subsidence caused by overly exploitation of ground water since 1921, which is a serious problem for this coastal city with altitude of 2.2–4.8 m above mean sea level. The largest cumulative land subsidence amounted to 2.6 m in the downtown area. Measures to decrease the ground water exploitation, change the pumping aquifers, and increase aquifer artificial recharge have been used to mitigate land subsidence since 1961. It is necessary to develop a proper numerical model to simulate and predict land subsidence. In this study, a decoupled three-dimensional (3-D) finite element land subsidence model including a 3-D ground water flow model and a 3-D geo-mechanical model was developed to simulate the 3-D deformation of the aquifer systems in the center area of Shanghai. The area of downtown Shanghai is 660 km, with 10 million inhabitants, dense high buildings, and 11 metro lines. The simulation spans the period from 1979 to 1995. Two different assumptions have been tested on the side boundary, i.e., precluding the three components of the displacement, or assuming a free-displacement condition. The distribution of calculated land subsidence and horizontal displacements in different aquifers was analyzed. The computed vertical displacement fitted well with the available observations. It has been verified that the two different assumptions on the lateral boundaries in the geo-mechanical model caused different results just limited on nodes close to boundary. The developed 3-D land subsidence model is reasonable and can be used to simulate and predict 3-D movement of aquifer systems in the center area of Shanghai, which could provide scientific support to local government in controlling land subsidence and differential movements of the land surface.


Introduction
Anthropogenic land subsidence in Shanghai, which is induced by excessive groundwater withdrawal, has been observed over the last 94 years, i.e. since the year 1921 (Zhang and Wei, 2002).Shanghai government started to set up extensive monitoring networks of groundwater levels and land subsidence from the middle 1980s.Some earlier developed models (Gu et al., 1990(Gu et al., , 1991(Gu et al., , 1993) ) were focused on the shallow sedimentary layers, i.e. the upper 100 m depth above the second confined aquifer, and did not consider the whole aquifer system.Soil deformation was assumed one dimensional (1-D).The researchers from the Nanjing University, China, and the Shanghai Institute of Geological Survey have conducted comprehensive studies on numerical modeling of land subsidence since 2000.A regional land subsidence model, composed of a three dimensional (3-D) groundwater flow module and a 1-D subsidence module embodying complex visco-elasto-plastic deformation characteristics with variable parameters, was developed and applied to simulate and predict land subsidence in Shanghai (Ye, 2004;Ye et al., 2005;Xue et al., 2008;Shi et al., 2008;Wu et al., 2010;Ye et al., 2011Ye et al., , 2012)).The 3-D groundwater flow model is used to calculate the time-dependent hydraulic heads in the sedimentary sequence, and the 1-D subsidence model is used to compute the time-dependent compaction of the sedimentary layers and the cumulative land subsidence.
All the land subsidence models mentioned above are based on the assumption of a 1-D deformation field in the vertical direction.However, the actual soil deformation caused by groundwater withdrawal is fully three dimensional.The deformation of an aquifer system caused by fluid withdrawal is theoretically described by the 3-D fully coupled poroelasticity model originally developed by Biot (1941).Because of ill-conditioning and large CPU-time consume of coupled models when applied to regional land subsidence problems, uncoupling between the flow and the strain fields is usually assumed in the classical groundwater hydrology.In an uncoupled formulation, also known as a two-step or explicitly coupled approach (Teatini et al., 2010), the groundwater flow equation is first solved and then land subsidence is computed solving the equilibrium equations with the groundwater changes specified as an external distributed source of strength within the porous medium.
A 3-D uncoupled land subsidence model is developed in this study to simulate 3-D displacements caused by ground water exploitation from the multiaquifer system in the center area of Shanghai.Downtown Shanghai is bounded by the outer ring road as shown in Fig. 1, with an area of 660 km 2 and the largest cumulative land subsidence of 2 m since 1921.This zone is the most important area in Shanghai, with the densest 10 million population and concentration in high-rise buildings, and crossed by 11 metro lines.Therefore, it is important to calculate the 3-D deformation of the aquifer system caused by ground water exploitation to better evaluate the possible impact on the underground infrastructures such as metro lines, tunnels etc.

Model setup
In geomechanical investigations at regional scale, uncoupled pressure -displacements solution can be safely used in predicting land subsidence in compacting sedimentary on any time scale of practical interest (Gambolati et al., 2000;Teatini et al., 2006Teatini et al., , 2010)).In this study, an uncoupled approach is thus followed, with the flow field initially computed at each time step and the deformation then calculated using the pore pressure gradient as a known source of strength in the geomechanical model.The 3-D uncoupled land subsidence model is expressed as (Verruijt, 1969): where E is the Young modulus, ν is the Poisson ratio, u x , u y and u z are the displacement components along the coordinate axes x, y and z, respectively, ε is the incremental volume strain ε = ∂u x ∂x + ∂u y ∂y + ∂u z ∂z , P is the incremental pore water pressure, which is equals to γ w • H with γ w the specific weight of water and H the incremental hydraulic head, and ∇ 2 is the Laplace operator.In Eq. (1d), K i,j is the hydraulic conductivity tensor, t is time, q is the source/sink term, and S s is the specific storage S s = γ w (α + nβ), with α the vertical oedometric compressibility, which is related to the Young modulus and Poisson ratio through α = , n is the medium porosity, and β the water volumetric compressibility.The parameter α is not constant and takes the elastic value α ke when the effective stress is less than the preconsolidation stress, otherwise the plastic value α kv (Ye et al., 2011(Ye et al., , 2012)).The specific storage is denoted by S ske and S skv , accordingly.Special care must be taken in the explicitly coupled approach to use S s values in the groundwater flow Eq. (1d) consistent with those used for the geomechanical parameters E and ν in the equilibrium Eqs.(1a)-(1c).
The simulated domain, which is the center area of Shanghai (Fig. 1), is comprised between the ground surface and a rigid 250-300 m deep basement (Fig. 2).A 3-D mesh consisting of 19 524 nodes and 34 309 triangular prisms is developed to accurately represent the multiaquifer system including six aquifers and five aquitards (Fig. 3).Zero flux on the basement and Dirichlet conditions on the top and outer boundaries are prescribed in the flow model.The known hydraulic heads on the top and outer boundaries are obtained from the ground flow model developed for the whole Shanghai city from 1961 to 2005 (Xue et al., 2008).A Dirichlet  The simulations span the interval from 30 December 1979 to 30 December 1995.There are many recharge wells and a few discharge wells in the second confined aquifer, and many discharge wells with high flow rates and some recharge wells with low flow rates in the fourth confined aquifer.The 16year long simulation period is divided into 64 time steps 3month long each.The initial groundwater levels in aquifer A2 ranged from −1.5 to 2 m above mean sea level (m.s.l.) with a flow generally from west to east, and that in aquifer A4 ranged from −22 to −10 m above msl with a flow generally from east to west.

Parameter subzones and parameter values
The whole studied aquifer system is finally divided into 103 subzones according to the difference in lithology, hydraulic, and mechanical properties of sediments.Due to lack of specific information, the Poisson ratio ν is set to 0.3 which is a typical values for recent sedimentary deposits.Hydraulic conductivities, elastic specific storage, and plastic specific storage of aquifers and aquitards are provided in Table 1.

Deformation and land displacements
Figure 6 compares the model results in term of vertical compaction/expansion and the observed deformation in a few depth intervals monitored at the F10 extensometric station.
Both the computed long-term compaction as well as the seasonal compaction/expansion are in good agreement with the measurements.
Figure 7 shows the vector maps of the cumulative horizontal displacement of the bottom of the fourth confined aquifer.The maximum values are calculated in the range between 4 and 8 mm.The largest horizontal displacements are obtained in the north-eastern part of the map, which corresponds to the zone with the largest changes in groundwater levels and also the highest compressibility.

Impacts of boundary conditions of 3-D geomechanical model
As reported above, two different outer boundary conditions are tested in the 3-D geomechanical model.One is a fixed boundary, and the other one is a traction-free boundary.ure 8 shows the simulated vertical and horizontal displacements as obtained implementing the two outer boundary conditions along a west-east profile through the center of the study area.The differences between the results are significant only in a relatively narrow strip along the domain boundary.groundwater flow and geomechanical models depending on the actual effective stress with respect to the preconsolidation stress.The model outcome satisfactorily agrees with the available observations.The cumulative land subsidence reaches 250 mm in the downtown area and the cumulative horizontal displacement peaks 11 mm on the ground surface and is smaller than 8.5 mm at depth.Fixed and traction-free conditions are tested on the outer boundary of the geomechanical model.The differences that remain limited to a relatively narrow strip close to the model frontier.The horizontal displacements appear more influenced by the boundary conditions than land subsidence.Although no measurement is available for the model validation in term of horizontal displacements, their magnitudes properly resemble the distribution of the hydraulic head changes and the compressibility of aquifers and aquitards.Hence, it can be concluded that the 3-D uncoupled geomechanical model can well simulate the three-dimensional displacement field of the aquifer system in the center area of Shanghai.

Figure 1 .
Figure 1.Location of the Shanghai City and the center area addressed in the presentmodelling study.

Figure 3 .
Figure 3. Perspective views of the 3-D finite-element mesh of the center area of Shanghai.The various geologic layers are distinguished in (a).In (b) the hydro-stratigraphic units are highlighted on the center part of the vertical section passing along the I-I' profile traced in Fig. 2.

Figure 4 Figure 4 .
Figure4compares the model results and the measured head in a few observation wells.The simulated heads agree generally well with the observed heads.Figures 5 shows the groundwater head maps in the second confined aquifer on 30 December 1995, in which the simulation outcomes agree well with the maps derived by interpolating the available records.

Figure 5 .
Figure 5. Computed (a) and measured (b) groundwater level in aquifer A2 on 30 December 1995, i.e. at the end of the simulation period.

Figure 7 .
Figure 7. Cumulative horizontal displacements at 30 December 1995, as provided by the 3-D geomechanical model in aquifer A4.

Figure 8 .
Figure 8.Comparison between the simulated vertical (a) and horizontal (b) displacements provided by the geomechanical model with different boundary conditions along a west-east profile crossing the downtown area.
one-way coupled geomechanical model is developed to simulate the groundwater level changes and 3-D displacements in the aquifer system at the center area of Shanghai due to groundwater withdrawal and recharge between 30 December 1979 and 30 December 1995.Elastic and plastic compressibility values are appropriately considered in the proc-iahs.net/372/443/2015/Proc.IAHS, 372, 443-448, 2015