# Vertical Seismic Effect on the Seismic Fragility of Large-Space Underground Structures

Authors: Zhiming He; Qingjun Chen
Journal: Advances in Civil Engineering
 (2019)

---

## Abstract

The measured vertical peak ground acceleration was larger than the horizontal peak ground acceleration. It is essential to consider the vertical seismic effect in seismic fragility evaluation of large-space underground structures. In this research, an approach is presented to construct fragility curves of large-space underground structures considering the vertical seismic effect. In seismic capacity, the soil-underground structure pushover analysis method which considers the vertical seismic loading is used to obtain the capacity curve of central columns. The thresholds of performance levels are quantified through a load-drift backbone curve model. In seismic demand, it is evaluated through incremental dynamic analysis (IDA) method under the excitation of horizontal and vertical acceleration, and the soil-structure-interaction and ground motion characteristics are also considered. The IDA results are compared in terms of peak ground acceleration and peak ground velocity. To construct the fragility curves, the evolutions of performance index versus the increasing earthquake intensity are performed, considering related uncertainties. The result indicates that if we ignore the vertical seismic effect to the fragility assessment of large-space underground structures, the exceedance probabilities of damage of large-space underground structures will be underestimated, which will result in an unfavorable assessment result.


---

## Body

## 1. Introduction

Large-space underground structures such as subway stations, commercial streets, and parking lots are wildly used urban construction measures. These structures can suffer severe damage under strong ground shaking [1, 2]. Especially for shallow embedded structures in soft soil, their susceptibility to damage can be increased, due to ground strain and velocity along with acceleration increase when approaching the ground surface [3, 4]. In addition, for some early built large-space underground structures, they are facing a high potential risk of seismic damage due to without properly considering seismic design. Consequently, to ensure the seismic safety of large-space underground structures, especially in seismic prone areas, the seismic fragility assessment is introduced to large-space underground structures by different researchers.So far, the seismic fragility studies related to large-space underground structures are conducted in different perspectives. Li et al. [5] used the pushover analysis method applicable for underground structures to investigate the vulnerability of Daikai subway station. In this process, the seismic demand calculated from pushover analysis method cannot consider the soil-structure-interaction and the frequency characters of ground motion. Liu et al. [6] investigated the vulnerability of Daikai subway station based on the full dynamic numerical analysis. The performance levels thresholds of aboveground structures were directly applied to underground subway station. Due to the geometric feature and the character which underground structures are buried and restrained by surrounding soil [7, 8], their seismic behavior and performance are highly distinct. Then, Huh et al. [9] investigated the vulnerability of a shallow double story underground RC box structure based on the ground response acceleration method. The performance levels thresholds of structures were defined through the pushover analysis method. In addition, the studies investigated by Castaldo et al. [10, 11] have demonstrated that the existing underground structure can cause a not neglectable increase in the seismic vulnerability of the nearby aboveground R.C structure due to the deep excavation-induced foundations’ displacement.The present fragility studies of large-space underground structures only consider the horizontal seismic loading. However, there are sufficient evidences to prove that the vertical seismic effect cannot be neglected. For instance, the measured vertical peak ground acceleration was larger than the horizontal peak ground acceleration in the 1995 Kobe earthquake [12]. In the 1994 Northridge earthquake [13], the peak vertical acceleration value of near-field motion reached 1.19 g, and the ratio of vertical-to-horizontal peak ground acceleration exceeded 1.5, significantly greater than 2/3. In addition, it has been proved that the contribution of vertical earthquake component is one of the major factors to the failure of the structure [7, 8, 14–16]. Both Parra-Montesinos [7] and Iida et al. [12] verified that the high axial load induced by the vertical component of ground motion can increase the axial compression ratio and reduce the ductility of central column of large-space underground structures, which is one of the most important factors of collapsing. Study of Nakamura et al. [16] demonstrated that the vertical ground motion caused the compression-bending failure of central columns. An et al. [8] stated that ductility of central column under high compression ratio which was induced by vertical vibration can decrease obviously. Although the vertical component of ground motion was introduced to seismic analysis of large-space underground structures, the existing studies still have not provide the solutions to consider the influence of the vertical seismic effect quantitatively in seismic fragility analysis.Along these lines, this paper presents a seismic fragility analysis approach for large-space underground structures, which account for the vertical seismic effect. In seismic capacity, the soil-underground structure pushover analysis method which considers the vertical seismic loading is applied to obtain the performance index thresholds. In seismic demand, it is evaluated through incremental dynamic analysis method under the excitation of horizontal and vertical acceleration, the soil-structure-interaction and ground motion characteristics are also considered. To construct the fragility curves, the evolutions of performance index versus the increasing earthquake intensity are performed, considering related uncertainties.


## 2. Methodology

### 2.1. Overview of the Proposed Method for Deriving Fragility Curve

The proposed method is based on pushover analysis and incremental dynamic analysis of large-space underground structures, considering the vertical seismic effect, which is described in Figure1. As the performance levels thresholds of large-space underground structures are not yet documented, one of the crucial tasks in this study is to obtain the thresholds. The existing seismic damage and failure mechanism of typical large-space underground structures demonstrate that central columns are the weakest position of the structure [12, 17]. Therefore, the seismic capacity of the entire structure can be represented by the seismic capacity of the central column. In order to obtain the capacity curve of the central column which can be used to quantify the thresholds of performance levels, the soil-underground structure pushover analysis method [18] is adopted. Because the seismic response of the soil-underground structure system is mainly controlled by the fundamental mode, the inverted triangular distribution that decreases linearly with depth is used for the distribution of body force [18]. This distribution is easy to obtain without ground response analysis, so is more practical than the other types of force distributions. But it should be noted that this pushover analysis method only considering the horizontal seismic loading. To consider the vertical seismic effect, the uniform vertical acceleration distribution [19] is added to the soil-underground structure system as vertical body force distribution, as shown in Figure 2. The influence of the vertical component of ground motion is often introduced by the av/ah ratio (av is the peak vertical acceleration and ah is the peak horizontal acceleration). As many design codes use an average av/ah ratio of 2/3 [20], the inputted vertical acceleration is scaled to 2/3 to the inputted horizontal acceleration, and both seismic accelerations are monotonically increasing simultaneously. The pushover analysis procedure is presented as the following sequence of steps:(1)
Establish the soil-underground structure analysis model(2)
Gravity response analysis: perform the static response of the soil-underground structure analysis model according to the gravity load(3)
Pushover analysis: based on the gravity response analysis result, conduct the pushover analysis by monotonically increasing forces using the inverted triangular distribution and uniform vertical acceleration distribution until structure collapse(4)
Record data of each analytical step and obtain the capacity curves of the central columnFigure 1
Proposed procedure for deriving analytical fragility curves of large-space underground structures.Figure 2
The computational model of pushover analysis.To quantify the thresholds of performance levels based on the obtained capacity curves, definition of performance index and performance levels is needed. As for performance index, the univariate performance index has been utilized for long time. Recently, the bivariate performance index has been applied widely in case-study such as in the probabilistic analysis of excavation-induced damages to existing structures [10] and in fragility estimation of reinforce concrete buildings [21], due to the advantage which can identify two damage criteria of structures. Since the earthquake damage investigations have demonstrated that the seismic failure of large-space underground structure is mainly caused by the horizontal deformation capacity of central columns [15], we select the univariate performance index in this study. Although different performance indexes and associated parameters have been proposed regarding to the vulnerability analysis of buildings [22] and bridge [23, 24], there is lack of relevant information for large-space underground structures. Therefore, in the present method, drift ratio is selected as the performance index due to the significant correlation between drift ratio and structural damage [25]. Moreover, this performance index is compatible with the measurement of seismic response and damage analysis of large-space underground structure in the previous studies [7, 8]. For aboveground structures, SEAOC [26] has proposed various performance levels and damage states according to previous experience of damages in structures. For large-space underground structures, earthquake damage and experimental data of structures or components are insufficient to define the damage state accurately. Therefore, as a preliminary conclusion, the relationships between the performance levels and damage states can be extended to large-space underground structures. Four performance levels are taken into account due to ground shaking in this study, which refer to fully operational (FO), operational (OP), life safety (LS), and near collapse (NC). The damage description of each performance level is described in Table 1.Table 1
Division of performance levels for large-space underground structures.
Performance levels
Damage states
Description of performance levels
Fully operational
No damage
Structure components do not suffer damage. Structures can be used normally.
Operational
Repairable
Strength and stiffness of structure maintain the same before and after the earthquake. Structures can be ensured its function with minor repair.
Life safety
Irreparable
Structure components experience serious damage and stiffness degeneration but are not collapsed.
Near collapse
Severe
Structure components are destroyed, and their strength and stiffness degenerate greatly. Structure nearly collapses.To obtain the seismic demand of soil-underground structure interaction system, the incremental dynamic analysis method [27] is adopted. A series of dynamic response analyses are conducted under the combination of horizontal and vertical seismic acceleration, for increasing levels of seismic intensity. Then, a number of curves describing the parameterized response versus the ground motion intensity level are produced (i.e., IDA curves). Utilizing the IDA results, the parameters deriving for fragility curves can be obtained by performing the linear regression of the logarithm of demand measure against the logarithm of intensity measure, which will introduce in detail in Section 2.2. To construct the fragility curves, the evolution of performance index versus the increasing earthquake intensity are performed, considering related uncertainties.The proposed method can be applied to evaluate new fragility curves which consider the different features of large-space underground structures geometries, the site properties and input motion characteristics.


