Structures 50 (2023) 1584–1597

Contents lists available at ScienceDirect

Structures
journal homepage: www.elsevier.com/locate/structures

Seismic damage evaluation of double-box metro tunnel accounting for soil

stiffness using three-dimensional finite element analysis
Van-Quang Nguyen a, b, Yong-Gook Lee a, Duhee Park a, *
a
Department of Civil and Environmental Engineering, Hanyang University, Seoul 04763, South Korea
b
Department of Civil Engineering, Vinh University, Vinh 461010, Viet Nam

A R T I C L E I N F O A B S T R A C T

Keywords: Development of a generalized damage state and index system for prediction of the seismic performance levels of
Double box tunnel box tunnels that accounts for coupled axial - flexural - shear response of the reinforced concrete lining and soil
Seismic damage analysis stiffness is a challenge. The objectives of this study are to develop physics-based definitions of damage states for
Finite element analysis
double-box metro tunnels and associated damage indices (DIs) accounting for the soil-tunnel relative stiffness
Damage state
Damage index
and inelastic behavior of concrete lining. A rigorous three-dimensional (3D) finite element (FE) model that
simulates concrete structure and steel rebar cage individually is used to monitor the damage accumulation. It is
proposed to utilize both the structural element response and seismic capacity curve to define four performance
levels. Because the drift ratio, a widely used DI, is demonstrated to be sensitive to soil stiffness, an alternative
index is linked with the damage states. The recommended DI is the ratio of the elastic moment (M) to the yield
moment (My ), because it less dependent on the soil stiffness compared with the drift ratio. Another advantage of
M/My is that the performance level can be determined from an elastic analysis, whereas an inelastic simulation
must be performed to calculate the drift ratio.

1. Introduction stiffness is needed to cover a wide range of different soil-tunnel

configurations.
Underground structures have been reported to be relatively safe The performance levels of box tunnels have been estimated using two
under earthquake excitations compared with above-ground structures approaches, which are 1) seismic capacity curve which relates hori­
[1–3]. However, extensive damage has been observed during recent zontal load with drift ratio and 2) structure element damage. The seismic
strong earthquake events. Typical examples of such events include the capacity curves have been widely used to predict the seismic perfor­
1995 Kobe earthquake in Japan [4], 1999 Chi-Chi earthquake in Taiwan mance of box tunnels. Zhuang, Ren [17] and Zhong, Shen [18] per­
[5], 1999 Duzce earthquake in Turkey [6], and 2008 Wenchuan earth­ formed a series of 2D FE numerical analyses of a large subway station
quake in China [7]. Therefore, there is a need to better understand the with two stories and three spans. The seismic performance levels and
seismic performance, structural damage patterns, and corresponding corresponding DIs related to the inter-story drift angle and ratio were
damage indices (DIs) of underground structures. presented. Du, Jiang [19] performed pushover analysis using 18 subway
Seismic damage assessment of box tunnels has been the subject of station structures embedded in inhomogeneous site profiles. The inter-
several studies [8–38]. A body of literature evaluated the damage evo­ story drift ratio limits associated with various performance levels were
lution of box tunnels, focusing on the collapse mechanism of the Daikai defined. Zou, Liu [40] proposed the damage states and corresponding
station structure during the 1995 Kobe earthquake [9,11,14,26,27,39]. DIs for three-story and two-span subway stations based on 2D pseudo-
It was concluded that the center column is the most sensitive structural static analyses. He, Xu [23] defined five damage states and associated
member, eventually leading to the collapse of the metro station. How­ inter-story drift ratios for the one-story and six-span subway station
ever, the observations are limited to the site-specific soil profile and the considering spatially varying soil properties. Zhong, Shen [22], Zhuang,
tunnel constructed in 1964, long before seismic design was applied to Yang [16], and Han, Xu [31] characterized damage states of the Daikai
underground structures. A typical tunnel located in different soil subway station based on the drift ratios. However, previous studies focus

* Corresponding author.
E-mail address: dpark@hanyang.ac.kr (D. Park).

https://doi.org/10.1016/j.istruc.2023.02.092
Received 2 September 2022; Received in revised form 19 February 2023; Accepted 20 February 2023
Available online 4 March 2023
2352-0124/© 2023 Published by Elsevier Ltd on behalf of Institution of Structural Engineers.
V.-Q. Nguyen et al. Structures 50 (2023) 1584–1597

on the site specific estimation of the tunnel damage and neglect the
potential influence of the soil stiffness on the tunnel damage charac­
teristics. Additionally, the applicability of the seismic capacity curves
developed for superstructures to tunnels has not been closely investi­
gated. The empirical procedures used to define the performance levels
from a seismic capacity curve need to be evaluated specifically for
tunnels.
In addition to the drift ratio, structure element level damage has been
used to define the damage states. Lee, Park [28], Nguyen, Lee [29],
Nguyen, Park [30], and Le, Huh [41] performed simplified 2D frame
analyses to simulate the structural damage accumulation. The number of
plastic hinges that develop in the reinforced concrete tunnel lining due
to flexural failure was associated with the damage states. Considering
the difficulty in performing inelastic analysis of the tunnel to simulate
plastic hinge formation, alternative indirect procedure to estimate the
damage states was presented. Each damage state was associated with the
moment ratio, which is defined as the ratio of the elastic moment (M)and
yield moment (My ) of the tunnel lining. Parallel elastic analyses were
Fig. 1. Geometry of the tunnel structure section (all dimensions are in performed to determine the elastic moments for each damage state deter­
millimeters). mined from inelastic analysis. Although utilization of the moment ratio

