Numerical simulation and analysis on the mechanical responses of the urban existing subway tunnel during the rising groundwater

Tunnelling and Underground Space Technology 98 (2020) 103297

Contents lists available at ScienceDirect

Tunnelling and Underground Space Technology journal homepage: www.elsevier.com/locate/tust

Numerical simulation and analysis on the mechanical responses of the urban existing subway tunnel during the rising groundwater

T

Dechun Lu , Xiaoqiang Li, Xiuli Du, Qingtao Lin, Qiuming Gong ⁎

Key Laboratory of Urban Security and Disaster Engineering of Ministry of Education, Beijing University of Technology, Beijing 100124, China

ARTICLE INFO

ABSTRACT

Keywords: Groundwater level rise Existing tunnel Numerical simulation Wetting-induced deformation Stress and pore pressure coupled analysis Mechanical behaviours

In this paper, a numerical method of the rising groundwater level is presented through an elastoplastic constitutive model for the unsaturated soil and a numerical realization method for the rising groundwater level. A model test of the rising groundwater level is conducted to validate the rationality of the presented numerical method. Taking an interval subway tunnel in Beijing as the background, a numerical simulation with a homogeneous soil layer, which contians the existing tunnel, is conducted according to the presented numerical method of the rising groundwater level. Based on the numerical simulation results, the variations of the ground saturation, the ground stress state, and the ground vertical deformation under the condition of the rising groundwater level are studied through a homogeneous soil layer, and the normal contact pressure on the outer surface of the lining, the radial deformation, the bending moment, the hoop thrust, and the shear force of the lining under the condition of the rising groundwater level are analyzed through a homogeneous soil layer with the existing tunnel.

1. Introduction Over recent years, an increasing trend in the groundwater level (GL) is observed all around the world, such as Beijing in China, London in the United Kingdom (Wilkinson, 1985), New York in the United States (Wilkinson, 1985), Tokyo in Japan(Kusaka et al., 2016), Aswan in Egypt (Selim et al., 2014), Lagos in Nigeria (Oyedele et al., 2009), Madinah in Saudi Arabia (Bob et al., 2016), Milan in Italy (Beretta et al., 2004; Gattinoni and Scesi, 2017) and other cities. Groundwater restorations can be triggered by one or more factors. These factors include the deliberate restoration of the past GL (Al-Rashed and Sherif, 2001), the reduced abstraction of the groundwater (Kreibich and Thieken, 2008; Shanahan, 2009), the seawater infiltration (Al-Sefry and Sen, 2006; Yasuhara et al., 2007), et al. The rising GL results in a host of engineering problems to the ground and the existing underground structures (Al-Sefry and Sen, 2006; Brassington, 1990; Chaudhary, 2012). For example, more than 10 cm rebound deformation is observed in Tokyo due to the rising GL (Kusaka et al., 2016). In addition, in Madinah of Saudi Arabia, a short-term remedial solution, which is used to pump groundwater through a network of nineteen wells, is implement to prevent the rise of the GL and protect the lager underground infrastructure (Bob et al., 2016). The Sabah Al-Salem area and Shamiyah-Kaifan area of Kuwait City install

the drainage system to control the rise of the GL and the rate of the water pumped out is beyond 7200 m3/day (Al-Rashed and Sherif, 2001). In London of the United Kingdom, since the year 2000, the rate of the water pumped out has been increased in order to control the rise of the GL and prevent the groundwater intrusion into subsurface structures (Shanahan, 2009). In Jeddah of Saudi Arabia, the underground drain system is constructed and more than 60000 m3 of groundwater is discharged into the Red Sea every day (Abu-Rizaiza et al., 1989). Groundwater is one of the main challenges associated with stability and safety issues in the operation of the existing tunnels (Fang et al., 2016; Shin et al., 2012, 2005). Problems with the rising GL are very serious. But up to now, the research on the effect of the rising GL on the existing underground structures is rare and the interaction mechanism between the rising GL and the response of the existing underground structures is still unclear (Kreibich and Thieken, 2008). This paper presents a numerical method of the rising GL, which adopts the constitutive model of the unsaturated soil to describe the mechanical behaviours of the soil and proposes a numerical realization method for the rise of the GL. Taking an interval tunnel of Beijing Metro Line 14 as the research background, a numerical simulation of the rising GL is conducted. Based on the numerical simulation results, the mechanical behaviors of the ground, the mechanical behaviours of the existing tunnel are studied in detail.

Corresponding author. E-mail addresses: [email protected] (D. Lu), [email protected] (X. Li), [email protected] (X. Du), [email protected] (Q. Lin), [email protected] (Q. Gong). ⁎

https://doi.org/10.1016/j.tust.2020.103297 Received 31 December 2018; Received in revised form 22 November 2019; Accepted 13 January 2020 0886-7798/ © 2020 Elsevier Ltd. All rights reserved.

Tunnelling and Underground Space Technology 98 (2020) 103297

D. Lu, et al.

