Soil Dynamics and Earthquake Engineering 162 (2022) 107489

Contents lists available at ScienceDirect

Soil Dynamics and Earthquake Engineering

journal homepage: [Link]/locate/soildyn

Seismic responses of rectangular tunnels in liquefiable soil considering

spatial variability of soil properties
Yubing Wang a, Jianjian He a, Shuang Shu b, Houle Zhang b, Yongxin Wu b, *
a
Center for Hypergravity Experiment and Interdisciplinary Research, College of Civil Engineering and Architecture, Zhejiang University, Hangzhou, 310058, China
b
Key Laboratory of Ministry of Education for Geomechanics and Embankment Engineering, College of Civil and Transportation Engineering, Hohai University, Nanjing,
210098, China

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

Keywords: The objective of this paper was to investigate the dynamic response of a tunnel-soil system in the soil with spatial
Tunnels variability of shear modulus with various coefficients of variation and scales of fluctuation. The soil shear
Spatial variability modulus was highlighted and modeled as an isotropic random field using the spectral representation method,
Liquefaction
then mapped into the finite difference model. The dynamic stochastic responses of the excess pore water pressure
Seismic behavior
Dynamic random finite difference analysis
ratio, the liquefied zone, the ground displacement, and the uplift displacement of the tunnel were discussed. It
was observed that the effect of coefficient of variation was pronounced in which an increase of the coefficient of
variation brought out a reduction in the corresponding dynamic stochastic responses. However, the scale of
fluctuation was found to be insignificant in the tunnel soil system.

1. Introduction shape of tunnel and input ground motions on the dynamic responses of
tunnels. Bao et al. [4] studied the seismic behavior of a metro tunnel in
Tunnels, serving as an essential component in the modern trans­ the liquefiable ground by using a fully coupling finite element-finite
portation network and infrastructural development, are being widely difference method. Miao et al. [5] explored the influence of spatially
constructed worldwide. A glance at the underground structures (e.g. varying seismic ground motion on seismic response of a subway shield
tunnels) design shows the seismic design are universally excluded since tunnel. Wang et al. [6] investigated the interaction between the tunnel
underground structures always suffer less damage than surface struc­ and near-surface structure during an earthquake based on the shaking
tures. However, in recent earthquake events, such as Kobe earthquake table tests. Yan et al. [7], Yu et al. [8] and Yu et al. [9] conducted a series
(1995, Japan), ChiChi earthquake (1999, Taiwan), and Wenchuan of shaking table tests on the long immersed tunnels subjected to
earthquake (2008, China), different levels of severe damages stimulated non-uniform seismic loadings. Both numerical and experimental studies
by seismic loads were observed in tunnels. These damages push the indicated that the seismic responses of tunnels under non-uniform
development of the seismic design for tunnels and other underground seismic loads are severer than those under uniform seismic loads.
structures. In this regard, dynamic response analyses underlines its In the numerical studies mentioned above, the soil property indexes
importance in tunnel structure under seismic loads. applied in the dynamic analyses are usually determinated from labora­
In the past decades, due to the deficit of knowledge in understanding tory tests or in-situ measurements. With these specific values, an anal­
the dynamic behavior of tunnels subjected to seismic loads, several re­ ysis in a deterministic manner yields to a solution that is dependent on
searchers have been motivated to carry out a series of investigations the quality of the input characteristic values. However, it is well known
through analytical, numerical, and experimental efforts. For example, that soil properties are aisotropic, that is, vary both vertically and hor­
Lee et al. [1] and Yu et al. [2] gave the analytical bending stiffness of izontally, due to sedimentational, physical or chemical changes in the
shield tunnels in the transverse and longitudinal directions, respectively. environment [10].
Recently, Patil et al. [3] systematically investigated the effects of tunnel In the past two decades, researchers have proved the importance of
embedment ratio, soil-tunnel interface conditions, lining thickness, spatial variability of soil properties for geotechnical designs. Notable

* Corresponding author.
E-mail addresses: wangyubing@[Link] (Y. Wang), hejianjian1992@[Link] (J. He), [Link]@[Link] (S. Shu), ceatczhl@[Link] (H. Zhang),
yxwuhhu@[Link] (Y. Wu).