Fig. 2. Details of reinforcement of the tunnel structure section (all dimensions are in millimeters).

1585
V.-Q. Nguyen et al. Structures 50 (2023) 1584–1597

Quantifiable definitions of each damage state are provided. Second

objective is to develop associated DIs that 1) accounts for the relative
soil-tunnel stiffness and 2) can be calculated from an elastic analysis for
easy application in practice. Because of the importance of capturing the
coupled axial-flexural-shear response of reinforced concrete lining, a
rigorous 3D FE numerical model of the tunnel and surrounding soil was
used. The concrete structure and rebar cage are modeled separately and
assembled. A range of soil stiffness was used to evaluate the effect of soil
stiffness. The unique damage accumulation pattern of a contemporary
double-box metro tunnel is investigated.
It should be noted that this study is an extension of previous works
conducted by the authors. Rigorous 2D and 3D finite element analyses
have been used to simulate the damage accumulation of box tunnels and
validated against the Daikai station collapse [14,15], but they have not
yet been used to develop generalized DI system. Although the index M/
My has been already proposed in previous studies as representative DIs,
simplified frame models capable of only simulating flexural plastic
hinges were used Lee, Park [28] Nguyen, Lee [29] Nguyen, Park [30] Le,
Huh [41]. As presented in the previous section, the coupled effect of
axial force - flexural moment - shear force on the damage of the rein­
forced lining needs to be simulated for prediction of the damage.
Additionally, the rationale for corelating the damage state with the
number of plastic hinges was not provided.

2. Numerical model

2.1. Underground tunnel model

A double-box tunnel structure, which is the most commonly used for

metro transportation in South Korea, was employed in this study. The
tunnel cross-section is 13 m wide and 7 m high. The sidewalls and slabs
Fig. 3. FE models of concrete and rebars used in the numerical simulation.
are 1 m thick. The center column has a clear height of 3.8 m with a cross-
section of 0.4 m × 1 m. The spacing between the columns is 3 m. The
Table 1 tunnel is buried 7 m below the ground surface. The details of the ge­
Concrete material properties. ometry and reinforcement are shown in Fig. 1 and Fig. 2.
A 3D FE model of the soil-tunnel structure system was developed
Parameter Value
using the commercial FE program LS-DYNA [42]. The FE model of the
Mass density (kg/m3) 2500 concrete and steel rebar in the tunnel lining are shown in Fig. 3. The
Elastic modulus, E (GPa) 24.8
Poisson’s ratio 0.18
depth of the model is equivalent to the spacing of the center column of 3
Uniaxial compressive strength, fc (MPa)
′
27.5 m. The concrete structure was modeled using eight-node solid elements,
Uniaxial tensile strength, ft (MPa)
′
3.7 whereas the steel rebars were modeled using beam elements. The mesh
Axial strain at compressive strength, εc1 (%) 0.22 element size is 0.1 m for both the concrete and rebars.
Axial strain at tensile strength, εck1 (%) 0.03 Concrete was simulated using the Winfrith concrete model (MAT
Ultimate strain value, εo (%) 0.14
084/085). This model has been used in numerous seismic analyses, as
reported by Nguyen, Nizamani [12] and Coleman [43]. The concrete
(M/My ) as DI allows quick and easy estimation of the damage, previous model properties used by Lee, Park [28] are listed in Table 1. The
studies which revealed several drawbacks. Coupled flexural and shear fracture energy, a fracture mechanics parameter used to describe the
damage modes likely to be induced in the tunnel lining was not simu­ concrete resistance to cracking, was set to 65 Nm/m2 for concrete grade
lated. Therefore, shear type of damage was not explicitly modeled. The C30 with a maximum aggregate size of 8 mm [44]. Fig. 4a illustrates the
number of plastic hinges cannot be directly associated with performance uniaxial stress–strain curve of the Winfrith concrete model. In LS-DYNA,
level. an erosion option can be incorporated with the Winfrith model to
The literature review reveals that the procedures currently used in simulate the damage and failure of concrete using the card *MAT_A­
practice and previous studies have numerous critical limitations. The DD_EROSION. The element is eroded (removed) when its principal
damage states have been defined qualitatively without considering the strain exceeds ± 20% [45], which is in line with the threshold value
unique structural performance of double-box tunnels. There is a need to utilized in previous studies [12,46–49].
define the damage states based on the quantifiable element level The response of the steel rebars was modeled using a bilinear kine­
structural response, thus removing any ambiguity in their definition. matic hardening material model (MAT 003). It is a commonly used
The effect of soil stiffness has been ignored in many studies, resulting in model in seismic analyses [12,50,51]. The failure strain for erosion was
non-generic DIs. Calculation of seismic capacity curves along the length set to 20%, as recommended by Thai and Kim [52], Omran and Mollaei
of the tunnel is not realistic in practice considering the difficulties in the [53], and Nguyen, Nizamani [12]. The material properties of the steel
simulation procedure and the need to develop multiple models for rebars were adopted from Lee, Park [28], as listed in Table 2. The
different site conditions. stress–strain relationship of the steel rebar model is depicted in Fig. 4b.
The first objective of the study is to first develop a physics-based Notably, although a particular set of structural properties is used for the
definition of five damage states of double-box metro tunnels based on tunnel structure, a wide range of the flexibility ratio (F) is considered by
both the structural element response and seismic capacity curve. using different soil stiffness (presented in Section 2.2), where F repre­
sents the relative stiffness between soil and tunnel [54].

