Seismic Fragility Analysis of Tunnels with Different Buried Depths in a Soft Soil

Seismic fragility of an engineering structure is the conditional probability that damage of a structure equals or exceeds a limit state under a specified intensity motion. It represents the seismic performance of structures and the correlation between ground motion and structural damage, playing an indispensable role in structural security assessment. A practical evaluation procedure of acquiring the fragility curves of tunnels in a soft soil has been proposed in this paper. Taking a typical metro tunnel in Shanghai as an example, two-dimensional finite element models of soil-tunnel cross-section were established. The nonlinear characteristics of soil layers were considered by the one-dimensional equivalent linear analysis in the equivalent-linear earthquake site response analyses (EERA) program. The ground motions were selected based on seismic station records. Comparing the analytical fragility curves with the empirical curves derived from American Lifelines Alliance (ALA) and HAZUS shows that the proposed method is reliable and feasible. Further study about the influence of buried depths on the fragility of the tunnel was performed. The results indicate that the failure probability of the tunnel is not monotonically decreasing with the increase of the buried depth for a given peak ground acceleration (PGA.)


Introduction
An earthquake is a serious natural disaster that constitutes a threat for urban infrastructure and people's life. People have insufficient understanding about seismic resistance of underground structures. It is generally believed that an underground structure is much safer than a ground structure because of smaller deformations under the condition of surrounding rock or soil constraints. Since the 1990s, many devastating earthquakes causing serious damage to subway stations and tunnels have occurred, such as the Kobe earthquake, Chi-Chi earthquake and Wenchuan earthquake, which means that underground structures are still susceptible to damage under severe earthquakes. According to actual cases of tunnel failure, it has been obtained that shallow tunnels and tunnels built in the unfavorable geological layer are more vulnerable [1]. Meanwhile, it is very difficult to repair when the tunnel is seriously damaged. Therefore, analyzing the seismic vulnerability of tunnels is necessary.
Seismic fragility of engineering structures is the conditional probability that damage of a structure equals or exceeds a limit state under a given seismic input ground motion. It not only quantitatively represents the seismic performance of engineering structures from the probability, but also describes the relationship between the intensity of ground motion and the degree of structural damage from the macro perspective. The vulnerability assessment of tunnels has been mainly based on empirical fragility curves or expert judgments [2][3][4]. American Lifelines Alliance (ALA) and the organizing committee of HAZUS system obtained empirical fragility curves of tunnels for different types of 1.
The suite of applicable ground motion is chosen for analysis according to the site information.

2.
Free site response analysis is conducted through a 1D equivalent linear analysis on the soil profiles, paired with input motions at different seismic intensity levels, to acquire the dynamic characteristics of soil layers. The nonlinear soil behavior is considered through equivalent shear modulus G and damping ratio ζ compatible with shear strain.

3.
Establish the finite element soil-tunnel model and evaluate seismic response by dynamic time history analysis in a finite element software such as ABAQUS.

4.
Selection of the intensity measure (IM) parameter of input ground motions for developing the fragility curves.
Sustainability 2020, 12, 892 3 of 17 5. Definition of damage states (DS) corresponding to the damage index (DI) for determining the damage of tunnel lining. 6.
The damage assessment can be performed based on lognormal approach with linear regression to determine the relationship of DI and IM.

7.
A lognormal standard deviation (β tot ) that describes the total variability associated with each fragility curve must be estimated. 8.
These outputs are obtained and used for developing the fragility curves of tunnels.
Sustainability 2020, 12, x FOR PEER REVIEW 3 of 17 6. The damage assessment can be performed based on lognormal approach with linear regression to determine the relationship of DI and IM. 7. A lognormal standard deviation (βtot) that describes the total variability associated with each fragility curve must be estimated. 8. These outputs are obtained and used for developing the fragility curves of tunnels.

Intensity Measure (IM) Selection
The peak ground acceleration (PGA), peak ground velocity (PGV), peak ground displacement (PGD) or response spectral acceleration (SA) are usually chosen as an intensity measure for seismic fragility in engineering applications. The PGA that is peak acceleration of input ground motion in the time history is a good index. ALA produced empirical fragility curves by PGA, and PGA is selected as the IM in many tunnel fragility studies [2,5,7,24]. Thus, PGA is adopted as the IM for calculating the conditional probability of tunnel lining failure in this paper.