### 2.2. Approach for Deriving Fragility Curves

The seismic fragilities express the conditional probabilities which the seismic demand (D) of the structure reach or exceed the capacity (C) of the structure. This exceedance probability is conditioned on a selected intensity measure (IM), which represents the levels of ground shaking. The generic expression of this conditional probability is(1)fragility=PD>CIM.The fragility function given in equation (1) can be evaluated using a probability distribution for the demand conditioned on the IM (i.e., probabilistic seismic demand model (PSDM)) and convolving it through a distribution for the capacity. Cornell et al. [28] suggested the estimation of the median demand (D^) follows a power model:(2)D^=aIMb,where a and b are regression coefficients.The Incremental Dynamic Analysis can be used to establish the PSDM. As a number of curves depicting the demandD versus the ground motion intensity are produced, then, by conducting the linear regression of the logarithm of D against the logarithm of IM, the PSDM parameters (a and b) can be determined.Furthermore, the distribution of the demand median is often assumed following a two-parameter lognormal probability distribution. Therefore, after estimating the dispersion (βDIM) of the demand median, which is conditioned on the IM, the fragility can be written as(3)PD>CIM=1−ΦlnC^−lna⋅IMbβDIM2+βC2+βDS2+βM2,where Φ is the standard cumulative probability function, C^ is the median of structural capacity associated with the performance level. βC, βDS, and βM are the uncertainties in capacity, the definition of damage states, and modeling, respectively [29]. Considering the lack of more rigorous estimation, the values of uncertainty βDS and βC are set as 0.4 following HAZUS in building [30] and 0.3 based on engineering judgment [31]. The uncertainty related to modeling is assumed to 0.2 [32]. The parameter βDIM can be calculated using(4)βDIM=1n−2∑i=1nlnDi−lna⋅IMb2,where n is the number of (D,IM) data points.


## 2.1. Overview of the Proposed Method for Deriving Fragility Curve

The proposed method is based on pushover analysis and incremental dynamic analysis of large-space underground structures, considering the vertical seismic effect, which is described in Figure1. As the performance levels thresholds of large-space underground structures are not yet documented, one of the crucial tasks in this study is to obtain the thresholds. The existing seismic damage and failure mechanism of typical large-space underground structures demonstrate that central columns are the weakest position of the structure [12, 17]. Therefore, the seismic capacity of the entire structure can be represented by the seismic capacity of the central column. In order to obtain the capacity curve of the central column which can be used to quantify the thresholds of performance levels, the soil-underground structure pushover analysis method [18] is adopted. Because the seismic response of the soil-underground structure system is mainly controlled by the fundamental mode, the inverted triangular distribution that decreases linearly with depth is used for the distribution of body force [18]. This distribution is easy to obtain without ground response analysis, so is more practical than the other types of force distributions. But it should be noted that this pushover analysis method only considering the horizontal seismic loading. To consider the vertical seismic effect, the uniform vertical acceleration distribution [19] is added to the soil-underground structure system as vertical body force distribution, as shown in Figure 2. The influence of the vertical component of ground motion is often introduced by the av/ah ratio (av is the peak vertical acceleration and ah is the peak horizontal acceleration). As many design codes use an average av/ah ratio of 2/3 [20], the inputted vertical acceleration is scaled to 2/3 to the inputted horizontal acceleration, and both seismic accelerations are monotonically increasing simultaneously. The pushover analysis procedure is presented as the following sequence of steps:(1)
Establish the soil-underground structure analysis model(2)
Gravity response analysis: perform the static response of the soil-underground structure analysis model according to the gravity load(3)
Pushover analysis: based on the gravity response analysis result, conduct the pushover analysis by monotonically increasing forces using the inverted triangular distribution and uniform vertical acceleration distribution until structure collapse(4)
Record data of each analytical step and obtain the capacity curves of the central columnFigure 1
Proposed procedure for deriving analytical fragility curves of large-space underground structures.Figure 2
The computational model of pushover analysis.To quantify the thresholds of performance levels based on the obtained capacity curves, definition of performance index and performance levels is needed. As for performance index, the univariate performance index has been utilized for long time. Recently, the bivariate performance index has been applied widely in case-study such as in the probabilistic analysis of excavation-induced damages to existing structures [10] and in fragility estimation of reinforce concrete buildings [21], due to the advantage which can identify two damage criteria of structures. Since the earthquake damage investigations have demonstrated that the seismic failure of large-space underground structure is mainly caused by the horizontal deformation capacity of central columns [15], we select the univariate performance index in this study. Although different performance indexes and associated parameters have been proposed regarding to the vulnerability analysis of buildings [22] and bridge [23, 24], there is lack of relevant information for large-space underground structures. Therefore, in the present method, drift ratio is selected as the performance index due to the significant correlation between drift ratio and structural damage [25]. Moreover, this performance index is compatible with the measurement of seismic response and damage analysis of large-space underground structure in the previous studies [7, 8]. For aboveground structures, SEAOC [26] has proposed various performance levels and damage states according to previous experience of damages in structures. For large-space underground structures, earthquake damage and experimental data of structures or components are insufficient to define the damage state accurately. Therefore, as a preliminary conclusion, the relationships between the performance levels and damage states can be extended to large-space underground structures. Four performance levels are taken into account due to ground shaking in this study, which refer to fully operational (FO), operational (OP), life safety (LS), and near collapse (NC). The damage description of each performance level is described in Table 1.Table 1
Division of performance levels for large-space underground structures.
Performance levels
Damage states
Description of performance levels
Fully operational
No damage
Structure components do not suffer damage. Structures can be used normally.
Operational
Repairable
Strength and stiffness of structure maintain the same before and after the earthquake. Structures can be ensured its function with minor repair.
Life safety
Irreparable
Structure components experience serious damage and stiffness degeneration but are not collapsed.
Near collapse
Severe
Structure components are destroyed, and their strength and stiffness degenerate greatly. Structure nearly collapses.To obtain the seismic demand of soil-underground structure interaction system, the incremental dynamic analysis method [27] is adopted. A series of dynamic response analyses are conducted under the combination of horizontal and vertical seismic acceleration, for increasing levels of seismic intensity. Then, a number of curves describing the parameterized response versus the ground motion intensity level are produced (i.e., IDA curves). Utilizing the IDA results, the parameters deriving for fragility curves can be obtained by performing the linear regression of the logarithm of demand measure against the logarithm of intensity measure, which will introduce in detail in Section 2.2. To construct the fragility curves, the evolution of performance index versus the increasing earthquake intensity are performed, considering related uncertainties.The proposed method can be applied to evaluate new fragility curves which consider the different features of large-space underground structures geometries, the site properties and input motion characteristics.