1586
V.-Q. Nguyen et al. Structures 50 (2023) 1584–1597

Fig. 4. Material models of the reinforced concrete structure: (a) concrete, (b) steel rebars.

recognized that the nonlinear soil behavior plays a critical role in the
Table 2 seismic performance of tunnels. The equivalent linear procedure is most
Reinforcement material properties. often used to represent the nonlinear soil response. In this method, the
Parameter Value effective shear strain, most often defined as 0.65 of the maximum strain,
Mass density (kg/m3) 7800 at the mid depth of the tunnel is initially calculated. A one-dimensional
Young’s modulus, E (GPa) 200 ground response analysis is typically used to extract the effective shear
Tangent modulus, Et (GPa) 0.4 strain. Then, the normalized shear modulus reduction curve is used to
Poisson’s ratio 0.3
determine the secant shear modulus at this effective strain. The equiv­
Yield strength, fy (MPa) 413
Yield strain, εy (%) 0.2
alent linear method is essentially a linear procedure which indirectly
Ultimate strength, fu (MPa) 620 accounts for the nonlinear stress versus strain behavior. This method has
Ultimate strain, εu (%) 20 been reported to provide favorable estimate of the nonlinear soil
response [57]. It is recommended to use this procedure in determining
the equivalent linear shear wave velocity. Therefore, the range of shear
2.2. Soil domain model wave velocities used in this study should be considered as equivalent
linear values, rather than elastic properties. The range of velocities used
Four uniform soil profiles were used, with shear wave velocities of in this study cover the wide spectrum of in-situ soil stiffness and motion
100, 200, 300, and 400 m/s. The density and Poisson’s ratio of the soil intensity encountered in the field.
profiles were set as 1800 kg/m3 and 0.3, respectively. The numerical
model of the soil-tunnel configuration is shown in Fig. 5. The 3D FE
2.3. Contact and boundary condition
model was 140 m in width, 30 m in height, and 3 m in depth. The width
and height of the soil domain were selected based on a sensitivity
The contact between the rebar and concrete was assumed to be fully
analysis conducted by the authors such that the boundaries of the FE
bonded by using the *CONSTRAINED_LARGRANGE_IN_SOLID option.
model does not influence the tunnel response. The depth of the nu­
The soil-tunnel interface contact was modeled using the *AUTOMA­
merical model was selected to be equal to the column spacing of the
TIC_SURFACE_TO_SURFACE_CONTACT option. This contact option al­
tunnel. The soil was modeled using an eight-node solid element with a
lows slipping and gapping between the slave (soil) and master (concrete)
uniform mesh size of 0.5 m.
segments during the loading [12].
The soil material was assumed to be elastic, although it is well
Dynamic analysis is broadly recognized to most closely represent the

Fig. 5. Finite element model for the numerical simulation of the soil-tunnel configuration.

1587
V.-Q. Nguyen et al. Structures 50 (2023) 1584–1597

Fig. 6. The damage patterns of the tunnel for the case of VS = 100 m/s.

seismic response of tunnels [29,55,56]. However, the pseudostatic 3. Results and discussions
method [57,58], where the transient seismic motion is simplified as an
equivalent static load, has also been reported to produce reasonable 3.1. Damage evolution of the box tunnel
estimates of the tunnel response [55,59,60]. The pseudo-static method
was applied in this study [57,58] to better capture the evolution of the Fig. 6 to Fig. 9 present the damage accumulation patterns of the
damage within the tunnel with an incremental increase in the seismic tunnel for VS = 100, 200, 300, and 400 m/s, respectively. Generally,
demand. Accordingly, the free-field racking deformation was assigned to cracks and eroded concrete elements initially occur at the two ends of
the side and top boundaries of the soil domain. A uniform horizontal the center column and then spread to the four corners of the tunnel.
displacement was imposed at the top boundary, and applied in the form These observations are in agreement with the findings of Sayed, Kwon
of an inverted triangle along the lateral boundaries to produce a uniform [11], Nguyen, Nizamani [12], Iida, Hiroto [4], and Ma, Lu [61], where
shear strain field within the soil domain. The bottom boundary of the 3D FE models were used to investigate the damage mechanism of a box
soil was fixed in both horizontal and vertical directions. The boundary tunnel (i.e., Daikai station). Rebar yielding initiates at the base of the
conditions are illustrated in Fig. 5. A gravitational load was applied center column, followed by the top of the center column, top-left corner,
before displacement assignment. base-right corner, base-left corner, and top-right corner (the 5th and 6th
rebar yielding is discussed in the next section). The sequence of rebar
2.4. Numerical model validation yielding is inconsistent with that of Lee, Park [28], who reported that it
first occurs at the four tunnel corners and then at the two ends of the
The 3D FE numerical model was validated against the observed center column. A possible explanation for this is that the center column
damage of the Daikai station during the 1995 Kobe earthquake in Japan was modeled as a continuous slab in the 2D frame model, leading to
[4]. The collapsed tunnel was well simulated using the numerical model. higher stiffness of the center column system. In addition, the shear
Additionally, the calculated surface settlement also matched favorably failure was not simulated in the study of Lee, Park [28].
with the measurement. Details of the numerical model validation are
presented in Nguyen, Nizamani [12]. Notably, dynamic analysis was 3.2. Development of seismic capacity curves
performed in the study of Nguyen, Nizamani [12], whereas the pseu­
dostatic analysis was used in this study. However, as previously pre­ The performance-based seismic design (PBSD) establishes perfor­
sented, the pseudostatic analysis is reported to provide reasonable mance targets for the seismic design of structures. Performance levels
estimates compared with the dynamic analysis [55,59,60]. are associated with the damage caused by different levels of seismic
hazards. PBSD procedures have been widely integrated into seismic
design codes for above-ground structures [62–69]. However, their
application to underground structures is limited. Furthermore, damage
states and corresponding damage indices have been developed for