2. Numerical method for the rising GL In the process of the rising GL, the void of the soil is gradually occupied by the water and the saturation of the soil gradually increases to 1. The rising GL in the ground is the wetting process of the unsaturated soils. The mechanical behaviours of the unsaturated soil are not only affected by the effective stress, but also related to the suction in soil. The numerical method of the rising GL is an effective way to research the problems related to the rising GL. In the numerical method of the rising GL presented in this paper, the changes of the pore water pressure for all ground elements are determined in advance and are taken as the load condition. The constitutive model of the unsaturated soil is implemented in the finite element code ABAQUS by combining the known changes of the pore water pressure. 2.1. Numerical model for the rising GL In the rapid development of Beijing, the long-term overexploitation of the groundwater causes the GL to drop with an average annual decrease of nearly 1 m from 1998 to 2014. Subsequently, with the opening of the South-to-North Water Transfer Project in 2014, the GL begins to rise year by year. The average rate of the rising GL is 0.64 m/ year. This paper takes the interval shield tunnel between the Dongfengbeiqiao Station and the Jiangtai Station of Beijing Metro Line 14 as the research background. A homogeneous soil layer is selected to study the action mechanism of the rising GL on the existing structures and the ground, and the interaction between the ground and the existing structures. The interval tunnel is a circular shield tunnel. The external and internal diameters of the lining are 10 m and 9 m, respectively. The distance from the ground surface to the crown is 15 m. In the numerical simulation of this paper, the numerical model of the ground is simplified as the homogeneous soil. Along with the direction of the tunnel excavation, the shape and the size of the tunnel are the same. According to the principle of symmetry, when the tunnel excavation is modelled through full section excavation with one step, the ground displacement along with the direction of the tunnel excavation is zero. The numerical simulation can adopt the plane strain model. Thus, the length of tunnel excavation direction (Y-axis) is taken as 4 m to saving computing time as shown in Fig. 1. In selecting the geometric dimensions of the other two directions, a sensitivity study is conducted. Based on the results of the sensitivity study, the size of the model ground is determined as 100 m (X-axis) × 4 m (Y-axis) × 60 m (Z-axis) as shown in Fig. 1. The distance between the tunnel center and the left and right lateral boundary is the same and is 5 times of the tunnel diameter. The numerical model of the ground is discretized using eight-node trilinear displacement and trilinear pore pressure elements with reduced integration and with the hourglass control. In order to improve the accuracy of numerical calculation and the convergence of the numerical simulation, the elements of the area near the existing tunnel are refined. For the interval shield tunnel, the width of each ring of the lining is 1.8 m. Every ring of the lining is assembled with nine precast reinforced concrete segments. The thickness of the concrete segment is 0.5 m. The waterproof materials with good performance are used in the outer surface of the lining. The lining can be considered impervious. For the mechanical behviours of the lining, the tunnel joints reduce the overall stiffness of the lining significantly (Hefny and Chua, 2006; Teachavorasinskun and Chub-uppakarn, 2010). This paper adopts the method of continuous solid elements with reduced stiffness (Liao et al., 2008; Zhang et al., 2015) to consider the effect of the tunnel joints on the mechanical behaviours of the lining. The lining is simplified as a homogeneous ring with an inner diameter of 9 m, an outer diameter of 10 m and a axial length of 4 m as shown in Fig. 2. The numerical model of the lining is discretized using the eight-node trilinear displacement elements with reduced integration. In the numerical model, the displacement boundary conditions are

Fig. 1. The element mesh of the numerical model of the ground.

Fig. 2. The element mesh of the numerical model of the lining.

applied to the lateral boundaries and the bottom boundary. At the left and right lateral boundaries, the displacements along with the X-axis are restrained. At the front and back lateral boundaries, the displacements along with the Y-axis are fixed. At the bottom boundary, the pin supports are applied to the bottom boundary and the displacements along with the X-axis, the Y-axis, and the Z-axis are fixed. The top boundary adopts the force boundary and the force at the top boundary is set as 0 kN. 2.2. Modelling the mechanical behaviours 2.2.1. The mechanical behaviours of the soil The elastoplastic constitutive model for the unsaturated soil, which is proposed by Sun et al. (2008), is adopted to describe the wettinginduced deformation of the soil and is implemented in the finite element code ABAQUS via the user material (UMAT) subroutine. The Sun model (Sun et al., 2008) is simple in form and adopts the suction and the effective stress as the stress-state variables. The parameters in the Sun model have the definite meanings (Sun et al., 2007). 2

Tunnelling and Underground Space Technology 98 (2020) 103297

D. Lu, et al.

In the Sun model (Sun et al., 2008), the strain increment is divided into the elastic strain increment and the plastic strain increment. The elastic strain increment d ije (i, j = 1, 2, 3) is calculated through the Hooke’s law and can be express as:

1+ E

=

d

ij

E

d

80

Sr / %

e ij

d

100

(1)

kk ij

where δij is the Kronecker delta, E is the Young’s elastic modulus. ν is the Poisson's ratio, and σ′ij is the effective stress tensor. σ′ij is defined by ij

=

ua

ij

ij

+ Sr s

0 10 0

where σij is the total stress tensor, and Sr is the saturation. s is the suction and is defined by

s = ua

(3)

uw

=

Sr u w

ij

E is defined by

E=

p (1 + e )

(5)

where e is the void ratio, and κ is the slope of the swelling line in the e – ln p′ plane under the isotropic compression condition. p′ is the effective mean stress and can be express as:

p =

1 3

(6)

ij ij

The plastic strain increment d is calculated through the condition of consistency and the associated flow rule. The yield function f is the same as the plastic potential function g, and they can be expressed as:

q2 +p M 2p

(0) (s)

p 0y

pn

pn

=0

(7)

where q is the deviatoric stress and can be expressed as

q=

3 ( 2

ij

p

ij )( ij

p

ij)

(8)

M is the slope of the critical state line in the p′ – q plane, p′0y are the intersection of the yield surface and the p′-axis when s = 0 kPa. p′n is a material parameter and is the value of the p′0y without the wettinginduced deformation. λ(s) is the slope of the normal consolidated line in the e – ln p′ plane when the suction s is constant. λ(0) is the value of λ(s) when s = 0 kPa. λ(s) has the following form:

(s ) = (0) +

ss pa + s

Table 1 Material parameters of the soil. λs

κ

ν

M

pa (kPa)

p′n (kPa)

16.5

0.01

0.01

0.003

0.3

1.8

101

4000

10 5

In the numerical simulation, the initial state plays a significant role to study the nonlinear problems. The numerical model in this study contains the ground and the existing tunnel. The initial state of the ground includes the initial hydraulic state and the initial mechanical state. The initial state of the existing tunnel only includes the initial mechanical state. Due to that the rate of the rising GL is very slow, the initial hydraulic state of the ground is a steady state. Thus, the initial pore water pressure distribution can be determined according to the position of the GL and is directly applied to the ground. Then, the initial saturation distribution is calculated on the basis of the adopted soil-water characteristic curve. The initial mechanical state of the ground includes the initial stress state and the initial displacement field. In the process of determining the initial mechanical state, the initial hydraulic state remains unchanged. The initial mechanical state of the ground is determined through two steps. The first step is conducted to obtain the initial mechanical state of the ground without the existing tunnel. According to the mechanical boundary conditions, the weight of the homogeneous soil, and the initial hydraulic state, the initial mechanical state is determined through the Geostatic analysis step. The Geostatic

where λs is a material parameter. pa represents one atmospheric pressure. In this paper, the model parameters of the soil are listed in Table 1. The values of these parameters are determined by reference to the parameters of the Pearl clay (Sun et al., 2008). The initial void ratio is taken as 1.0. The soil-water characteristic curve of the unsaturated soil is used to describe the relationship between the saturation and the suction in the

λ

10 4

2.3. Determination of the initial state for the numerical model

(9)

γ (kN/m3)

10 3 s / kPa