## 2.2. Approach for Deriving Fragility Curves

The seismic fragilities express the conditional probabilities which the seismic demand (D) of the structure reach or exceed the capacity (C) of the structure. This exceedance probability is conditioned on a selected intensity measure (IM), which represents the levels of ground shaking. The generic expression of this conditional probability is(1)fragility=PD>CIM.The fragility function given in equation (1) can be evaluated using a probability distribution for the demand conditioned on the IM (i.e., probabilistic seismic demand model (PSDM)) and convolving it through a distribution for the capacity. Cornell et al. [28] suggested the estimation of the median demand (D^) follows a power model:(2)D^=aIMb,where a and b are regression coefficients.The Incremental Dynamic Analysis can be used to establish the PSDM. As a number of curves depicting the demandD versus the ground motion intensity are produced, then, by conducting the linear regression of the logarithm of D against the logarithm of IM, the PSDM parameters (a and b) can be determined.Furthermore, the distribution of the demand median is often assumed following a two-parameter lognormal probability distribution. Therefore, after estimating the dispersion (βDIM) of the demand median, which is conditioned on the IM, the fragility can be written as(3)PD>CIM=1−ΦlnC^−lna⋅IMbβDIM2+βC2+βDS2+βM2,where Φ is the standard cumulative probability function, C^ is the median of structural capacity associated with the performance level. βC, βDS, and βM are the uncertainties in capacity, the definition of damage states, and modeling, respectively [29]. Considering the lack of more rigorous estimation, the values of uncertainty βDS and βC are set as 0.4 following HAZUS in building [30] and 0.3 based on engineering judgment [31]. The uncertainty related to modeling is assumed to 0.2 [32]. The parameter βDIM can be calculated using(4)βDIM=1n−2∑i=1nlnDi−lna⋅IMb2,where n is the number of (D,IM) data points.


## 3. Case Study

### 3.1. Numerical Models