Definition of Damage States (DS)
The definition of damage states, which are described by a damage index expressing relationship between the actual response and the capacity of the lining, and selection of corresponding damage indicator, are critical for vulnerability assessment of engineering structures. Generally, responses of structures can be evaluated by engineering demand parameters (EDP), which are useful to predict the damage to structures. Although previous researches have defined various damage indexes employed in seismic fragility analysis of buildings and bridges, such as maximum displacement, the inter-story drift ratio, and stress for evaluating the EDP [25,26], no such information is available for tunnels. In the process of developing empirical fragility curves, the damage states of tunnels are often divided based on the damage degree of tunnels in past earthquakes. The process of deriving numerical fragility curves is inconsistent with the empirical curves.

Intensity Measure (IM) Selection
The peak ground acceleration (PGA), peak ground velocity (PGV), peak ground displacement (PGD) or response spectral acceleration (SA) are usually chosen as an intensity measure for seismic fragility in engineering applications. The PGA that is peak acceleration of input ground motion in the time history is a good index. ALA produced empirical fragility curves by PGA, and PGA is selected as the IM in many tunnel fragility studies [2,5,7,24]. Thus, PGA is adopted as the IM for calculating the conditional probability of tunnel lining failure in this paper.

Definition of Damage States (DS)
The definition of damage states, which are described by a damage index expressing relationship between the actual response and the capacity of the lining, and selection of corresponding damage indicator, are critical for vulnerability assessment of engineering structures. Generally, responses of structures can be evaluated by engineering demand parameters (EDP), which are useful to predict the damage to structures. Although previous researches have defined various damage indexes employed in seismic fragility analysis of buildings and bridges, such as maximum displacement, the inter-story drift ratio, and stress for evaluating the EDP [25,26], no such information is available for tunnels. In the process of developing empirical fragility curves, the damage states of tunnels are often divided based on the damage degree of tunnels in past earthquakes. The process of deriving numerical fragility curves is inconsistent with the empirical curves.
Based on the common failure modes of tunnels, the damage index is defined as the ratio of the actual (M) and the capacity (M RD ) bending moment of a tunnel cross-section, according to some existing studies [15,27]. The mechanical behavior of a tunnel lining is simulated by an elastic beam subjected to deformations imposed by the oscillating surrounding ground due to seismic waves propagating perpendicular to the tunnel axis. The actual bending moment is calculated by the dynamic time history analysis. The capacity of the bending moment can be obtained by considering the materials and geometry properties of the tunnel section. Five different damage states are determined, namely, none, minor, moderate, extensive and collapse. The damage states and corresponding damage index are shown in Table 1.

Derivation of Fragility Function
Many engineering structural fragility studies show that the following relationship is satisfied between the structural engineering demand parameter (EDP) and the ground motion parameter (IM): where the α and β are regression parameters. Median value of damage index EDP m is a special structural engineering demand parameter. The following exponential relationship between EDP m and IM can also work: Taking the logarithm on both sides of the above Equation (2), thus: The above parameters can be calculated by the data of finite element analysis and obtained by linear regression.
The currently applied approach consists in response function of structures obeying lognormal distribution, that is: [ln(EDP) − ln(EDP m )] 2 are the logarithm median and logarithm standard deviation of EDP, respectively.
Similarly, the probability function of the capacity (C) of structure can also be expressed by lognormal distribution function. Therefore, the function can also be expressed by two parameters µ c and σ c as follows: where µ c and σ c are the logarithm median and logarithm standard deviation of C, respectively. Seismic fragility curve of structures is a statistical tool representing the conditional probability that structural response EDP exceeds the capacity (C), defined by the damage states given a specified intensity level (IM), and it is expressed in the following equation [4,5,28]: The Equation (6) can be rewritten as P f = P(C − EDP < 0|IM). Setting Z = C − EDP, with statistical parameter Z that represents the state function of structures which also obeys normal distribution. It is easy to calculate the average value and standard deviation of Z, which are µ z = µ c − µ d , σ z = σ 2 c + σ 2 d , separately. Therefore, the failure probability of the structure conditioned on ground motion parameter can be expressed as: For the convenience of calculation, the N(µ Z , σ Z ) should be transformed into a standard normal distribution N(0, 1). Setting γ = (Z − µ z )/σ Z , and then the new equation of the failure probability is: where the P f (DS i IM) is the probability of being equal to or exceeding a particular damage state for a given seismic intensity level; Φ(·) is the standard cumulative probability function; C m is the median threshold value of the capacity required to reach the particular damage state that can be obtained by looking DI m up in Table 1; the parameter σ c describes the uncertainty caused by the capacity of tunnels, and σ c = 0.4 based on HAZUS for building [4]; the parameter σ d is caused by the EDP, and σ d = 0.3 according to analyses for bored tunnels of BART system [29]. Some experts suggest that σ Z be replaced by the total variability σ tot to calculate the probability of a tunnel's failure, which is reasonable in the vulnerability assessment of a tunnel [4, 15,16]. The modified fragility function is: where σ tot describes the total variability associated with each fragility curve, and where σ IM is the standard deviation caused by the uncertainty of other factors related to vulnerability assessment, being usually determined by calculating natural logarithmic standard deviation of damage index (DI) for the different input motions at each level of PGA [4].