1588
V.-Q. Nguyen et al. Structures 50 (2023) 1584–1597

Fig. 7. Damage patterns of the tunnel for the case of VS = 200 m/s.

specific soil-tunnel configurations [19,22,23], thereby failing to account the crack and yield performance levels, the drift ratios are within narrow
for the effect of soil stiffness. In this study, the performance levels are ranges. However, the drift ratios at the peak and failure are significantly
defined from both the structural element response and seismic capacity different. This is likely to be highly sensitive to the inelastic constitutive
curve, which presents the relationship between the lateral force and drift model employed for concrete and steel. Moreover, the drift ratios at the
ratio, as illustrated in Fig. 10 [18,19,70]. Four performance levels are peak state are significantly higher than those determined in previous
defined: crack, yield, peak, and collapse. The onset of the crack level is studies.
determined using the concrete model, whereas the yielding level is
initiated as the steel rebar yields. The peak and collapse performance
levels are determined using a seismic capacity curve. Accordingly, they 3.3. Evaluation of the empirical procedure for estimating the yield
are defined as the levels at which the ultimate shear capacity (peak) of performance level
the tunnel (Vmax ) and 0.85 Vmax are reached, respectively [19,22,23,71].
Five damage states are defined based on the four performance levels, as This section compares the drift ratios at which the rebars start
shown in Fig. 10. The no-damage state represents the condition before yielding calculated from the numerical simulation and the yield drift
the development of concrete cracks, whereas minor damage state is ratio determined from the empirical methods. When calculating the
defined as the condition after concrete cracks occur, but the steel rebars ductility factors, the yield deformation (displacement, rotation or cur­
do not yield. Moderate damage state ranges from the yield-up to the vature, and drift ratio) must be defined. Different empirical procedures
peak damage level, after which extensive damage state follows. The used to estimate the yield drift ratio are shown in Fig. 13. These are as
collapse damage state begins at the collapse performance level. follows:
The seismic capacity curves of the tunnel enclosed within the four
soil profiles based on the numerical simulation results are illustrated in • Procedure 1 [22,23,72]: The yield drift ratio is the drift ratio of the
Fig. 11. The lateral load and drift ratio pairs at which the concrete point where the extension line of the initial stiffness intersects the
cracks, rebar yields, Vmax and finally 0.85 Vmax are reached are marked horizontal line corresponding to the peak shear force (Fig. 13a).
on the curves. Vmax is shown to increase with an increase in VS . Both the • Procedure 2 [71–74]: The yield drift ratio is the drift ratio of the
drift ratios at which cracks develop and rebar yields also increase with point of intersection between the horizontal line at Vmax and the
an increase in VS . However, the drift ratios at Vmax and 0.85 Vmax are secant line at a load of 0.75 Vmax (Fig. 13b).
negatively correlated with VS . For VS = 100 and 200 m/s, ductile re­ • Procedure 3 [18,19]: The yield drift ratio is the drift ratio of the
sponses are observed, whereas brittle curves are observed for VS = 300 intersection point between the horizontal line at Vmax and the secant
and 400 m/s. Overall, it is demonstrated that the drift ratios at which line at the yield point determined in Procedure 1. (Fig. 13c).
damage occurs are sensitive to VS of the surrounding soil. Therefore, the
drift ratio is not an ideal damage index. Fig. 14 compares the calculated drift ratios at which the yield per­
Fig. 12 presents a comparison of the drift ratios calculated for the formance level is reached for all soil profiles and those determined
four VS with those determined in previous studies [16,18,19,22,23]. For empirically. It can be observed from the figure that Procedure 1 (P1)
produces the best estimate of the initiation of the yield compared to the

1589
V.-Q. Nguyen et al. Structures 50 (2023) 1584–1597

Fig. 8. Damage patterns of the tunnel for the case of VS = 300 m/s.