The underground commercial street structure which was designed and built in China is taken as the research background. Its cross-sectional dimensions are depicted in Figure3. The reinforcement ratios of structural components are 2% (columns), 1% (wall), and 0.8% (ceiling and bottom slab). Sites I–IV are selected according to practical engineering condition, and, respectively, correspond to stiff site (site I), medium stiff site (site II), medium soft site (site III), and soft site (site IV), based on the code for seismic design of building in China [33]. The shear wave velocity of each site is listed in Figure 4.Figure 3
Cross section of the selected large-space underground structure.Figure 4
Shear velocity. (a) Site I. (b) Site II. (c) Site III. (d) Site IV.
(a)
(b)
(c)
(d)To obtain the seismic capacity and seismic demand of large-space underground structures, respectively, a detailed soil-underground structure pushover analysis model and soil-underground structure interaction analysis model are developed using the finite element code ABAQUS. Studies have demonstrated that the results obtained through finite element method are always affected by the epistemic uncertainties (i.e., the definition of structural model) due to the analytical model become increasingly simple or increasingly complex [34–36]. Therefore, to consider the epistemic uncertainties during the finite element analysis, different methods have been proposed under the framework of Bayesian approach [34–36]. In these methods, the experimental data are one of the importance factors in evaluating the posterior distribution of the uncertainty. As for large-space underground structures, the epistemic uncertainty is not considered in the present analytical model due to the limited experimental data. Depths of the analytical domain are 22.65 m (site I), 48 m (site II), 60 m (site III), and 70 m (site IV) respectively, which is the thickness of the assumed soil deposit. In dynamic analysis, the width of the analytical model mainly depends on two factors: the distance between the underground structure and free-field zone, and the reflection of the lateral boundary. Because lateral boundary does not absorb energy, the distance between underground structure and tied boundary must be sufficiently large so that the reflected vibration will not affect the dynamic behavior of underground structures. A large domain of 226.5 m (site I), 480 m (site II), 600 m (site III), and 700 m (site IV) are used in this analysis, which has been validated as sufficiently large. The calculation domain and element mesh of the models are depicted in Figures 5 and 6.Figure 5
Soil-underground structure pushover analysis model with the additional free field.Figure 6
Soil-underground structure interaction analysis model.The inelastic response characteristics of the surrounding soil are simulated through an elastoplastic model [37]. The bone curve of this model is established using the isotropic hardening law and characterized by the Mohr–Coulomb yield criterion [38]. The Mohr–Coulomb yield criterion assumes that yield happens when the shear stress in a material achieves a value that depends linearly on the normal stress in the same plane. The corresponding mechanical properties (i.e., the shear velocity (Vs), the unit weight (γ), the Poisson ratio (υ), the cohesion (c), and the friction angle (φ)) of different soils in each site are listed in Tables 2–5. During the analyses, the elasticity of soil can be calculated through E=G×2×1+υ, while G=ρ×Vs2. The plasticity of soil is characterized by the cohesion and the friction angle. The plastic-damage model [39] is used for underground structure; it assumes that the uniaxial tensile and compressive behavior of concrete are characterized through damage plasticity. During the mathematic calculation, the crucial parameters of the plastic-damage model (Tables 6 and 7) are adopted to describe the material properties and deformation behavior of concrete. Moreover, the idealized elastoplastic model is selected for the steel bars, with the density of 7800 kg/m3 and yield stress of 300 MPa.Table 2
Soil physical properties for site I.
Soil layer
Thickness (m)
Soil type
Unit weight (kg·m−3)
Shear wave velocity (m·s−1)
Poisson ratio
Cohesion (kPa)
Friction angle (°)
1
1.75
Miscellaneous fill
1750
142
0.40
10.0
23.8
2
1.90
Silty clay
1900
348
0.35
13.4
12.0
3
9.00
Strong weathered granite
2300
695
0.35
—
32
4
6.00
Medium weathered granite
2450
1246
0.25
—
32
5
4.00
Slightly weathered rock
2500
1578
0.22
—
32Table 3
Soil physical properties for site II.
Soil layer
Thickness (m)
Soil type
Unit weight (kg·m−3)
Shear wave velocity (m·s−1)
Poisson ratio
Cohesion (kPa)
Friction angle (°)
1
3.00
Fill
2000
152
0.30
10.0
23.8
2
4.00
Silty clay
2000
170
0.30
13.4
12.0
3
9.00
Clay
2000
258
0.30
10.0
25.0
4
6.00
Silt
2000
300
0.30
10.0
25.0
5
4.00
Silt
2000
320
0.30
0.0
32.0
6
3.00
Clay
2000
400
0.30
10.0
25.0
7
2.00
Round gravel
2000
350
0.30
0.0
35.0
8
8.00
Medium-fine sand
2000
420
0.30
0.0
34.0
9
9.00
Clay
2000
400
0.30
10.0
25.0
10
18
Rock
2000
550
0.30
—
32Table 4
Soil physical properties for site III.
Soil layer
Thickness (m)
Soil type
Unit weight (kg·m−3)
Shear wave velocity (m·s−1)
Poisson ratio
Cohesion (kPa)
Friction angle (°)
1
2.00
Clay
1940
173
0.33
10.0
25.0
2
3.00
Silty clay
1980
200
0.32
13.4
12.0
3
13.30
Silty clay
1990
212
0.32
10.0
23.8
4
7.10
Silty clay
2000
244
0.40
13.4
12.0
5
19.60
Sandy clay
2000
273
0.30
10.0
25.0
6
15.00
Silty clay
2050
333
0.26
10.0
25.0Table 5
Soil physical properties for site IV.
Soil layer
Thickness (m)
Soil type
Unit weight (kg·m−3)
Shear wave velocity (m·s−1)
Poisson ratio
Cohesion (kPa)
Friction angle (°)
1
1.00
Fill
1890
74
0.40
10.0
23.8
2
5.40
Silty clay
1850
87
0.35
13.4
12.0
3
11.26
Silty clay
1830
110
0.38
10.0
23.8
4
26.40
Sandy clay
1820
220
0.35
10.0
25.0
5
7.40
Clay
2040
195
0.35
10.0
25.0
6
18.54
Fine sand
1935
225
0.30
0.0
32.0Table 6
Material properties of concrete.
Parameters
Value
Density
ρ
kg
⋅
m
−
3
2400
Dilation angle
ψ
°
30
Elastic modulus
E
MPa
32500
Poisson’s ratio
υ
0.2
Initial compressive yield stress
σ
c
0
MPa
9.55
Limited compressive yield stress
σ
cu
MPa
26.8
Initial tensile yield stress
σ
t
0
MPa
2.39
Tensile stiffness recovery parameter
ω
t
0
Compressive stiffness recovery parameter
ω
c
1
Damage variable
(Table3)Table 7
Stress and damage factor versus plastic strain of concrete.
Plastic strain
Stress (MPa)
Damage variable
Compression
0.0000E + 00
18.026
0.00000
2.9233E − 05
21.576
0.01439
1.2438E − 04
29.746
0.04864
2.0932E − 04
32.617
0.07492
3.2177E − 04
34.594
0.10661
6.2346E − 04
36.052
0.18124
8.3638E − 04
35.380
0.22836
1.0796E − 03
33.723
0.27834
1.3385E − 03
31.557
0.32727
1.8656E − 03
26.920
0.41454
2.6240E − 03
21.018
0.51496
3.7809E − 03
14.657
0.62447
6.2920E − 03
9.087
0.73452
9.2686E − 03
5.002
0.83782
Tensile
0.0000E + 00
1.608
0.00000
1.3586E − 05
3.131
0.03406
2.3198E − 05
3.215
0.05538
4.0877E − 05
3.037
0.09857
6.2300E − 05
2.738
0.15603
1.0436E − 04
2.165
0.28146
1.7538E − 04
1.444
0.49663
2.6383E − 04
0.951
0.69271
3.4435E − 04
0.716
0.79632
4.7258E − 04
0.517
0.88134
5.9794E − 04
0.412
0.92187
7.2209E − 04
0.346
0.94431
9.0721E − 04
0.283
0.96303
1.0917E − 03
0.242
0.97345The interaction between underground structure and the surrounding soil is related to their interface properties. In the normal direction of the interface, we define a hard contact. The normal contact compressive stress can mutually transfer via the contact constraint. The element nodes on the interface satisfy Hooke’s Law and the Harmonized Condition of Displacement. In the tangential direction of the interface, tangential contact shear stress is also transferred, and we assume that tangential mechanics behavior of the interface follows the Coulomb friction law with a friction coefficient,µ, of 0.4, which corresponds to a frictional angle of 22° for the soil-underground structure interface. This contact relationship was commonly used by different researchers [14, 40, 41].In pushover analysis, to consider the gravity-load effect, an additional free-field model is established. The mixed lateral boundary condition is used for the lateral sides of the model (i.e., the displacement boundary condition is fixed in the vertical direction and a forced boundary condition is adopted in the horizontal direction). The bottom boundary is fixed in both directions. Therefore, the vertical displacements and horizontal forces of the lateral boundary are obtained through self-gravity response analysis of the additional free-field model and then subjected to the corresponding lateral boundary of the soil-underground structure model. The dynamic analyses are conducted in two steps. First step is a static step, which the geostatic stresses are introduced, considering the underground structure in place. During this step, the base of the analysis model is fixed in the horizontal and vertical directions. The subsequent step is a dynamic analyses step, the horizontal and vertical seismic acceleration is applied uniformly along the base of the model, and the bottom and lateral boundary are free in the horizontal and vertical direction.The viscous damping of the soil-underground structure system is modeled in the form of Rayleigh type. Because this damping type is a linear combination of mass as well as stiffness matrix, it can incorporate into the analysis procedure efficiently. The material damping for soil was selected in different values during dynamic response analysis by different researchers [42–44], mainly ranging from 0.05 to 0.1. In [43], the modeling issues of Rayleigh damping matrix in soil layers are discussed which demonstrate the damping value of 0.08 for soil can better represent a realistic value of soil damping during strong ground motions. Therefore, we utilize a 0.08 damping in the dynamic analysis for soil. Similarly, a widely used 0.05 damping is selected for underground structures.The soil and structure are discretized with 8-node reduced integration linear solid element (C3D8R), while the reinforcement is modeled with three-node linear beam elements (B31). In dynamic analysis, the element mesh for soil is based on Liao’s study [45], which can ensure the efficient reproduction of all the waveforms of whole frequency range under study, and the maximum height of element hmax in soil is determined as(5)hmax=1/75−1/160Vsfmax,where Vs is the shear and compression motion velocity, which can be determined through G0=ρVs2, ρ is the density of soil, and fmax is the maximum vibration frequency of the inputted motion. Therefore, the maximum heights of the element are 1 m to 3 m from the surface to bottom. A finer discretization is adopted near the underground structure.


### 3.2. Selection of Input Motions

The seismic IMs can be expressed in terms of PGA, PGV, orSa (T1, 5%) (i.e., 5% damped first-mode spectral acceleration). As for underground tunnel structures, American Lifelines Alliance [46] produced empirical fragility curves by PGA, while Corigliano [47] proposed empirical fragility curves by PGV. Due to no consensus on the IMs required for seismic analysis of underground structures, PGA and PGV are adopted as IMs to investigate the effects of each parameter in the scattering of the engineering demand parameter (EDP) values. A total of 15 strong motion records are selected from strong motion database of the United States (PEER), Japan (NRIESDP), and China and Taiwan (SMART-I) for each site, respectively, exhibiting different spectral acceleration amplitudes, frequency content, significant duration, and seismotectonic environment. This is considered sufficient to capture the randomness of ground motion and provide sufficient accuracy in estimating seismic demand [48]. These records are all selected from different earthquakes, with source distances larger than 10 km, relatively large interpolated earthquakes (MW ≥ 6), and peak ground acceleration of 0.1 g or more. The horizontal and vertical acceleration response spectra (damping ratio = 0.08) of the selected ground motions are plotted in Figures 7 and 8.Figure 7
Acceleration response spectra of the horizontal component. (a) Site I. (b) Site II. (c) Site III. (d) Site IV.
(a)
(b)
(c)
(d)Figure 8
Acceleration response spectra of the vertical component. (a) Site I. (b) Site II. (c) Site III. (d) Site IV.
(a)
(b)
(c)
(d)The selected motion records are ground motions. Therefore, the horizontal and vertical components are inverted to bedrock to obtain the input motions. According to the established empirical curves which decrease in shear modulus and increase in material damping with the increasing shear strain amplitude for each soil layer, these processes were accomplished through the shake91 code based on the one-dimensional equivalent-linear wave-propagation analyses. For a gradually increasing level of seismic intensity, the horizontal earthquake components are scaled from lower to higher levels of peak acceleration on bedrock (i.e., 0.035 g, 0.1 g, 0.2 g, 0.4 g, 0.6 g, 0.8 g, 1.0 g, and 1.2 g). Then, a series of nonlinear time history analyses of soil-underground structure interaction analysis model are conducted.