Specification of Tunnel Lining
This research takes a prototype of a subway tunnel in Shanghai whose essential parameters are known to conduct the seismic vulnerability assessment [30]. The tunnel with a double-hole of a 13 m center distance is considered. The shape of the tunnel lining section and its parameters are displayed in Figure 2. The section of the hole is circular, and its outer diameter is 6.2 m. The thickness of the tunnel lining is 0.35 m. The tunnel is made of concrete, and the material properties are as follows: Young's modulus E = 34.5 GPa, Poisson ratio v = 0.2. The buried depths chosen are 10 m, 15 m, 20 m, 35 m, respectively, to explore the variation of the tunnel lining vulnerability with buried depth. To express the key points of the lining clearly later, the section has been simplified, as shown in Figure 2.
Sustainability 2020, 12, x FOR PEER REVIEW 6 of 17 in Figure 2. The section of the hole is circular, and its outer diameter is 6.2 m. The thickness of the tunnel lining is 0.35 m. The tunnel is made of concrete, and the material properties are as follows: Young's modulus E = 34.5 GPa, Poisson ratio v = 0.2. The buried depths chosen are 10 m, 15 m, 20 m, 35 m, respectively, to explore the variation of the tunnel lining vulnerability with buried depth. To express the key points of the lining clearly later, the section has been simplified, as shown in Figure  2.

Soil Profile
The selected site is a deep soft soil site with rock formations at about 280 m depth [31]. The entire soil profile is mainly composed of clay and sand, which could be divided into 14 soil layers and 6 soil materials (see Table 2) in analysis. The shear wave velocity and density of soil layers are shown in Figure 3. The nonlinear characteristics of soil layers are considered by one-dimensional equivalent linearization method, to obtain equivalent shear modulus ratio G/Gmax and damping ratio ζ compatible with equivalent shear strain. Before applying the method, the soil dynamic parameters must be defined, collected from lots of experimental results that are the variations of shear modulus ratio G/Gmax and damping ratios ζ with shear strain, as listed in Table 2 [32]. The maximum shear

Soil Profile
The selected site is a deep soft soil site with rock formations at about 280 m depth [31]. The entire soil profile is mainly composed of clay and sand, which could be divided into 14 soil layers and 6 soil materials (see Table 2) in analysis. The shear wave velocity v s and density ρ of soil layers are shown in Figure 3.  in Figure 2. The section of the hole is circular, and its outer diameter is 6.