2.2.2. The mechanical behaviours of the lining The linear elastic model is adopted to describe the mechanical behaviours of the lining. The Young’s elastic modulus of the lining is determined through the Young’s elastic modulus of the concrete segment and the reduction factor. The Young’s elastic modulus of the concrete segment is 35 GPa. The reduction factor is taken as 0.7 (Liao et al., 2008; Zhang et al., 2015). The Young’s elastic modulus of the lining is determined as 24.5 GPa. The Poisson's ratio ν is taken as 0.17. The unit dry weight γ is 25 kN/m3. On the interface between the ground and the lining, the normal mechanical behaviours are modelled by the hard contact model, and the tangential mechanical behaviours obey the Coulomb law. The friction coefficient of the Coulomb law is taken as 0.3.

p ij

f=g=

10 2

unsaturated soil. In addition, the initial state and the mechanical properties of the unsaturated soil also affect the soil-water characteristic curve. Determining the soil-water characteristic curve is not trivial. For simplify, in the numerical simulation of the rising GL, the soil-water characteristic curve is used to grasp the main rules between the saturation and the suction. For the model ground above the GL, with the increase of the distance from the GL, the saturation decreases and the suction increases. Therefore, as shown in Fig. 3, the adopted soil-water characteristic curve (Vanapalli et al., 1999) can represent the saturation variation rules of the model ground with the suction. In ABAQUS, the soil-water characteristic curve is applied to the model ground through inputting the negative pore water pressure and the corresponding saturation in form of a list.

(4)

ij

10 1

Fig. 3. The soil-water characteristic curves of the soil (Vanapalli et al., 1999).

where ua is the pore air pressure, and uw is the pore water pressure. ua adopts the atmospheric pressure as the zero pressure and remains constant throughout the domain being modeled. Thus, for the unsaturated soils, the suction is the absolute value of the pore water pressure. Eq. (2) can be simplified as: ij

40 20

(2)

ij

60

3

Tunnelling and Underground Space Technology 98 (2020) 103297

D. Lu, et al.

analysis step has the ability to determine the stress state which satisfies the initial boundary conditions and the maximum allowable displacement. After the Geostatic analysis step, the initial state of the ground without the existing tunnel has been determined. Based on the initial state obtained through the first step, the second step is conducted to obtain the initial mechanical state of the ground with the existing tunnel through simulating the tunnel excavation. The tunnel excavation is modelled through full section excavation with one step. The stress state and the deformation field, which are caused by the tunnel excavation, are calculated through the steady-state stress and pore pressure coupling analysis. The equilibrium state after the tunnel excavation is taken as the initial state of the ground with the existing tunnel. In the process of the second step to determine the initial state of the ground with the existing tunnel, the initial state of the existing tunnel is obtained at the same time.

Barrel 2

Barrel 1

Permeable plate

Hose

2.4. Implement method of the rising GL

Fig. 5. The test device of the rising GL.

According to the rate of the rising GL in different cities of the world (Bob et al., 2016; Brassington, 1990; Colombo et al., 2018; Oyedele et al., 2009), it can be found that the rate of the rising GL in most cities is very slow and is less than 1 m/year. The groundwater in the process of the rising GL is considered under the steady sate and the penetration of groundwater is ignored. For the steady state, the total hydraulic head is the same throughout the domain being modeled. The total hydraulic head H is the sum of the position hydraulic head h and the pressure head uw/γw. γw is the unit weight of the water, and is taken as 10 kN/m3. Fig. 4 shows a steady state for the saturated and unsaturated soils. For the convenience of analysis, the position, where h = 0 m, is defined on the GL. At the GL, uw is equal to 0 kPa. Thus, the total hydraulic head for the soil is equal to zero. uw is equal to -γwh. At the unsaturated sate as shown in point A, hA is positive and uwA is negative. At the saturated sate as shown in point B, hB is negative and uwB is positive. Whether in the saturated or unsaturated soil, the pore water pressure increases linearly with the increase of the depth, and the slope is γw. Because the lining is impervious, the effect of the existing shield tunnel on the pore water pressure distribution is ignored. The height of the numerical model of the ground is 60 m and the position of the initial GL is at a depth of 40 m. Combining with the features of the pore water pressure distribution at the steady state, the initial pore water pressure above the GL increases linearly from −400 kPa to 0 kPa and the initial pore water pressure below the GL increases linearly from 0 kPa to 200 kPa. In the process of the rising GL, the pore water pressure throughout the domain being modeled increases linearly along with the GL. When

h

2.5. Verification of the numerical model A test device of the rising GL is designed to study the influence of the rising GL as shown in Fig. 5. The test device is mainly composed of a permeable plate and two cylindrical transparent plexiglass barrels with the wall thickness of 10 mm. There are holes with a diameter of 3 mm and a space of 2 cm on the permeable plate. The Barrel 1 is used to place the soil sample. The Barrel 2 is used to place the water. The connection between the Barrel 1 and the Barrel 2 is achieved through a hose at the bottom of the barrels. Thus, at the steady state, the GL in the Barrel 1 can be expressed through the position of the water level in the Barrel 2. In the model test of the rising GL, the soil sample is prepared through sprinkling fine sand into the Barrel 1. The soil sample is a cylinder with a diameter of 0.35 m and a height of 0.5 m. The position, where the tick label is 0.0 m, is set at the bottom of the sample. Correspondingly, the tick label at the top of the sample is 0.5 m. A displacement sensor is placed on the top surface of the specimen to monitor the displacement of the top surface. The displacement sensor has a measuring range of 50 mm and an accuracy of 0.05 mm. The test procedure can be expressed through the following steps. (1) Before starting the test, the water level in the Barrel 2 is maintained at the tick label of 0.0 m by adding water. (2) After an hour, water in the soil sample is considered to be in a steady state and the data of the displacement sensor is recorded. (3) The water level in the Barrel 2 is raised by 10 cm and remains unchanged by adding water. (4) Repeat steps (2) and (3) until the water level reaches the top of the specimen. During the process of test, the error of the water level height in the Barrel 2 is less than 0.005 m. The test data are depicted by the blue hollow circle in Fig. 6. The positive value of the surface displacement represents the settlement of top surface. According to the presented numerical method of the rising GL, a corresponding numerical simulation is conducted. The model parameters of the fine sand are shown in Table 2. The soil–water characteristic curve, as shown in Fig. 3, is adopted. The simulation results are shown by the black solid line in Fig. 6. It can be seen that the numerical simulation results agree well with the experimental data and the presented numerical method of rising GL is reasonable to simulate the process of the rising GL.

uw=-γwh1

h1

hA