[Link]
Received 20 May 2022; Received in revised form 5 August 2022; Accepted 5 August 2022
Available online 13 August 2022
0267-7261/© 2022 Elsevier Ltd. All rights reserved.
Y. Wang et al. Soil Dynamics and Earthquake Engineering 162 (2022) 107489

Table 1
Summary of soil parameters used.
Soil property Unit Value
2
Unit weight of soil, γsoil kN/m 15
Shear modulus, G MPa 20
Bulk modulus, B MPa 30
Friction angle, φ degree 25
Cohesion, C kPa 0

Fig. 1. Finite difference model.

achievements cover a wide range of geotechnical problems, such as the

bearing capacity of foundations [11–15], settlements of foundations
[16,17], the stability of slopes [18–23], etc.
Up to now, only a few works have been undertaken to elucidate the
importance of the spatially varying feature of soils in tunneling engi­
neering. For example, Song et al. [24] employed the finite difference
method and explored the effects of spatially varied weathered rock on
various tunnel behaviors, including deformation, elastic-plastic inter­
face, ground reaction curve, and failure mechanism. Mollon et al. [25]
investigated the effects of spatially varied shear strength parameters on
the critical collapse pressure of a pressurized tunnel face. Ali et al. [26]
studied the influence of spatial variabilities on the undrained stability of
an unlined circular tunnel using Random Adaptive Finite Element Limit Fig. 2. Kobe earthquake excitation.
Analysis (RAFELA). Huang et al. [27] presented a detailed numerical
analysis to investigate the probabilistic response of tunnel convergence surface.
in spatially varied soft soils. Chen et al. [28] conducted a simplified The Finn model was used to simulate the nonlinear soil behavior.
procedure to evaluate the failure probability of crossing tunnels Table 1 summarized the soil property indexes in the reference model.
considering the spatial variabilities of rock mass properties. Yue and Ang The relationship between the variations of cyclic shear strain (γ) and
[29,30] studied the stochastic response of tunnels based on the proba­ volumetric strain increment (Δεvd ), proposed by Martin et al. [31], was
bility density evolution method with the spatially varied soil properties shown in Eq. (1). This relationship has been widely utilized in the
in the nonlinear finite element model with 2-D and 3-D, respectively. All analysis of seismic behavior of underground structures [32,33]. The
the studies mentioned above established that the spatially varied soil variations of cyclic shear strain (γ) and volumetric strain increment
properties played a pivotal role in the tunnel behavior. However, the (Δεvd ) was defined as
dynamic stochastic response of tunnel-soil systems could be better
C3 ε2vd
convinced if the spatially varied soil properties included. Δεvd = C1 (γ − C2 εvd ) + (1)
γ + C4 εvd
Thus, this paper mainly focused on the dynamic response of soil-
tunnel behaviors considering the spatial variability of soils, i.e., the
where εvd is the total cumulative volumetric strain; C1, C2, C3 and C4 are
excess pore water pressure ratio, the liquefied zone, the ground
the constant coefficients. Based on the results of cyclic triaxial tests for
displacement, as well as the uplift displacement of the tunnel. The
the crystalline silica sand with the relative density of 45%, and the co­
remainder of this paper is organized as follows: (1) first briefly intro­
efficients were set to be 0.79, 0.52, 0.2 and 0.5, respectively.
duced the non-Gaussian simulation method based on the spectral rep­
The Kobe earthquake’s (1940) time history with scaled to 0.3g
resentation method to generate the random fields; (2) next presented the
(Fig. 2) was applied at the base boundary of the analysis domain. The
numerical modeling using the nonlinear finite-element software com­
length of the earthquake excitation was truncated as 30 s in order to
bined with the random fields; (3) finally investigated the influences of
reduce the computational efforts. The numerical simulation included
the coefficient of variation and the scale of fluctuation on the dynamic
static and dynamic stages. In the static analysis, the bottom boundary
responses.
was constrained in all directions, and the lateral boundaries were fixed
in the horizontal direction. Besides, the back and former boundaries
2. Computational framework
were fixed in the vertical direction to satisfy the plane-strain assumption
in the finite-difference computations. In the dynamic analysis, the
2.1. Finite difference modeling
boundaries were utilized to absorb the reflected waves and enforce the
discrete half-space conditions of the numerical model. The free-field
The nonlinear dynamic behavior of the rectangular tunnel under a
boundaries were adopted for the lateral boundaries to simulate the
seismic load was analyzed using the finite-difference program FLAC3D.
half-space condition. The earthquake loading in the horizontal direction
The computational model is shown in Fig. 1. The tunnel was with a
was applied at the base boundary that assumed to be rigid [34].
dimension of 6 m × 6 m and its center was located 10 m below the
ground surface. The analysis domain had a length of 50 m and a height
of 30 m, with a total of 1500 elements. The water table was at the ground

