Shui-Long Shen
email: slshen@sjtu.edu.cn

Lei Ma

Ye-Shuang Xu
email: xuyeshuang@sjtu.edu.cn

Zhen-Yu Yin

S.-L Shen

Y.-S Xu

Interpretation of increased deformation rate in aquifer IV due to groundwater pumping in Shanghai

Keywords: groundwater pumping, hydraulic gradient, aquifer IV deformation, shear stress, volumetric contraction pompage de l'eau souterraine, gradient hydraulique, déformation de l'aquifère IV, contrainte de cisaillement, contraction volumétrique

Shanghai is a typical example of a region in China suffering from severe land subsidence. Previous investigations have shown that land subsidence is strongly related to groundwater pumping in both the temporal and spatial domains. A series of field investigations shows that the deformation rate of aquifer IV (AqIV) has increased, contributing significantly to the total subsidence in Shanghai. To predict the deformation behavior of AqIV in Shanghai due to groundwater pumping, four different models have been applied in this study: (i) a model based on groundwater flow theory, (ii) an elastic model based on Cauchy mechanics, (iii) a creep model incorporated into Cauchy mechanics, and (iv) a recently proposed land subsidence model by Budhu based on Cosserat mechanics. Based on a number of assumptions, Budhu's concept was extended to consider the time-dependent behavior of groundwater potential in aquifers and the consolidation behavior of overburdened aquitards. The calculated results of AqIV's deformation were compared with the field-measured data. The results show that when using groundwater flow theory and Cauchy mechanics, it is difficult to predict the increased deformation rate of the sandy layers in aquifers even when the creep behavior is incorporated into the calculation. However, Budhu's method based on Cosserat mechanics can predict this increased deformation rate well. The predicted result indicates that increased deformation of aquifer IV is due to volumetric contraction caused by shear stress rather than the creep of aquifer sand. The phenomenon of deformation delay in aquifers can also be interpreted with the proposed model.

Introduction

Land subsidence is a geo-hazard characterized by sinking ground surfaces due to the compression of geological strata. In the coastal areas of China, land subsidence occurs in many places, including the Shanghai, Tianjin, Taiyuan, and Su-Xi-Chang areas [START_REF] Chai | Land subsidence due to groundwater drawdown in Shanghai[END_REF][START_REF] Shi | Characterization of land subsidence induced by groundwater withdrawals in Su-Xi-Chang area, China[END_REF][START_REF] Zhang | Stress-strain measurements of deforming aquifer systems that underlie Shanghai, China[END_REF][START_REF] Xu | The state of land subsidence and prediction approaches due to groundwater withdrawal in China[END_REF].

Recent research shows that groundwater pumping is still a primary contributor to land subsidence in Shanghai [START_REF] Shen | Numerical evaluation of land subsidence induced by groundwater pumping in Shanghai[END_REF]Xu et al. 2012aXu et al. , 2012b)). To mitigate the subsidence there, groundwater pumping has been limited both in its volume and by aquifer layers since 1966. The groundwater basin of Shanghai is a multi-aquifer-aquitard system (MAAS), which consists of a phreatic aquifer group (labeled as Aq0) and five confined aquifers (labeled as AqI to AqV in the following text) separated by six aquitards (labeled as AdI to AdVI). Since the 1970s, groundwater withdrawal has been conducted mainly in AqIII and AqIV (SGEAEB 2002;[START_REF] Shen | Numerical evaluation of land subsidence induced by groundwater pumping in Shanghai[END_REF].

According to the principles of soil mechanics, the consolidation of soft clay in aquitards occurs with increased effective stresses. Because the compressibility of clayey soil in overburdened aquitards is much greater than that of the sandy soil in the lower aquifer, it is generally agreed that the consolidation of aquitards is the main reason for land subsidence [START_REF] Thu | Modelling subsidence in the Hanoi City area, Vietnam[END_REF][START_REF] Phien-Wej | Land subsidence in Bangkok, Fig. 21. Relation between shear stress and parameter a in aquifer[END_REF][START_REF] Horpibulsuk | Assessment of engineering properties of Bangkok clay[END_REF][START_REF] Ortiz | Evolution of long-term land subsidence near Mexico City: review, field investigations, and predictive simulations[END_REF][START_REF] Shen | Numerical evaluation of land subsidence induced by groundwater pumping in Shanghai[END_REF]. Field observations also show that in most cases the compression of aquitards is much greater than that of aquifers with the same scale of thickness [START_REF] Gambolati | Mathematical simulation of the subsidence of Venice: 1[END_REF][START_REF] Su | Investigation on the deformation characteristics of soft deposit in Shanghai under pumping and recharge[END_REF][START_REF] Qian | Computation of land subsidence in Shanghai[END_REF][START_REF] Shen | Estimation of land subsidence based on groundwater flow model[END_REF][START_REF] Teatini | Groundwater pumping and land subsidence the Emilia-Romagna coastland, Italy: modeling the past occurrence and the future trend[END_REF][START_REF] Bell | Permanent scatterer InSAR reveals seasonal and long-term aquifer-system response to groundwater pumping and artificial recharge[END_REF][START_REF] Ortiz | Evolution of long-term land subsidence near Mexico City: review, field investigations, and predictive simulations[END_REF]. In Shanghai, it was found that the stress-strain of an aquifer occurs when it is under a state of elasticity and most of the deformation of the aquifer can be recovered after the groundwater level returns to its previous state [START_REF] Su | Investigation on the deformation characteristics of soft deposit in Shanghai under pumping and recharge[END_REF][START_REF] Qian | Computation of land subsidence in Shanghai[END_REF][START_REF] Shen | Numerical evaluation of land subsidence induced by groundwater pumping in Shanghai[END_REF]. For this reason, several prediction models on land subsidence were established to calculate the consolidation of soft clayey soil [START_REF] Shen | Estimation of land subsidence based on groundwater flow model[END_REF][START_REF] Teatini | Groundwater pumping and land subsidence the Emilia-Romagna coastland, Italy: modeling the past occurrence and the future trend[END_REF][START_REF] Bell | Permanent scatterer InSAR reveals seasonal and long-term aquifer-system response to groundwater pumping and artificial recharge[END_REF][START_REF] Xu | The state of land subsidence and prediction approaches due to groundwater withdrawal in China[END_REF][START_REF] Ortiz | Evolution of long-term land subsidence near Mexico City: review, field investigations, and predictive simulations[END_REF].

A series of recent investigations has found that the deformation rate in AqIV has increased since 1996 [START_REF] Zhang | Stress-strain measurements of deforming aquifer systems that underlie Shanghai, China[END_REF]); however, the groundwater level did not experience a drawdown. This new deformation phenomenon indicates that the aquifers, rather than the aquitards, have become the main compression layer, which is a phenomenon contrary to that found in earlier research [START_REF] Su | Investigation on the deformation characteristics of soft deposit in Shanghai under pumping and recharge[END_REF][START_REF] Qian | Computation of land subsidence in Shanghai[END_REF]. The groundwater flow model incorporating Terzaghi's one-dimensional (1D) consolidation theory was used to predict subsidence in Shanghai [START_REF] Chai | 1D analysis of land subsidence in Shanghai[END_REF][START_REF] Xue | Simulation of land subsidence in the south of Yangtze River Delta[END_REF][START_REF] Shen | Numerical evaluation of land subsidence induced by groundwater pumping in Shanghai[END_REF]. The predicted deformation of the aquifer was much lower than that of field-measured results. This phenomenon of increased deformation rate in AqIV was difficult to predict using the existing models because they were based on Cauchy continuum mechanics, in which only the isotropic consolidation of the unit cell using the volumetric compressibility of soil was considered. [START_REF] Budhu | Mechanics of land subsidence due to groundwater pumping[END_REF] proposed land-subsidence mechanics based on [START_REF] Cosserat | Theorie des corps deformables[END_REF] continuum mechanics to simulate the behavior of shear and torsion effect on sand due to the hydraulic gradients in unconfined aquifers (in the following text, this method is referred to as Budhu's approach). This model analyzes the occurrence of shear stress and torsion in aquifers due to large hydraulic gradients. In Shanghai, groundwater pumping wells in aquifers are not distributed uniformly, and a large hydraulic gradient exists in many places [START_REF] Shen | Numerical evaluation of land subsidence induced by groundwater pumping in Shanghai[END_REF][START_REF] Tan | Measured performance of a 26 m deep top-down excavation in downtown Shanghai[END_REF]Xu et al. 2012aXu et al. , 2012bXu et al. , 2013aXu et al. , 2013b)). This hydraulic gradient may cause unbalanced shear stress in confined or unconfined aquifers and the shear stress increases with the increase of horizontal distance. The deformation caused by the imbalance of shear stress constitutes a large part of the total deformation of aquifers [START_REF] Budhu | Mechanics of land subsidence due to groundwater pumping[END_REF][START_REF] Budhu | Earth fissure formation from the mechanics of groundwater pumping[END_REF].