### 3.3. Pushover Analysis Results

The capacity curves of the central column considering the vertical seismic effect are obtained (Figure9). To investigate the influence of the vertical seismic effect, we also present the capacity curves without considering the vertical seismic effect (Figure 10). The characteristic points (i.e., θy (the yield drift ratio), θp (the drift ratio of peak loading), and θu (the ultimate drift ratio)) are shown within the figures. We can know the values of θy and θp under considering the vertical seismic loading are almost the same to the values without considering the vertical seismic effect. That is mainly because during the initial loading period, the subjected acceleration is small. Therefore, the vertical seismic effect is not significant compared to self-gravity. With both horizontal and vertical seismic acceleration increasing, the structure gradually enters to the nonlinear stage. Under the vertical seismic effect, the capacity curve decreasing rate is faster compared to only consider the horizontal seismic loading. Therefore, we can find the value of θu under the vertical seismic effect is nearly 50% less compared to without considering the vertical seismic effect. It demonstrates the vertical seismic effect can seriously reduce the ductility of the central column. This finding is in accordance with the previous studies [7, 8, 14–16]. Thus, the seismic capacity of the central column will be overestimated if the vertical seismic effect is not considered.Figure 9
Capacity curves considering the vertical seismic effect. (a) Site I. (b) Site II. (c) Site III. (d) Site IV.
(a)
(b)
(c)
(d)Figure 10
Capacity curves without considering the vertical seismic effect. (a) Site I. (b) Site II. (c) Site III. (d) Site IV.
(a)
(b)
(c)
(d)


### 3.4. Thresholds of Performance Levels

In order to determine the thresholds of different performance levels based on pushover analysis results, in this study, lateral load-drift curve (Figure11) is adopted as the backbone curve model based on the experimental study [49]. This model consists of four performance levels for the component, which are in accordance with the definition of performance levels in this study. It includes four characteristic points: origin point (A), yield point (B), peak strength point (C), and ultimate point (D). The associated drift ratios are θy, θp, and θu. The ultimate point is defined as the force dropped to 85% of the peak force. The values of θOP and θLS are defined as(6)θOP=θy+0.5θp−θy,θLS=θp+0.5θu−θp.Figure 11
Backbone load-drift quantization model.The average values of characteristic points of the capacity curves under different site are used to quantify the thresholds of different performance levels. The quantified thresholds are listed in Table8. For FO and OP performance level, the thresholds(H+V) (i.e., thresholds considering the vertical seismic effect) are the same to thresholds(H) (i.e., thresholds without considering the vertical seismic effect). However, for LS and NC performance level, thresholds(H+V) are obviously less than thresholds(H). This means that the underground structure is easier to reach these two performance levels for the cases considering the vertical seismic effect. If we ignore the vertical seismic effect to the thresholds of performance levels, the seismic capacity of underground structure can be overestimated.Table 8
Thresholds of performance levels.
Performance levels (10−3 rad)
FO
OP
LS
NC
Thresholds(H+V)
1.5
2.2
5.8
8.6
Thresholds(H)
1.5
2.2
7.9
12.8


### 3.5. Incremental Dynamic Analysis Results

In incremental dynamic analysis, the failure condition of the structure has been considered recently within the numerical analysis by the mean defining the collapse condition of the structure [50]. The failure condition can be reflected in the seismic fragility assessment through the mean of the total probability theorem. As for large-space underground structures, they almost reach completely failure condition during earthquake due to the surrounding soil. Therefore, the failure condition is not considered in the present incremental dynamic analysis. The EDP (i.e., drift ratio) and IMs (i.e., PGA and PGV) are set as horizontal and vertical axis, respectively. Based on the spline interpolation between discrete points, then IDA curves of different seismic motions under different sites are obtained (Figures 12 and 13). Due to various values of the material properties of soil such as modulus elasticity, Poisson’s ratio, cohesion, and friction angle, as well as the different ground motion characteristics, each IDA curve shows different developing tendency. Unlike aboveground structures, the IDA curves do not exhibit a noticeable plateau after a certain seismic intensity. This phenomenon, which is due to the soil, prevents structures from reaching the nonlinear stage. Moreover, the restraint of the surrounding soil results in the underground structure not reflecting its self-vibration characteristics. Consequently, the dynamic response between underground structure and aboveground structures exhibits clear differences.Figure 12
IDA curves in terms of PGA as IM. (a) Site I. (b) Site II. (c) Site III. (d) Site IV.
(a)
(b)
(c)
(d)Figure 13
IDA curves in terms of PGV as IM. (a) Site I. (b) Site II. (c) Site III. (d) Site IV.
(a)
(b)
(c)
(d)The dispersion of the IDA curves is mainly due to different input motions (i.e., same PGA, but different ground motion records). To quantify this dispersion, the average of the In (EDP)’s standard deviation is computed and compared. The smaller value of In (EDP)’s standard deviation represents lower dispersion so the associated IM is more appropriate for IDA curves. Table9 listed the average of In (EDP)’s standard deviation of the drift ratio with respective PGV and PGA as the IM of different sites. It is found that PGA is more suitable than PGV as IM in stiff site and has a relatively large scatter in soft site, while PGV shows less scatter as IM in soft site. In order to yield the regression coefficient in equation (3), we select PGV as the representative parameter of seismic IM in the following fragility curves because its mean In (EDP)’s standard deviation is smaller than PGA.Table 9
Average ofσlnEDP.
IM/site
Site I
Site II
Site III
Site IV
PGA
0.403
0.436
0.451
0.503
PGV
0.444
0.441
0.432
0.469


### 3.6. Fragility Analysis

The derivations of fragility curves are based on the constructions of diagrams of the obtained performance index thresholds (Table8) versus PGV at the ground surface. The diagrams are estimated through linear regression analyses. In these analyses, the natural logarithm of the performance index is set as the dependent variable, and the natural logarithm of PGV is set as the independent variable (Figure 14).Figure 14
Linear regression analyses of IDA results. (a) Site I. (b) Site II. (c) Site III. (d) Site IV.
(a)
(b)
(c)
(d)The sets of fragility curves of large-space underground structures derived for each site are given in Figure15. The comparison reveals that the vulnerability of the large-space underground structure is increased when it is embedded in softer site conditions, i.e., from site I to site IV. Therefore, in this context, the important effect of sites conditions on the vulnerability of large-space underground structure is highlighted. It is observed that the probability of the NC performance level is very low in site I, site II, and site III. These phenomena are in accordance with the obtained observations from recent strong earthquakes. The site which Daikai subway station located is similar to site II in this study. However, the obtained fragility curves of large-space underground structure embedded in site II are less satisfactory in the case of Daikai subway station. This difference may be caused by the difference of the design and construction standard of the large-space underground structure. The current obtained fragility curves represent the underground structures designed and constructed based on the improved modern seismic code. However, Daikai station is designed without considering the modern seismic codes. In addition, these fragility curves are developed for square section, which is not the same as the vulnerable central columns in Daikai subway station.Figure 15
Comparison of fragility curves of different sites under different performance levels: (a) FO, (b) OP, (c) LS, and (d) NC.
(a)
(b)
(c)
(d)In Figure16, the fragility curves yielding from thresholds(H+V) and thresholds(H) are compared. As the thresholds(H+V) and thresholds(H) are the same in FO and OP performance levels, we only compare the fragility curves derived from LS and NC performance levels. The exceedance probability of damage in the case of considering the vertical seismic effect is increased from 1% to 13% for LS performance level in both site conditions. For the CP performance level, this increase ranges between 1% and 15%. This demonstrates, when considering the vertical earthquake motion in seismic fragility analyses of large-space underground structures, the exceedance probabilities will be underestimated if thresholds(H) are adopted, resulting in an unfavorable assessment result.Figure 16
Comparison of fragility curves obtained from thresholds(H+V) and thresholds(H). (a) Site I. (b) Site II. (c) Site III. (d) Site IV.
(a)
(b)
(c)
(d)