2
Y. Wang et al. Soil Dynamics and Earthquake Engineering 162 (2022) 107489

Fig. 3. Typical random field realizations of shear modulus.

2.2. Random field simulation power spectral density (PSD) function, which is obtained by:
∫ +∞ ∫ +∞
1 ( )
The random field of the elastic modulus can be defined by the two- SG (κx , κz ) = 2
ρG ξx , ξz e− i(κx ξx +κz ξz ) dξx dξz (7)
dimensional correlation function as following: (2π) − ∞ − ∞
⎡ √̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ ⎤
( )2 ( )2 where ξx and ξz are the distances in the horizontal and vertical di­
2τx 2τz ⎦
ρ = exp⎣ − + (2) rections, respectively. For the correlation function defined in Eq. (5), the
δx δz
corresponding PSD function can be given as:

where τx and τz are the lag distance (τ) in the horizontal and vertical S(κx , κz ) = (
2δx
) (
2δz
) (8)
directions; δx and δz are the scale of fluctuation (δ) in the horizontal and π 4 + δ2x κ2x π 4 + δ2z κ2z
vertical directions.
The random field of the elastic modulus is defined as lognormal
The elastic modulus of soil followed a lognormal distribution with a distribution, which is a typical non-Gaussian random field, while the
mean value at 20 MPa. To investigate the effects of the coefficient of
random field generated by Eq. (3) follows Gaussian distribution.
variation (COVG) on the dynamic responses, three different COVG (0.1, Therefore, the generated Gaussian distribution should be mapped to the
0.3 and 0.5) were selected, while δ stayed constant at 5 m; for the effects
corresponding non-Gaussian distribution, i.e. lognormal distribution in
of scale of fluctuation, δ was assumed at 5, 50 and 100 m, and COVG was this study. A general way to map is through the following translation:
equal to 0.3. Except the parameters mentioned above, the random field
of G was fully defined, with a certain method to generate the realizations y(x, z) = F y − 1 [Φ(g(x, z))] (9)
of the random fields. To simulate random field, several methods have
been developed, such as the local averaging subdivision method [35], where Φ( ⋅) is the standard unit Gaussian cumulative probability func­
the Karhunen-Loève expansion [36], the midpoint method [37] and the tion (CDF) and Fy ( ⋅) is the marginal non-Gaussian CDF of the standard
spectral representation method [38]. In this study, the spectral repre­ unit non-Gaussian y(x, z), with inverse function Fy − 1 ( ⋅). Because of a
sentation method was adopted. nonlinear transformation, the PSD function or the correlation function
For a Gaussian two-dimensional random field g(x, z) with zero mean will be changed. Therefore, a proper underlying Gaussian PSD function
and unit deviation, it can be generated by: contributed that the translated non-Gaussian function capable to match
2 − 1 √̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅[
( the target one shown in Eq. (7). The sample-free iterative scheme is
∑
N ∑
1− 1 N

g(x, z) = SG (κxl1 , κzl2 )Δκx Δκz cos κxl1 x + κzl2 z utilized to find the underlying Gaussian PSD function in this study.
l1 =0 l2 =0 Once the underlying Gaussian PSD function is determined, a stan­
) ( )]
+ φ1l1 l2 + cos κxl1 x − κzl2 z + φ2l1 l2 (3) dard unit Gaussian random field g(x, z) can be generated by Eq. (3), and
then translated to a standard unit Gaussian random field y(x, z) through
Eq. (9). Finally, the desired non-Gaussian random field y(x, z) can be
κxli = l1i Δκx , l1i = 0, 1, …, N1 − 1 (4)
obtained by
κzli = l2i Δκz , l2i = 0, 1, …, N2 − 1 (5) y(x, z) = y(x, z) × COVG × μG + μG (10)