other two procedures. However, it should be noted that all empirical shear strains at which each of the four performance levels decrease with
procedures result in pronounced differences compared with the refer­ an increase in soil stiffness. This is because, as the soil stiffness increases,
ence numerical outputs. Therefore, it is recommended that a numerical the flexibility of the tunnel also increases, thereby making it more
analysis be used to estimate the onset of yield. vulnerable. It is demonstrated that the rebars tend to yield in pairs (i.e.,
1st-2nd, 3rd-4th, 5th-6th). This trend is consistent with the results re­
ported by Lee, Park [28]. The rebar yielding sequences are identical for
3.4. Damage indices of tunnel structures
all the soil profiles, as shown in Fig. 15. Moreover, it is also revealed that
the number of yields is associated with the damage state. The first four
As presented in the previous section, the drift ratio is not an ideal
rebar yielding sequences occur within the moderate damage state,
index for predicting the damage state of a box tunnel owing to its
whereas the last two rebar yields occur in the collapse damage state.
dependence on VS of the soil. In this section, the applicability of M/My as
Therefore, the number of yields can be used as an approximate index of
a potential damage index is investigated, where M denotes the elastic
the damage state.
moment demand, and My denotes the moment capacity (yield moment)
Fig. 16 displays DIs plotted against the free-field shear strain for
at the critical section. It should be noted that M/My was used as the
different soil conditions. The markers in each line denote the four per­
damage index in a study by Lee, Park [28]. Parallel 3D FE elastic ana­
formance levels defined previously. The minor damage state initiates at
lyses were performed to determine the elastic moment at the critical
M/My = 0.8, whereas the moderate damage state is triggered at M/My =
section at the bottom of the center column.
2.35. The crack and yield performance levels are shown to be inde­
Fig. 15 shows the plots of the elastic and inelastic bending moments
pendent of VS . The extensive damage state initiates at M/My = 2.35 to
against the free-field shear strain at the base of the center column, where
2.65, whereas the collapse damage state begins at M/My from 2.8 to
the yielding initiates. Six square markers denote the strain and bending
3.11. For conservative estimates, it is recommended to useM/My = 2.35
moment pairs at which the rebar yields at each of the four outer corners
of the tunnel, in addition to the top and bottom of the center column. and 2.80 for extensive and collapse damage states, respectively.
Evidently, the slopes of bending moment curves in elastic analyses are Table 3 summarizes the damage states, damage descriptions, and
identical to those of inelastic analyses from the beginning of the loading corresponding DIs proposed in this study. Moreover, the DIs proposed by
until the concrete cracks. This slightly differs from the result of Lee, Park Lee, Park [28] are presented for comparison purposes. Notably, the
[28], where it was reported that the curves have the same slopes until collapse damage state was not defined by Lee, Park [28]. The thresholds
the first rebar yielding. The difference is primarily caused by the limi­ for the minor, moderate, extensive, and collapse damage states were 1.0,
tation of the frame element model used in Lee, Park [28], which cannot 1.40, 2.35, and 2.80, respectively. The DIs proposed in this study differ
simulate the crack behavior of the concrete material. from those proposed by Lee, Park [28]. As previously explained, the
After the cracks develop, the bending moment curves determined possible reasons for this difference are 1) the capability to model shear
from the inelastic analyses are nonlinear, whereas they are linear re­ failure and 2) a 3D representation of the tunnel. In other words, a 2D
sponses in the elastic analyses. It is also demonstrated that the free-field beam-spring model utilizing a fiber element for the structural element

1590
V.-Q. Nguyen et al. Structures 50 (2023) 1584–1597

Fig. 9. Damage patterns of the tunnel for the case of. VS = 400 m/s.

Fig. 10. Damage state definition.

cannot capture the critical characteristics of double-box tunnel damage. 4. Conclusions

It should be noted that although the tunnel section used in this study
is the most common type used for metro transportation in South Korea, A box tunnel undergoes complex structural response induced by in­
but it may differ for other countries. A complete analysis of the effect of crements in the axial and shear forces, along with moment variations
the tunnel section is required to assess the variability of the DI. The when subjected to strong earthquakes. A generalized damage state and
proposed DIs should only be used for box tunnels with similar size and associated index system for prediction of the seismic performance levels
overburden depth to those employed in this research. of box tunnels that accounts for this coupled axial - flexural - shear
response of the reinforced concrete lining and soil stiffness has not yet
been developed. In this study, physics-based definitions of damage states
are proposed for a double-box metro tunnel. Then, associated DIs ac

1591
V.-Q. Nguyen et al. Structures 50 (2023) 1584–1597

Fig. 11. Seismic capacity curves of the tunnel located in different soil stiffness.

2.5
Average (this study)

2 Vs = 100 m/s
Vs = 200 m/s
Drift ratio (%)

1.5 Vs = 300 m/s

Vs = 400 m/s
1 Zhong et al. (2020)
Zhuang et al. (2021)
0.5 Du et al. (2021)
Zhong et al. (2021)
0 He et al. (2022)
Crack Yield Peak Failure

Fig. 12. Comparison of the drift ratios at four performance levels determined in this study and previous studies.

counting for the soil-tunnel relative stiffness and inelastic behavior of determine the performance levels, it is proposed to use the ma­
concrete lining are presented. Rigorous 3D FE seismic damage analyses terial model outputs to define the crack and yield performance
were performed for a range of soil-tunnel configurations to capture the levels. The peak and collapse performance levels are character­
damage accumulation. A bilinear model with the ability to map cracks ized from the seismic capacity curve, in line with previous
was used to simulate the concrete material behavior, whereas an elas­ studies. All empirical procedures that estimate the onset of yield
tic–plastic model was employed to imitate the steel rebar response. An from the seismic capacity curves are demonstrated to provide
elastic model was used to model the soil. A single-story double-box poor predictions of the actual steel rebar yielding. Therefore, it is
tunnel, the most widely designed and built for metro subway systems in recommended to use the physics-based definitions for the per­
South Korea, was used. A range of stiffness was used to capture the effect formance levels. Five damage states are defined based on the four
of soil stiffness. Based on the numerical results, the following conclu­ performance levels.
sions are drawn. (3) The drift ratios are shown to be sensitive to the stiffness of the
surrounding soil. Therefore, using the drift ratio as the DI is not
(1) The damage accumulation sequence of a contemporary double- recommended. An alternative DI, defined as the ratio of the
box metro tunnel is investigated. Similar to the widely studied elastic moment (M) to the yield moment (My ), is less dependent
Daikai station, which collapsed during the 1995 Kobe earth­ on the soil stiffness compared with the drift ratio, and therefore
quake, the damage initiates in the center column. After the steel recommended. Another advantage of using M/My is its capability
rebars yield in the center column, the damage is propagated to to determine the damage level from an elastic analysis, whereas
the four corners of the tunnel. However, the metro tunnel does an inelastic analysis using a rigorous 3D model needs to be per­
not fully collapse, whereas the Daikai station suffered a structural formed to calculate the drift ratio.
collapse. The advanced 3D numerical model that models the
concrete structure and steel cage separately is capable of
capturing the structural element level damage including cracks Declaration of Competing Interest
and yielding of steel rebars, as well as shear damage.
(2) A new procedure for defining the four performance levels is The authors declare that they have no known competing financial
proposed. Instead of using only the seismic capacity curve to interests or personal relationships that could have appeared to influence
the work reported in this paper.