Soil Profile
The selected site is a deep soft soil site with rock formations at about 280 m depth [31]. The entire soil profile is mainly composed of clay and sand, which could be divided into 14 soil layers and 6 soil materials (see Table 2) in analysis. The shear wave velocity and density of soil layers are shown in Figure 3. The nonlinear characteristics of soil layers are considered by one-dimensional equivalent linearization method, to obtain equivalent shear modulus ratio G/Gmax and damping ratio ζ compatible with equivalent shear strain. Before applying the method, the soil dynamic parameters must be defined, collected from lots of experimental results that are the variations of shear modulus The nonlinear characteristics of soil layers are considered by one-dimensional equivalent linearization method, to obtain equivalent shear modulus ratio G/G max and damping ratio ζ compatible with equivalent shear strain. Before applying the method, the soil dynamic parameters must be defined, collected from lots of experimental results that are the variations of shear modulus ratio G/G max and damping ratios ζ with shear strain, as listed in Table 2 [32]. The maximum shear modulus G max of soil layers should be estimated by shear wave velocity v s and density ρ, namely the following equation: The equivalent shear modulus ratio G/G max and damping ratio ζ are calculated by EERA which is an iterative computational program based on SHAKE91, paired with input motions at each of the seismic intensity levels. In the iterative procedure, the ratio of the effective and maximum shear strain is assumed equal to 0.65. The results can be converted into elastic modulus and Rayleigh damping of the soil-tunnel model, and then the dynamic time history analysis can be carried out.
The free field response analysis is a key step in the analysis of soil-structure interaction for underground structures [33]. Figure 4 shows a typical example of the computed ground response in terms of equivalent shear modulus ratio G/G max , equivalent damping ratio ζ, peak acceleration a max and peak shear strain γ max . The computed a max value at the soil surface (PGA) is selected as the seismic intensity parameter in the fragility curves.
The equivalent shear modulus ratio G/Gmax and damping ratio ζ are calculated by EERA which is an iterative computational program based on SHAKE91, paired with input motions at each of the seismic intensity levels. In the iterative procedure, the ratio of the effective and maximum shear strain is assumed equal to 0.65. The results can be converted into elastic modulus and Rayleigh damping of the soil-tunnel model, and then the dynamic time history analysis can be carried out.
The free field response analysis is a key step in the analysis of soil-structure interaction for underground structures [33]. Figure 4 shows a typical example of the computed ground response in terms of equivalent shear modulus ratio G/Gmax, equivalent damping ratio ζ, peak acceleration and peak shear strain . The computed value at the soil surface (PGA) is selected as the seismic intensity parameter in the fragility curves.

Establishment of Soil-Tunnel Model and Boundary Condition
The finite element models (FEM) of the two-dimensional soil-tunnel system is generated using the software ABAQUS to calculate the actual internal forces (including the moment and the axis force) of the lining under the excitations of vertically propagating shear wave. The soil-tunnel model with a 2800 m width and a 280 m depth is displayed in Figure 5. Considering that the longitudinal length of the tunnel is much larger than the cross-section size, the soil and tunnel lining are simulated by

Establishment of Soil-Tunnel Model and Boundary Condition
The finite element models (FEM) of the two-dimensional soil-tunnel system is generated using the software ABAQUS to calculate the actual internal forces (including the moment and the axis force) of the lining under the excitations of vertically propagating shear wave. The soil-tunnel model with a 2800 m width and a 280 m depth is displayed in Figure 5. Considering that the longitudinal length of the tunnel is much larger than the cross-section size, the soil and tunnel lining are simulated by plane-strain element and beam element, respectively. The grid size satisfies the size control conditions proposed by Kuhlemyer et al. [34]; that is, the maximum side of each element is less than 1/8 wavelength corresponding to the maximum frequency. To obtain an accurate model, the model should be composed closely with the real behavior of the soil-tunnel interface. Therefore, the area around the tunnel is more boundaries could dissipate the arriving scattered waves. Building a sufficiently large soil domain (2800 m × 280 m) could dissipate the waves. Constraining the boundary nodes to move in pure shear is to simulate a free-field condition. The pure shear constraint can be achieved by ensuring the nodes at the lateral boundaries at the same elevation to move together horizontally [35,36], also called the periodic boundary. The periodic boundary is applied to the model and keeps compatible deformation.

Selection of Input Ground Motions
Selection of ground motions is significant for seismic fragility analysis; it is not only to ensure adequate quantity, but also to meet site characteristics. Although the strong motion observation network has been set up in Shanghai, the strong motion records have not been obtained so far. Thus, the required seismic waves can only be selected from the existing records. Fifteen ground motions with different intensities are selected from the Pacific Earthquake Engineering Research Center (PEER) by the method based on the station/earthquake information [37,38] listed in Table 3, and the response spectrum of the input motions is shown in Figure 6. The criteria for selecting seismic waves are as follows: 1. In reference to Chinese codes for seismic design of buildings, the seismic wave whose shear wave velocity Vs_30 m is less than 200 m/s is selected firstly. 2. To reduce the effect of near-field effects, the epicentral distance of the selected seismic waves must be more than 10 km. 3. To avoid the influence of the same focal mechanism, only one horizontal seismic wave is selected for each station and each seismic event. 4. To avoid the dependence of the selected seismic waves on the recording direction, only one of the larger PGA is selected for the two horizontal time history recording curves of the same seismic event.  The bottom boundary condition for the soil domain is fully fixed at the base. The lateral boundaries could dissipate the arriving scattered waves. Building a sufficiently large soil domain (2800 m × 280 m) could dissipate the waves. Constraining the boundary nodes to move in pure shear is to simulate a free-field condition. The pure shear constraint can be achieved by ensuring the nodes at the lateral boundaries at the same elevation to move together horizontally [35,36], also called the periodic boundary. The periodic boundary is applied to the model and keeps compatible deformation.