Δκx = κxu /N1 , Δκz = κzu /N2 (6) where μG and COVE are the mean value and COV of the elastic modulus,
respectively.
where Δκx and Δκz are the discretization steps in the wave number Compared with the conventional static random field analysis, the
domain along horizontal direction and vertical direction, respectively; dynamic random field analysis required more additional computational
κxu and κzu are the corresponding upper cut-off wave numbers; φ1l1 l2 and efforts. It required approximately 3 h completing one finite difference
φ2l1 l2 are two sets of N1 N2 independent random phase angles uniformly realization on the authors’ computer (3.8 Hz CPU and 64G memory).
distributed over the interval 0 to 2π. SG (κx , κz ) is the target Gaussian Hence, to gather statistics of the dynamic responses, 200 realizations of

3
Y. Wang et al. Soil Dynamics and Earthquake Engineering 162 (2022) 107489

ratio, as depicted in Fig. 4. The EPWP ratio was calculated as the excess
pore water pressure divided by the corresponding initial effective stress
in each selected point or layer.
Fig. 5 plots the EPWP ratio varying with time history under different
COVG at Points A and B. Different variation trends of EPWP ratio were
observed at the points above and below the tunnel. As evident in Fig. 5
(a), the stochastic time history curves of the EPWP ratio were almost
identical with the deterministic curve at Point A, indicating that COVG
had an insignificant effect on the EPWP ratio above the tunnel. How­
ever, as presented in Fig. 5(b), a noticeable difference was found in the
range of 6–22s where a higher COVG corresponded to a lower EPWP
ratio, and the EPWP ratios calculated in the random fields with different
COVG at the end of the analysis were similar compared to the deter­
ministic result.
The influence of δ on the EPWP ratio is illustrated in Fig. 6. It clearly
Fig. 4. Display of the selected points and layers. demonstrated that the time history curves under different scales of
fluctuation were similar, indicating that the scale of fluctuation had no
significant effect on the EPWP ratio. This was in agreement with the
random field were generated for each case. A typical series of re­
findings of Popescu et al. [39].
alizations for different value of δ with COVG = 0.3 are shown in Fig. 3.
Fig. 7 shows the time history curve of EPWP ratio with different
COVG. Similar to the results shown in Fig. 5(a), the influence of COVG on
3. Results and discussion
the EPWP ratio was negligible in the Layer A (Fig. 7(a)). It was specu­
lated that the shallow soil layers above the tunnel were closer to the free
The effect of COVG and δ on the excess pore water pressure ratio, the
water table, with more conducive to drainage and a quicker dissipation
liquefied zone, the ground displacement, as well as the uplift displace­
of pore pressure, which leads to a marginal influence of soil spatial
ment of tunnel, were investigated in order to gain a comprehensive
variability. In Layers B and C, there was a negative correlation between
insight into the dynamic behavior of the liquefiable soil-tunnel structure
the EPWP ratio and COVG. The behaviors could be explained as follows.
system.
The spatially distributed G became highly variable with the increasing
COVG, producing a larger volume of soil with larger or smaller G. A
3.1. Assessment of the excess pore water pressure (EPWP) ratio smaller G tended to accelerate the dissipation of excess pore water
pressure and further led to a reduction in the EPWP ratio. It is worth
In this section, we selected two representative points and three noting that the green curve in Fig. 7(c) corresponding with the case of
representative soil layers to estimate the stochastic response of EPWP

Fig. 5. Influence of COVG on the EPWP ratio at (a) Point A and (b) Point B.

Fig. 6. Influence of δ on the EPWP ratio at (a) Point A and (b) Point B.

4
Y. Wang et al. Soil Dynamics and Earthquake Engineering 162 (2022) 107489