1592
V.-Q. Nguyen et al. Structures 50 (2023) 1584–1597

Fig. 13. Empirical procedures to define the yield performance level.

1593
V.-Q. Nguyen et al. Structures 50 (2023) 1584–1597

Fig. 14. Comparison of the yield performance levels determined from the numerical simulations and empirical procedures.

1594
V.-Q. Nguyen et al. Structures 50 (2023) 1584–1597

Fig. 15. Comparison of the bending moment calculated at the base of the center column from elastic and inelastic analyses.

Fig. 16. M/My ratios corresponding to the damage states.

1595
V.-Q. Nguyen et al. Structures 50 (2023) 1584–1597

Table 3 [21] Chen Z, Chen W, Li Y, Yuan Y. Shaking table test of a multi-story subway station
Proposed damage states and corresponding damage indices. under pulse-like ground motions. Soil Dyn Earthq Eng 2016;82:111–22.
[22] Zhong Z, Shen Y, Zhao M, Li L, Du X, Hao H. Seismic fragility assessment of the
Damage Damage description Damage Previous DI Daikai subway station in layered soil. Soil Dyn Earthq Eng 2020;132:106044.
state index, DI (M/My ) of Lee, [23] He Z, Xu H, Gardoni P, Zhou Y, Wang Y, Zhao Z. Seismic demand and capacity
(M/My ) Park [28] models, and fragility estimates for underground structures considering spatially
varying soil properties. Tunn Undergr Space Technol 2022;119:104231.
None No damage DI < 0.7 DI < 1.0 [24] Van Kien D. Numerical simulation of the stability of rock mass around large
Minor Crack development 0.7 ≤ DI < 1.0 ≤ DI < 1.2 underground cavern. Civ Eng J 2022;8:81–91.
1.40 [25] Qiu H, Zhang X, Daddow M. Prediction of ground settlement induced by slurry
Moderate Rebars yielding in the center 1.40 ≤ DI 1.2 ≤ DI < 2.0 shield tunnelling in granular soils. Civ Eng J 2020;6:2273–89.
column < 2.35 [26] Liu J, Liu X. Pushover analysis of Daikai subway station during the Osaka-Kobe
Extensive Concrete crushing in the center 2.35 ≤ DI 2.0 ≤ DI earthquake in 1995. In 14th world conference on earthquake engineering, Beijing,
column < 2.8 China, 2008.
Collapse Rebars yield at all four tunnel 2.8 ≤ DI N.A. [27] Kishi N, Sonoda K, Komuro M, Kawarai T. Numerical simulation of the Daikai
station subway structure collapse due to sudden uplift during earthquake. J Eng
conners and two ends of the
Mech 2021;147:04020152.
center column, concrete crushing
[28] Lee T-H, Park D, Nguyen DD, Park J-S. Damage analysis of cut-and-cover tunnel
in the tunnel conners. structures under seismic loading. Bull Earthq Eng 2016;14:413–31.
[29] Nguyen D-D, Lee T-H, Nguyen V-Q, Park D. Seismic damage analysis of box metro
tunnels accounting for aspect ratio and shear failure. Appl Sci 2019;9:3207.
Acknowledgment [30] Nguyen D-D, Park D, Shamsher S, Nguyen V-Q, Lee T-H. Seismic vulnerability
assessment of rectangular cut-and-cover subway tunnels. Tunn Undergr Space
Technol 2019;86:247–61.
This work was supported by the National Research Foundation of [31] Han R, Xu C, Liu D, Chen J, Du X. Static pushover test of spring-underground
Korea (NRF) grant funded by the Korea government (MSIT) (No. structure system for seismic performance analysis of underground structure. Eng
2022R1A2C3003245). Struct 2022;271:114936.
[32] Ding P, Tao L, Yang X, Zhao J, Shi C. Three-dimensional dynamic response analysis
of a single-ring structure in a prefabricated subway station. Sustain Cities Soc 2019;
References 45:271–86.
[33] Tao L, Ding P, Shi C, Wu X, Wu S, Li S. Shaking table test on seismic response
characteristics of prefabricated subway station structure. Tunn Undergr Space
[1] Dowding CH, Rozan A. Damage to rock tunnels from earthquake shaking.
Technol 2019;91:102994.
J Geotech Eng Div 1978;104:175–91.
[34] Tao L, Ding P, Yang X, Lin P, Shi C, Bao Y, et al. Comparative study of the seismic
[2] Arango I. Theme Paper: Earthquake engineering for tunnels and underground
performance of prefabricated and cast-in-place subway station structures by
structures. A case history. Geotechnical Earthquake Engineering and Soil Dynamics
shaking table test. Tunn Undergr Space Technol 2020;105:103583.
IV2008. p. 1-34.
[35] Tao L, Shi C, Ding P, Li S, Wu S, Bao Y. A study on bearing characteristic and failure
[3] Hashash YM, Hook JJ, Schmidt B, John I, Yao C. Seismic design and analysis of
mechanism of thin-walled structure of a prefabricated subway station. Front Struct
underground structures. Tunn Undergr Space Technol 2001;16:247–93.
Civ Eng 2022;16:359–77.
[4] Iida H, Hiroto T, Yoshida N, Iwafuji M. Damage to Daikai subway station. Soils
[36] Tao L, Shi C, Ding P, Yang X, Bao Y, Wang Z. Shaking table test of the effect of an
Found 1996;36:283–300.
enclosure structure on the seismic performance of a prefabricated subway station.
[5] Wang W, Wang T, Su J, Lin C, Seng C, Huang T. Assessment of damage in mountain
Tunn Undergr Space Technol 2022;125:104533.
tunnels due to the Taiwan Chi-Chi earthquake. Tunn Undergr Space Technol 2001;
[37] Nguyen V-Q, Tran V-L, Nguyen D-D, Sadiq S, Park D. Novel hybrid MFO-XGBoost
16:133–50.
model for predicting the racking ratio of the rectangular tunnels subjected to
[6] Ghasemi H, Cooper JD, Imbsen R, Piskin H, Inal F, Tiras A. The November 1999
seismic loading. Transp Geotech 2022:100878.
Duzce Earthquake: post-earthquake investigation of the structures on the TEM. Rep
[38] Sadiq S, Nguyen VQ, Jung H, Park D. Effect of flexibility ratio on seismic response
No FHWA-RD-00. 2000;146:00-146.
of cut-and-cover box tunnel. Adv Civ Eng Mater 2019;2019.
[7] Wang Z, Gao B, Jiang Y, Yuan S. Investigation and assessment on mountain tunnels
[39] Li W, Chen Q. Seismic performance and failure mechanism of a subway station
and geotechnical damage after the Wenchuan earthquake. Sci China Technol Sci
based on nonlinear finite element analysis. KSCE J Civ Eng 2017:1–12.
2009;52:546–58.
[40] Zou Y, Liu H, Jing L, Cui J. A pseudo-static method for seismic responses of
[8] Zhuang H, Zhao C, Chen S, Fu J, Zhao K, Chen G. Seismic performance of
underground frame structures subjected to increasing excitations. Tunn Undergr
underground subway station with sliding between column and longitudinal beam.
Space Technol 2017;65:106–20.
Tunn Undergr Space Technol 2020;102:103439.
[41] Le TS, Huh J, Park J-H. Earthquake fragility assessment of the underground tunnel
[9] Ma C, Lu D, Du X, Qi C. Effect of buried depth on seismic response of rectangular
using an efficient SSI analysis approach. J Appl Math Phys 2014;2:1073–8.
underground structures considering the influence of ground loss. Soil Dyn Earthq
[42] LSTC. LS-DYNA keyword user’s manual version 971. 2007.
Eng 2018;106:278–97.
[43] Coleman DK. Evaluation of concrete modeling in LS-DYNA for seismic application
[10] Ma C, Lu D-C, Du X-L, Qi C-Z, Zhang X-Y. Structural components functionalities
[Master thesis]. Austin, USA: The University of Texas at; 2016.
and failure mechanism of rectangular underground structures during earthquakes.
[44] Ceb-Fip MC. Design code. Comite Euro International du Beton 1990:51–9.
Soil Dyn Earthq Eng 2019;119:265–80.
[45] Sagals G, Orbovic N, Blahoianu A. Sensitivity studies of reinforced concrete slabs
[11] Sayed MA, Kwon O-S, Park D, Van Nguyen Q. Multi-platform soil-structure
under impact loading. Transactions of the 21st SMiRT. 2011.
interaction simulation of Daikai subway tunnel during the 1995 Kobe earthquake.
[46] Wilt T, Chowdhury A, Cox P. Response of reinforced concrete structures to aircraft
Soil Dyn Earthq Eng 2019;125:105643.
crash impact. Prepared for US Nuclear Regulatory Commission Contract NRC-02-
[12] Nguyen V-Q, Nizamani ZA, Park D, Kwon O-S. Numerical simulation of damage
07-006. 2011.
evolution of Daikai station during the 1995 Kobe earthquake. Eng Struct 2020;206:
[47] Thai D-K, Kim S-E. Numerical simulation of pre-stressed concrete slab subjected to
110180.
moderate velocity impact loading. Eng Fail Anal 2017;79:820–35.
[13] Zhuang H, Wang R, Shi P, Chen G. Seismic response and damage analysis of
[48] Luccioni BM, Aráoz GF, Labanda NA. Defining erosion limit for concrete. Int J Prot
underground structures considering the effect of concrete diaphragm wall. Soil Dyn
Struct 2013;4:315–40.
Earthq Eng 2019;116:278–88.
[49] Bi K, Ren W-X, Cheng P-F, Hao H. Domino-type progressive collapse analysis of a
[14] Lu C-C, Hwang J-H. Nonlinear collapse simulation of Daikai Subway in the 1995
multi-span simply-supported bridge: A case study. Eng Struct 2015;90:172–82.
Kobe earthquake: Necessity of dynamic analysis for a shallow tunnel. Tunn
[50] Ding J-H, Jin X-L, Guo Y-Z, Li G-G. Numerical simulation for large-scale seismic
Undergr Space Technol 2019;87:78–90.
response analysis of immersed tunnel. Eng Struct 2006;28:1367–77.
[15] Li W, Chen Q. Effect of vertical ground motions and overburden depth on the
[51] Luat N-V, Shin J, Park J-H, Lee K. Experimental and numerical investigation on
seismic responses of large underground structures. Eng Struct 2020;205:110073.
seismic performance of retrofitted RC columns with web direct/indirect bonding
[16] Zhuang H, Yang J, Chen S, Dong Z, Chen G. Statistical numerical method for
external H-section. J Build Eng 2021;44:103404.
determining seismic performance and fragility of shallow-buried underground
[52] Thai D-K, Kim S-E. Failure analysis of reinforced concrete walls under impact
structure. Tunn Undergr Space Technol 2021;116:104090.
loading using the finite element approach. Eng Fail Anal 2014;45:252–77.
[17] Zhuang H, Ren J, Miao Y, Jing L, Yao E, Xu C. Seismic performance levels of a large
[53] Omran ME, Mollaei S. Assessment of empirical formulas for estimating residual
underground subway station in different soil foundations. J Earthq Eng 2019:1–26.
axial capacity of blast damaged RC columns. Eur J Sustain Dev 2017;6:383–96.
[18] Zhong Z, Shen Y, Zhao M, Li L, Du X. Seismic performance evaluation of two-story
[54] Wang J-N. Seismic Design of Tunnels: A Simple State-of-the-art Design Approach.
and three-span subway station in different engineering sites. J Earthq Eng 2021:
Parsons Brinckerhoff; 1993.
1–31.
[55] Tsinidis G, de Silva F, Anastasopoulos I, Bilotta E, Bobet A, Hashash YM, et al.
[19] Du X, Jiang J, El Naggar MH, Xu C, Xu Z. Interstory drift ratio associated with
Seismic behaviour of tunnels: From experiments to analysis. Tunn Undergr Space
performance objectives for shallow-buried multistory and span subway stations in
Technol 2020;99:103334.
inhomogeneous soil profiles. Earthq Eng Struct Dyn 2021;50:655–72.
[56] Zhu T, Wang R, Zhang J-M. Evaluation of various seismic response analysis
[20] Chen Z, Chen W, Zhang W, Lou M. Effects of axial compression ratio of central
methods for underground structures in saturated sand. Tunn Undergr Space
columns on seismic performance of a multi-story underground structure. Int J
Technol 2021;110:103803.
Comput Methods 2016;13:1641014.