the GL is located at the depth of hGL, the pore water pressure at the depth of D is γw(D − hGL). In this paper, the variations of the pore water pressure along with the GL are applied to all elements of the ground as the load condition via the displacement (DISP) subroutine. The mechanical behaviours of the soil and the existing tunnel are calculated through the steady-state stress and pore pressure coupled analysis.

A HA= hA +uwA/γw

GL

O hB h2

O B HB= hB +uwB/γw uw=-γw h2

Fig. 4. The total hydraulic head at the steady state. 4

Tunnelling and Underground Space Technology 98 (2020) 103297

1.5

0

1.2

10 20

0.9

Depth / m

Surface displacement / mm

D. Lu, et al.

0.6 Test data Simulation results

0.3 0.0 0.0

0.1 0.2 0.3 0.4 0.5 Position of the water level / m

60

0.6

γ (kN/m3)

λ

λs

κ

ν

M

pa (kPa)

p′n (kPa)

17

0.012

0.01

0.002

0.3

1.8

101

4000

0

20

40

Sr / %

60

80

100

Fig. 8. The distribution of the saturation at different GL.

linearly with the rise of the GL. The suction exists in the unsaturated soils and is the absolute value of the negative pore water pressure. In the process of the rising GL, the suction gradually decreases to zero.

Table 2 Mechanical parameters of the fine sand.

3.1.1. Saturation of the ground Based on the variations of the pore water pressure along with the GL, the saturation distribution of the numerical model, when the GL is at the depths of 10 m, 20 m, 30 m, and 40 m, respectively, is illustrated in Fig. 8. Below the GL, the ground is regarded as saturated and the saturation is 1. Above the GL, the ground is unsaturated and the saturation is determined through the soil-water characterisitic curve. The shape of the pore water pressure distribution is the same at different depths. The curves of the the saturation distribution can be regarded to move upward gradually along with the GL.

3. Numerical simulation and analysis of the rising GL Taking the rise of the GL in Beijing as the background, a numerical simulation with a homogeneous soil layer, which contians the existing tunnel, is conducted according to the presented numerical method of the rising GL. Based on the numerical simulation results, the mechanical behaviours of the ground and the lining are studied under the condition of the rising GL.

3.1.2. Stress state of the ground In the process of the rising GL, the void of the soil is gradually occupied by the water and the total unit weight of the soil increases. The total vertical stress is adopted to show the the total stress of the ground. Fig. 9 displays the variation of the total vertical stress at the depths of 40 m, 30 m, 20 m, and 10 m, respectively, along with the rise of the GL. Below the GL, for every 1 m rise in the GL, the total vertical stress increases by an average of 5 kPa. Above the GL, the total vertical stress is basically unchanged. This is due to that the saturation above the GL is basically unchanged as shown in Fig. 8. Baesed on the pore water pressure distribution in Fig. 7, the saturation distribution in Fig. 8, and the total vertical stress distribution in Fig. 9, the effective vertical stress distribution can be determined through Eq. (4). The variations of the effective stress is the key factor leading to the ground deforamtion. Fig. 10 depicts the variations of the effective vertical stress at different depths along with the GL. The decreasing rate of the effective vertical stress above the GL is larger than that below the GL. Above the GL, the effective vertical stress is greater

3.1. Mechanical behaviours of the ground with a homogeneous soil layer The mechanical behaviours of the ground, which contain the force and deformation, significantly affect the safe operation of the existing underground structures. In the process of the rising GL, the saturation of the soil gradually increases to 1. The mechanical behaviours of the unsaturated soil are controlled by the suction and the effective stress. The effective stress can be calculated through the total stress, the pore water pressure and the saturation according to Eq. (4). During the rising GL, the main reason for the changes of the total stress and the saturation of soil is the variation of the pore water pressure. In the numerical simulation of the rising GL, the distribution of the pore water pressure is applied to all elements of the model ground as the load condition. Combined with the pore water pressure distribution at the steady state as shown in Fig. 4, the variations of the pore water pressure at the depths of 40 m, 30 m, 20 m, and 10 m, respectively, are illustrated in Fig. 7. The pore water pressure at a certain depth increases

1000 Total vertical stress / kPa

600 400 uw / kPa

GL 10 m 20 m 30 m 40 m

40 50

Fig. 6. Test data and the numerical simulation results of the surface displacement.

200 Depth 40 m 30 m 20 m 10 m

0 -200 -400 40

30

30

20

10

800 600 400 200 0 40

0

Depth 40 m 30 m 20 m 10 m 30

20

10

0

The depth of the GL / m

The depth of the GL / m Fig. 7. The variation of the pore water pressure along with the GL.

Fig. 9. The distribution of the total vertical stress in the process of the rising GL. 5

Tunnelling and Underground Space Technology 98 (2020) 103297

D. Lu, et al.

Effective vertcal stress / kPa

800

Normal contact pressure

600

Depth 40 m 30 m 20 m 10 m

200 0 40

Bending moment

Bending moment

400

30

20

10

Hoop thrust

0

Shear force

Hoop thrust Shear force

Fig. 12. Parameters to define the lining behaviours.

The depth of the GL / m Fig. 10. The distribution of the effective vertical stress in the process of the rising GL.

which includes the hoop thrust, the bending moment and the shear force, to describe the lining behaviours. The normal contact pressure is the normal stress acting on the outer surface of the lining and contains the effective stress of the ground and the pore water pressure. For the internal force of the lining, the length of the section along the excavation direction is taken as 1 m. The radial deformation is the relative displacement between a position of the outer surface of lining and the center of the tunnel. For convenience, ω is defined as the angle between the line, that passing through the research point and the center of the tunnel, and the line, that passing the crown and the tunnel center. ω is measured in the clockwise. According to the relative position of the lining and the GL, the process of the rising GL is divided into three stages: (1) the GL is below the lining and rises from the depth of 40 m to 25 m, (2) the GL is on the lining and rises from the depth of 25 m to 15 m, and (3) the GL is above the lining and rises from the depth of 15 m to the ground surface.

