Seismic responses of rectangular subway tunnels in a clayey ground

Observations from past earthquakes have shown that subway tunnels can suffer severe damage or excessive deformation due to seismic shaking. This study presents the results of finite element analyses on subway tunnels installed in normally consolidated clay deposits subjected to far-field ground motions. The clay strata were modelled by a hyperbolic-hysteretic constitutive model. The influences of three factors on the seismic response of the clay-tunnel systems were examined, namely ground motion intensity, tunnel wall thickness and clay stiffness. Furthermore, the computed racking deformations of the tunnel were compared with several analytical estimations from the literature, and the relationship between racking ratio and flexibility ratio for rectangular tunnels installed in normally consolidated clay deposits was proposed. The findings may provide a useful reference for practical seismic design of tunnels.


Introduction
To ensure sustainable development, much attention has been paid to constructing more tunnels and other underground structures in view of the relatively limited land resources and the increasing number of immigrants encountered in most of metropolitan cities. The performance of these tunnels and underground structures against natural hazards (e.g. earthquakes) deserves to be examined fully. Many previous publications have reported that tunnels and underground structures can suffer severe damage or excessive deformation due to seismic shaking [1][2][3][4]. To date, the guidelines for seismic design of tunnels predominantly rely on simplified methods [3,[5][6][7][8][9][10][11][12][13]. At the absence of the dynamic soil-structure interaction mechanism, those methods may result in a considerable overestimation or underestimation of the seismic response of a tunnel [14][15].
Numerous studies have been performed to investigate the seismic response of tunnels considering the dynamic soil-structure interaction, using a 1-g shaking table or centrifuge tests [14,[16][17][18][19][20][21][22][23][24][25][26] and numerical analyses [4,[14][15][27][28][29][30][31][32][33][34][35][36][37][38]. Most of these studies mainly focused on the tunnels installed in dry sands or liquefiable soils. Only a few exceptions can be found. For instance, Amorosi and Boldini [30] performed 2D FE simulations to investigate the influences of soil constitutive models (namely visco-elastic and visco-elasto-plasic effective stress models) a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 on the seismic response of a circular tunnel embedded in clayey soils; Moss and Crosariol [22] conducted 1-g shaking table tests to examine the seismic behaviour of rectangular tunnels installed in soft clay deposits, focusing on responses of the ground surface acceleration and lateral distortion of the tunnel structure; Chen et al. [24] carried out both 1-g shaking table tests and 3D FE analyses to investigate the damage mechanism of a 3-storey and 3-span subway structure installed in soft clays subjected to strong ground motions, with the focuses placed on the acceleration and lateral displacement responses of the subway structure. In general, compared to a sand deposit, the clay deposit's stiffness is comparatively smaller and the hysteresis damping effect is more significant. Besides, when subjected to strong ground motions, the excess pore pressures generated within clay deposits may remain at a relatively low level, and liquefaction of the clay deposits is not likely to happen; this seismic behaviour of the clay deposit is significantly different from that of a saturated sandy deposit. For these reasons, the clay deposits will be considered in this study.
In coastal cities such as Singapore and Shanghai, where large portions of the land are underlain by thick layers of normally consolidated clays, many superstructures and underground structures are built on/in normally consolidated clay deposits. A distinct feature of normally consolidated clay is its significant amplification effect on earthquake-induced ground motions for small to moderate ground motion intensities [39][40][41][42][43]. As a result of this amplification, the structures built on or installed in this type of clay strata may suffer amplified loading arising from the surrounding soils. Furthermore, given that the relevant studies involving underground structures in a normally consolidated clay stratum are still relatively few, especially for those involving far-field ground motions, the seismic interaction mechanism of rectangular tunnel-clay system has yet to be fully understood. Hence, it will be beneficial to carry out such research to advance the understanding of their seismic performance.
In this study, a series of two-dimensional (2D) plane strain finite element (FE) analyses was performed to investigate the seismic response of rectangular tunnels that are installed in normally consolidated clay deposits and subjected to far-field ground motions. A rectangular tunnel is widely encountered in real practice; for instance, the cross-river tunnel (i.e. Hejiangtao Tunnel) in Hengyang city (China) is a rectangular shape with 11.3m in tunnel height and 10.3m in tunnel width. The clay was modelled following a hyperbolic-hysteretic soil model, and the basic properties are listed in Table 1. This hyperbolic-hysteretic model was developed based on results from a series of resonant column and cyclic triaxial tests on kaolin and soft marine clays [43][44], which can favourably well replicate the dynamic behaviour of soft clay under seismic loadings, namely hysteresis damping, strain level-and loading cycle-dependent stiffness reduction, and significant amplification effect on seismic-induced motion. The current study focused on the ground acceleration response and tunnel structure responses, with the influences of ground motion intensity, tunnel wall thickness and clay stiffness taken into consideration; seismic tunnel structure responses comprise the shear force, bending moment and lateral deformation of the left and right tunnel walls during the seismic shaking events.

Input far-field ground motions
As Fig 1 shows, three types of base motions were adopted in this study. The base motion with a peak acceleration of 0.06g shown in Fig 1(a) represents a type of far-field shaking event that may be felt in Singapore. It is a synthesized signal with reference to the past seismic-induced motions measured in Singapore, which was also used by Zhang [45] and Zhang et al. [46] for both the centrifuge tests and FE analyses. In addition, the ground motions shown in Fig 1(b) and 1(c) were compiled based on two horizontal earthquake records (east-west direction) measured at Stations El Monte-Fairview Av and TTN005 during the 1992 Landers and 1999 Chi-Chi earthquakes, respectively. These two earthquake records were selected from the Ground Motion Database provided by the Pacific Earthquake Engineering Research (PEER) Centre. They have site-to-source distances of approximately 136 km and 83 km, respectively, which fall within the scope of far-field ground motion according to Dadashi and Nasserasadi [47] and Bhandari et al. [48]. For simplicity, these ground motions are termed Type-1, Type-2 and Type-3 base motions. The response spectra corresponding to a peak acceleration of 0.06g are plotted in Fig 1(d). The displacement plots of the input base motions are also provided in Fig  2. As can be seen, due to the differences in time duration and frequency content, the peak displacements are generally different even though the peak accelerations are same. In addition, to investigate the influence of ground motion intensity, each of these base motions was scaled into four ground motions with peak values of 0.03g, 0.06g, 0.12g and 0.24g.

Numerical modelling procedure
The present numerical analysis consists of two phases, namely gravity loading and earthquake loading phases. Initial stress condition is assigned to the soil domain to account for the gravitational stress field of the clay-tunnel system in the gravity loading phase; subsequently, far-field ground motions are excited at the model base in the earthquake loading phase. Some other information on the numerical modelling procedure is given below.

Basic information
As Fig 3 shows, the 2D plane strain FE model of the clay-tunnel systems employed in this study was set up using ABAQUS, which contains 7900 linear quadrilateral elements and 140 linear beam elements to model the clay and rectangular tunnels, respectively. Although there are many different suggestions (e.g. [49][50]) concerning the element size in seismic numerical analyses, it is generally accepted that the element number per seismic wavelength of 10 is sufficient for accurately modelling the seismic wave propagation. This criterion is strictly followed in the present study (refer to reference [51] for further relevant information). In addition, although not included herein, a pure clay FE model was also employed to obtain the free-field ground acceleration response for comparison purposes. As Fig 3 shows, the rectangular tunnel is divided by a central column into two square parts and embedded at a depth of 10 m below the ground surface, which was adopted to represent a subway station. As Table 2

Boundary conditions
To minimize the possible boundary reflection effect arising from the truncated lateral boundaries, with the reference to the previous studies by Tsinidis et al. [25], Bao et al. [37] and Tsinidis [36], the lateral boundaries were adopted to have a dimension 12 times the tunnel width. In addition, the nodes at the same depths on the two lateral boundaries were constrained to have the same horizontal and vertical displacements. This approach has been validated to be effective in replicating a free-field ground motion, and was adopted in many previous studies [15,33,42,51]. In addition, a rigid boundary condition was applied at the bottom of the FE model so that the model base simply moves along the direction of the shaking.

Clay-tunnel interface
Like the studies by Debiasi et al. [33] and Tsinidis et al. [15], hard-contact and penalty friction algorithms were used to model the normal and tangential behaviours of the soil-tunnel interface, respectively. For tangential behaviour, the friction coefficient was set equal to 0.4663, which was obtained by assuming a frictional angle the same as that of clay. The effect of friction coefficient is not included in this study, while preliminary separate numerical analyses suggest that the effect of friction coefficient is minimal for input base motions with peak base accelerations (PBA) less than or equal to 0.12 g and considerable for PBA = 0.24 g. Following the in-flight T-bar test results during centrifuge tests on normally consolidated clays reported by Zhang [45] and Zhang et al. [46], the shear stress limit was set equal to 0.24 σ' v in kPa, where σ' v is the effective overburden pressure of the clay at the depth of interest.

Clay constitutive model
As Fig 4 shows, the normally consolidated clay in this study was modelled by a hyperbolic-hysteretic soil model, which was developed based on a series of test results from cyclic triaxial and resonant column tests on soft clays [44]. The basic shear stress-strain relationship of this model can be expressed by: where q and ε s are the current deviatoric stress and generalized shear strain, respectively; q r1 and q r2 are the respective deviatoric stresses at the reversal points; ε r1 and ε r2 are the respective generalized shear strains at the reversal points; G max is the small-strain shear modulus; q f is the deviatoric stress at failure.  The deviatoric stress q and generalized shear strain ε s can be respectively represented by the following equations: where σ', ε and τ respectively denote the effective normal stress, strain and shear stress, with subscript showing the direction. For normally consolidated clay, following Viggiani and Atkinson [52] and Banerjee [44], the small-strain (or maximum) shear modulus can be expressed by the following formula: where p' 0 is the initial mean effective normal stress, and the units for both G max and p' 0 are kPa; A is the calibration constant or stiffness parameter. For kaolin clay, the parameter A equals 2060 [44]. This dynamic soil constitutive model was coded into a VUMAT subroutine in ABAQUS for explicit dynamic analysis, which can significantly save the computational time and resource involving seismic soil-structure interaction as compared to the normally used implicit dynamic analysis [51,53]. For example, compared to the implicit dynamic analysis, the computational time involving the present explicit analysis can be significantly reduced by approximately 4 times for the 2D clay-tunnel model considered in this study. More detailed information on this soil model can be found in references [42,44].

Analysis results
Three influencing factors on the seismic response of clay-tunnel systems were considered in the present study; that is, tunnel wall thickness, ground motion intensity and clay stiffness. The instantaneous profiles of the maximum bending moment and shear force of the tunnel wall correspond to the instant when the tunnel wall attained its respective maximum values. Unless otherwise stated, the parameter A is set as 2060. base motions with PBAs ranging from 0.03 g to 0.24 g, where the deviatoric stress was computed using Eq (2). As can be seen, the hysteresis loop and variation in the effective normal stress become more significant with the increasing ground motion intensity. Fig 6 also plots   the corresponding shear stress versus shear strain curves, which indicates that the stiffness degradation with increasing shear strain and hysteresis damping can be reasonably captured. This study focuses mainly on the responses of acceleration at different depths and tunnel structure, and in the following subsections the three influencing factors on them are to be investigated.  Table 3, subjected to the Type-1 base motion with PBA of 0.06 g. These profiles were obtained by selecting the peak values experienced by the nodes along the symmetrical axis shown in Fig 3. As also reported by many other researchers [39][40][41][42] , Fig 7 suggests that the clay stratum can significantly amplify the input base motion with a general increasing trend from clay base to the ground surface, regardless of the existence of the tunnel structure or the tunnel wall thickness. In addition, as Fig 7 suggests, the ground acceleration response can be influenced considerably by the dynamic clay-tunnel interaction, especially for the clay nodes located above and underneath (depths equal to 10 m and 15 m, respectively) the tunnel structure, respectively. The influence of dynamic clay-tunnel interaction is also clearly demonstrated in Fig 8, plotting the influence of tunnel thickness on the ground acceleration amplification factors experienced at different depths, where the amplification factor is defined as the ratio of the peak acceleration experienced by a soil node to the peak acceleration of input base motion (i.e. PBA). The acceleration amplification factor against tunnel wall thickness can be represented approximately by the

Influence of tunnel wall thickness
where f a is the acceleration amplification factor; L and H are the thickness and height of the tunnel wall, respectively; L/H is the normalized tunnel wall thickness. The ground acceleration amplification factor at the node underneath tunnel structure against the tunnel thickness shows a general increasing trend, which is reasonable in that a larger portion of the seismic wave tends to be reflected into the soil domain by a stiffer tunnel structure as the wall becomes thicker. However, the opposite slight decreasing trend is found for both the nodes above the tunnel structure and at the ground surface because more ground motion energy is entrapped in the soil domain underneath the tunnel with stiffer walls. The clear discrepancies in the acceleration results shown on Figs 7 and 8 indicate the importance of considering the dynamic soil-tunnel interaction that is also suggested by Ulgen et al. [14] and Tsinidis et al. [15]. Besides, ground amplification effect is generally more significant for input base motion with PBA = 0.03g than that with PBA = 0.24g, while the influence of tunnel wall thickness seems to be more evident associated with PBA = 0.24g. This difference is likely Seismic responses of rectangular subway tunnels due to that the stiffness and damping ratio associated with the hyperbolic-hysteretic soil model is strain-dependent.
where M max and F max are the respective maximum bending moment and shear force of the tunnel wall; EI and GA are the flexural rigidity and shear rigidity of the tunnel wall, respectively; M max L/(2EI) and F max /(GA) are the maximum bending strain and maximum shear strain experienced by the tunnel wall, in percentage terms. The instantaneous maximum lateral deformation profiles of the tunnel walls, which generally follow a relatively simple trend and can be approximately represented by a liner line, are not presented herein. Instead, the results of the maximum lateral deformation of the tunnel wall, which is defined as the maximum relative lateral displacement between the two ends of the tunnel wall throughout the respective seismic shaking event, are presented in this study. Fig 13 shows the plots of the maximum normalized lateral deformation of the lateral wall Seismic responses of rectangular subway tunnels against the normalized tunnel wall thickness for clay-tunnel systems subjected to the three types of base motions with PBAs of 0.03g and 0.24g, which suggests that the relationship between the maximum lateral deformation of the tunnel wall and the tunnel wall thickness can be approximately represented by the following equations: where U max is the maximum lateral deformation of the tunnel wall; U max /H is maximum normalized lateral deformation of the tunnel wall, in terms of percentage.
Furthermore,  show that the influences of tunnel wall thickness on the seismic response of tunnel structure are generally different for base motions with varying PBAs, indicating that the seismic clay-tunnel interaction depends on the ground motion intensity. In the next section, more computed results pertaining to the effect of ground motion intensity on the clay-tunnel systems are to be reported.

Influence of ground motion intensity
Four different-intensity ground motions are used to study the effect of ground motion intensity, the PBAs ranging from 0.03g to 0.24g, as shown in both Fig 1 and Table 3 Seismic responses of rectangular subway tunnels shows, ground acceleration response can also be significantly influenced by the PBA, with a clearly increasing trend. In particular, accelerations experienced at the location underneath the tunnel structure seem to be more sensitive to the variation in PBA, which are generally smaller than those experienced above the tunnel for PBAs less than or equal to 0.06g and then become comparatively larger for PBAs of 0.12g and 0.24g. Fig 15 shows the plots of the acceleration amplification factor against PBA at the different depths, which suggest a general decreasing trend. In addition, the acceleration amplification factors against PBA can be where g is the gravitational acceleration, and PBA/g is the normalized PBA. Seismic responses of rectangular subway tunnels wall at the clay-tunnel interface and a crossover point of the distribution of dynamic stress arising from the surrounding clays. The earthquake-induced kinematic force upon the tunnel structure is highly dependent on the motion intensities of the surrounding soils. In this study, the free-field peak acceleration at depth of the tunnel centre (depth = 12.5 m), termed PCA, is used to represent the average motion intensity of the soils at the tunnel level. As Figs 18 and 19 show, the maximum bending moment and the shear force of the tunnel wall can be well correlated with the PCA by the following two equations: where PCA/g is the normalized PCA.
Furthermore, Fig 20 also plots the maximum normalized lateral deformation of the tunnel wall against the PCA, which suggests a good relationship between them, as shown below.
As Figs 14 and 16-20 show, the ground motion intensity generally has a clear increasing trend of influence on the seismic response of the tunnel structure and the ground acceleration response. However, the ground acceleration amplification factors (Fig 15) experienced at different depths tend to gradually decrease with the increasing ground motion intensity, which is  likely due to the nonlinear stress-strain relationship, stiffness reduction and damping augmentation with the increasing strain level associated with the adopted hyperbolic-hysteretic clay model.

Influence of clay stiffness
The influence of the clay stiffness was considered by varying the parameter A in Eq (4), the values of which are listed in Table 3. Due to space constraints, the profiles of the peak ground acceleration, the maximum bending moment and the shear force of the tunnel wall associated with different clay shear moduli are not presented herein. Fig 21 shows the trends of the acceleration amplification factor and maximum shear force versus the maximum shear modulus of clay. Both the acceleration amplification factor and the maximum shear force tend to increase steadily with the increasing clay stiffness up to a critical value, beyond which they tend to have the opposite trends or remain relatively unchanged with the increasing clay stiffness. Similar findings are also observed in Fig 22, showing the trends of the maximum lateral deformation and bending moment against the clay stiffness. For the plots shown in Figs 21(a) and 22(a) involving the Type-1 base motion, the critical value of the maximum shear modulus of clay can be expressed approximately by the following equation. According to Banerjee [44] and Banerjee et al. [42], Eq (13) can well represent the maximum shear modulus of normally consolidated Malaysian kaolin clay.
To better understand the trends demonstrated in Figs 21 and 22, it is helpful to compare the fundamental period of the clay stratum with the dominant period of the input base motion. Fig 23 plots the typical peak shear strain profile along the clay depth for pure clay bed subjected to the Type-1 base motion with PBA of 0.03 g. As can be seen, the peak shear strains are generally very small, with the maximum value less than 0.08%. For such small shear strain levels, the stiffness degradation effect is insignificant and can be neglected, which is also demonstrated in Fig 6(a). Hence, the fundamental period calculated based on the maximum soil shear modulus is approximately the same as that of clay stratum subjected to small base motions with PBA of 0.03g. The depth-averaged maximum soil shear modulus G soil may be represented by the following equation: where H soil is the thickness of the clay from the ground surface to the clay base, with the unit of m; G max is the depth-dependent maximum and can be calculated based on Eq (4) and the stiffness parameters listed in Table 3. Substituting the critical clay stiffness represented by Eq (13) into Eq (14) leads to: According to Kramer [54], the fundamental period of a soil stratum can be estimated using the following equation: where T soil is the fundamental period of a soil stratum, with the units of s; V s is the equivalent shear wave velocity in a soil stratum, with the units of m/s, which can be obtained following the equation: where ρ soil is the density of soil, with the units of kg/m 3 . By substituting the soil thickness and stiffness parameters into Eq (16), the fundamental periods of the clay deposits corresponding to different maximum shear moduli can be obtained, as shown in Fig 24.  Fig 24 shows that, for maximum shear modulus of clay less than the critical shear modulus corresponding to the Type-1 base motion, the fundamental period of the soil stratum decreases from approximately 2 s to 1 s, which gradually becomes closer to the dominant period (approximately 0.9 s, as Fig 1 shows) of the Type -1 base motion. Beyond the critical shear modulus, the fundamental period of the soil stratum continues to decrease gradually until down to approximately 0.5 s corresponding to the stiffness parameter A equal to 8240, which has a trend of augmenting the difference with the dominant period of the input base motion. In general, smaller differences between the fundamental period of the soil stratum and the dominant period of the input base motion lead to more intense seismic response, and vice versa. Hence, the trends shown in Figs 21 and 22 of both the ground acceleration and the tunnel structure responses against the clay stiffness are highly frequency-dependent and can be reasonably well explained by comparing the fundamental period of the clay deposit with the dominant period of the input base motion.

Comparison with the existing analytical approaches
As mentioned in Section Introduction, the experimental and numerical data of the seismic behavior of rectangular tunnel installed in soft clay deposits is very limited. In this study, the computed racking deformation results of the tunnel structure are compared with the estimations of analytical approaches developed by Wang [5], Penzien [6] and Anderson et al. [11].
For relatively stiff tunnels where the rotation of the tunnel structure can be neglected, the maximum lateral deformation of the tunnel wall is also termed racking deformation or distortion [5][6]. However, as suggested by Debiasi et al. [33] and Ulgen et al. [26], for relatively flexible tunnel structures, the lateral deformation component due to rocking rotation should be subtracted from the maximum lateral deformation to obtain the racking deformation of tunnel structure. In this study, the effect of tunnel rotation on the maximum lateral deformation of the tunnel wall was found to be insignificant and hence neglected. The racking ratio of a tunnel can be represented by the following equation: where R is the racking ratio of a tunnel, Δ str is the racking distortion of a tunnel, and Δ str is the free-field ground distortion between the depths, respectively equal to tunnel top and bottom depths. The racking ratio is usually correlated with the flexibility ratio F, which is expressed by the following equation: where W is the tunnel width, H is the tunnel height, S is the force required to cause a unit racking deflection of the tunnel, computed through a simple static elastic frame analysis, and G s is the soil shear modulus.
For the hyperbolic-hysteretic soil model used in this study, a degraded soil shear modulus should be used for G s , which can be estimated using the following equation [44]: Fig 25 presents plots of the racking ratio against the flexibility ratio obtained from the present FE analyses and using analytical approaches developed by Wang [5], Penzien [6] and Anderson et al. [11]. As can be seen, for flexibility ratios less than about 1.5, notwithstanding the local discrepancies, the overall trend of the racking ratio against the flexibility ratio obtained from this study compares favorably well with these suggested by the analytical solutions. On the other hand, for flexibility ratios larger than 2, all the three analytical approaches tend to considerably overestimate the racking ratio of tunnel structure embedded in soft clay deposit, with the closest prediction given by Anderson et al. [11]. Furthermore, as Fig 25(d) shows, the computed racking ratio versus the flexibility ratio in this study can be well represented by the following expression: Eq (21) is useful in estimating the racking distortion of a rectangular tunnel embedded in soft deposits subjected to seismic ground motions. It can be estimated that the predictions of tunnel racking ratio using Eq (21) are about 12.5% less than that suggested by Anderson et al. [11].

Conclusions
A series of 2D FE parametric studies incorporated with a hyperbolic-hysteretic soil constitutive model were performed to investigate the seismic response of rectangular tunnels installed in clay strata subjected to three types of far-field ground motions, in which three crucial factors, namely tunnel wall thickness, peak base acceleration and clay stiffness, were considered. All the three factors have significant influence on both the responses of ground acceleration and tunnel structure. Some specific conclusions are summarized below. 1. Due to the effect of clay-tunnel interaction, the ground acceleration response with the presence of tunnel structure is significantly different than the free-field ground acceleration response, especially for the locations at the depths near the tunnel level.
2. The soft clay deposits can significantly amplify the earthquake-induced ground motions, regardless of the values of tunnel wall thickness or PBA, with the maximum acceleration amplification factor found to be larger than 3 (as shown in Figs 8a and 15a).
3. The influence of clay stiffness on the seismic response of the clay-tunnel system is highly frequency-dependent. In the future study, random FE analyses (e.g. Zhang et al. [55], Liu and Zhang [56]) may be employed to better investigate the influences of soil properties on the seismic response of the clay-tunnel system, which may facilitate the reliability-based design for tunnels installed in clay deposits. Seismic responses of rectangular subway tunnels 4. Some useful correlations are formulated using dimensionless terms to relate the maximum seismic responses of clay-tunnel systems to each of the influencing factors considered, which may provide a useful reference for the practical seismic design of tunnels installed in soft clay deposits.
5. The existing analytical approaches developed by Wang [5], Penzien [6] and Anderson et al. [11] can favourably predict the racking deformation response of rectangular tunnel installed in soft clay deposits for flexibility ratios less than 1.5, beyond which these approaches tend to significantly overestimate the racking deformation response. In view of this, a new semi-empirical equation to correlate the racking ratio with flexibility ratio for rectangular tunnels embedded in soft clay deposits subjected to seismic shakings is derived, with the aid of regression analysis using the present FE analysis results.
The present study focuses on the seismic performance of clay-tunnel systems subjected to far-field ground motions. Hence, the findings and correlations presented in this study may be valid only for rectangular tunnels installed in soft clay deposits subjected to the earthquakes with the frequency contents comparable to that adopted in this study, and care should be exercised when applying the findings drawn in this study to other types of earthquakes. Future studies can be performed to further examine the seismic behaviour of the clay-tunnel systems subjected to different earthquakes with distinct dominant periods.