## 3.1. Numerical Models

The underground commercial street structure which was designed and built in China is taken as the research background. Its cross-sectional dimensions are depicted in Figure3. The reinforcement ratios of structural components are 2% (columns), 1% (wall), and 0.8% (ceiling and bottom slab). Sites I–IV are selected according to practical engineering condition, and, respectively, correspond to stiff site (site I), medium stiff site (site II), medium soft site (site III), and soft site (site IV), based on the code for seismic design of building in China [33]. The shear wave velocity of each site is listed in Figure 4.Figure 3
Cross section of the selected large-space underground structure.Figure 4
Shear velocity. (a) Site I. (b) Site II. (c) Site III. (d) Site IV.
(a)
(b)
(c)
(d)To obtain the seismic capacity and seismic demand of large-space underground structures, respectively, a detailed soil-underground structure pushover analysis model and soil-underground structure interaction analysis model are developed using the finite element code ABAQUS. Studies have demonstrated that the results obtained through finite element method are always affected by the epistemic uncertainties (i.e., the definition of structural model) due to the analytical model become increasingly simple or increasingly complex [34–36]. Therefore, to consider the epistemic uncertainties during the finite element analysis, different methods have been proposed under the framework of Bayesian approach [34–36]. In these methods, the experimental data are one of the importance factors in evaluating the posterior distribution of the uncertainty. As for large-space underground structures, the epistemic uncertainty is not considered in the present analytical model due to the limited experimental data. Depths of the analytical domain are 22.65 m (site I), 48 m (site II), 60 m (site III), and 70 m (site IV) respectively, which is the thickness of the assumed soil deposit. In dynamic analysis, the width of the analytical model mainly depends on two factors: the distance between the underground structure and free-field zone, and the reflection of the lateral boundary. Because lateral boundary does not absorb energy, the distance between underground structure and tied boundary must be sufficiently large so that the reflected vibration will not affect the dynamic behavior of underground structures. A large domain of 226.5 m (site I), 480 m (site II), 600 m (site III), and 700 m (site IV) are used in this analysis, which has been validated as sufficiently large. The calculation domain and element mesh of the models are depicted in Figures 5 and 6.Figure 5
Soil-underground structure pushover analysis model with the additional free field.Figure 6
Soil-underground structure interaction analysis model.The inelastic response characteristics of the surrounding soil are simulated through an elastoplastic model [37]. The bone curve of this model is established using the isotropic hardening law and characterized by the Mohr–Coulomb yield criterion [38]. The Mohr–Coulomb yield criterion assumes that yield happens when the shear stress in a material achieves a value that depends linearly on the normal stress in the same plane. The corresponding mechanical properties (i.e., the shear velocity (Vs), the unit weight (γ), the Poisson ratio (υ), the cohesion (c), and the friction angle (φ)) of different soils in each site are listed in Tables 2–5. During the analyses, the elasticity of soil can be calculated through E=G×2×1+υ, while G=ρ×Vs2. The plasticity of soil is characterized by the cohesion and the friction angle. The plastic-damage model [39] is used for underground structure; it assumes that the uniaxial tensile and compressive behavior of concrete are characterized through damage plasticity. During the mathematic calculation, the crucial parameters of the plastic-damage model (Tables 6 and 7) are adopted to describe the material properties and deformation behavior of concrete. Moreover, the idealized elastoplastic model is selected for the steel bars, with the density of 7800 kg/m3 and yield stress of 300 MPa.Table 2
Soil physical properties for site I.
Soil layer
Thickness (m)
Soil type
Unit weight (kg·m−3)
Shear wave velocity (m·s−1)
Poisson ratio
Cohesion (kPa)
Friction angle (°)
1
1.75
Miscellaneous fill
1750
142
0.40
10.0
23.8
2
1.90
Silty clay
1900
348
0.35
13.4
12.0
3
9.00
Strong weathered granite
2300
695
0.35
—
32
4
6.00
Medium weathered granite
2450
1246
0.25
—
32
5
4.00
Slightly weathered rock
2500
1578
0.22
—
32Table 3
Soil physical properties for site II.
Soil layer
Thickness (m)
Soil type
Unit weight (kg·m−3)
Shear wave velocity (m·s−1)
Poisson ratio
Cohesion (kPa)
Friction angle (°)
1
3.00
Fill
2000
152
0.30
10.0
23.8
2
4.00
Silty clay
2000
170
0.30
13.4
12.0
3
9.00
Clay
2000
258
0.30
10.0
25.0
4
6.00
Silt
2000
300
0.30
10.0
25.0
5
4.00
Silt
2000
320
0.30
0.0
32.0
6
3.00
Clay
2000
400
0.30
10.0
25.0
7
2.00
Round gravel
2000
350
0.30
0.0
35.0
8
8.00
Medium-fine sand
2000
420
0.30
0.0
34.0
9
9.00
Clay
2000
400
0.30
10.0
25.0
10
18
Rock
2000
550
0.30
—
32Table 4
Soil physical properties for site III.
Soil layer
Thickness (m)
Soil type
Unit weight (kg·m−3)
Shear wave velocity (m·s−1)
Poisson ratio
Cohesion (kPa)
Friction angle (°)
1
2.00
Clay
1940
173
0.33
10.0
25.0
2
3.00
Silty clay
1980
200
0.32
13.4
12.0
3
13.30
Silty clay
1990
212
0.32
10.0
23.8
4
7.10
Silty clay
2000
244
0.40
13.4
12.0
5
19.60
Sandy clay
2000
273
0.30
10.0
25.0
6
15.00
Silty clay
2050
333
0.26
10.0
25.0Table 5
Soil physical properties for site IV.
Soil layer
Thickness (m)
Soil type
Unit weight (kg·m−3)
Shear wave velocity (m·s−1)
Poisson ratio
Cohesion (kPa)
Friction angle (°)
1
1.00
Fill
1890
74
0.40
10.0
23.8
2
5.40
Silty clay
1850
87
0.35
13.4
12.0
3
11.26
Silty clay
1830
110
0.38
10.0
23.8
4
26.40
Sandy clay
1820
220
0.35
10.0
25.0
5
7.40
Clay
2040
195
0.35
10.0
25.0
6
18.54
Fine sand
1935
225
0.30
0.0
32.0Table 6
Material properties of concrete.
Parameters
Value
Density
ρ
kg
⋅
m
−
3
2400
Dilation angle
ψ
°
30
Elastic modulus
E
MPa
32500
Poisson’s ratio
υ
0.2
Initial compressive yield stress
σ
c
0
MPa
9.55
Limited compressive yield stress
σ
cu
MPa
26.8
Initial tensile yield stress
σ
t
0
MPa
2.39
Tensile stiffness recovery parameter
ω
t
0
Compressive stiffness recovery parameter
ω
c
1
Damage variable
(Table3)Table 7
Stress and damage factor versus plastic strain of concrete.
Plastic strain
Stress (MPa)
Damage variable
Compression
0.0000E + 00
18.026
0.00000
2.9233E − 05
21.576
0.01439
1.2438E − 04
29.746
0.04864
2.0932E − 04
32.617
0.07492
3.2177E − 04
34.594
0.10661
6.2346E − 04
36.052
0.18124
8.3638E − 04
35.380
0.22836
1.0796E − 03
33.723
0.27834
1.3385E − 03
31.557
0.32727
1.8656E − 03
26.920
0.41454
2.6240E − 03
21.018
0.51496
3.7809E − 03
14.657
0.62447
6.2920E − 03
9.087
0.73452
9.2686E − 03
5.002
0.83782
Tensile
0.0000E + 00
1.608
0.00000
1.3586E − 05
3.131
0.03406
2.3198E − 05
3.215
0.05538
4.0877E − 05
3.037
0.09857
6.2300E − 05
2.738
0.15603
1.0436E − 04
2.165
0.28146
1.7538E − 04
1.444
0.49663
2.6383E − 04
0.951
0.69271
3.4435E − 04
0.716
0.79632
4.7258E − 04
0.517
0.88134
5.9794E − 04
0.412
0.92187
7.2209E − 04
0.346
0.94431
9.0721E − 04
0.283
0.96303
1.0917E − 03
0.242
0.97345The interaction between underground structure and the surrounding soil is related to their interface properties. In the normal direction of the interface, we define a hard contact. The normal contact compressive stress can mutually transfer via the contact constraint. The element nodes on the interface satisfy Hooke’s Law and the Harmonized Condition of Displacement. In the tangential direction of the interface, tangential contact shear stress is also transferred, and we assume that tangential mechanics behavior of the interface follows the Coulomb friction law with a friction coefficient,µ, of 0.4, which corresponds to a frictional angle of 22° for the soil-underground structure interface. This contact relationship was commonly used by different researchers [14, 40, 41].In pushover analysis, to consider the gravity-load effect, an additional free-field model is established. The mixed lateral boundary condition is used for the lateral sides of the model (i.e., the displacement boundary condition is fixed in the vertical direction and a forced boundary condition is adopted in the horizontal direction). The bottom boundary is fixed in both directions. Therefore, the vertical displacements and horizontal forces of the lateral boundary are obtained through self-gravity response analysis of the additional free-field model and then subjected to the corresponding lateral boundary of the soil-underground structure model. The dynamic analyses are conducted in two steps. First step is a static step, which the geostatic stresses are introduced, considering the underground structure in place. During this step, the base of the analysis model is fixed in the horizontal and vertical directions. The subsequent step is a dynamic analyses step, the horizontal and vertical seismic acceleration is applied uniformly along the base of the model, and the bottom and lateral boundary are free in the horizontal and vertical direction.The viscous damping of the soil-underground structure system is modeled in the form of Rayleigh type. Because this damping type is a linear combination of mass as well as stiffness matrix, it can incorporate into the analysis procedure efficiently. The material damping for soil was selected in different values during dynamic response analysis by different researchers [42–44], mainly ranging from 0.05 to 0.1. In [43], the modeling issues of Rayleigh damping matrix in soil layers are discussed which demonstrate the damping value of 0.08 for soil can better represent a realistic value of soil damping during strong ground motions. Therefore, we utilize a 0.08 damping in the dynamic analysis for soil. Similarly, a widely used 0.05 damping is selected for underground structures.The soil and structure are discretized with 8-node reduced integration linear solid element (C3D8R), while the reinforcement is modeled with three-node linear beam elements (B31). In dynamic analysis, the element mesh for soil is based on Liao’s study [45], which can ensure the efficient reproduction of all the waveforms of whole frequency range under study, and the maximum height of element hmax in soil is determined as(5)hmax=1/75−1/160Vsfmax,where Vs is the shear and compression motion velocity, which can be determined through G0=ρVs2, ρ is the density of soil, and fmax is the maximum vibration frequency of the inputted motion. Therefore, the maximum heights of the element are 1 m to 3 m from the surface to bottom. A finer discretization is adopted near the underground structure.