The objectives of this study are the following: (i) to predict the deformation of confined aquifers in Shanghai using different existing models (e.g., a groundwater flow model, classic Cauchy continuum mechanics considering both the elastic and creep behavior of soil, and Budhu's approach), (ii) to interpret the mechanism that increases the deformation rate of the aquifers, and (iii) to investigate the reasons for deformation delay with the drawdown of groundwater potential.

Background

Geological profile and soil properties

Shanghai is located on the deltaic deposits of the Yangtze River on the eastern coast of China. The soil stratum in Shanghai is a typical deltaic soft deposit with MAAS. The groundwater basin is composed of six confined aquifers and six aquitards. Figure 1 shows the locations of the investigated areas: Laodong Park and Shuangyang School, which are located near the Yangtze River in eastern Shanghai. In recent decades, a large amount of groundwater was pumped continuously from these areas for domestic and industrial use. As depicted in the figure, there was a drawdown cone of groundwater from a depth of -35 to -50 m in the western part of Laodong Park. This area is analyzed in this study. The soil profile within a 400 m depth from Laodong Park to Baoshan is depicted in Fig. 2 (for cross section see Fig. 1). The subsoil within 400 m is separated into 10 layers: 5 aquifers and five aquitards. Beneath the soft soil is the base rock, which extends to great depths.

The soil profile and properties of the soft subsoil at this site, within a 300 m depth of the borehole near Laodong Park, are depicted in Fig. 3. As shown in Fig. 3, the water content is close to the liquid limit in the upper layers and decreases to the plastic limit in lower layers. The unit weight is approximately 18 kN/m 2 in the uppermost layer and increases to 20 kN/m 2 in the lowest layer [START_REF] Tan | Measured performance of a 26 m deep top-down excavation in downtown Shanghai[END_REF][START_REF] Tan | Observed behaviors of a long and deep excavation constructed by cut-and-cover technique in Shanghai soft clay[END_REF]Tan andWang 2013a, 2013b). With the exception of the top soft soil, the void ratios of soft soil within a 300 m depth are less than 1.0, and less than 0.7 within a 100 m depth. The variation pattern of the compression index is similar to the void ratio. Also shown in Fig. 3 is the hydraulic conductivity, which is greater than 10 -5 m/s in the aquifer, which is three orders of magnitude greater than that of the aquitards [START_REF] Chai | Performance of embankments with and without reinforcement on soft subsoil[END_REF][START_REF] Shen | Analysis of field performance of embankments on soft clay deposit with and without PVDimprovement[END_REF]Huang and Han 2009;Huang et al. 2009;[START_REF] Xu | Geological and hydrogeological environment in Shanghai with geohazards to construction and maintenance of infrastructures[END_REF].

Groundwater level and deformation of AqIV

Groundwater pumping for domestic and industrial use has taken place in Shanghai since the 1890s. It has been strictly limited since 1966 [START_REF] Chai | Land subsidence due to groundwater drawdown in Shanghai[END_REF][START_REF] Gong | Analysis on land subsidence due to construction engineering in soft soil region of Shanghai[END_REF][START_REF] Shen | Numerical evaluation of land subsidence induced by groundwater pumping in Shanghai[END_REF]Xu et al. 2012a;Xu et al. 2013aXu et al. , 2013b)). Figure 4 shows the amount of groundwater withdrawn in recent decades. Each aquifer is subject to different rates of groundwater pumping. The volume of groundwater withdrawn reached 90 × 10 6 m 3 /year in 1995 in AqIV. The uncontrolled groundwater pumping caused soil deformation in aquitards as well as in aquifers.

Based on field data presented by the Shanghai government in 2002, the percentage of average compression for each aquifer and the sum of the compressions of the aquitards contributing to the total subsidence at Laodong Park is shown in Fig. 5. From 1980 to 1995, the compression of all aquitards was responsible for approximately 50% of the total land subsidence and the compression of AqIV was responsible for approximately 32%. However, during the period from 1996 to 2000, the contributing effect of the compression of AqIV increased to approximately 50% of the total subsidence. Therefore, there was a phenomenon of increased deformation rate in AqIV.

Figure 6 shows the variation in groundwater level and deformation in AqIV from 1980 to 2000 at Laodong Park. The groundwater level has lowered from -15 to -35 m in the past two decades. There was a dramatic decline in the groundwater level at this site beginning in 1985, and reached its lowest point in 1997. Figure 6 also shows deformation of AqIV at Laodong Park. In 2000, the cumulative deformation of AqIV reached 300 mm, and it continues to increase today. The variation pattern of deformation was different from the variation pattern of the groundwater level. The deformation from 1980 to 1990 was very small and fluctuated within 10 mm. However, the deformation increased dramatically after 1990 and continued increasing for the next 10 years. This demonstrates that there was a significant delay of deformation after reduction in groundwater levels.

Figure 7 shows the variation in groundwater level and deformation in AqIV from 1980 to 2010 at Shuangyang School. As shown in this figure, the groundwater level lowered from -15 to -35 m in 20 years. After 2000, the groundwater level began to recover by the process of groundwater recharging and reached -20 m by 2010. The deformation rate at this site was similar to that at Laodong Park. By 2010, the cumulative deformation of AqIV reached 160 mm at Shuangyang School. The deformation rate increased with the decline of the groundwater level before 1997. Since 1997, although the groundwater level has begun to recover, the defor- [START_REF] Chai | Land subsidence due to groundwater drawdown in Shanghai[END_REF]). Aq, aquifer; Ad, aquitard; w n , water content ratio; w l , liquid limit; w p , plastic limit; ␥ t , unit weight; e 0 , void ratio; C c , compression index; k v , vertical permeability. mation continues to increase, demonstrating a clear delay in deformation.

Based on the field-observed data on groundwater pumping, groundwater level, and deformation in AqIV, observations can be summarized as follows:

1. Although the volume of groundwater withdrawn from AqIV has been dramatically reduced, the deformation of AqIV continues to increase; 2. After 1995, although the groundwater level had begun to recover with the reduction of withdrawn volume, the deformation of AqIV continued to increase; and

3. The groundwater level is related to the withdrawn volume; however, deformation shows a delayed behavior.

Existing interpretations

Most prediction models of land subsidence were established by incorporating a consolidation model into a groundwater flow model, in which ground compaction is calculated by considering the compressibility of soil in both aquifers and aquitards [START_REF] Gambolati | Mathematical simulation of the subsidence of Venice: 1[END_REF][START_REF] Shen | Estimation of land subsidence based on groundwater flow model[END_REF][START_REF] Teatini | Groundwater pumping and land subsidence the Emilia-Romagna coastland, Italy: modeling the past occurrence and the future trend[END_REF][START_REF] Bell | Permanent scatterer InSAR reveals seasonal and long-term aquifer-system response to groundwater pumping and artificial recharge[END_REF][START_REF] Ortiz | Evolution of long-term land subsidence near Mexico City: review, field investigations, and predictive simulations[END_REF]. The consolidation process can be calculated using both Terzaghi's 1D consolidation theory [START_REF] Gambolati | Mathematical simulation of the subsidence of Venice: 1[END_REF][START_REF] Shen | Estimation of land subsidence based on groundwater flow model[END_REF]) and Biot's threedimensional (3D) consolidation theory [START_REF] Kim | Three-dimensional numerical simulation of fully coupled groundwater flow and land deformation in unsaturated true anisotropic aquifers due to groundwater pumping[END_REF][START_REF] Teatini | Groundwater pumping and land subsidence the Emilia-Romagna coastland, Italy: modeling the past occurrence and the future trend[END_REF]. The compaction of soil can be calculated by using the coefficient of compressibility, which considers the deformation behavior of poroelastic [START_REF] Shen | Estimation of land subsidence based on groundwater flow model[END_REF]) and plastic soils [START_REF] Gambolati | Mathematical simulation of the subsidence of Venice: 1[END_REF], and soil creep [START_REF] Yin | Equivalent times and one-dimensional elastic viscoplastic modelling of time-dependent stress-strain behaviour of clays[END_REF][START_REF] Yin | Elastic visco-plastic modelling of onedimensional consolidation[END_REF][START_REF] Yin | Non-linear creep of soils in oedometer tests[END_REF][START_REF] Gu | A 3-D coupled model with consideration of rheological properties[END_REF][START_REF] Yin | A new elastic viscoplastic model for time-dependent behaviour of normally and overconsolidated clays: theory and verification[END_REF][START_REF] Xue | Simulation of land subsidence in the south of Yangtze River Delta[END_REF][START_REF] Yin | Identifying parameters controlling soil delayed behaviour from laboratory and in situ pressuremeter testing[END_REF]Yin et al. 2010aYin et al. , 2010b[START_REF] Yin | Modeling time-dependent behavior of soft sensitive clay[END_REF]. Therefore, based on traditional soil mechanics, the deformation of aquitards is much greater than that of aquifers [START_REF] Su | Investigation on the deformation characteristics of soft deposit in Shanghai under pumping and recharge[END_REF][START_REF] Qian | Computation of land subsidence in Shanghai[END_REF][START_REF] Chai | Land subsidence due to groundwater drawdown in Shanghai[END_REF][START_REF] Chai | 1D analysis of land subsidence in Shanghai[END_REF][START_REF] Shen | Estimation of land subsidence based on groundwater flow model[END_REF][START_REF] Shen | Numerical evaluation of land subsidence induced by groundwater pumping in Shanghai[END_REF].