than the total vertical stress due to the existence of the suction. Below the GL, the effective vertical stress is less than the total vertical stress. The difference value between the total vertical stress and the effective vertical stress is the pore water pressure. In the process of the rising GL, the effective vertical stress decreases due to the increase of the pore water pressure. 3.1.3. Deformation of the ground The deforamtion behaviours of the ground have a obvious impact on the safe operation of existing underground structures. As shown in Fig. 11, the distribution of the vertical displacement at different GL is selected to express the features of the ground deformation during the rising GL. Below the GL, the ground is saturated and the mechanical behaviours of ground are only controlled by the effective stress. Thus, the ground is swelling due to that the pore water pressure increases and the effective stress decreases. Above the GL, the ground is unsaturated and the mechanical behaviours of the ground are controlled by the effective stress and the suction. In the process of the rising GL, the saturation increases and the suction decreases. Thus, the bonding among the soil particles gradually weakens and the bearing capacity of the soil reduces. The wetting-induced settlement deformation is produced and the value of it incresases along with the GL.

3.2.1. Rising GL below the lining The rise of the GL below the lining is the process of the rising GL from the depth of 40 m to 25 m. The initial state, the state when the GL rises from the depth of 40 m to 30 m, and the state when the GL rises from the depth of 40 m to 25 m are selected to exhibit the mechanical behaviours of the ground and the lining. (1) Vertical displacement of the ground Fig. 13 displays the nephogram of the vertical displacement. The positive value of the vertical displacement represents the settlement deformation. The negative value of the vertical displacement represents the rebound deformation. As shown in Fig. 13, the wetting-induced deformation mainly occurs near the GL. At the same depth, the wettinginduced deformation of soil around the lining is smaller than that away from the lining. Thus, the lining bears the additive loads due to the uneven deformation of the ground.

3.2. Mechanical responses of the ground with the exisitng tunnel A numerical simulation of the rising GL is conducted through the numerical model of the ground with the existing tunnel. The mechanical response of the ground is shown through the vertical displacement distribution. The mechanical behaviours of the lining are expressed though the normal contact pressure on the outer surface of the lining, the radial deformation of the lining, and the internal force in the lining. Fig. 12 defines the normal contact pressure and the internal force,

(2) Normal contact pressure on the outer surface of the lining Corresponding to the nephogram of the vertical displacement as shown in Fig. 13, the distribution of the normal contact pressure on the outer surface of the lining is illustrated in Fig. 14. The normal contact pressure on the lower part of the lining is larger than that on the supper part of the lining due to the gravity of the lining, and the position of the maximum normal contact pressure is located at the invert. Because of the uneven deformation of the ground as shown in Fig. 13, the normal contact pressure increases obviously. The position with the largest additive normal contact pressure is the invert (ω = 180°). The position with the smallest additive normal contact pressure is the crown (ω = 0°). When the GL rises to 25 m, the additive normal contact pressure at the invert is 144.2 kPa and is 0.7 times of the initial normal contact pressure at the invert.

Fig. 11. The distribution of the vertical displacement along with GL.

6

Tunnelling and Underground Space Technology 98 (2020) 103297

D. Lu, et al.

14 7 0 -7 -14 -14 -7 0 7 14

0

330

30

300

60

270

GL 40 m 30 m 25m

90

240

120 210

180

150

Radial deformation /mm Fig. 15. The radial deformation during the rising GL below the lining.

maximum negative radial deformation increases to 17.0 mm. The deformation mode of the lining is the compression of the crown and the invert.

(a) GL from 40 m to 30 m

(4) Internal force of the lining Fig. 16 illustrates the distribution of the internal force. As shown in Fig. 16(a), the position of the maximum positive bending moment gradually moves from the spring line to the position with the ω of 99° and 261°. The maximum positive bending moment increases from 137.8 kN⋅m to 414.3 kN⋅m and the increment is 2.0 times of the maximum positive bending moment at the initial state. The position of the maximum negative bending moment gradually moves from the crown to the invert. The maximum negative bending moment increases from 126.1 kN⋅m to 451.0 kN⋅m and the increment is 2.6 times of the maximum negative bending moment at the initial state. Corresponding to the distribution of normal contact pressure, as illustrated in Fig. 16(b), the position of the maximum hoop thrust gradually moves from the spring line to the position with the ω of 99° and 261°. The maximum hoop thrust increases from 1037.3 kN to 1741.0 kN and the increment is 0.7 times of the maximum hoop thrust increases at the initial state. The minimum hoop thrust is located on the crown. The hoop thrust at the crown increases from 812.0 kN to 979.0 kN and the increment is 0.2 times of the initial hoop thrust at the crown. Fig. 16 (c) displays the distribution of the shear force. The distribution curve of shear force is smooth. The position of maximum positive shear force moves from the position with the ω of 45° to the position with the ω of 225°. The maximum positive shear force increases from 41.5 kN to 183.3 kN and the increment is 3.4 times of the maximum positive shear force in the initial state. The position of maximum negative shear force moves from the position with the ω of 315° to the position with the ω of 135°. The maximum negative shear force increases from 90.8 kN to 259.8 kN and the increment is 1.9 times of the maximum negative shear force at the initial state. The variations of the bending moment, the hoop thrust and the shear force along with the GL indicate that the rising GL increases the non-uniformity of the internal force distribution obviously.

(b) GL from 30 m to 25 m Fig. 13. The nephogram of the vertical displacement during the rising GL below the lining.

400 330 300 300 200 100 0 270 100 200 240 300 210 400

0

30 60

GL 40 m 30 m 25 m

90 120

150 180 Normal contact pressure / kPa

3.2.2. Rising GL on the lining The rise of the GL on the lining is the process of the rising GL from the depth of 25 m to 15 m. The state when the GL rises from the depth of 40 m to 25 m, the state when the GL rises from the depth of 40 m to 20 m, and the state when the GL rises from the depth of 40 m to 15 m are adopted to present the mechanical responses of the lining and the ground.

Fig. 14. The normal contact pressure during the rising GL below the lining.

(3) Deformation of the lining Fig. 15 shows the the distribution of the radial deformation of the lining. The position of the maximum negative radial deformation is located on the invert. The position of the maximum positive radial deformation is located on the spring line (ω = 90°, 270°). With the rise of the GL, the absolute values of the maximum negative radial deformation and maximum positive radial deformation increase, and the difference between the maximum positive radial deformation and the

(1) Vertical displacement of the ground The nephogram of the vertical displacement is shown in Fig. 17. When the GL rises from the invert (the depth D = 25 m) to the spring line (D = 20 m), as shown in Fig. 17(a), the wetting-induced 7

Tunnelling and Underground Space Technology 98 (2020) 103297

D. Lu, et al.

330 600 300 300 0 -300 -600 270 -300 0 240 300 210 600

0

2100

30 60 90

1400

0

330

30 60

300

700 90

0 270 700

120 180

150

1400

120

240

2100

210

180

40 m

30 m

30 60 90 120 150

180

Shear force / kN