## 3.2. Selection of Input Motions

The seismic IMs can be expressed in terms of PGA, PGV, orSa (T1, 5%) (i.e., 5% damped first-mode spectral acceleration). As for underground tunnel structures, American Lifelines Alliance [46] produced empirical fragility curves by PGA, while Corigliano [47] proposed empirical fragility curves by PGV. Due to no consensus on the IMs required for seismic analysis of underground structures, PGA and PGV are adopted as IMs to investigate the effects of each parameter in the scattering of the engineering demand parameter (EDP) values. A total of 15 strong motion records are selected from strong motion database of the United States (PEER), Japan (NRIESDP), and China and Taiwan (SMART-I) for each site, respectively, exhibiting different spectral acceleration amplitudes, frequency content, significant duration, and seismotectonic environment. This is considered sufficient to capture the randomness of ground motion and provide sufficient accuracy in estimating seismic demand [48]. These records are all selected from different earthquakes, with source distances larger than 10 km, relatively large interpolated earthquakes (MW ≥ 6), and peak ground acceleration of 0.1 g or more. The horizontal and vertical acceleration response spectra (damping ratio = 0.08) of the selected ground motions are plotted in Figures 7 and 8.Figure 7
Acceleration response spectra of the horizontal component. (a) Site I. (b) Site II. (c) Site III. (d) Site IV.
(a)
(b)
(c)
(d)Figure 8
Acceleration response spectra of the vertical component. (a) Site I. (b) Site II. (c) Site III. (d) Site IV.
(a)
(b)
(c)
(d)The selected motion records are ground motions. Therefore, the horizontal and vertical components are inverted to bedrock to obtain the input motions. According to the established empirical curves which decrease in shear modulus and increase in material damping with the increasing shear strain amplitude for each soil layer, these processes were accomplished through the shake91 code based on the one-dimensional equivalent-linear wave-propagation analyses. For a gradually increasing level of seismic intensity, the horizontal earthquake components are scaled from lower to higher levels of peak acceleration on bedrock (i.e., 0.035 g, 0.1 g, 0.2 g, 0.4 g, 0.6 g, 0.8 g, 1.0 g, and 1.2 g). Then, a series of nonlinear time history analyses of soil-underground structure interaction analysis model are conducted.


## 3.3. Pushover Analysis Results

The capacity curves of the central column considering the vertical seismic effect are obtained (Figure9). To investigate the influence of the vertical seismic effect, we also present the capacity curves without considering the vertical seismic effect (Figure 10). The characteristic points (i.e., θy (the yield drift ratio), θp (the drift ratio of peak loading), and θu (the ultimate drift ratio)) are shown within the figures. We can know the values of θy and θp under considering the vertical seismic loading are almost the same to the values without considering the vertical seismic effect. That is mainly because during the initial loading period, the subjected acceleration is small. Therefore, the vertical seismic effect is not significant compared to self-gravity. With both horizontal and vertical seismic acceleration increasing, the structure gradually enters to the nonlinear stage. Under the vertical seismic effect, the capacity curve decreasing rate is faster compared to only consider the horizontal seismic loading. Therefore, we can find the value of θu under the vertical seismic effect is nearly 50% less compared to without considering the vertical seismic effect. It demonstrates the vertical seismic effect can seriously reduce the ductility of the central column. This finding is in accordance with the previous studies [7, 8, 14–16]. Thus, the seismic capacity of the central column will be overestimated if the vertical seismic effect is not considered.Figure 9
Capacity curves considering the vertical seismic effect. (a) Site I. (b) Site II. (c) Site III. (d) Site IV.
(a)
(b)
(c)
(d)Figure 10
Capacity curves without considering the vertical seismic effect. (a) Site I. (b) Site II. (c) Site III. (d) Site IV.
(a)
(b)
(c)
(d)


## 3.4. Thresholds of Performance Levels