Selection of Input Ground Motions
Selection of ground motions is significant for seismic fragility analysis; it is not only to ensure adequate quantity, but also to meet site characteristics. Although the strong motion observation network has been set up in Shanghai, the strong motion records have not been obtained so far. Thus, the required seismic waves can only be selected from the existing records. Fifteen ground motions with different intensities are selected from the Pacific Earthquake Engineering Research Center (PEER) by the method based on the station/earthquake information [37,38] Table 3, and the response spectrum of the input motions is shown in Figure 6. The criteria for selecting seismic waves are as follows:

listed in
1.
In reference to Chinese codes for seismic design of buildings, the seismic wave whose shear wave velocity Vs_30 m is less than 200 m/s is selected firstly.

2.
To reduce the effect of near-field effects, the epicentral distance of the selected seismic waves must be more than 10 km.

3.
To avoid the influence of the same focal mechanism, only one horizontal seismic wave is selected for each station and each seismic event. 4.
To avoid the dependence of the selected seismic waves on the recording direction, only one of the larger PGA is selected for the two horizontal time history recording curves of the same seismic event.

5.
Fifteen different seismic waves are selected, and the seismic waves with larger magnitudes are prioritized in the wave selection process to eliminate the earthquake that is unlikely to affect structural safety. be performed by an inverse analysis to obtain the bedrock waves as the input. Meanwhile, to acquire the internal forces of the tunnel lining for gradually increasing level of seismic intensity, the peak of time histories is, in turn, scaled into 0.03 g, 0.075 g, 0.15 g, 0.3 g, 0.45 g, 0.6 g, 0.75 g, 0.9 g, 1.05 g and 1.2 g. The principle is that the frequency characteristics before and after amplitude modulation is consistent.

Internal Forces and Deformation of the Lining
The weak position can be estimated by investigating the response distribution along the circumferential direction of the lining structure under earthquake excitations. Meanwhile, the bending moment is an important index of the tunnel damage state in the fragility analysis. Linear The waves are input at the bottom of the model, so the selected surface ground motions need to be performed by an inverse analysis to obtain the bedrock waves as the input. Meanwhile, to acquire the internal forces of the tunnel lining for gradually increasing level of seismic intensity, the peak of time histories is, in turn, scaled into 0.03 g, 0.075 g, 0.15 g, 0.3 g, 0.45 g, 0.6 g, 0.75 g, 0.9 g, 1.05 g and 1.2 g. The principle is that the frequency characteristics before and after amplitude modulation is consistent.

Internal Forces and Deformation of the Lining
The weak position can be estimated by investigating the response distribution along the circumferential direction of the lining structure under earthquake excitations. Meanwhile, the bending moment M is an important index of the tunnel damage state in the fragility analysis. Linear elastic analysis was carried out by using ABAQUS software. The internal force distribution was studied with changing the buried depth of tunnels and the peak acceleration of ground motion. It was found that most of the results showed the similar trend. The variation trend of bending moment and axial force with different buried depths under the excitation of Northridge-01 wave (0.15 g and 0.60 g) is illustrated in Figure 7. The peak values of the bending moment M and the axial force N of the circular tunnel lining always appear near the four positions having a 45 • central angle with the connecting line of the dome (point A) and the arch bottom (point E), corresponding to θ = 45 • , θ = 135 • , θ = 225 • , and θ = 315 • . It also can be seen that the bending moment increases first and then decreases with the increase of burial depth, and that the peak occurs when the depth is 15 m. However, the axial force of the tunnel increases with the increase of the buried depth.
Sustainability 2020, 12, x FOR PEER REVIEW 10 of 17 elastic analysis was carried out by using ABAQUS software. The internal force distribution was studied with changing the buried depth of tunnels and the peak acceleration of ground motion. It was found that most of the results showed the similar trend. The variation trend of bending moment and axial force with different buried depths under the excitation of Northridge-01 wave (0.15 g and 0.60 g) is illustrated in Figure 7. The peak values of the bending moment and the axial force of the circular tunnel lining always appear near the four positions having a 45° central angle with the connecting line of the dome (point A) and the arch bottom (point E), corresponding to θ = 45°, θ = 135°, θ = 225°, and θ = 315°. It also can be seen that the bending moment increases first and then decreases with the increase of burial depth, and that the peak occurs when the depth is 15 m. However, the axial force of the tunnel increases with the increase of the buried depth.   PGA under the same wave, and the increasing trend starts being slow or no longer exists when PGA is greater than 0.9 g. Besides, the dispersion of axial force reaction is less than that of moment reaction.
The failure of the tunnel is mainly caused by the overlarge deformation based on the statistical analysis of seismic damage data of lining. To further explore the relationship between structural deformation and forces, the correlation of the relative displacement between point A and point E in Figure 2 and bending moment were fitted, as shown in Figure 9. There is a linear correlation between them according to Figure 9. So, selecting bending moment as the damage index is compatible with the use of the deformation. Therefore, it is feasible to take the moment as the damage index when the deformation index is difficult to determine in seismic fragility analysis.