Fig. 7. Variations of EPWP ratio under different COVG at: (a) Layer A, (b) Layer Fig. 8. Variations of EPWP ratio under different δ at: (a) Layer A, (b) Layer B
B and (c) Layer C. and (c) Layer C.

5
Y. Wang et al. Soil Dynamics and Earthquake Engineering 162 (2022) 107489

Fig. 9. Variation of A80 with time history under different COVG.

Fig. 10. Cumulative distribution function (CDF) of A80 under different (a) COVG and (b) δ.

COVG = 0.1 and δ = 5 m was higher than the deterministic result. The
Table 2
result highlighted that attention should be paid to the combination of
Summary of the tunnel uplift displacement from Monte Carlo simulations.
relatively small COVG and δ, which might result in a higher EPWP ratio.
The time history of EPWP ratio for the three soil layers is illustrated COVG 0.1 0.3 0.5 0.3
in Fig. 8. Similar to the results in Fig. 6, nearly idnetical stochastic re­ δ (m) 5 5 50 100
sponses were observed for soil layers in the vicinity of the tunnel and at Mean of d (m) 24.67 23.67 23.49 23.67 24.04 24.59
other depths. Standard deviation of 2.88 6.07 10.52 6.07 4.44 3.34
d (m)
Deterministic case (m) 26.05
3.2. Assessment of the liquefied zone

Referring to the research of Popescu et al. (2006), the index of liq­ addition, it can be seen that the influence of δ on the liquefied zone was
uefied zone was defined as not as significant as that by COVG.
( ) Fig. 10 provides the cumulative distribution functions (CDFs) of A80
A u(t)
′ > 0.8 varied with COVG and δ. The cumulative distributions of A80 from three
σv0
A80 (t) = (12) typical COVG showed that all the stochastic responses of A80 were
A
smaller than the deterministic result for COVG = 0.5. Contrast with the
where u(t) is the EPWP at a time instant t after the seismic loading is case for COVG = 0.3 and 0.1, 97.3% and 77.4% of the stochastic results
applied, σv0 is the initial effective stress of a spatial position,
′
were less than the deterministic result (Fig. 10(a)). On the basis of
( ) Figs. 10(b), 2.8% of the stochastic responses of A80 exceeded the
A u(t)
′ > 0.8 and A represent the area in which EPWP ratio is greater
σv0 deterministic result with δ = 5 m and this percentage increased to 7.9%
than 0.8 and the total analysis domain, respectively. and 10.4% with δ = 50 m and 100 m, respectively. It can be claimed that
Fig. 9 presents the time history curve of the index of liquefied zone a smaller COVG along with a greater δ contributed to a larger area of
under different COVG and δ. The results derived from the stochastic liquefaction based on the stochastic analyses.
analysis presents a similar trend, but smaller than those obtained by the
deterministic analysis. Hence, the liquefied area could be overestimated
without considering the spatial variability of soil properties. With the 3.3. Assessment of the tunnel uplift displacement (Up)
increase in COVG, the A80(t) decreased accordingly. It was mainly due to
that a higher COVG enhanced the dispersion of soil shear modulus and Under the seismic loading, one of the most severe problems associ­
increased the occurrence possibility of soil with small shear modulus. ated with the underground structures, such as tunnel, in liquefiable soils
Fig. 9(b) shows the effect of scale of fluctuation on A80(t). There was a is the floatation response due to the soil liquefaction. The mean and
negative correlation between A80 and δ with the time history. Because standard deviation of the uplift displacement of tunnel at the end of
the soil profile became more uniform with a larger δ, the stochastic analysis for all the cases are tabulated in Table 2. The mean of
response of A80(t) therefore tended to the deterministic result. In d decreased with the increase of COVG for a given δ, while the mean of

6
Y. Wang et al. Soil Dynamics and Earthquake Engineering 162 (2022) 107489

Fig. 11. Cumulative distribution function (CDF) of Up under different (a) COVG and (b) δ.