In order to determine the thresholds of different performance levels based on pushover analysis results, in this study, lateral load-drift curve (Figure11) is adopted as the backbone curve model based on the experimental study [49]. This model consists of four performance levels for the component, which are in accordance with the definition of performance levels in this study. It includes four characteristic points: origin point (A), yield point (B), peak strength point (C), and ultimate point (D). The associated drift ratios are θy, θp, and θu. The ultimate point is defined as the force dropped to 85% of the peak force. The values of θOP and θLS are defined as(6)θOP=θy+0.5θp−θy,θLS=θp+0.5θu−θp.Figure 11
Backbone load-drift quantization model.The average values of characteristic points of the capacity curves under different site are used to quantify the thresholds of different performance levels. The quantified thresholds are listed in Table8. For FO and OP performance level, the thresholds(H+V) (i.e., thresholds considering the vertical seismic effect) are the same to thresholds(H) (i.e., thresholds without considering the vertical seismic effect). However, for LS and NC performance level, thresholds(H+V) are obviously less than thresholds(H). This means that the underground structure is easier to reach these two performance levels for the cases considering the vertical seismic effect. If we ignore the vertical seismic effect to the thresholds of performance levels, the seismic capacity of underground structure can be overestimated.Table 8
Thresholds of performance levels.
Performance levels (10−3 rad)
FO
OP
LS
NC
Thresholds(H+V)
1.5
2.2
5.8
8.6
Thresholds(H)
1.5
2.2
7.9
12.8


## 3.5. Incremental Dynamic Analysis Results

In incremental dynamic analysis, the failure condition of the structure has been considered recently within the numerical analysis by the mean defining the collapse condition of the structure [50]. The failure condition can be reflected in the seismic fragility assessment through the mean of the total probability theorem. As for large-space underground structures, they almost reach completely failure condition during earthquake due to the surrounding soil. Therefore, the failure condition is not considered in the present incremental dynamic analysis. The EDP (i.e., drift ratio) and IMs (i.e., PGA and PGV) are set as horizontal and vertical axis, respectively. Based on the spline interpolation between discrete points, then IDA curves of different seismic motions under different sites are obtained (Figures 12 and 13). Due to various values of the material properties of soil such as modulus elasticity, Poisson’s ratio, cohesion, and friction angle, as well as the different ground motion characteristics, each IDA curve shows different developing tendency. Unlike aboveground structures, the IDA curves do not exhibit a noticeable plateau after a certain seismic intensity. This phenomenon, which is due to the soil, prevents structures from reaching the nonlinear stage. Moreover, the restraint of the surrounding soil results in the underground structure not reflecting its self-vibration characteristics. Consequently, the dynamic response between underground structure and aboveground structures exhibits clear differences.Figure 12
IDA curves in terms of PGA as IM. (a) Site I. (b) Site II. (c) Site III. (d) Site IV.
(a)
(b)
(c)
(d)Figure 13
IDA curves in terms of PGV as IM. (a) Site I. (b) Site II. (c) Site III. (d) Site IV.
(a)
(b)
(c)
(d)The dispersion of the IDA curves is mainly due to different input motions (i.e., same PGA, but different ground motion records). To quantify this dispersion, the average of the In (EDP)’s standard deviation is computed and compared. The smaller value of In (EDP)’s standard deviation represents lower dispersion so the associated IM is more appropriate for IDA curves. Table9 listed the average of In (EDP)’s standard deviation of the drift ratio with respective PGV and PGA as the IM of different sites. It is found that PGA is more suitable than PGV as IM in stiff site and has a relatively large scatter in soft site, while PGV shows less scatter as IM in soft site. In order to yield the regression coefficient in equation (3), we select PGV as the representative parameter of seismic IM in the following fragility curves because its mean In (EDP)’s standard deviation is smaller than PGA.Table 9
Average ofσlnEDP.
IM/site
Site I
Site II
Site III
Site IV
PGA
0.403
0.436
0.451
0.503
PGV
0.444
0.441
0.432
0.469


## 3.6. Fragility Analysis

The derivations of fragility curves are based on the constructions of diagrams of the obtained performance index thresholds (Table8) versus PGV at the ground surface. The diagrams are estimated through linear regression analyses. In these analyses, the natural logarithm of the performance index is set as the dependent variable, and the natural logarithm of PGV is set as the independent variable (Figure 14).Figure 14
Linear regression analyses of IDA results. (a) Site I. (b) Site II. (c) Site III. (d) Site IV.
(a)
(b)
(c)
(d)The sets of fragility curves of large-space underground structures derived for each site are given in Figure15. The comparison reveals that the vulnerability of the large-space underground structure is increased when it is embedded in softer site conditions, i.e., from site I to site IV. Therefore, in this context, the important effect of sites conditions on the vulnerability of large-space underground structure is highlighted. It is observed that the probability of the NC performance level is very low in site I, site II, and site III. These phenomena are in accordance with the obtained observations from recent strong earthquakes. The site which Daikai subway station located is similar to site II in this study. However, the obtained fragility curves of large-space underground structure embedded in site II are less satisfactory in the case of Daikai subway station. This difference may be caused by the difference of the design and construction standard of the large-space underground structure. The current obtained fragility curves represent the underground structures designed and constructed based on the improved modern seismic code. However, Daikai station is designed without considering the modern seismic codes. In addition, these fragility curves are developed for square section, which is not the same as the vulnerable central columns in Daikai subway station.Figure 15
Comparison of fragility curves of different sites under different performance levels: (a) FO, (b) OP, (c) LS, and (d) NC.
(a)
(b)
(c)
(d)In Figure16, the fragility curves yielding from thresholds(H+V) and thresholds(H) are compared. As the thresholds(H+V) and thresholds(H) are the same in FO and OP performance levels, we only compare the fragility curves derived from LS and NC performance levels. The exceedance probability of damage in the case of considering the vertical seismic effect is increased from 1% to 13% for LS performance level in both site conditions. For the CP performance level, this increase ranges between 1% and 15%. This demonstrates, when considering the vertical earthquake motion in seismic fragility analyses of large-space underground structures, the exceedance probabilities will be underestimated if thresholds(H) are adopted, resulting in an unfavorable assessment result.Figure 16
Comparison of fragility curves obtained from thresholds(H+V) and thresholds(H). (a) Site I. (b) Site II. (c) Site III. (d) Site IV.
(a)
(b)
(c)
(d)


## 4. Conclusions

An approach is proposed to construct fragility curves for large-space underground structures considering the vertical seismic effect. Thresholds of performance levels are obtained based on the soil-underground structure pushover analysis method, which considers the contribution of vertical seismic effect to seismic capacity. The seismic demand of the large-space underground structure is evaluated through 3D incremental dynamic analyses, accounting for soil inelasticity and ground motion characteristics. Based on a selected typical large-space underground structure, this approach is applied to the derivation of seismic fragility curves of different sites, accounting for soil-underground structure interaction. The fragility curves are compared, highlighting the important role of site and the vertical seismic effect in the vulnerability of large-space underground structures. The main findings of this study are as follows:(1)
For fully operational and operational performance level, the thresholds are the same whether considering the vertical seismic effect. For life safety and near collapse performance level, the thresholds considering the vertical seismic effect are obviously less compare to without considering the vertical seismic effect. If we ignore the vertical seismic effect to the thresholds of performance levels, the seismic capacity of the large-space underground structure will be overestimated.(2)
The vulnerability of large-space underground structures is increased when the site is getting softer. It also reveals that the exceedance probabilities of damage for the thresholds considering the vertical seismic effects are larger than the exceedance probabilities of damage for the thresholds without considering the vertical seismic effect in life safety and near collapse performance levels.(3)
When considering the vertical earthquake motion in a seismic fragility analysis of large-space underground structures, the use of performance level thresholds without considering the vertical seismic effect will underestimate the exceedance probabilities of damage, which can result in an unfavorable assessment result. With the proposed fragility analysis approach, the vertical seismic effect can be considered quantitatively on obtaining the seismic fragility curves of large-space underground structures.



---

*Source: 9650294-2019-04-07.xml*