Seismic Fragility Curves with Different Buried Depths
The correlation between DI and PGA must be determined before the fragility curve could be obtained. Through a series of dynamic time-history analysis of the soil-tunnel system, the maximum bending moment of the lining with four buried depths can be acquired. Thus, the DI (M/M RD ) can also be calculated. DI is a special structural engineering demand parameter (EDP), and PGA is selected as the intensity measure (IM) for the analysis. Therefore, DI and PGA satisfy the following formula: ln(DI) = α + βln(PGA).  analysis of seismic damage data of lining. To further explore the relationship between structural deformation and forces, the correlation of the relative displacement between point A and point E in Figure 2 and bending moment were fitted, as shown in Figure 9. There is a linear correlation between them according to Figure 9. So, selecting bending moment as the damage index is compatible with the use of the deformation. Therefore, it is feasible to take the moment as the damage index when the deformation index is difficult to determine in seismic fragility analysis.

Seismic Fragility Curves with Different Buried Depths
The correlation between DI and PGA must be determined before the fragility curve could be obtained. Through a series of dynamic time-history analysis of the soil-tunnel system, the maximum bending moment of the lining with four buried depths can be acquired. Thus, the DI ( / ) can also be calculated. DI is a special structural engineering demand parameter (EDP), and PGA is selected as the intensity measure (IM) for the analysis. Therefore, DI and PGA satisfy the following formula: The correlation between DI and PGA is determined by linear regression analysis, as shown in Figure 10. The goodness-of-fit (R 2 ) of all results are better than 0.9, which means the discreteness of results is relatively small, and that fitting results are satisfactory. This may be related to the wave selection method based on the station/earthquake information that will reduce the dispersion of calculation results according to Reference [37]. The correlation between DI and PGA is determined by linear regression analysis, as shown in Figure 10. The goodness-of-fit (R 2 ) of all results are better than 0.9, which means the discreteness of results is relatively small, and that fitting results are satisfactory. This may be related to the wave selection method based on the station/earthquake information that will reduce the dispersion of calculation results according to Reference [37]. The logarithmic standard deviation of damage index (DI) for the different input motions at each level of PGA can be calculated from the regression results, which means that the σ could be estimated, as illustrated in Table 4.  The logarithmic standard deviation of damage index (DI) for the different input motions at each level of PGA can be calculated from the regression results, which means that the σ IM could be estimated, as illustrated in Table 4. After acquiring the correlation between DI and PGA and the total uncertainty σ tot , they could be taken into the fragility Equation (9) to calculate the failure probability of the tunnel lining for a given seismic intensity level. Thus, the formulas of the tunnel fragility for the different buried depths are: where Equations (13)- (16) refers to the buried depth of 10 m, 15 m, 20 m and 35m, respectively.
The fragility curves of the tunnel are shown in Figure 11. It can be seen that, for the tunnels with the four depths in Shanghai soft soil, the tunnel with 15 m buried depth is most vulnerable to earthquake and 35 m buried depth is least vulnerable at a given PGA. Therefore, the failure probability of the tunnel is not monotonically decreasing with the increase of the buried depth. The maximum deformation of the soil layer appears in the 15 m depth case, shown in Figure 4. In general, the reason for this phenomenon is related to the deformation due to a small equivalent shear modulus and a large damping ratio. The greater the deformation of the soil layer, the greater the internal force generated.