stochastic behaviors.
Table 3 Fig. 12 plotted the variations of the ground displacements from the
Summary of the ground displacements from Monte Carlo simulations.
tunnel centerline at each side. On the basis of the figure, the stochastic
COVG 0.1 0.3 0.5 0.3 responses of the heave displacement were similar, and the stochastic
δ (m) 5 5 50 100 heave displacements were smaller than the deterministic result within
Mean of d (m) 19.46 18.95 18.74 18.95 19.44 19.69
the range of about double the tunnel dimension (i.e., ±12 m away from
Standard deviation of 2.17 6.02 8.01 3.94 3.93 3.01 the centerline). It is reasonable that the presence of tunnel diminished
d (m) the spatial variability of soil modulus in the vicinity of the tunnel due to
Deterministic case (m) 21.01 the local averaging effect. Nevertheless, more than ±12 m away from
the centerline, an increased COVG was responsible with a larger settle­
ment, as displayed in Fig. 12(a). This phenomenon can be explained by
d was positively correlated with δ under a fixed COVG.
that the impact of tunnel-soil interaction on the soil variability became
The cumulative distribution functions (CDFs) of Up with different
gradually weakened, hence the COVG retrieved to affect the ground
COVG and δ is demonstrated in Fig. 11. The distributions of CDF under a
displacement. Fig. 12(b) indicated that the ground displacements had
smaller δ and larger COVG were wider. In Fig. 11(a), at COVG = 0.1, the
correlation with the scale of fluctuation, except that the slight differ­
distribution range of Up was from 18 to 40 cm, in which 84.8% of the
ences in tunnel settlements from − 25 to − 15 m.
stochastic results were less than the deterministic result; at COVG = 0.3
and 0.5, there were only 10.2% and 12% of the stochastic results are
greater than the deterministic result, respectively. As shown in Fig. 11 4. Conclusion
(b), there were 80.5%, 89.0% and 90.1% of the stochastic results less
than the deterministic result when δ = 5, 50 and 100 m, respectively. In this study, the effects of coefficient of variation and scale of
Moreover, the curves in Fig. 11(b) were quite consistent for different δ, fluctuation of soil shear modulus on the dynamic reliability of liquefi­
indicating that the influence of δ on the tunnel uplift displacement was able soil and the dynamic response of rectangular tunnel were investi­
marginal. gated. For this purpose, different parameters have been studied and
assessed. The main conclusions can be summarized as follows:

3.4. Assessment of the ground displacement (d) (1). The excess pore water pressure ratio was insensitive to the scale
of fluctuation in the analyzed soil domain, while a corresponding
The uplift of the underground structure also interacts with the soil, increase was observed with the coefficient of variation for the soil
resulting in vertical deformation such as uplift or settlement at different above the tunnel. However, the influence of coefficient of varia­
locations of the ground surface. The mean and standard deviation of d at tion became less pronounced in a deeper soil layer.
the centerline of tunnel is summarized in Table 3. Compared to the re­ (2). Both the coefficient of variation and the scale of fluctuation of soil
sults in Table 2, the ground displacements at the centerline of tunnel shear modulus significantly affected the liquefied zone. With a
were less than the tunnel uplift displacements, but without difference in larger coefficient of variation, the maixmum liquefied zone is

Fig. 12. Variations of the ground displacement.

7
Y. Wang et al. Soil Dynamics and Earthquake Engineering 162 (2022) 107489