1596
V.-Q. Nguyen et al. Structures 50 (2023) 1584–1597

[57] Hashash Y, Karina K, Koutsoftas D, O’Riordan N. Seismic design considerations for [66] ASCE41-13. Seismic evaluation and upgrade of existing buildings. Reston, Virginia:
underground box structures. Earth Retention Conference. Bellevue, WA, USA, American Society of Civil Engineers; 2013.
2010. p. 620-37. [67] CEN. Eurocode 8: Design of structures for earthquake resistance-part 1: general
[58] Tsinidis G, Pitilakis K, Madabhushi G, Heron C. Dynamic response of flexible rules, seismic actions and rules for buildings. Bruss Eur Comm Stand 2005.
square tunnels: centrifuge testing and validation of existing design methodologies. [68] Yamanouchi H, Hiraishi H, Ohashi Y, Fujitani H, Okada T. Performance-based
Geotechnique 2015;65:401–17. engineering for structural design of buildings: Building Research Institute, Ministry
[59] Sun QQ, Dias D. Assessment of stress relief during excavation on the seismic tunnel of Construction; 2000.
response by the pseudo-static method. Soil Dyn Earthq Eng 2019;117:384–97. [69] MOHURD. Code for Seismic Design of Buildings. Ministry of Housing and Urban-
[60] Argyroudis S, Pitilakis K. Seismic fragility curves of shallow tunnels in alluvial Rural Development of the People’s Republic of China, Beijing. 2010.
deposits. Soil Dyn Earthq Eng 2012;35:1–12. [70] ATC-40. Seismic evaluation and retrofit of existing concrete building. Applied
[61] Ma C, Lu D, Du X. Seismic performance upgrading for underground structures by Technology Council, RedWood City, CA. 1996.
introducing sliding isolation bearings. Tunn Undergr Space Technol 2018;74:1–9. [71] Jiang H, Fu B, Chen L. Damage-control seismic design of moment-resisting RC
[62] FEMA273. NEHRP Commentary on the guidelines for the rehabilitation of frame buildings. J Asian Archit Build Eng 2013;12:49–56.
building. Washington DC: Federal Emergency Management Agency; 1996. [72] Park R. Evaluation of ductility of structures and structural assemblages from
[63] FEMA274. NEHRP Commentary on the guidelines for the rehabilitation of laboratory testing. Bull N Z Soc Earthq Eng 1989;22:155–66.
building. Washington DC: Federal Emergency Management Agency; 1996. [73] Nguyen HD, Shin M. Effects of soil–structure interaction on seismic performance of
[64] FEMA356. Prestandard and commentary for the seismic rehabilitation of buildings. a low-rise R/C moment frame considering material uncertainties. J Build Eng 2021;
Washington DC: Federal Emergency Management Agency; 2000. 44:102713.
[65] FEMA-P58-1. Seismic Performance Assessment of Buildings (Volume 1-Methodol­ [74] Nguyen HD, LaFave JM, Lee Y-J, Shin M. Rapid seismic damage-state assessment of
ogy). Federal Emergency Management Agency, Washington DC. 2018. steel moment frames using machine learning. Eng Struct 2022;252:113737.

1597