Hoop thrust / kN

Bending moment / kN⋅m

(a)

150

0

300 330 150 300 0 -150 -300 270 -150 0 240 150 210 300

25m

(b)

(c)

Fig. 16. The internal force during the rising GL below the lining.

400 330 300 300 200 100 0 270 100 200 240 300 210 400

0

30 60 90

GL 25 m 20 m 15 m

120

150 180 Normal contact pressure / kPa

Fig. 18. The distribution of the normal contact pressure during the rising GL on the lining.

(2) Normal contact pressure on the outer surface of the lining

(a) GL from 25 m to 20 m

The variations of the normal contact pressure are shown in Fig. 18. The position of the maximum normal contact pressure is the invert and the normal contact pressure near the invert is basically unchanged. The reason is that the ground deformation below the invert is almost the same as shown in Fig. 17. In the process of the rising GL from the invert to the spring line, the normal contact pressure near the spring line reduces slightly due to that the wetting induced deformation mainly occurs at the spring line. Except for area near the invert and the spring line, the normal contact pressure of the lining increases. The maximum increment of the normal contact pressure, which is located at the position with the ω of 135° and 225°, is 50.4 kPa and is 0.2 times of the initial normal contact pressure at same position. In the process of the rising GL from the spring line to the crown, the normal contact pressure below the spring line is basically unchanged, the normal contact pressure near the crown increases, and the normal contact pressure at the remaining area decreases. The maximum increment of the normal contact pressure is located on the crown, is 25.0 kPa, and is 0.1 times of the maximum normal pressure at the invert in the initial state.

(b) GL from 20 m to 15 m

(3) Deformation of the lining

Fig. 17. The nephogram of the vertical displacement during the rising GL on the lining.

Fig. 19 illustrates the distribution of the radial deformation of the lining. When the GL rises from the invert to the spring line, the maximum positive radial deformation is located on the spring line. The position of the maximum negative radial deformation moves from the invert to the crown. The difference between the positive radial deformation and the negative radial deformation increases from 17.0 mm to 24.0 mm. In the process of the rising GL from the spring line to the crown, the radial deformation below the spring line is basically unchanged and the radial deformation above the spring line changes slightly. The position of the maximum negative radial deformation

deformation mainly occurs at the spring line. Below the invert, the vertical displacement is almost the same. The uneven deformation area is located only above the invert. When the GL rises from the spring line to the crown (D = 15 m), as shown Fig. 17(b), the wetting induced deformation mainly occurs at the crown. The uneven deformation of the ground occurs above the spring line. According to Figs. 13(b), 17(a), and (b), the area of the uneven deformation moves upwards along with GL.

8

Tunnelling and Underground Space Technology 98 (2020) 103297

D. Lu, et al.

14 330 7 300 0 -7 -14 270 -7 0 240 7 210 14

0

GL 25 m 20 m 15 m

30 60 90 120

180

150

Radial deformation /mm Fig. 19. The distribution of the radial deformation during the rising GL on the lining.

(a) GL from 15 m to 10 m

moves to the crown and the difference between the positive radial deformation and the negative radial deformation increases from 24.0 mm to 25.4 mm. This feature agrees with the variations of the normal contact pressure. (4) Internal force of the lining The distribution of the internal force is shown in Fig. 20. As shown in Fig. 20(a), when the GL rises from the invert to the spring line, the position of the maximum positive bending moment moves from the position with the ω of 99° and 261° to the spring line. The maximum positive bending moment increases from 414.3 kN⋅m to 591.8 kN⋅m and the increment is 1.3 times of the maximum positive bending moment in the initial state. The maximum negative bending moment, which is located on the invert, increases from 451.0 kN⋅m to 514.5 kN⋅m and the increment is 0.6 times of the initial bending moment at the invert. When the GL rises from the spring line to the crown, the bending moment changes slightly. For the hoop thrust, as shown in Fig. 20(b), the hoop thrust of the invert and the crown is basically unchanged. The maximum hoop thrust, which is located at the position with the ω of 99° and 261°. When the GL rises from the invert to the spring line, the hoop thrust in the lining increases except for the invert and the crown. The maximum hoop thrust increases from 1741.0 kN to 1917.0 kN and the increment is 0.2 times of the initial hoop thrust at same position. When the GL rises from the spring line to the crown, the hoop thrust in the lining decreases except for the invert and the crown. The maximum hoop thrust decreases from 1917.0 kN to 1822.3 kN and the reduction is 0.1 times of the initial hoop thrust at same position. As shown in Fig. 20(c), the distribution of shear force is similar to that of the bending moment. When the GL rises from the invert to the spring line, the maximum positive shear force, which is located at the position with the ω of 225°, increases from 183.3 kN to 220.1 kN and the

0 600 330 300 300 0 -300 -600 270 -300 0 240 300 210 600 180

2100

30 60 90 120 150

1400

(b) GL from 5 m to 0 m Fig. 21. The nephogram of the vertical displacement during the rising GL above the lining.

increment is 1.4 times of the initial shear force at the same position. The maximum negative shear force, which is located at the ω of 135°, increases from 259.8 kN to 309.8 kN and the increment is 0.6 times of that the initial shear force at the same position. When the GL rises from the spring line to the crown, the position of the maximum positive shear force moves to the position with the ω of 45°. The maximum positive shear force increases from 220.1 kN to 251.8 kN and the increment is 0.8 times of the maximum positive shear force in the initial state. The

0

330

30

300

60

700 0 270 700 1400 2100

90

240

120 210

180

150

300 330 150 300 0 -150 -300 270 -150 0 240 150 210 300

25m

20m

(a)

60 90 120 180

15m

(b)

Fig. 20. The distribution of the internal force during the rising GL on the lining. 9

30

150

Shear force / kN

Hoop thrust / kN

Bending moment /kN⋅m

0

(c)

Tunnelling and Underground Space Technology 98 (2020) 103297

D. Lu, et al.

400 330 300 300 200 100 0 270 100 200 240 300 210 400

0

between the maximum positive radial deformation and the maximum negative radial deformation decrease from 25.4 mm to 21.9 mm.

GL 15 m 10 m 0m

30 60

(4) Internal force of the lining

90

Fig. 22 shows the variation of the normal contact pressure. When the GL rises from the crown to the ground surface, the maximum normal contact pressure is located on the invert and increases from 340.0 kPa to 395.5 kPa. The increase of the maximum normal contact pressure is 0.3 times of the maximum normal contact pressure at initial state.