smaller and the decreasing rate of liquefied zone is slower in the [9] Yu HT, Yan X, Bobet A, Yuan Y, Xu GP, Su QK. Multi-point shaking table test of a
long tunnel subjected to non-uniform seismic loadings. Bull Earthq Eng 2018;16
post-earthquake period. With the increase of the scale of fluctu­
(2):1041–59.
ation, the liquefied zone became larger and approached to the [10] Cao ZJ, Wang Y, Li DQ. Quantification of prior knowledge in geotechnical site
deterministic result. characterization. Eng Geol 2016;203:107–16.
(3). Due to the interaction between the tunnel and the surrounding [11] Griffiths DV, Fenton GA, Manoharan N. Bearing capacity of rough rigid strip
footing on cohesive soil: probabilistic study. J Geotech Geoenviron Eng 2002;128
soils, the influence of scale of fluctuation was not as significant as (9):743–55.
the coefficient of variation. However, the increase of scale of [12] Popescu R, Deodatis G, Nobahar A. Effects of random heterogeneity of soil
fluctuation still showed the same influence trend as the decrease properties on bearing capacity. Probabilist Eng Mech 2005;20(4):324–41.
[13] Srivastava A, Sivakumar Babu GL. Effect of soil variability on the bearing capacity
of coefficient of variation. of clay and in slope stability problems. Eng Geol 2009;108(1–2):142–52.
(4). The spatial variability of soil shear modulus showed no influence [14] Li JH, Tian Y, Cassidy M. Failure mechanism and bearing capacity of footings
within the range of double the tunnel dimension in terms of the buried at various depths in spatially random soil. J Geotech Geoenviron Eng 2015;
141(2):04014099.
ground displacement. As it exceeded double the tunnel dimen­ [15] Li L, Li JH, Huang JS, Liu HJ, Cassidy MJ. The bearing capacity of spudcan
sion, the increased coefficient of variation could decrease the foundations under combined loading in spatially variable soils. Eng Geol 2017;227:
ground settlement. 139–48.
[16] Brzakala W, Puła W. A probabilistic analysis of foundation settlements. Comput
Geotech 1996;18(4):291–309.
Author statement [17] Al-Bittar T, Soubra AH. Probabilistic analysis of strip footings resting on spatially
varying soils and subjected to vertical or inclined loads. J Geotech Geoenviron Eng
2014;140(4):04013043.
Yubing Wang: Conceptutation; Funding acquisition; Writing- [18] Jiang SH, Huang J. Efficient slope reliability analysis at low-probability levels in
original draft. spatially variable soils. Comput Geotech 2016;75:18–27.
Jianjian He: Formal analysis; Investigation; Software; Writing- [19] Jiang SH, Huang J, Yao C, Yang JH. Quantitative risk assessment of slope failure in
2-D spatially variable soils by limit equilibrium method. Appl Math Model 2017;
original draft. 47:710–25.
Shuang Shu: Data curation; Methodology; Writing-original draft; [20] Javankhoshdel S, Bathurst RJ. Influence of cross correlation between soil
Writing-review & editing. parameters on probability of failure of simple cohesive and c-φ slopes. Can Geotech
J 2015;53(5):839–53.
Houle Zhang: Formal analysis; Software; Visualization.
[21] Zhang J, Wang H, Huang HW, Chen LH. System reliability analysis of soil slopes
Yongxin Wu: Supervision; Validation; Writing-review & editing. stabilized with piles. Eng Geol 2017;229(7):45–52.
[22] Jiang SH, Papaioannou I, Straub D. Bayesian updating of slope reliability in
spatially variable soils with in-situ measurements. Eng Geol 2018;239:310–20.
Declaration of competing interest [23] Jiang SH, Huang J, Qi XH, Zhou CB. Efficient probabilistic back analysis of
spatially varying soil parameters for slope reliability assessment. Eng Geol 2020;
271:105597.
The authors declare that they have no known competing financial [24] Song KI, Cho GC, Lee SW. Effects of spatially variable weathered rock properties on
interests or personal relationships that could have appeared to influence tunnel behavior. Probabilist Eng Mech 2011;26(3):413–26.
the work reported in this paper. [25] Mollon G, Phoon KK, Dias D, Soubra AH. Validation of a new 2D failure mechanism
for the stability analysis of a pressurized tunnel face in a spatially varying sand.
J Eng Mech 2010;137(1):8–21.
Data availability [26] Ali A, Lyamin AV, Huang JS, Sloan SW, Cassidy MJ. Undrained stability of a single
circular tunnel in spatially variable soil subjected to surcharge loading. Comput
Geotech 2017;84:16–27.
Data will be made available on request.
[27] Huang HW, Xiao L, Zhang DM, Zhang J. Influence of spatial variability of soil
young’s modulus on tunnel convergence in soft soils. Eng Geol 2017;228:357–70.
Acknowledgments [28] Chen FY, Wang L, Zhang WG. Reliability assessment on stability of tunnelling
perpendicularly beneath an existing tunnel considering spatial variabilities of rock
mass properties. Tunn Undergr Space Technol 2019;88:276–89.
The authors acknowledge the support of the National Natural Sci­ [29] Yue QX, Ang HS. Nonlinear response and reliability analysis of tunnels under
ence Foundation of China under Grant No. 52078186. strong earthquakes. Struct Infrastruct E 2016;12(5):618–30.
[30] Yue QX, Ang HS. 3D reliability evaluation of tunnels under strong-motion
earthquakes considering spatial randomness. Struct Infrastruct E 2017;13(7):
References 882–93.
[31] Martin GR, Finn W, Seed HB. Fundamentals of liquefaction under cyclic loading.
[1] Lee KM, Hou XY, Ge XW, Tang Y. An analytical solution for a jointed shield-driven J geotech engrg div 1975;101(5):423–38.
tunnel lining. Int J Numer Anal Met 2001;25:365–90. [32] Azadi M, Hosseini SMMM. Analyses of the effect of seismic behavior of shallow
[2] Yu HT, Cai C, Bobet A, Zhao X, Yuan Y. Analytical solution for longitudinal bending tunnels in liquefiable grounds. Tunn Undergr Space Technol 2010;25(5):543–52.
stiffness of shield tunnels. Tunn Undergr Space Technol 2019;83:27–34. [33] Zheng G, Yang PB, Zhou HZ, Zeng CF, Yang XY, He XP, Yu XX. Evaluation of the
[3] Patil M, Choudhury D, Ranjith PG, Zhao J. Behavior of shallow tunnel in soft soil earthquake induced uplift displacement of tunnels using multivariate adaptive
under seismic conditions. Tunn Undergr Space Technol 2018;82:30–8. regression splines. Comput Geotech 2019;113:103099.
[4] Bao XH, Xia ZF, Ye GL, Fu YB, Su D. Numerical analysis on the seismic behavior of [34] Koutsourelakis S, Prevost JH, Deodatis G. Risk assessment of an interesting
a large metro subway tunnel in liquefiable ground. Tunn Undergr Space Technol structure-soil system due to liquefaction. Earthq Eng Struct Dynam 2002;31:
2017;66:91–106. 851–79.
[5] Miao Y, Yao EL, Ruan B, Zhuang HY, Chen GX, Long XH. Improved hilbert spectral [35] Griffiths DV, Fenton GA. Probabilistic slope stability analysis by finite elements.
representation method and its application to seismic analysis of shield tunnel J Geotech Geoenviron Eng 2004;130(5):507–18.
subjected to spatially correlated ground motions. Soil Dynam Earthq Eng 2018; [36] Phoon KK, Huang SP, Quek ST. Implementation of Karhunen-Loeve expansion for
111:119–30. simulation using a wavelet-Galerkin scheme. Probabilist Eng Mech 2002;17(3):
[6] Wang GB, Yuan MZ, Miao Y, Wu J, Wang YX. Experimental study on seismic 293–303.
response of underground tunnel-soil-surface structure interaction system. Tunn [37] Jiang SH, Huang J. Efficient slope reliability analysis at low-probability levels in
Undergr Space Technol 2018;76:145–59. spatially variable soils. Comput Geotech 2016;75:18–27.
[7] Yan X, Yuan JY, Yu HT, Bobet A, Yuan Y. Multi-point shaking table test design for [38] Shu S, Gao YF, Wu YX. Probabilistic bearing capacity analysis of spudcan
long tunnels under non-uniform seismic loading. Tunn Undergr Space Technol foundation in soil with linearly increasing mean undrained shear strength. Ocean
2016;59:114–26. Eng 2020;204:106800.
[8] Yu HT, Yuan Y, Xu GP, Su QK, Yan X, Li C. Multi-point shaking table test for long [39] Popescu R, Prevost JH, Deodatis G, Chakrabortty P. Dynamics of nonlinear porous
tunnels subjected to non-uniform seismic loadings - part ii: application to the HZM media with applications to soil liquefaction. Soil Dynam Earthq Eng 2006;26(6–7):
immersed tunnel. Soil Dynam Earthq Eng 2016;108:187–95. 648–65.