To help clarify the contradiction between the aforementioned deformation phenomenon of AqIV and the concept of traditional soil mechanics, the following three mechanisms are presented: (i) additional load from highrise buildings [START_REF] Tang | Model test study of land subsidence caused by high-rise building group in Shanghai[END_REF]), (ii) creep behavior of sandy soil [START_REF] Xue | Simulation of land subsidence in the south of Yangtze River Delta[END_REF], and (iii) impact of seepage-blocking underground structures in an aquifer (Xu et al. 2012a(Xu et al. , 2012b(Xu et al. , 2013b)). Generally, high-rise buildings have induced settlement concentrated in a narrow region (e.g., a width of the building's foundation). At the investigated sites (Laodong Park and Shuangyang School), the maximum height of the building was 24 stories with a distance of 200 m from the subsidence measuring point. This can be corroborated using Google Maps (http:// maps.google.com.hk/). Moreover, pile foundations were founded on top of AqI at a depth of approximately 45 m [START_REF] Chen | Field tests, modification, and application of deep soil mixing method in soft clay[END_REF]Shanghai Statistics 2011;Tan andWang 2013a, 2013b), which is more than 130 m above the top of AqIV. The compression of the soil layers between the pile tip and AqIV was less than that of AqIV. Moreover, the compression of AqII and AqIII was much lower than that of AqIV from 1996 to 2010. As shown in Fig. 8, the compression of AqII and AqIII correlates to the groundwater level. Therefore, we can confirm that high-rise buildings have had no effect on the deformation of AqIV since 1995. As for the impact of seepage-blocking underground structures in an aquifer, Xu et al. (2012b) proved that seepage blocking occurred mostly in AqI and much less in AqII. There were no cases of seepage blocking in AqIII and AqIV. Deformation caused by the creep behavior of sand soil required further investigation in this study.

Prediction methods

To interpret the mechanism of the aforementioned phenomenon, simulation of the subsidence was conducted using the fol-

Groundwater flow model

In a groundwater flow model, a traditional two-step approach (Gambolati and Freeze 1973) is adopted. In the first step, the hydrodynamics of the pumped multi-aquifer system is simulated by a 3D groundwater flow model. Then, the subsidence (compaction) of soil strata is computed with the aid of 1D consolidation theory with the pore pressure field specified as an externally distributed source of strength within the porous medium [START_REF] Shen | Estimation of land subsidence based on groundwater flow model[END_REF][START_REF] Shen | Numerical evaluation of land subsidence induced by groundwater pumping in Shanghai[END_REF].

The aquifer hydrodynamics are governed by the classic groundwater flow equation (e.g., [START_REF] Bear | Dynamics of fluids in porous media[END_REF][START_REF] Gambolati | Mathematical simulation of the subsidence of Venice: 1[END_REF] (1)

ٌ(kٌh) Ϫ q ϭ S s Ѩh Ѩt
where k is hydraulic conductivity, h is hydraulic head, q is the external source-sink flux, S s is the coefficient of specific storage, and t is time.

The release of water in soil is normally due to compression of the soil and expansion of the water. Therefore, the coefficient of specific storage, S s , is expressed in the following way:

(2)

S s ϭ ␥ w (m v ϩ n␤ w )
where ␥ w is the unit weight of water, m v is the soil coefficient of volume compressibility, n is porosity, and ␤ w is the coefficient of volume compressibility of water. Therefore, we can calculate the compression of a soil layer with varying hydraulic head from a specific storage, S s . The total amount of land subsidence can be calculated by accumulating the compression of soil layers. For details on use of this method, refer to [START_REF] Shen | Estimation of land subsidence based on groundwater flow model[END_REF] and [START_REF] Teatini | Groundwater pumping and land subsidence the Emilia-Romagna coastland, Italy: modeling the past occurrence and the future trend[END_REF].

Figure 9 shows the 3D finite element method (FEM) calculation model and the calculated range in the groundwater flow model, which were set as 2 km horizontally and 400 m vertically. The displacement and drainage boundary conditions were set as the following: at the bottom both vertical and horizontal displacements were restrained, and the bottom line was set as undrained. The left and right vertical boundaries were set as undrained, and their horizontal displacements were fixed. The ground surface was set as drained. The parameters of soil used in the water flow model are listed in Table 1 [START_REF] Shen | Numerical evaluation of land subsidence induced by groundwater pumping in Shanghai[END_REF]. The drawdown of the groundwater level was set at the pumping well locations with the groundwater level and time relation shown in Figs. 67. That is, in each time step of the calculation the groundwater level was lowered at the well location as a boundary condition and the redistribution of groundwater level in the surrounding computational domain was calculated by eq. ( 1).

Elastoplastic model based on Cauchy mechanics

The classical Cauchy mechanics were also adopted. For simplicity, this analysis was conducted with the FEM under axisymmetric stress conditions using the software package Plaxis [START_REF] Brinkgreve | PLAXIS-2D Version 8 user's manual[END_REF]). The calculated range and boundary conditions used in Plaxis were the same as those used in the groundwater flow model. The aquitards were represented by the soft soil model in Plaxis, and the aquifers were represented by an elastoplastic model based on the Mohr-Coulomb criterion [START_REF] Brinkgreve | PLAXIS-2D Version 8 user's manual[END_REF]). The drawdown of groundwater level was set at the pumping well locations, with the groundwater level and time relation (Figs. 6 and7) set in the same way as for the groundwater flow model. The redistribution of groundwater level in the surrounding computational domain was calculated by the same procedure as shown in eq. ( 1) [START_REF] Brinkgreve | PLAXIS-2D Version 8 user's manual[END_REF]), but in a two-dimensional (2D) condition.

The parameters of soft clay used in Plaxis are listed in Table 2 (SGEAEB 2002;[START_REF] Shen | Numerical evaluation of land subsidence induced by groundwater pumping in Shanghai[END_REF]. The Young's modulus E used in Plaxis is based on the following equation:

(3) where E s is compression modulus. In solid mechanics, the relation between E and E s is

E s ϭ 1/m v
(4) E ϭ E s [1 Ϫ 2 2 /(1 Ϫ )]
where v is Poisson's ratio.

Model considering the creep behavior of sandy soil

In this study, the creep model was adopted for comparison. The creep behavior of sandy soil was modeled by incorporating the elastic viscoplastic behavior proposed by [START_REF] Yin | Elastic viscoplastic modelling of the timedependent stress-strain behaviour of soils[END_REF], [START_REF] Yin | Elastic viscoplastic consolidation modelling and interpretation of pore-water pressure responses in clay underneath Tarsiut Island[END_REF], and Yin et al. (2010aYin et al. ( , 2010c[START_REF] Yin | Modeling time-dependent behavior of soft sensitive clay[END_REF]. According to the model of Yin et al. (2010a[START_REF] Yin | Modeling time-dependent behavior of soft sensitive clay[END_REF], the total strain rate is additively composed of the elastic strain rates and viscoplastic strain rates. The elastic behavior is assumed to be isotropic, similar to the Modified Cam Clay model. The viscoplastic strain rate, ˙ij vp , is governed by the following equation:

(5) ˙ij

vp ϭ ͩ p m d p m r ͪ N Ѩf d Ѩ ij
where is referred to as the fluidity parameter, p m d is the size of the dynamic loading surface, p m r is size of the static yield surface, N is the strain-rate coefficient, f d is treated as a viscoplastic potential function, and ij is the effective stress tensor. In this model, an elliptical surface is adopted to describe the dynamic loading surface and the reference surface (see Yin et al. 2010a). In addition, and N can be obtained from the secondary compression coefficient, C ␣e , defined in e-lnt space:

(6) ϭ C ␣e (1 ϩ e 0 ) c ; N ϭ Ϫ C ␣e
where c is the reference time taken as equal to 1 day. By considering creep strain in this model, the viscous volumetric strain rate is obtained from the flow rule of eq. ( 5) using the secondary compression coefficient. The model parameters include parameters of the Modified Cam Clay model (Poisson's ratio, ; compression index, ; swelling index, ; preconsolidation pressure, p0 , slope of critical-state line, M) and the secondary compression coefficient, C ␣e .

Budhu's approach based on Cosserat mechanics

To predict the deformation of AqIV in Shanghai, a new approach proposed by [START_REF] Budhu | Mechanics of land subsidence due to groundwater pumping[END_REF] based on Cosserat continuum mechanics was introduced in this study. According to this approach, when there is a hydraulic gradient in the spatial domain, unbalanced forces will occur, which cause shear deformation in a vertical direction. Figure 10 shows the two elements of deformation caused by increasing effective stress. Based on Budhu's approach, the changes in effective stress due to the decreasing groundwater level are decomposed into two parts. Part I is an isotropic effective stress, which leads to isotropic compression. Part II is a vertical effective stress varying with lateral distance, which causes simple shear and rotation. In Budhu's theory, as it is assumed that the shear deformation occurs under the elastic condition, the shear strain has its maximum value at the top of the aquifer layer and decreases linearly with depth to zero at the bottom of the aquifer layer. However, this new approach was derived under the conditions of an unconfined aquifer. Because the stress condition in a confined aquifer is different from that in an unconfined aquifer, some hypotheses are assumed to extend the applicable range of Budhu's approach.

Figure 11a shows the depression cone of the groundwater level and the soil consolidation condition in the confined aquifer with overburdened aquitards. When groundwater in the confined aquifer is pumped out, a cone of depression of groundwater level will form around the pumping well in the aquifer. Unlike in an unconfined phreatic aquifer, groundwater pumping in a confined aquifer will cause groundwater replenishment from both the aquifer itself and the overburdened aquitards, leading to the consolidation of the overburdened aquitards. Because the hydraulic conductivity of overburdened aquitards is significantly less than that of an aquifer (by more than three orders), the amount of replenishment of groundwater from aquitard to aquifer is very low during the times when groundwater is pumped from the aquifer. Consequently, the replenishment of groundwater occurs mainly from the aquifer over a large horizontal distance from the pumping well, which causes a hydraulic gradient around the well.

Figure 11b depicts the stress state at the boundary line between aquifer and aquitard. Because the total overburden stress () remains constant before and after the decline of groundwater level, there is reaction stress in the aquifer to balance the overburdened total stresses. According to the principle of effective stress, the resistance to in an aquifer is provided by a combination of the stresses transferring through soil particles, i.e., the effective stress (=) and pore-water pressure caused by water in the pores (u). During groundwater pumping in the confined aquifer, the load from overburdened aquitards remains unchanged, which means that the imposed total stress on the confined aquifer remains the same. Because of groundwater pumping, the drawdown of groundwater head decreases the pore pressure in an aquifer so that the effective stress transferring through soil particles increases to keep the total stress in balance. These increases in effective stress in soil particles vary across the spatial domain with the drawdown cone of groundwater, and cause shear stresses in the aquifer [START_REF] Budhu | Mechanics of land subsidence due to groundwater pumping[END_REF][START_REF] Budhu | Earth fissure formation from the mechanics of groundwater pumping[END_REF], which may lead to an increased tendency toward deformation.

To apply Budhu's approach in a confined aquifer, the following additional hypotheses for confined aquifers are presented:

1. A confined, isotropic, homogeneous aquifer overburdened by an aquitard with much lower hydraulic conductivity than that of the underlying confined aquifer (e.g., 1/10th); 2. Variation of groundwater replenishment from an aquitard that depends on the consolidation rate of overburdened aquitards; 3. An aquifer that is a continuous medium without cracks; 4. An overburdened load from aquitards above an aquifer that remains unchanged during groundwater pumping in confined aquifers; 5. Groundwater potential at its initial state in a confined aquifer that is the same in the horizontal direction; 6. Groundwater level that declines slowly with pumping time.

To calculate the time-dependent stress state in a confined aquifer using Budhu's concept and the aforementioned assumptions, the distribution of groundwater level is assumed to be characterized by the following equation:

(7) h w (t) ϭ Ϫ͓␣ a x a ϩ h 0 (t)͔
where h w (t) is the drawdown of the groundwater head with time from the initial head surface at the pumping well, ␣ a is a soil layer constant, x is the distance from the origin, a is the time factor, and h 0 (t) is the drawdown of the groundwater head with time from the initial head surface at an arbitrary selected origin.

If parameter a equals 2, the groundwater level equation becomes the original equation proposed by [START_REF] Budhu | Mechanics of land subsidence due to groundwater pumping[END_REF], which denotes a stable state. In the field, the groundwater level gradually decreased with groundwater pumping. There was a series of groundwater level functions, and this gradual change was expressed by the time factor a. For simplicity, the relationship between the time factor, a, and t is assumed to be characterized by the following equation:

(8) a ϭ Ϫ lnt where = is a parameter representing the time period of drawdown from initial to stable state and t is time (unit: year). The determination method of = can be based on Fig. 12 from a time period of stable groundwater levels.

Based on this assumption, the time factor, a, can correspond directly to time, t. Figure 13 shows the relationship between the time factor, a, and t with the variation of the parameter = when the groundwater level reaches a stable state from its initial state. As shown in Fig. 13, the time factor, a, increases with parameter =. For example, when keeping the time period ⌬t (= = 4.303) equal to 10, a drops quickly at the beginning of pumping, which also means the groundwater level near the pumping well drops sharply. After 5 years, the reduction ratio of a decreases, signifying that the depression cone of the hydraulic head is enlarged.

Equation ( 9) shows the effective stress equation used in the study of [START_REF] Budhu | Mechanics of land subsidence due to groundwater pumping[END_REF]. Unlike the situation in an unconfined aquifer, the sand in a confined aquifer must be saturated. The value for ␤, which represents the degree of saturation, 8) from stable time period of groundwater level. 9) is always 1 in confined aquifers. To simulate the variation rate of groundwater replenishment from overburdened aquitards, b is introduced in eq. ( 10). Therefore, eq. ( 10) represents the situation of confined aquifers with a similar form to eq. ( 9) to represent the consolidation of overburdened aquitards and the time-dependent drawdown of the groundwater level in the confined aquifer ( 9)

¯ij ϭ ␤ ͑ 1 V ͒ ͵ V ij dV (10) ¯ij ϭ b ͑ 1 V ͒ ͵ V ij dV
where = is the effective stress, ¯ij is the average effective stress tensor, ␤ is a factor related to the degree of saturation and the porosity, b is the groundwater replenishment factor with value between 0 and 1, and V is the characteristic volume of the element.

In this study, b is assumed to increase linearly with time. At the beginning of pumping, when there is a large hydraulic gradient between aquifer and aquitard, the groundwater replenishment is high and the effective stress in the aquifer increases slowly so that b equals 0. After a period of time, once groundwater motion becomes a stable seepage and the soil completes consolidation in the aquitard, the effective stress in the aquifer reaches its maximum value; in this case, parameter b equals 1.

The groundwater level drawdown with time h 0 (t) at the well location is shown in Figs. 6 and7. To follow a similar condition of redistribution of groundwater level in the calculated domain as shown in the previous three methods, the groundwater level h w (t) was calculated by eqs. ( 7) and ( 8) at each time step. The detailed procedure for each case is described as follows.

The groundwater level h 0 (t) with time used in the calculation for Laodong Park is depicted in Fig. 6. This shows that from 1980 to 1985, there was no significant change in the groundwater level. Therefore, this value was selected as the initial groundwater level; after 1985, the groundwater level decreased linearly until the end of 2000, reaching its lowest point at -35 m. Therefore, the time period of drawdown from initial to stable state was approximately 15 years with the corresponding parameter = = 4.708 (which was used in this case).

Figure 7 presents the measured and simulated values of the groundwater level h 0 (t) with time at Shuangyang School. In a manner different from the case at Laodong Park, the groundwater level at Shuangyang School has recovered since 1997. For simplicity in the simulation model, the groundwater level after 1997 was set as a stable groundwater level at the average value of -30 m. Therefore, the time period of drawdown from the initial to the stable state was approximately 12.5 years with the corresponding parameter = = 4.35 (which was used in this case).

The other parameters used in the extended Budhu model implemented at both Laodong Park and Shuangyang School are the following: (i) the character range of the model is 1 km; (ii) ␣ a = 0.5 × 10 -5 m -1 (Laodong Park) and ␣ a = 0.1 × 10 -5 m -1 (Shuangyang School); (iii) the thickness of AqIV at these two sites, l z , is 60 m; (iv) the aquifer soil is elastic and isotropic with E = 144 MPa, G = 70 MPa, and G c = 35 MPa; and (v) the unit weight of water, ␥ w , is 10 kN/m 3 .

Results

Laodong Park

The deformation and groundwater level of AqIV has been continually measured at the site of Laodong Park throughout the past 20 years. Figure 14 plots the measured and calculated deformation curve of AqIV around Laodong Park, using four methods. As shown in Fig. 14, the extended Budhu model could predict defor-mation of AqIV well. Of the predicted values, the groundwater flow model produced the lowest one, the elastic Cauchy mechanics model using the Plaxis package produced the second-lowest, and the creep model using the Plaxis package produced the thirdlowest. All three methods can model the deformation before 1993 well. However, after 1995, there are large discrepancies between the measured values and the predicted values produced by these three methods. The solid line calculated using Budhu's method based on Cosserat mechanics can model the measured value well from 1980 to 2000.

Figure 15 shows the cross-sectional view of deformation over 20 years calculated from Budhu's method. With a nonuniform arrangement of wells, a large hydraulic gradient was found in AqIV. This hydraulic gradient caused uneven deformation in AqIV. As shown in Fig. 15, the deformation of AqIV developed with time. At the beginning of groundwater pumping, the deformation was very small and distributed uniformly across the seepage domain. With the development of a hydraulic gradient in groundwater level, the deformation across the seepage domain became unevenly distributed.

Shuangyang School

The deformation and groundwater level of AqIV has also been measured continually at the site of Shuangyang School throughout the past 30 years. Figure 16 plots the measured and calculated deformation curves of AqIV using four methods, around Shuangyang School. As shown in Fig. 16, the variation tendency at this site using the four different methods was the same as that at Laodong Park. Because the decline in the groundwater level at this site was less than that at Laodong Park, the deformation in AqIV at Shuangyang School was initially low. Therefore, all four methods can model the deformation before 1997 well. However, after 2000, there are large discrepancies between the measured and predicted values produced by three of the methods (with the exception being Budhu's method). The solid line calculated using Budhu's concept based on Cosserat mechanics could predict the measured value well from 1980 to 2010.

Figure 17 shows the cross-sectional view of deformation throughout the last 30 years calculated using Budhu's method. As shown in Fig. 17, the increased deformation rate in AqIV at this site is the same as that at Laodong Park. When groundwater pumping began, the deformation was very small and distributed uniformly across the seepage domain. With the development of a hydraulic gradient of groundwater level, the deformation across the seepage domain became unevenly distributed.

Variation of shear stress and volumetric contraction

From the aforementioned analyses of the cases at Laodong Park and Shuangyang School, it can be deduced that the classic models (e.g., a groundwater flow model, elastic Cauchy mechanics, and a creep model) would underestimate deformation of the confined AqIV since 1995. To understand the reason for this underestimation with the classic models, the effect of shear stress caused by the hydraulic gradient of groundwater level required investigation using Budhu's method. Figure 18 shows the variation of calculated shear stress in AqIV at Laodong Park. From this figure, we find that the shear stress, xz , equals zx at a distance of ͙2l z according to the results proposed by [START_REF] Budhu | Mechanics of land subsidence due to groundwater pumping[END_REF]. Beyond ͙2l z , xz is larger than zx . As the distance from the pumping well increases, xz increases at a faster rate and reaches 200 kPa one kilometre (1000 m) away. This large shear stress may be the primary contributor to the increased tendency toward deformation in AqIV.

The existence of shear stress can induce not only a change of shape, but also the volumetric strain of granular material [START_REF] Rowe | The stress-dilatancy relation for static equilibrium of an assembly of particles in contact[END_REF][START_REF] Chang | Modeling stress-dilatancy for sand under compression and extension loading conditions[END_REF][START_REF] Yin | Stress-dilatancy behavior for sand under loading and unloading conditions[END_REF]. Many laboratory investigations [START_REF] Rowe | The stress-dilatancy relation for static equilibrium of an assembly of particles in contact[END_REF][START_REF] Yin | Stress-dilatancy behavior for sand under loading and unloading conditions[END_REF] show that when the shear strain is small (less than 0.5%), the shear stress induces a volumetric contraction in soil [START_REF] Ishihara | Fundamentals of soil dynamics[END_REF]). For the two cases investigated in this study the total strain is about 0.52%, thus the shear strain is less than 0.5% and the shear stress induces contraction of aquifer soil. Thus, in the field the deformation of AqIV in Shanghai may be considered due to the contraction induced by shear stress rather than due to the creep effect.

Discussion on the deformation delay

As presented in Figs. 6 and7, there was a deformation delay phenomenon with the variation of groundwater level at both Laodong Park and Shuangyang School. To help clarify this phenomenon, the data on the groundwater level and deformation in AqIV at both sites are presented again in Fig. 19. The time period for the measurements was from 1981 to 2010. As shown in this figure, the behavior of aquifer deformation with the decline of the groundwater level occurred in the same way at these two sites. During the initial 10 years, the groundwater level declined from -7.5 m to approximately -30 m. Concurrently, the total deformation changed from 20 to 40 mm. After 1990, groundwater pumping in AqIV had been strictly limited by reducing the number of pumping wells, and the groundwater level fluctuated around -35 m. However, the deformation quickly increased, and even after 2000, the rate of deformation was not yet convergent. Although the groundwater level began to recover after 1997 at Shuangyang School, deformation of the aquifer continued to increase. Therefore, we have found that significant delays in deformation with respect to variations in the groundwater level in AqIV have occurred and are continuing today.

This deformation delay with groundwater level in the aquifer implies that the deformation of sand does not instantly change with the groundwater level. Because the shear deformation and macrorotation may constitute a major part of the total deformation according to Budhu's concept, the reason for deformation delay may be due to the progressive accumulation of shear stress with time. [START_REF] Xue | Simulation of land subsidence in the south of Yangtze River Delta[END_REF] attributed the deformation delay of aquifers to the effect of soil creep. However, sand creep is much smaller than that of soft clay from a short-term perspective [START_REF] Yin | Modeling time-dependent behavior of soft sensitive clay[END_REF]. The calculated results of this study show that creepinduced deformation is not the primary contributor to AqIV's total deformation. Based on the laboratory investigation [START_REF] Nawir | Experimental evaluation of the viscous properties of sand in shear[END_REF], the occurrence of creep for sand requires a high shearstress level, which is difficult to interpret using a classic Cauchy continuum model. Moreover, in the studies proposed by [START_REF] Xue | Simulation of land subsidence in the south of Yangtze River Delta[END_REF] and several other researchers, the reason for the presence of a high shear-stress level is not investigated. In Shanghai, groundwater-pumping wells in aquifers are not distributed uniformly [START_REF] Shen | Numerical evaluation of land subsidence induced by groundwater pumping in Shanghai[END_REF]Xu et al. 2012a). After limiting pumping since 1990, many pumping wells were closed. The reduction of pumping wells may have caused a large hydraulic gradient in many places [START_REF] Shen | Numerical evaluation of land subsidence induced by groundwater pumping in Shanghai[END_REF]Xu et al. 2012a). Based on the calculated results of the extended Budhu's model, a hydraulic gradient can induce a very high unbalanced shear stress in AqIV. This shear stress can induce contraction (i.e., negative dilatancy), which is the primary contributor to the significant deformation of AqIV rather than creep of aquifer sand.

At the beginning of groundwater pumping, the groundwater level decreased rapidly around the well, as shown in Fig. 20. However, the area affected by groundwater level reduction was not large, as there was only a small depression cone. At that moment, the isotropic compression stress was large in the aquifer layer, and the shear stress was relatively low. With the exception of particle crushing, deformation of sand under isotropic compression stress is generally relatively low. According to the investigation by [START_REF] Miura | Stress-strain characteristics of sand in a particle-crushing region[END_REF], sand particle crushing requires a very high compression stress level (e.g., 10 MPa). The compression stress level induced by groundwater pumping in AqIV was much lower than that necessary to cause particle crushing so that the state of the sandy soil in AqIV was also far from yielding. Therefore, sand compression was minimal at this initial stage. After a period of groundwater level drawdown, the drawdown cone became enlarged. With the gradient of groundwater level covering a large range, shear stress in the aquifer layer increased according to the proposed model, which may have caused shear stress in the aquifer sand layer. This shear stress can induce volumetric contraction in the small strain, e.g., less than 0.5%. Moreover, this volumetric contraction of sand may have accumulated with the accumulation of shear stresses. This cumulative high shear stress level may have caused the creep effect of the sand. Although the groundwater level began recovering, the contraction kept increasing for a period of time due to the residual shear stress in the sand layer. The existence of shear stress in sand may be one of the significant reasons for the phenomenon of deformation delay in the aquifers. The variation process of shear stress is discussed in detail in the following text based on the proposed model.

To clarify the cumulative process of shear stress in aquifers, a hypothetical case is presented. The parameters used are the following: (i) the length, l x , is 1 km; (ii) h 0 = 0 m, ␣ a = 1.5 × 10 -5 m -1 , and b = 1; (iii) thickness of aquifer, l z , is 100 m. Figure 21 depicts the variation of shear stress with the a value and the distance ratio in this case. As shown in Fig. 21, zx decreases linearly as the distance from the pumping well increases and remains the same as the groundwater level decreases. However, xz changes with a, especially at a location near the pumping well. xz increases as time increases, but decreases as the time factor, a, increases. Therefore, xz accumulates with the formation of a drawdown cone created due to groundwater pumping. This process of shear stress accumulation in an aquifer may also cause the creep of sand. The high shear stress level might induce cracks and (or) earth fissures in the aquifer layers [START_REF] Budhu | Earth fissure formation from the mechanics of groundwater pumping[END_REF]) so that the cumulative shear stress might not reach 400 kPa, as shown in Fig. 21. The occurrence of cracks changes the hydraulic conductivity of both the aquitard and the aquifer, which may change the deformation pattern underground.

Conclusions

This paper presents a case history of the phenomenon of increased deformation rate in AqIV in Shanghai. To interpret this deformation phenomenon, four different modeling approaches were implemented: a groundwater flow model, an elastic model that have exhibited increased deformation of their confined aquifers have been analyzed using the aforementioned four models. The predicted results show that using the classic models (groundwater flow model, elastic model based on Cauchy mechanics, and creep model), it was difficult to predict the field-measured data after 1995. However, Budhu's method based on Cosserat mechanics can reasonably predict the increased rate of deformation of aquifer IV in Shanghai. 4. The explanation for the reasonable prediction of aquifer deformation is as follows. Pumping in an aquifer can cause a high hydraulic gradient within it. This can lead to unbalanced shear stress in the aquifer. The shear stress can induce volumetric contraction under a small strain condition, e.g., less than 0.5%. Therefore, the deformation of AqIV in Shanghai can be considered from the shear stress-induced volumetric contraction as predicted by the proposed method rather than due to the creep effect and (or) yielding of the sand in aquifer IV. 5. The deformation delay with drawdown of the groundwater level of the aquifer can also be interpreted by considering the time factor, which was incorporated into Budhu's method. That is, the progressive increase in shear stress is one of the primary reasons for the deformation delay with decline of groundwater levels in aquifers, which causes a gradual change in the hydraulic gradient. Another significant reason is the consolidation of overburdened aquitards. Overall, both reasons are related to the increasing rate of effective stress in aquifers.

List of symbols

Fig. 1 .

 1 Fig. 1. Contour map of groundwater level in AqIV, Shanghai.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Sectional view of geological and hydrogeological profile to 400 m deep at the investigated site, Shanghai.

Fig. 4 .

 4 Fig. 4. History of groundwater pumping from aquifers in Shanghai from 1965 to 2001 (modified from Shen and Xu 2011).

Fig. 5 .

 5 Fig. 5. Distribution of total deformation in past two decades in Shanghai (data from SGEAEB 2002).

Fig. 6 .

 6 Fig. 6. Variation of groundwater level and deformation of AqIV at Laodong Park, Shanghai.

Fig. 7 .

 7 Fig. 7. Variation of groundwater level and deformation of AqIV at Shuangyang School, Shanghai.

Fig. 8 .

 8 Fig. 8. Variation of groundwater level and deformation in (a) AqII and (b) AqIII at Laodong Park, Shanghai (data from SGEAEB 2002).

Fig. 9 .

 9 Fig. 9. Finite element analysis mesh in groundwater flow model (y-axis is enlarged 3 times).

Fig. 10 .

 10 Fig. 10. Variation of effective stress and resulting deformation (created from the concept in Budhu and Adiyaman 2010).

Fig. 11 .

 11 Fig. 11. Cross-sectional view of groundwater pumping in aquifer: (a) aquitard over confined aquifer with potential; (b) stress distribution at the interface between aquifer and aquitard.

Fig. 12 .

 12 Fig. 12. Determination of time parameter = in eq. (8) from stable time period of groundwater level.

Fig. 13 .

 13 Fig. 13. Relation between time factor, a, and t.

Fig. 14 .

 14 Fig. 14. Comparison of the deformation of AqIV at Laodong Park from 1980 to 2000.

Fig. 15 .

 15 Fig. 15. Cross-sectional view of accumulated deformation at Laodong Park calculated based on the extended Budhu model.

Fig. 17 .

 17 Fig. 17. Cross-sectional view of accumulated deformation of AqIV at Shuangyang School calculated based on the extended Budhu model.

Fig. 18 .

 18 Fig. 18. Comparison of shear stress of Laodong Park case.

Fig. 19 .

 19 Fig. 19. Relationship between groundwater level and deformation of AqIV.

Fig. 20 .

 20 Fig. 20. Schematic diagram of variation of the pumping process.

a

  time factor C c compression index C ␣e secondary compression coefficient c cohesion e void ratio e 0 initial void ratio E Young's modulus [M/L/T 2 ] E s compression modulus [M/L/T 2 ] f d viscoplastic potential function G classical shear modulus [M/L/T 2 ] G c Cosserat continua shear modulus [M/L/T 2 ] h hydraulic head [L] h 0 (t) drawdown of groundwater head with time from the initial head surface at an arbitrary selected origin [L] h w (t) drawdown of groundwater head with time from the initial head surface at pumping well [L]k hydraulic conductivity [L/T] k h horizontal hydraulic conductivity [L/T] k v vertical hydraulic conductivity l z thickness of aquifer IV [L] M slope of critical-state line M c (M e )slope of the critical-state line m v soil coefficient of volume compressibility [L•T 2 /M] N strain-rate coefficient n porosity p m d size of the dynamic loading surface p m f size of the static yield surface q external source-sink flux [T -1 ] S s coefficient of specific storage [L -1 ] t time [T] u pore pressure Poisson's ratio V characteristic volume of element [L 3 ] w l liquid limit w n water content ratio w p plastic limit x distance from the origin [L] ␣ a soil layer constant ␤ factor related to the degree of saturation and the porosity ␤ w coefficient of volume compressibility of water [L•T 2 /M] ␥ t unit weight ␥ w unit weight of water [M/L 2 /T 2 ] ˙ij vp viscoplastic strain ˙v vp viscoplastic volumetric strain vl vp limit of viscoplastic volumetric strain slope of the swelling line slope of the compression line = constant fluidity [T -1 ] b water replenishment factor total stress p0 preconsolidation pressure = effective stress [M/L/T 2 ] ij effective stress tensor [M/L/T 2 ] ¯ij average effective stress tensor [M/L/T 2 ] shear stress [M/L/T 2 ] c reference time angle of internal friction

Table 1 .

 1 Parameter used in groundwater flow model.

	Layer m Aquitard I and phreatic water Depth (m) e 0ϳ-35 1.05 3.52E-04	1.60E-03	4.01E-03	3.45E-03
	Aquifer I	-35ϳ-50	0.78	1.95E-05	1.04E+00	2.59E+00	1.91E-04
	Aquitard II	-50ϳ-65	0.85	3.28E-04	2.34E-03	5.85E-03	3.21E-03
	Aquifer II	-65ϳ-100	0.82	9.87E-06	8.64E+00	2.16E+01	9.67E-05
	Aquitard III	-100ϳ-110	0.78	1.95E-04	1.65E-03	4.13E-03	1.91E-03
	Aquifer III	-110ϳ-150	0.71	8.58E-06	8.64E+00	2.16E+01	8.41E-05
	Aquitard IV	-150ϳ-165	0.68	2.41E-04	9.85E-04	2.46E-03	2.36E-03
	Aquifer IV	-165ϳ-225	0.66	6.27E-06	8.64E+00	2.16E+01	6.14E-05
	Aquitard V	-225ϳ-270	0.62	1.61E-04	7.26E-04	1.81E-03	1.58E-03
	Aquifer V	-270ϳ-300	0.64	4.97E-06	8.64E+00	2.16E+01	4.87E-05
	Aquitard VI	-300ϳ-350	0.62	1.61E-04	7.26E-04	1.81E-03	6.90E-05
	Rock	-350ϳ-400	0.10	-	1.00E-07	1.00E-07	1.50E-09

v (kPa -1 ) k v (m/day) k h (m/day) S s (m -1 )

Note: e, void ratio; m v , volume compression coefficient; k h and k v , horizontal and vertical hydraulic conductivities, respectively; S s , coefficient of specific storage.

Table 2 .

 2 Parameter used in classical Cauchy continuum model.

	Layer	Depth (m)		E (MPa) ␥ t (kN/m 3 ) *	*	c (kPa) (°) M
	Aquitard I and phreatic water 0ϳ-35	0.3 -	18.1	0.098	0.0098	17	29	1.2
	Aquifer I	-35ϳ-50	0.2	46	19.2	-	-	-	-	-
	Aquitard II	-50ϳ-65	0.3 -	19.1	0.049 0.0049	4	31	1.2
	Aquifer II	-65ϳ-100	0.2	91	19.1	-	-	-	-	-
	Aquitard III	-100ϳ-110	0.3 -	19.7	0.034	0.0034	5	34	1.2
	Aquifer III	-110ϳ-150	0.2	105	20.0	-	-	-	-	-
	Aquitard IV	-150ϳ-165	0.3 -	20.2	0.036	0.0036	13	27	1.2
	Aquifer IV	-165ϳ-225	0.2	144	20.2	-	-	-	-	-
	Aquitard V	-225ϳ-270	0.3 -	20.3	0.025	0.0025	30	20	1.2
	Aquifer V	-270ϳ-300	0.2	181	20.3	-	-	-	-	-
	Aquitard VI	-300ϳ-350 0.3 -	20.2	0.025	0.0025	30	20	1.2
	Rock	-350ϳ-400 0.2	2000	22.4	-	-	-	-	-

Note: , Poisson's ratio; E, Young's modulus; ␥ t , unit weight; *, modified compression index in Plaxis (*, /(1 + e), where is compression index); *, modified swelling index in Plaxis (k*, /(1 + e), where is swelling index); c, cohesion; , angle of internal friction; M, slope of critical-state line.

Acknowledgements

The research work described herein was funded by the National Nature Science Foundation of China (NSFC) (Grant No. 41072209 & 41102175). It is also supported by the starting funding of young faculty members and Chenxing award for young scholar of Shanghai Jiao Tong University. These financial supports are gratefully acknowledged. The writers are also grateful to anonymous reviewers, whose constructive comments have improved the quality of this paper.


Lirmm M Iago Bonnici

Dr Abdelkader Gouaich

Dr Fabien Michel

Formalisation of metamorph Reinforcement Learning

This technical report describes the formalisation of a particular Reinforcement Learning (RL) situation that we call "metamorph" (mRL). In this situation, the signature of the learner agent, i.e. its set of inputs, outputs and feedback slots, can change over the course of learning.

RL can be viewed as signal processing, because the learner agent transforms the inputs/feedbacks signals it is continuously fed with into output signals. The following formalisation is therefore concerned with signals description and the transformation from one signal to another. Also, since the signature of the agent is expected to change, we get concerned in the definition of what is a "signature" and a "signature change".

In the first part, we describe mRL learning context, or how the metamorph agent is embedded into its environment and interacts with it. In the second part, we describe one generic example of a metamorph learner agent: a dynamical computational graph that could theoretically be used in controlling the agent. In the last part, we reformulate the classical problem of RL, a.k.a. "maximizing feedback" in terms of this formalised mRL.

Learning Context 1 Values and their domains

We call domain a tuple d = (d.i, d.σ ), where:

• d.σ is a set • i.σ is an index, so that d.σ = d .σ =⇒ d = d
Note: within this whole document, the dot notation a.b does not represent a product between a and b, but the b attribute of object a.

Let D be the set of all domains. Let D be its flattened version containing all values from all domains:

D = d∈D d.σ (1)
Let ⊥ be a joker, null value belonging to none of the domains:

⊥ / ∈ D (2)
2 Flows

A flow is a function associating a value to each real date t:

h : R + → D t → h(t) (3)
3 Sequence approximation

For any precision ε ∈ R + * , a flow h may be approximated by a sequence ( ε h n ) n∈N such that:

∀n ∈ N, ε h n = h(nε) (4)
Note that the initial value is the same for all precisions:

∀ε ∈ R + * , ε h 0 = h(0) = h 0 (5)

In the next, we will often work with sequences and discrete time n ∈ N, with the idea that they approximate continuous flows as ε → 0.

Flow determination

Let (u n ) n∈N be a sequence in D. For any n, its "past" values live in D [N] with:

D [N] = i∈N D i (6)
so that ∀n ∈ N, (u 0 , . . . , u n ) ∈ D [N] (7)

Let f be a function associating a value in D to each precision and each past.

f : R + * × D [N] → D (ε, (u 0 , . . . , u n )) → f ε (u 0 , . . . , u n ) (8)
Let h, k be two flows in D, D : we say that f determines k from h when all sequences that f generates from approximations of h are approximations of k:

∀ε ∈ R + * , ∀n ∈ N * , ε k n = f ε ( ε h 0 , . . . , ε h n ) (9)
This may be written another, graphical way:

h ( f ) k (10)
A slightly different case is needed to bootstrap circular situations: we say that f initially determines k from h when the sequence determination is shifted towards zero:

∀ε ∈ R + , ε k 0 = f ε (∅) ∀n ∈ N * , ε k n = f ε ε h 0 , . . . , ε h n-1 (11)
Which is written:

h ( f * ) k (12)
5 Graphical aliases

Fleshing this graphical notation, we will use the following alias when k is determined by the past of h and h combined:

h h ( f ) k (13) ∀ε ∈ R + * , ∀n ∈ N * , ε k n = f ε ε h 0 , ε h 0 , . . . , ε h n , ε h n (14)
And the following alias when f determines both k and k from h:

k k ( f ) h (15) ∀ε ∈ R + * , ∀n ∈ N * , ε k n , ε k n = f ε ( ε h 0 , . . . , ε h n ) (16)

Flow determination determination

When the function f determining one flow from another is itself determined by another flow, we will use this construct:

h h (F) ( f ) k (17)
Meaning that there exists a flow of determining functions f : R + → F D [N] → D , approximated for each precision ε by a sequence ( ε f n ) n∈N so that:

∀ε ∈ R + * , ∀n ∈ N, ε f n = F ε h 0 , . . . , ε h n ε k n = ε f n (h 0 , . . . , h n ) (18)

States and Multiple flows

We call state a tuple s = (s.∆, s.ν), where:

• s.∆ is a subset of D that we call the state signature

• (s.ν d ) d∈s.
∆ is a family of values indexed by s.∆, such that each value belongs to the set of the corresponding domain:

∀d ∈ s.∆, s.ν d ∈ d.σ (19)
Let S be the set of all states.

A flow of states h : R + → S is called a multiple flow. It carries two interesting informations:

• a flow of domains that we call the signature of h. We write it h.∆:

h.∆ : R + → 2 D t → h.∆(t) = h(t).∆ (20)
• a flow of values h.ν, at any time consistent with the signature:

h.ν : R + → D [N] t → h.ν(t) = h(t).ν (
∀d ∈ D, ∀t ∈ R + , h d (t) = h.ν(t) d if d ∈ h.∆(t) ⊥ if not (22)
In other words, h can represents many flows and h d "selects" the subflow corresponding to one domain d.

Partial flow determination

Let k be a multiple flow, one such construct will be used:

h ( f ) k.∆ k.ν ( f ) h k (23)
It states that the signature of k and its values are determined by two different sources (h, f ) and (h , f ):

∀ε ∈ R + * , ∀n ∈ N, ε k.∆ n = f ε ( ε h 0 , . . . , ε h n ) ε k.ν n = f ε ε h 0 , . . . , ε h n (24)
10 Bit flows and Contiguous multiple flows

Let h d : R + → d.σ ∪ {⊥} be a subflow in domain d.
h d is said to be a bit flow if and only if its value is always ⊥ except on one contiguous window in R + . I.e., there exists t 0 , t 1 ∈ R + (t 1 may be infinite) such that:

∀t ∈ R + , if t ∈ ]t 0 , t 1 [ , h d (t) ∈ d.σ else h d (t) = ⊥ (25)
The domain d is said to be active at all dates t ∈ ]t 0 , t 1 [, regarding to the subflow h d , or else it is inactive.

When approximated with precision ε, the activation window translates to the corresponding discrete window ε n 0 , ε n 1 ⊂ N, ε n 0 being the date of first non-⊥ value in ε h, and ε n 1 the date after the last non-⊥ value or ∞. Now let h : R + → S be a multiple flow. h is said to be a contiguous multiple flow if and only if all its subflows h d are bit flows.

A contigous multiple flow can be represented as figure 1. Notice how each bit subflow is only active on a contiguous window of R + . In this example, the first four domains describe the following set of values for subflows:

[0, 1] Z r g b R t a h.ν(t a )=(9.75, 14) h.∆(t a )={d 1 , d 3 } t b h.ν(t b )=(-4.87, b, 12, 0.61) h.∆(t b )={d 1 , d 2 , d 3 , d 4 } t c h.ν(t c )=(g) h.∆(t c )={d 2 } : h d 4 : h d 3 : h d 2 : h d 1 ⊥ ⊥ ⊥ ⊥ ⊥ t∈R + h : s∈S subflows \ evaluated in t a t b t c h d 4 : R + → [0, 1] ∪ {⊥} h d 4 (t a ) = ⊥ h d 4 (t b ) = 0.61 h d 4 (t c ) = ⊥ h d 3 : R + → Z ∪ {⊥} h d 3 (t a ) = 14 h d 3 (t b ) = 12 h d 3 (t c ) = ⊥ h d 2 : R + → {r, g, b} ∪ {⊥} h d 2 (t a ) = ⊥ h d 2 (t b ) = r h d 2 (t c ) = g h d 1 : R + → R ∪ {⊥} h d 1 (t a ) = 9.75 h d 1 (t b ) = -4.87 h d 1 (t c ) = ⊥
d 1 .σ = R, d 2 .σ = {r, g, b} , d 3 .σ = Z, d 4 .σ = [0, 1].
The gray dotted lines correspond to three samples t a , t b , t c where the multiple flow h and its subflows are evaluated. Dashed areas mark the end of each bit subflow active window.

High-level agent representation

At the highest level, our learner agent can be represented by 4 contiguous multiple flows (i, o, ϕ, P), with the following determination scheme:

(E * ) o.∆ ϕ (A) o.ν (P) i o (26)
Here:

• E represents the environment in which the agent is immersed

• i represents the agent's inputs or sensors: their nature may change in time as the signature i.∆ evolves.

• o represents the agent's outputs or actuators: their nature may also change in time as o.∆ evolves.

• ϕ represents the agent's objectives, a continuously fed evaluation of the actions it undertakes. They also may change in nature as ϕ.∆ evolves.

Note that the environment determines i, ϕ and o.∆, so the agent cannot decide the data it is fed with or its output signature.

• P represents the agent instant behaviour. It is an inner computational procedure determining the output values based on all input history.

• A is the abstract agent strategy, constantly adapting its behaviour based on environmental information.

Only two objects are not depending on time here: the environment E and the inner agent strategy A. E will be defined by user, and it is our job to design A (see objectives later).

Computational graph: an example of metamorph agent

In this chapter, we describe a graph process that can be translated into the inner reaction P of our learner agent.

Computational graphs are elements of Γ, that we will progressively describe hereafter.

12 Abuse of notation

In the next, three objects will essentially represent the same thing: the inner behaviour of the agent:

• One flow P : R + → Γ, assigning one computational graph to each date.

• One flow P : R + → F S [N] → D [N] of determining functions, determining output values o.ν from input pasts.

• One contiguous multiple flow P : R + → S describing the evolution of every intermediate value in the flowing computational graph.

Since they are conceptually related but still different in nature, we think that using the same symbol for all of them will cause less confusion than using one different symbol for each. As such, we will refer to each of them with the same letter P, and let the reader rely on the context to figure them out. Also, considering P as a multiple flow and P d its subflow in domain d, we will refer to the sequence approximation ( ε P d n ) n∈N with the simpler expression (d n ) n∈N so that double subscript is avoided.

As a consequence, in such a context, the symbol d refers both to the domain d ∈ D and to the corresponding subflow sequence approximation d = ε P d .

Operations

A computational graph in Γ is a structure that typically contains "operations".

We call operation an object of the form:

a, (a 1 , ρ 1 ), (a 2 , ρ 2 ), . . . , (a p , ρ p ) , (a I , a II , . . . , a q ) (27)

Where:

• All a i ∈ D are domains, and a P a i subflow corresponds to each of them.

• (a 1 , ρ 1 ), . . . , (a p , ρ p ) are called the operation inputs

• a I , . . . , a q are called the operation outputs.

• All ρ i ∈ a i .σ are called rest values, and belong to their corresponding input domain.

• a : D p [N] → D q is a determining function called the operation formula. Its exact signature is a :

(a 1 .σ × • • • × a p .σ ) [N] → a I .σ × • • • × a q .σ .
The computational graph process P, like any flow, can be approximated by a sequence ε P. At any date n ∈ N, the instant computational graph ε P n ∈ Γ may contain such an operation, so it may read somewhere:

a 2 :ρ 2 (a) a 1 :ρ 1 a 3 :ρ 3 a I a II ∈ ε P n = P(nε) (28)
This means that all domains a i are being active within the multiple flow P at this date. In other words: the bit subflows P a i are in their non-⊥ window:

∀n ∈ N, operation "a" present in ε P n ⇐⇒ ∀i ∈ {1, . . . , p, I, . . . , q} , a i n = ⊥ =⇒ a i n ∈ a i .σ (29)
Also, the outputs of the operation are determined by its inputs histories: a I n , . . . , a q n = a (a 1 0 , . . . , a p 0 ), . . . , (a 1 n , . . . , a p n )

Operations are solid in that no input (a 1 , ρ 1 ), . . . , (a p , ρ p ), no output a I , . . . , a q and no formula a can be present in the graph if not all other elements of the operation are present.

Note: since the multiple flow P is contiguous, P a i are bit subflows. Thus, as a consequence of operation solidity, every operation can only appear once in the graph process P during a ]t 0 , t 1 [ activity window. And this activity window is the same for all its related P a i bit subflows.

As a consequence, for each operation, there exists one maximal "arrival date" in the graph process, t 0 ∈ R + , and a corresponding maximal "approximated arrival date" ε n 0 ∈ N, such that:

∀t < t 0 , P a i (t) = ⊥ (31) ∀n < ε n 0 , a i n = P a i (nε) = ⊥ (32)
Operations are autonomous (in the sense of dynamical systems) in that their determining function a ignores the arrival date of the operation. In other words, their result is independent on the number of initial ⊥ values, so ∀n ε n 0 : a (⊥, . . . , ⊥), . . . , (⊥, . . . , ⊥), (a 1 n 0 , . . . , a p n 0 ), . . . , (a 1 n , . . . , a p n ) = a (a 1 n 0 , . . . , a p n 0 ), . . . , (a 1 n , . . . , a p n )

(33)

This makes the operation capable of processing flowing information, gifted with a memory, but unaware of the absolute date.

Fluxes and delays

Operations in a computational graph in Γ are typically connected together via "fluxes". Operation formulae determine outputs subflows, and fluxes determine input subflows.

We call flux an object of the form:

ω = (δ , o, (ι, ρ)) (34)
Where

• δ ∈ R + * is a delay value • o is an output • (ι, ρ) is an input
At any date, the computational graph may contain such a flux, so it reads:

o ι:ρ δ ∈ ε P n or ω ∈ ε P n (35)
This means that the subflow P ι copies the subflow P o with a delay δ as long as ω is part of the graph. However, it cannot copy values of P o anterior to the apparition of ω in the graph, and will use its rest value ρ instead.

More formally, let t 0 ∈ R + be the arrival date of ω into the graph:

∀t > t 0 , P ι (t) = if t -δ t 0 : ρ else : P o (t -δ ) (36)
Inputs cannot have more that one flux attached at any given date.

For any other input (ι, ρ) present in the graph, but with no flux attached, their corresponding subflow has the rest value ρ:

∀t ∈ R + , ∅ ι:ρ ∈ P(t) =⇒ P ι (t) = ρ (37)
With fluxes, operations can be connected to each other so that information gets processed by their formulae from one subflow to another.

Since each flux transfers information with a positive delay δ , no circular configuration of operations yields undefined subflows, and we can always find a precision ε < δ yielding approximated subflows ε P d that are defined for all n ∈ N (see example section 15).

15 Embedding into environment: one example At any time, the computational graph is connected to the agent input and output flows.

Each active subflow in the agent inputs corresponds to one free output in the graph, so that operations inputs can feed from them.

Each active subflow in the agent outputs corresponds to one free input in the graph, so that operations outputs can feed them.

For instance consider the following initial situation:

         (i 0 , ϕ 0 , o.∆ 0 ) = E(∅) i.∆ 0 = {A, B, C} ⊂ D i.ν 0 = (A 0 , B 0 , C 0 ) ∈ A.σ × B.σ ×C.σ o.∆ 0 = {D, E, F} ⊂ D (38)
Then, the minimal computational graph P 0 ∈ Γ would be:

P 0 = A B C D:ρ D E:ρ E F:ρ F ∅ ( 39
)
With no operations, only one free output per agent input, one free input per agent output and one rest value for each free input.

In this situation, only the rest values would be used, so there would be, D 0 = ρ D , E 0 = ρ E , F 0 = ρ F . So, with any precision ε ∈ R + * , the first outputs values would be determined as:

ε o.ν 0 = (ρ D , ρ E , ρ F ) ∈ D.σ × E.σ × F.σ (40)
However, the initial computational graph could also be more fleshed. In the following example, it contains three operations with their inner formulae a, b, c, connected to each other and to the agent environment via 8 fluxes:

P 0 = A B C D:ρ D E:ρ E F:ρ F (a) a 1 :ρ a 1 a 2 :ρ a 2 a 3 :ρ a 3 a I a II b 1 :ρ b 1 (b) b I c 1 :ρ c 1 c I (c) c II c 2 :ρ c 2 δ 1 δ 2 δ 3 δ 4 δ 5 δ 6 δ 7 δ 8 (41)
Then, assuming that the input signature stays constant (i.∆ 0 = ε i.∆ 1 = ε i.∆ 2 ) then evolves into ε i.∆ 3 = {B, H}, and the output signature stays constant (o.∆ 0 = ε o.∆ 1 = ε o.∆ 2 ) then evolves into ε o.∆ 3 = {F, K}, the graph could also stay con-stant (P 0 = ε P 1 = ε P 2 ) then become, for example: Let's assume that the ε-approximation of continuous time yields that the delays above correspond to these discrete delays: delay (R + * ) δ 1 δ 2 δ 3 δ 4 δ 5 δ 6 δ 7 δ 8 δ 9 δ 10 δ 11 δ 12 ε-delay(N * ) 2 1 1 2 2 3 1 1 1 2 2 3 (43)

ε P 3 = B H F:ρ F K:ρ K (a) a 1 :ρ a 1 a 2 :ρ a 2
Then, the approximated multiple flows ε i, ε P, ε o can be computed like in the following tables:

n A B C H 0 A 0 B 0 C 0 ⊥ 1 A 1 B 1 C 1 ⊥ 2 A 2 B 2 C 2 ⊥ 3 ⊥ B 3 ⊥ H 3 4 ⊥ B 4 ⊥ H 4 5 ⊥ B 5 ⊥ H 5 n D E F K 0 ρ D ρ E ρ F ⊥ 1 ρ D ρ E c I 0 ⊥ 2 b I 0 ρ E c I 1 ⊥ 3 ⊥ ⊥ ρ F ρ K 4 ⊥ ⊥ ρ F ρ K 5 ⊥ ⊥ ρ F e I

inputs outputs

Table 1: First few values of the example run: each column represents the εapproximation of one input or output domain. The global agent signature has changed between approximation steps 2 and 3. Notice how input subflows are either constant rest flows or delayed copies of other subflows. For instance, a 1 contains a copy of B, while c 1 , whose delays δ 6 then δ 11 are long, is mostly resting to its rest value ρ c 1 .

In other words, no matter the experienced environment or the current signature, we want the agent A to be able to find which output signals will make all feedback maximal. When this is achieved, we say that our agent keeps "learning", even though its objectives, inputs and outputs change over time.

Figure 1 :

 1 Figure 1: An example contiguous multiple flow h and its bit subflows h d 1 , . . . , h d 4 .Notice how each bit subflow is only active on a contiguous window of R + . In this example, the first four domains describe the following set of values for subflows:d 1 .σ = R, d 2 .σ = {r, g, b} , d 3 .σ = Z, d 4 .σ = [0, 1].The gray dotted lines correspond to three samples t a , t b , t c where the multiple flow h and its subflows are evaluated. Dashed areas mark the end of each bit subflow active window.

n a 1

 1 a 2 a 3 a I a II b 1 b I c 1 c 2 c I c II e 1 e 2 e I 0 ρ a 1 ρ a 2 ρ a 3 a I 0 a II 0 ρ b 1 b I 0 ρ c 1 ρ c 2 c IFirst few values of the example inner multiple flow P. Each column represents the evolution of an operation input or output. The structure of the computational graph has changed between approximation steps 2 and 3.

  ) d∈D , where each flow h d : R + → d.σ ∪ {⊥} associates the value taken by h in the domain d on each date, and ⊥ if d is not part of the signature of h on this date:

	21)
	8 Subflows
	A multiple flow h : R + → S can also be seen as / is isomorphic to a family of
	subflows (h d

Objective

In this chapter, we state the last constraints on the system so it can be simulated with ε-approximations and our objective can be expressed.

Finite multiple flows

A multiple flow h : R + → S is finite if and only if its signature is finite at any date:

All constitutive multiple flows (i, o, ϕ, P) of the agents we are interested in are finite.

Orderable objectives

In the agents we are interested in, all sets visited by the objective multiple flow ϕ are numerical, i.e. subsets of R.

18 Increasing objectives A collection of E i is given as problem data (it may be reduced to only one E).

Our goal is to design the inner agent strategy A such that all objective subflows ϕ d are maximal.

More formally, for any bit subflow ϕ d , in any environment E i , let ]t 0 , t 1 [ be its activation window. Then, we want: (46)