Fig. 24 displays the variations of the internal force. The lining is in the saturated soil and the increase of the pore water pressure is the main reason for the variation of the internal force. Therefore, as shown in Fig. 24(a), the maximum positive bending moment, which is located on the spring line, decreases from 571.5 kN⋅m to 520.5 kN⋅m and the reduction is 0.4 times of the maximum positive bending moment at the initial state. The maximum negative bending moment, which is located on the crown, decreases from 596.3 kN⋅m to 516.0 kN⋅m and the reduction is 0.6 times of the maximum negative bending moment, at the initial state. The variations of the hoop thrust along with the GL are depicted in Fig. 24(b). The maximum hoop thrust, which is located on the invert, increases from 1822.3 kN to 2071.3 kN and the increment of the hoop thrust is 0.2 times of the maximum hoop thrust at the initial state. The minimum hoop thrust, which is located on the crown, increases from 953.5 kN to 1206.0 kN and the increment of the hoop thrust is 0.3 time of the minimum hoop thrust at the initial state. The hoop thrust of the lining increases uniformly with the rise of the GL. The distribution of the shear force is similar to the distribution of the bending moment. The maximum positive shear force, which is located at about the ω of 45°, decreases from 251.8 kN to 215.0 kN and the reduction is 0.9 times of the maximum positive shear force at the initial state. The maximum negative shear force, which is located at about the ω of 315°, decreases from 320.5 kN to 297.8 kN and the reduction of the shear force is 0.3 times of that at the initial state. The variation of the internal force in the lining along with the GL indicates that the rising GL decreases the non-uniformity of the internal force distribution and the load acting on the lining increases. When the GL rises from the initial GL to the invert, the load acting on the lining increases obviously, and the non-uniformity of internal force distribution increases significantly. When the GL rises from the invert to the crown, the load acting on the lining changes slightly. The non-uniformity of internal force distribution increases with the rising GL from the invert to the spring line. The non-uniformity of internal force distribution decreases with the rising GL from the spring line to the crown. When the GL rises from the crown to the ground surface, the load acting on the lining increases uniformly and the non-uniformity of internal force distribution decreases.

(3) Deformation of the lining

4. Conclusion

The variations of the radial deformation are shown in Fig. 23. The ground below the lining is saturated and the pore water pressure increases along with the GL. The increase of the pore water pressure reduces the positive radial deformation on the spring line and the negative radial deformation on the invert and crown. The difference

This paper presents a numerical method to simulate the rise of the GL. In the presented numerical method, a numerical realization method for the rise of the GL is proposed by applying the known pore pressure distribution to the numerical model as the load condition, and an elastoplastic constitutive model of the unsaturated soil is implemented into the finite element code by incorporating the known variation of the pore water pressure. The rationality of the numerical method is validated through a model test. Taking an interval subway tunnel in Beijing as the background, a numerical simulation of rising GL is conducted. Based on the numerical model and the model parameters adopted in this paper, the following conclusions are drawn:

120

150 180 Normal contact pressure / kPa

Fig. 22. The distribution of the normal contact presssure during the rising GL above the lining.

maximum negative shear force is basically unchanged. The shape of the internal force distribution has been changed and the non-uniformity of the internal force distribution increases. 3.2.3. Rising GL above the lining The rise of the GL above the lining is the process of the rising GL from the depth of 15 m to ground surface. The state when the GL rises from the depth of 40 m to 15 m, the state when the GL rises from the depth of 40 m to 10 m, and the state when the GL rises from the depth of 40 m to the ground surface are selected to present the mechanical behaviours of the lining and the ground. (1) Vertical displacement of the ground Fig. 21 displays the nephogram of the vertical displacement. The soil around the lining occurs the rebound deformation and the area of the uneven deformation gradually disappears. (2) Normal contact pressure on the outer surface of the lining

14 330 7 300 0 -7 -14 270 -7 0 240 7 210 14

0

30 60

GL 15 m 10 m 0m

90

(1) When the GL rises from the initial GL to the invert, the normal contact pressure at the invert is the largest and the increment of the maximum normal contact pressure is 0.7 times of the initial normal contact pressure at the invert. In the process of the rising GL from the invert to the crown, the maximum normal contact pressure is basically unchanged. When the GL rises from the crown to the ground surface, the normal contact pressure increases uniformly and the increment of the maximum normal contact pressure is 0.3 times of the maximum normal contact pressure at initial state.

120

150 180 Radial deforamtion /mm

Fig. 23. The distribution of the radial deformation during the rising GL above the lining. 10

Tunnelling and Underground Space Technology 98 (2020) 103297

D. Lu, et al.

330 600 300 300 0 -300 -600 270 -300 0 240 300 210 600

0

2100

30 60 90 120

180

150

1400

330 300

700 0 270 700 1400

240 210

2100

180

15m

0

30

300 330 150 60 300 0 -150 90 -300 270 -150 0 240 120 150 150 210 300

Hoop thrust / kN

Bending moment / kN⋅m

(a)

0

10m

30 60 90 120

180

150

Shear force / kN 0m

(b)

(c)

Fig. 24. The distribution of the internal force during the rising GL above the lining.

interests or personal relationships that could have appeared to influence the work reported in this paper.

(2) The deformation mode of the lining is the compression of the crown and the invert in the process of the rising GL from the initial state to the spring line, and is the compression of the spring line in the process of the rising GL from the spring line to the ground surface. (3) Compared with the initial maximum positive bending moment, the maximum positive bending moment increases by 2.0 times when the GL rises from the initial GL to the invert, and by 1.3 times when the GL rises from the invert to the spring line. In the process of the rising GL from the spring line to the crown, the bending moment changes slightly. When the GL rises from the crown to the ground surface, the reduction of the maximum positive bending moment is 0.4 times of the initial maximum positive bending moment. The variation of the shear force is similar to the variation of the bending moment. (4) Compared with the initial maximum hoop thrust, the maximum hoop thrust increases by 0.7 times when the GL rises from the initial GL to the invert, and by 0.2 times when the GL rises from the invert to the spring line. In the process of the rising GL from the spring line to the crown, the reduction of the maximum hoop thrust is 0.1 times of the initial maximum hoop thrust. In the process of the rising GL from the crown to the ground surface, the increment of the maximum hoop thrust is 0.2 times of the initial maximum hoop thrust.