Comparison Between Numerical and Empirical Fragility Curves
Empirical fragility curves could be obtained based on the history seismic damage data, which can reflect the actual damage of structures. However, empirical fragility curves do not always perform well for a city lacking strong seismic records. Meanwhile, the acquisition of the empirical fragility curves does not take into account specific soil characteristics, tunnel depth, tunnel material information, and so on. However, it is still instructive. The analytical fragility curves obtained in this study and the empirical fragility curves derived from ALA and HAZUS are compared in Figure 12 [2,4]. It can be seen that only the analytical fragility curves of the tunnel with 35 m depth are close to the empirical fragility curves of HAZUS. The empirical fragility curves are often statistically averaged. The shallow metro tunnel in Shanghai soft soil is easy to cause unexpected responses. Then, the curve of the tunnel with a larger buried depth is consistent with the empirical one. The intensity level corresponding to the failure probability of 50% is defined as the median PGA, and the median PGA is listed in Table 5. earthquake and 35 m buried depth is least vulnerable at a given PGA. Therefore, the failure probability of the tunnel is not monotonically decreasing with the increase of the buried depth. The maximum deformation of the soil layer appears in the 15 m depth case, shown in Figure 4. In general, the reason for this phenomenon is related to the deformation due to a small equivalent shear modulus and a large damping ratio. The greater the deformation of the soil layer, the greater the internal force generated. Figure 11. Fragility curves of the tunnel in soft soils with different buried depths (PGA: peak ground acceleration).

Comparison Between Numerical and Empirical Fragility Curves
Empirical fragility curves could be obtained based on the history seismic damage data, which can reflect the actual damage of structures. However, empirical fragility curves do not always perform well for a city lacking strong seismic records. Meanwhile, the acquisition of the empirical fragility curves does not take into account specific soil characteristics, tunnel depth, tunnel material information, and so on. However, it is still instructive. The analytical fragility curves obtained in this study and the empirical fragility curves derived from ALA and HAZUS are compared in Figure 12 [2,4]. It can be seen that only the analytical fragility curves of the tunnel with 35 m depth are close to the empirical fragility curves of HAZUS. The empirical fragility curves are often statistically averaged. The shallow metro tunnel in Shanghai soft soil is easy to cause unexpected responses. Then, the curve of the tunnel with a larger buried depth is consistent with the empirical one. The intensity level corresponding to the failure probability of 50% is defined as the median PGA, and the median PGA is listed in Table 5.

Conclusions
The empirical fragility curves always have the statistical significance, ignoring specific soil characteristics, buried depth and tunnel material. A numerical method of performing the fragility assessment of the metro tunnel which can effectively avoid the defects has been proposed in this paper. The benefit of the proposed method is that it can provide a basis for seismic fragility analysis

Conclusions
The empirical fragility curves always have the statistical significance, ignoring specific soil characteristics, buried depth and tunnel material. A numerical method of performing the fragility assessment of the metro tunnel which can effectively avoid the defects has been proposed in this paper. The benefit of the proposed method is that it can provide a basis for seismic fragility analysis of tunnels without strong seismic records. Two-dimensional finite element models of a typical metro tunnel with four buried depths were built by using ABAQUS. The nonlinear characteristics of soil profiles were considered by one-dimensional equivalent linear method. The input ground motions were selected by the method based on station/earthquake information. A series of dynamic time-history analysis were carried out.
The bending moment M is an important index of describing the tunnel damage state. According to the distribution of M, the maximum M of the lining always appears near four positions that are θ = 45 • , θ = 135 • , θ = 225 • and θ = 315 • . There is a linear correlation between M and tunnel relative deformation. Selecting bending moment as the damage index is compatible with the use of deformation. Therefore, it is valid to take the moment as the damage index when the deformation index is difficult to determine in seismic fragility analysis. The correlation between damage index and PGA is determined by linear regression. The goodness-of-fit results are satisfactory, which means the discreteness of results is relatively small. This may be related to the method of wave selection based on the station/earthquake information that will reduce the dispersion of calculation results.
Comparing the analytical fragility curves and the empirical curves derived from ALA and HAZUS indicates that the proposed method is reliable. Further analysis about the effect of buried depth on the fragility of the tunnel has been performed. The results show that the failure probability of the tunnel is not monotonically decreasing with the increase of the buried depth for a given PGA, and that the tunnel with 15 m buried depth among the four depths is most vulnerable to earthquake, which is related to the overlarge deformation of soil layers at 15 m depth due to a low equivalent shear modulus and a high equivalent damping ratio.