Acknowledgement This study was supported by the National Natural Science Foundation of China (Grant Nos. 51421005 and 51778026). Appendix A. Supplementary material Supplementary data to this article can be found online at https:// doi.org/10.1016/j.tust.2020.103297. References Abu-Rizaiza, O.S., Sarikaya, H.Z. Ali, Khan, M.Z., 1989. Urban groundwater rise control: case study. J. Irrig. Drain. Eng. 115 (4), 588–607. Al-Rashed, M.F., Sherif, M.M., 2001. Hydrogeological aspects of groundwater drainage of the urban areas in Kuwait City. Hydrol. Process. 15 (5), 777–795. Al-Sefry, S.A., Sen, Z., 2006. Groundwater rise problem and risk evaluation in major cities of Arid Lands-Jedddah Case in Kingdom of Saudi Arabia. Water Resour. Manage. 20 (1), 91–108. Beretta, G.P., Avanzini, M., Pagotto, A., 2004. Managing groundwater rise: Experimental results and modelling of water pumping from a quarry lake in Milan urban area (Italy). Environ. Geol. 45 (5), 600–608. Bob, M., Rahman, N., Elamin, A., Taher, S., 2016. Rising groundwater levels problem in urban areas: A case study from the central area of Madinah city, Saudi Arabia. Arab. J. Sci. Eng. 41 (4), 1461–1472. Brassington, F.C., 1990. Rising groundwater levels in the United Kingdom. Proc. Inst. Civ. Eng. 88 (6), 1037–1057. Chaudhary, M.T.A., 2012. Implications of rising groundwater level on structural integrity of underground structures-investigations and retrofit of a large building complex. Struct. Surv. 30 (2), 111–129. Colombo, L., Gattinoni, P., Scesi, L., 2018. Stochastic modelling of groundwater flow for hazard assessment along the underground infrastructures in Milan (northern Italy). Tunn. Undergr. Space Technol. 79, 110–120. Fang, Y., Guo, J., Grasmick, J., Mooney, M., 2016. The effect of external water pressure on the liner behavior of large cross-section tunnels. Tunn. Undergr. Space Technol. 60, 80–95. Gattinoni, P., Scesi, L., 2017. The groundwater rise in the urban area of Milan (Italy) and its interactions with underground structures and infrastructures. Tunn. Undergr. Space Technol. 62, 103–114. Hefny, A.M., Chua, H.C., 2006. An investigation into the behaviour of jointed tunnel lining. Tunn. Undergr. Space Technol. 21 (3–4), 428. Kreibich, H., Thieken, A.H., 2008. Assessment of damage caused by high groundwater inundation. Water Resour. Res. 44 (9). Kusaka, T., Sreng, S., Tanaka, H., Sugiyama, H., Ito, T., Kobayashi, K., 2016. Experimental study on influence of ground rebound on tunnels caused by groundwater restoration. Jpn. Geotech. Soc. Spec. Publ. 2 (45), 1578–1582. Liao, S., Peng, F., Shen, S., 2008. Analysis of shearing effect on tunnel induced by load transfer along longitudinal direction. Tunn. Undergr. Space Technol. 23 (4), 421–430. Oyedele, K.F., Ayolabi, E.A., Adeoti, L., Adegbola, R.B., 2009. Geophysical and hydrogeological evaluation of rising groundwater level in the coastal areas of Lagos, Nigeria. Bull. Eng. Geol. Environ. 68 (1), 137–143. Selim, S.A., Hamdan, A.M., Rady, A.A., 2014. Groundwater rising as environmental problem, causes and solutions: case study from Aswan city, Upper Egypt. Open J. Geol. 4 (7), 324. Shanahan, P., 2009. Groundwater in the urban environment. In: Baker, L.A. (Ed.), The

The rising GL below the lining has the biggest impact on the mechanical response of the existing tunnel due to the uneven deformation of the ground. Not only the load on the existing subway tunnel is increased, but also the non-uniformity of internal force distribution becomes apparent. The rising GL on the lining has the medium influence on the mechanical response of the existing tunnel. The rising GL above the lining has the relatively small impact influence on the mechanical behaviours of the existing tunnel. The numerical results presented in this paper could be used to explain the damage mechanism of the existing tunnel due to the rise of the GL. In addition, the numerical method of the rising GL can be used to research the response of the underground structure during the rising GL in the future. CRediT authorship contribution statement Dechun Lu: Methodology, Writing - review & editing, Funding acquisition, Resources. Xiaoqiang Li: Conceptualization, Software, Data curation, Writing - original draft. Xiuli Du: Project administration, Supervision. Qingtao Lin: Software, Visualization, Data curation. Qiuming Gong: Supervision. Declaration of Competing Interest The authors declare that they have no known competing financial 11

Tunnelling and Underground Space Technology 98 (2020) 103297

D. Lu, et al. Water Environment of Cities. Springer, Boston, MA, pp. 29–48. Shin, J., Kim, S., Shin, Y., 2012. Long-term mechanical and hydraulic interaction and leakage evaluation of segmented tunnels. Soils Found. 52 (1), 38–48. Shin, J.H., Potts, D.M., Zdravkovic, L., 2005. The effect of pore-water pressure on NATM tunnel linings in decomposed granite soil. Can. Geotech. J. 42 (6), 1585–1599. Sun, D.A., Sheng, D., Sloan, S.W., 2007. Elastoplastic modelling of hydraulic and stress–strain behaviour of unsaturated soils. Mech. Mater. 39 (3), 212–221. Sun, D.A., Sheng, D., Xiang, L., Sloan, S.W., 2008. Elastoplastic prediction of hydro-mechanical behaviour of unsaturated soils under undrained conditions. Comput. Geotech. 35 (6), 845–852. Teachavorasinskun, S., Chub-uppakarn, T., 2010. Influence of segmental joints on tunnel lining. Tunn. Undergr. Space Technol. 25 (4), 490–494.

Vanapalli, S.K., Fredlund, D.G., Pufahl, D.E., 1999. The influence of soil structure and stress history on the soil–water characteristics of a compacted till. Géotechnique 49 (2), 143–159. Wilkinson, B., 1985. Rising groundwater levels in London and possible effects on engineering structures. In: Proceedings of the 18th Congress of the International Association of Hydrogeologists, pp. 145–157. Yasuhara, K., Murakami, S., Mimura, N., Komine, H., Recio, J., 2007. Influence of global warming on coastal infrastructural instability. Sustain. Sci. 2 (1), 13–25. Zhang, D.M., Ma, L.X., Zhang, J., Hicher, P.Y., Juang, C.H., 2015. Ground and tunnel responses induced by partial leakage in saturated clay with anisotropic permeability. Eng. Geol. 189, 104–115.

12