sustainability

Article
Deformation-Based Basal Heave Reliability Analysis and
Selection on Monitoring Points for General Braced Excavations
Kaiqi Meng 1 , Guangming Yu 1 , Liang Li 1, *, Zhen Xu 2 , Jun Lei 3 , Yanxiang Fan 4 , Hongbiao Yu 5 and Liang Xu 1

1 School of Civil Engineering, Qingdao University of Technology, Qingdao 266033, China

2 Qingdao Metro Line 1 Co., Ltd., Qingdao 266045, China
3 Engineering Bureau Co., Ltd., China State Construction, Changsha 410007, China
4 The Fifth Engineering Co., Ltd., China Railway 25th Bureau Group, Qingdao 266101, China
5 Tunnel Engineering Co., Ltd., China Communications Construction, Beijing 100024, China
* Correspondence: [email protected]

Abstract: A framework for evaluating deformation-based basal heave stability is proposed in order
to distinguish between the different responses under freely developed and prohibited basal heave
failures. In the case of freely developed basal heave failure, the maximum deformation values occur
at the center point of pit bottom, whereas this is not the case for the prohibited basal heave failure.
The critical thickness of soft soil layer between the end of supporting structures and the top of hard
stratum is about 0.3B (B = excavation width), beyond which the freely developed basal heave failure
arises. In situations otherwise, the prohibited basal heave failure occurs. The failure probability of
basal heave failure at the center point increases significantly as B ranges within a limited value; then,
it begins to decrease or to vary slightly at a certain value under a given thickness of soft soil layer. If
the thickness of soft soil layer is so sufficiently large that freely developed basal heave failure occurs
for any of B, the failure probability of basal heave failure at the center point increases as B increases.
The selection of the optimum monitoring points for basal heave stability is recommended to account
for the weights in the contribution to the basal heave deformations of the influencing factors such
as excavation width and thickness of soft soil layer. The proposed framework is applicable to basal
Citation: Meng, K.; Yu, G.; Li, L.; Xu,
heave reliability analysis for braced excavations where deformation values are focused.
Z.; Lei, J.; Fan, Y.; Yu, H.; Xu, L.
Deformation-Based Basal Heave Keywords: basal heave stability; basal heave deformation; monitoring points design; response
Reliability Analysis and Selection on surface method; Monte Carlo simulation
Monitoring Points for General Braced
Excavations. Sustainability 2023, 15,
8985. https://doi.org/10.3390/
su15118985 1. Introduction
Academic Editor: Maxim Basal heave stability has seen witnessed its rise importance within the large number of
A. Dulebenets braced excavations [1–4]. Since the failure to properly address the basal heave stability has
great potential to induce pit failures, basal heave stability has attracted extensive attention
Received: 12 May 2023
from worldwide researchers. Traditionally, basal heave stability is evaluated via factor of
Revised: 28 May 2023
safety against basal heave. The methods for calculating the factor of safety against basal
Accepted: 30 May 2023
heave are classified into three categories, i.e., (1) the Terzaghi method based on bearing
Published: 2 June 2023
capacity [5–7], (2) the slip circle method [8–10], and (3) the numerical methods combining
strength reduction techniques [11]. Since the former two methods cannot properly take
the effect of the supporting structures into account, the numerical method with strength
Copyright: © 2023 by the authors. reduction technique is widely adopted owing to its advantage in the consideration of
Licensee MDPI, Basel, Switzerland. complicated braces and construction stages [12–16]. In numerical methods, the factor of
This article is an open access article safety against basal heave is similar to that used in slope stability analysis. It is a factor by
distributed under the terms and which the shear strength is reduced to bring the pit system to a limit equilibrium state. It has
conditions of the Creative Commons been found that, before the pit system reaches the limit equilibrium state, fairly substantial
Attribution (CC BY) license (https:// deformation occurs at the locations adjacent to the pit walls. This results in the disfunction
creativecommons.org/licenses/by/ of pit from the perspective of serviceability limit state. Therefore, the factor of safety by
4.0/).

Sustainability 2023, 15, 8985. https://doi.org/10.3390/su15118985 https://www.mdpi.com/journal/sustainability

Sustainability 2023, 15, 8985 2 of 19

numerical methods overestimates the basal heave stability and the deformation-based
index may be relevant.
The basal heave deformation is a visual representation that can be observed within
the construction of a pit. The effects of excavation width, depth, and support structures on
basal heave deformation were investigated and summarized in [17,18]. Previous studies
have demonstrated the deformation rules. Considering the fact that deformation values
are readily limited to a specific value in terms of serviceability limit state design, the
deformation values at characteristic points provide much insight into basal heave stability
and they are also the important monitoring items within the construction of a pit. To the
knowledge of the authors, no research has been reported regarding the deformation-based
basal heave stability.
In addition, the deterministic basal heave stability analysis has to shift towards a
probabilistic approach where the variability of soil parameters can be properly taken into
account. Following the probabilistic approach, the reliability of basal heave stability is
evaluated by analyzing the distribution of factors of safety against basal heave under
different soil parameters or via the calculation of the reliability index using simplified
equations [5,19]. To enhance the computational efficiency when the Monte Carlo simulation
is adopted, machine learning-aided surrogate models or automatically driven FEM models
have been used in basal heave reliability analysis [20–23]. The previous studies have
facilitated the deformation-based basal heave reliability analysis.
This paper aims to develop a methodology with which to conduct basal heave reliabil-
ity in terms of deformation values and finally to select optimum monitoring points in order
to dominate the response of a braced excavated pit for sustainable underground space de-
velopment. The paper starts with the definition of limit state function in deformation-based
basal heave reliability analysis, followed by the calibration of response surface function
to predict the deformation values at characteristic points. Then, the proposed methodol-
ogy is fully described in Section 4. The proposed method is illustrated through a braced
excavation example where the influences of excavation width and the thickness of the soft
layer on the distributions of deformation values at characteristic points are investigated.
The extent to which the deformation value at the midpoint can represent the response of a
pit system is also discussed.

2. Deformation-Based Basal Heave Failures

The braced excavations are mainly conducted in the urban area with populated high-
rise buildings. The deformations induced by the excavations have significant influence
on the serviceability state of the surrounding structures. For this reason, the basal heave
deformation is adopted to construct the limit state function of basal heave failure. The
determination methods of basal heave deformation mainly consist of a summation method,
a residual stress method, a self-weight stress offset method, and so on. However, the
above methods cannot reasonably consider the influence of factors [24] such as pit shape,
size, groundwater and support structure on the basal heave deformation. Therefore, the
numerical analysis method [25,26], which properly addresses those factors, has been widely
used to model the construction stages in a pit system and to obtain the deformation value
of specific locations.

2.1. Hardening Soil Model

The hardening soil model of Schanz et al. [27] is adopted in the current study and,
therefore, its fundamentals are briefly reviewed herein for favorable readability. The
hardening soil model uses stress-dependent stiffnesses to calculate the elastic strains and to
calculate the plastic strains via multiple-surface yield criterion. In essence, the hardening
soil model is the combination of Duncan–Chang hyperbolic model and the isotropic strain
hardening model. The distinguished characteristics of the hardening soil model is the use of
different moduli for different loading conditions such as virgin loading or un-/reloading. In
total, three stress-dependent moduli are present in the hardening soil model: Eoed (tangent
Sustainability 2023, 15, 8985 3 of 19

stiffness modulus for primary loading), E50 (secant stiffness modulus for primary loading),
and Eur (stiffness modulus for unloading and reloading). To facilitate the determination of
the stress-dependent modulus, three respective reference moduli are used Eoed ref (reference

tangent stiffness modulus for primary oedometer loading), E50 ref (reference secant stiffness
ref
modulus for standard drained triaxial test), and Eur (the reference unloading and reloading
stiffness modulus). The stress-dependent modulus is related to the reference modulus
using the following equations:
σ3 !m
ref
ccosϕ − k0 sinϕ
Eoed = Eoed (1)
ccosϕ + pref sinϕ

σ3 !m
ref
ccosϕ − k0 sinϕ
Eur(or 50) = Eur (or 50) (2)
ccosϕ + pref sinϕ
where σ3 is the initial effective horizontal stress on site. Because tension is assumed to be
positive and compression is negative in the PLAXIS 2D CE V20 Update 4 software package,
a minus is added in Equations (1) and (2). k0 is the coefficient of lateral pressure, σ3/k0 is ef-
fective vertical stress. pref is the reference stress. c and ϕ are effective cohesion and effective
internal friction angle. Previous studies [28–32] have highlighted the empirical relation-
ships among E50 ref , Eref , and Eref . Eref
ur oed (1∼3),50
= 1.13Eref ref ref ref
(1∼3),oed , E4,50 = E4,oed , E(1∼2),ur =
9.13E(ref ref = 4Eref , Eref = 30Eref
, E3,ur
1∼2),oed 3,oed 4,ur 4,oed are adopted in this study, where the single
subscript or that before the coma is the number of the soil layer. For more details about the
hardening soil model, the readers are referred to [27] and the user manual for PLAXIS.

2.2. Limit State Function in Terms of Basal Heave Deformation

From the perspective of serviceability limit design, the deformation values at character-
istic points are confined to specific thresholds. The locations of characteristic points can be
derived from the specifications in the Chinese Code for Acceptance of Construction Quality
of Building Foundation (GB-50202-2001) [33]. For the specific thresholds, researchers can
refer to the Chinese Technical Code for Monitoring of Building Foundation Excavation
Engineering (GB50497-2009) [34]. Therefore, the limit state function of basal heave failure
is defined as:
G (X) = dmax −dc (X) = 0 (3)
where X denotes the vector of random variables considered in basal heave reliability
analysis, such as cohesion, internal friction angle, and unloading–reloading modulus of
soil, etc. dmax is the threshold of deformation value. G(X) is the limit state function of basal
heave failure, G(X) > 0 indicates that the basal heave failure does not occur, and vice versa.
dc (X) is the predicted deformation value obtained using finite element simulation or by
using response surface function with sufficient accuracy.

2.3. Identification of Basal Heave Failure Using MCS

Once the limit state function of basal heave failure is available, the Monte Carlo
simulation (MCS) is commonly used in reliability analysis [35]. The basic idea of MCS is
to emulate the system response (e.g., the basal heave deformation at characteristic point)
through repeated executions of finite element simulations under random samples. Consider
for example, that a finite number (e.g., N) of random samples (one sample corresponds to
one vector of (X)) are generated, and the limit state function value under each of N random
samples is calculated by Equation (3), leading to N deformation values as dc (Xi ), i = 1,
2, 3, . . . , N. The sample resulting in a negative G(X) value is defined as a failure sample.
The ratio of counted number of failure samples to N approximates the basal heave failure
probability Pf as:
m
Pf ≈ (4)
N
Sustainability 2023, 15, 8985 4 of 19

in which m is the number of failure samples, and N is the total number of samples. As
N approaches a sufficiently large value, the simulated Pf is equal to the ‘true’ value.
Traditionally, the choice of N is dependent on the accuracy of simulated Pf . Let δPf denote
the coefficient of variation in simulated Pf , and the following equation relates the target Pf ,
N and δPf [36]:
v 
u
u 1 − Pf
δP f = (5)
t
NP f
where δPf is the coefficient of variation of failure probability, Pf is the target failure proba-
bility, and N is the total number of samples. A sensitivity study of N on Pf is recommended
to select a proper N before the MCS is applied.

3. The Prediction of Basal Heave Deformation Value Using Response Surface Method
Although the wide applications of MCS in reliability analysis have demonstrated its
advantages of simplicity and easy implementation, the disadvantage of MCS is its com-
putational inefficiency, especially for the cases where repeated finite element simulations
(e.g., PLAXIS 2D simulation in this study) are unavoidable. To address this issue, a simple
surrogate model using a response surface function is calibrated to model the relationship
between X and dc . If the calibrated response surface function has a fairly good accuracy, it
is adopted instead of finite element simulation to obtain dc , thereby reducing the computa-
tional cost. The second-order polynomial response surface function without cross items has
been widely used in the reliability analysis of geotechnical engineering problems [37–41].
It is adopted to establish the relationship between the basal heave deformation values
and soil parameters in this study. The approximate basal heave deformation value can be
expressed as
n n
drc (X)= a+ ∑i=1 bi xi + ∑i=1 di xi2 (6)
where drc (X) is the predicted basal heave deformation value by response surface function,
X = (x1 , x2 , . . . , xn ) is the vector of n random variables, and a, bi , and di (i = 1, 2, . . . , n) are
constants to be determined. To determine the 2n + 1 coefficients, 2n + 1 experimental points
should be selected using the center point composite design method. In the set of mean
values µxi , the center point is taken as the first experimental point, i.e., X1 = (µ x1 , µ x2 , . . . ,
µxn ). For the ith random variable xi , it takes the value of µxi ± ωσxi , respectively, while the
random variable reset o all their remaining mean values to form two separate experiment
points. σxi is the standard deviation of xi , and ω is a coefficient for generating experiment
points. The effect of ω on the accuracy of the response surface function is insignificant and
ω = 1 is usually selected [42]. The generated 2n + 1 experiment points X1 , X2 , . . . , X2n +1 are
input into the finite element model in PLAXIS 2D to obtain the corresponding basal heave
deformation values (dc (X1 ), dc (X2 ), . . . , dc (X2n+1 )). By equating dc (Xi ) with drc (Xi ), 2n + 1
linear equations regarding coefficients (a, bi , di ), i = 1, 2, . . . , n can be obtained. By solving
these linear equations, 2n + 1 coefficients can be determined and hence the response surface
function in Equation (6) is constructed. It should be noted that, before the response surface
function can be used to make predictions, the accuracy of the response surface must be
verified. The necessary explanations for this will be given later in this paper.

4. Basal Heave Stability and Reliability Analysis Based on Deformation Analysis

The proposed methodology for deformation-based basal heave reliability analysis is
illustrated in Figure 1. It consists of two parts, i.e., part I (deterministic analysis using com-
mercial software) and part II (reliability analysis using validated response surface functions
and MCS). The function of part I is to depict the deformation values and principles under a
number of excavation plans and to provide benchmark values for the experimental points
in the calibration of response surface functions. Part II aims to identify the failure samples
in accordance with the specified threshold via MCS and to count the number of failure
samples. The braced excavation models with different excavation and reinforcing plans are
 The proposed methodology for deformation-based basal heave reliability analysis is
illustrated in Figure 1. It consists of two parts, i.e., part I (deterministic analysis using
commercial software) and part II (reliability analysis using validated response surface
functions and MCS). The function of part I is to depict the deformation values and princi-
ples under a number of excavation plans and to provide benchmark values for the exper-
Sustainability 2023, 15, 8985 imental points in the calibration of response surface functions. Part II aims to identify 5 ofthe
19

failure samples in accordance with the specified threshold via MCS and to count the num-
ber of failure samples. The braced excavation models with different excavation and rein-
forcing
firstly plans are
designed to firstly designed
facilitate to facilitate the
the implementation of implementation of part
part I. The statistics I. The
of soil statistics
parameters
involved in the excavations
of soil parameters involvedare in determined
the excavationsaccording to the engineering
are determined according geological survey
to the engineer-
reports and the survey
ing geological relevantreports
literature
andtotheexecute the second
relevant literaturepart. Part I is conducted
to execute the secondby PLAXIS
part. Part I
2D CE V20 Update
is conducted 4, as shown
by PLAXIS 2D CEinV20 the Update
left half 4,
of as
Figure
shown1. This
in theincludes
left halfthe
of establishment
Figure 1. This
of FEM model,
includes the selection of
the establishment ofthe
FEMcharacteristic
model, thepoints, the of
selection execution of the model,
the characteristic and the
points, the
use
execution of the model, and the use of outputs for the deformation values. Part IIdesign
of outputs for the deformation values. Part II starts with the experimental points starts
for
withthethe
calibration
experimentalof thepoints
response surface
design functions,
for the followed
calibration of theby the repeated
response surfaceexecution
functions,of
deterministic
followed by the analysis
repeatedto obtain the benchmark
execution deformation
of deterministic analysis values at experimental
to obtain the benchmark points.
de-
The response
formation surface
values functions are calibrated
at experimental points. The and validated
response usingfunctions
surface self-developed code in
are calibrated
Matlab. The final step of part II is to conduct MCS to obtain the number
and validated using self-developed code in Matlab. The final step of part II is to conduct of failure samples
and
MCS totocalculate
obtain thethenumber
failure probability.
of failure samples and to calculate the failure probability.

Designing different working conditions with

Determining the statistical values
different excavation widths and the thicknesses of
of the mechanical parameters of
the weak soil layer and determining the support
the foundation soil layer
design as well as construction plan

PLAXIS 2D CE V20 Update 4 Input PLAXIS 2D CE V20 Update 4 Output

For modeling and performing For viewing and outputting
computational simulations calculation results

Designing calculation software

input parameters based on
response surface method

Selecting characteristic points

Establishing Plaxis2D
(The red points can all be selected
numerical models (B=3H)
to output their deformations)
Repeating the deterministic analysis process to
obtain the basal heave deformation value of
each characteristic point

For building response

Running parameters and simulations surfaces and performing
related calculations

Outputting Counting the basal heave

deformation diagram deformation value of each Matlab code
For building response surfaces Matlab code
characteristic point
For calculating the probability of failure

Analyzing the basal Performing reliability

heave deformation law analysis

Part Ⅰ Deterministic analysis Part Ⅱ Reliability analysis

Figure1.1.Flow
Figure Flowchart
chartfor
forbasal
basalheave
heavestability
stabilityand
andreliability
reliabilityanalysis.
analysis.

5. Example
5.1. Finite Element Modeling
Consider a schematic model of braced excavation shown in Figure 2. The soil lay-
ers uncovered from top to bottom comprise fill layer (0~−3 m), powdered clay layer
(−3~−15 m), powdered layer (−15~−37.2 m) and moderately weathered muddy silt layer
(−37.2~−50 m). Let B, H, D, T, and L represent the excavation width, excavation depth,
penetration depth of support structure, vertical distance between pit bottom and top of
hard stratum, and horizontal distance between the support structure edge to the model
edge, respectively. A width coefficient k is introduced to study the effect of B on the basal
heave deformations:
B = kH (7)
where B is the excavation width, H is the excavation depth, and k = 0.5, 0.8, 1, 2, . . . , 10 in
this study. The braced works consist of 16 m deep and 0.35 m thick diaphragm walls and
two horizontal struts, installed at −2 m, and −6 m, respectively. The strut stiffness, denoted
 tively. A width coefficient k is introduced to study the effect of B on the basal heave de-
formations:

B = kH (7)

Sustainability 2023, 15, 8985 where B is the excavation width, H is the excavation depth, and k = 0.5, 0.8, 1, 2, …10 6 of in
19
this study. The braced works consist of 16 m deep and 0.35 m thick diaphragm walls and
two horizontal struts, installed at −2 m, and −6 m, respectively. The strut stiffness, denoted
6 6
by EA, is
by EA, is 4.8
4.8 ×
× 10106kN ⁄m. The
kN/m. The EA
EA and
and EI of diaphragm
EI of diaphragm walls
walls are
are1212×
× 10
106 kN ⁄m and
kN/m and
3 3 2 2
120
120 × 10 kN·m , respectively. The groundwater table is at a depth of −3 m. Table 11 lists
× 10 kN∙m , respectively. The groundwater table is at a depth of −3 m. Table lists
the
the five construction stages
five construction stages of the braced
of the braced excavations.
excavations. The
The means
means ofof soil
soil parameters
parameters are
are
shown
shown inin Table 2.
Table 2.

Figure 2. The schematic geometry and mesh of the FE model (k == 3).

3).

Table 1. The
Table 1. The construction
construction stages
stages of
of braced
braced excavations.
excavations.

Stage
Stage No.
No. Excavation
Excavation Depth/m
Depth/m Items
Items
0 0.0 Initial ground stress analysis
0 0.0 Initial ground stress analysis
1 1 0.00.0 Diaphragm
Diaphragm wall construction
wall construction
2 2 −3.0
−3.0 Excavation,
Excavation, strut installedatat−−2
strut installed 2 mm
3 3 −7.0
−7.0 Excavation, strut
Excavation, installedatat−−6
strut installed 6 mm
4 −10.0 Excavation
4 −10.0 Excavation

2. Means
Table 2.
Table Means for
for the
the soil
soil parameters.
parameters.

Soil Properties
Soil Properties γγ Eref
50 EEref
ref
oed
oed Eref
ur cc’ 0 φϕ’ 0
Soiltype
Soil type (kN/m
(kN/m3 )3) (kN/m
(kN/m22)) (kN/m
(kN/m22) ) (kN/m
(kN/m22)) (kN/m
(kN/m22)) (°)
(◦ )
Filllayer
Fill layer1 1 17
17 6780
6780 6000
6000 54,780
54,780 55 2222
Powdered
Powdered clay
claylayer
18
18 9040
9040 8000
8000 73,040
73,040 15
15 2020
layer
2 2
Powdered layer
Powdered layer3 3 17
17 11,300
11,300 10,000
10,000 40,000
40,000 11 3131
Moderately
Moderately
Weathered muddy
Weath-
24.8
24.8 700,700
700,700 700,700
700,700 2,100,000
2,100,000 200
200 2828
ered
silt muddy
layer 4/ hardsilt
stratum

5.2. The Determination of Model Boundary and Mesh Density

In this section, the plane strain finite element analysis for the excavation is carried out
using commercial software PLAXIS 2D CE V20 Update 4. The soil is modeled by 15-node
triangular elements and the hardening soil model is selected in PLAXIS 2D to model the
soil behavior [32,43]. The diaphragm walls are modeled by plat elements and the struts
are modeled by anchor elements. The interaction between structure and soil is modeled
by interface elements, which follow the Mohr–Coulomb criterion, and we use a strength
reduction factor (Rinter ) to consider the strength parameters between structure and soil. The
default set Rinter = 0.7 is used in this paper. Based on sensitivity analysis and experience
from previous studies [44], L should not be lower than 2H. Hence, L = 3H is adopted in
order to eliminate the truncated effects of the slope model on the calculation results. The
thickness of the hard stratum is maintained at 8 m or more in order to ensure stable results.
The side boundaries are fixed on the horizontal direction and the bottom boundary Is
constrained both horizontally and vertically.
Sustainability 2023, 15, 8985 7 of 19

The mesh density has a significant effect on the calculation results and the computa-
tional effort [45]. In order to determine the optimal mesh density in a way that reaches a
balance between precision and computational efficiency, FEM models with B = H, 3H and
6H are used to study the influence of mesh density on the calculation results and efforts. In
this section, D = 6 m and T = 13 m remain unchanged. For each of three FEM models, five
mesh densities, ranging between highly rough, rough, medium, fine, and ultra-fine, are
adopted to discretize the model domain. Under each of the five mesh densities, the defor-
mation value at the center point of pit bottom is extracted from the FEM simulation result.
Let the deformation value calculated from the ultra-fine mesh density be the benchmark.
The relative discrepancy, denoted by δ, between each of the five mesh densities and the
ultra-fine one can be determined. Table 3 summarizes the relative discrepancies at different
FEM models. It is noted that the FEM simulation is conducted on a desktop computer
with Intel (R) Core (TM) i7-9700K [email protected] GHz and 16.0 GB RAM. It can be noticed
from Table 3 that deformation values at the center point of the pit bottom under ultra-fine
mesh density are 16.53, 40.70, and 64.00 mm regarding the FEM model with B = H, 3H,
and 6H, respectively. In the case of B = H, taking 16.53 mm as benchmark, the respective
deformation values for highly rough, rough, medium, and fine mesh densities are 17.26,
16.89, 16.62, and 16.57 mm. The calculated δ values are 4.42%, 2.18%, 0.54%, and 0.24%,
respectively.

Table 3. Sensitivity analysis of mesh density.

Excavation Width
B=H B = 3H B = 6H
Mesh Density Calculation Calculation Calculation
δ/% Time δ/% Time δ/% Time
/s /s /s
Highly rough mesh 4.42 14.40 0.49 17.19 0.16 16.90
Rough mesh 2.18 23.94 0.02 23.72 0.05 26.93
medium mesh 0.54 24.28 0.02 27.78 0.05 46.70
Fine mesh 0.24 68.35 0.05 83.04 0.05 72.38
Ultra-fine mesh 134.63 129.72 157.14

The smaller the excavation width of the pit, the greater the influence of the mesh
division on the calculation results. Considering these three pits, the calculation results
obtained by using a fine and medium degree of mesh division are closer to the results
obtained by using ultra-fine mesh division. However, using medium mesh grants a signifi-
cant time advantage over fine and ultra-fine mesh. Therefore, considering the accuracy of
the calculation results as well as the calculation cost, medium mesh was used in this paper.

5.3. Deterministic Analysis

5.3.1. The Influence of B on the Basal Heave Deformation
A total of 12 Bs are selected to investigate the influence of B on the basal heave
deformation. To eliminate the influence of T on the basal heave deformation,
√ the critical
depth Tc proposed in the literature [44] is used in this study. Tc = B/ 2 is the minimum
depth to ensure a freely developed sliding basal heave failure with D = 0 [44]. In this
section, T = Tc + D is adopted.
A series of characteristic points is selected at the pit bottom at which the basal heave
deformation values are extracted from the completed FEM simulation. Let the x axis
represent the horizontal distance of each of the characteristic points from the center point
of the pit bottom and the y axis represent the corresponding heave deformation value.
Figure 3 plots the variation in the basal heave deformation values under different Bs. The
basal heave deformation is symmetrically distributed, with the center point of the pit
bottom. When B is narrow, i.e., B = 0.5H, 0.8H, H, the basal heave deformation at both sides
of the pit is larger, and that at the central part is smaller (Trend I).
 deformation values are extracted from the completed FEM simulation. Let the x axis rep-
resent the horizontal distance of each of the characteristic points from the center point of
the pit bottom and the y axis represent the corresponding heave deformation value. Figure
3 plots the variation in the basal heave deformation values under different Bs. The basal
heave deformation is symmetrically distributed, with the center point of the pit bottom.
Sustainability 2023, 15, 8985 8 of 19
When B is narrow, i.e., B = 0.5H, 0.8H, H, the basal heave deformation at both sides of the
pit is larger, and that at the central part is smaller (Trend I).

Figure 3.
Figure The variations
3. The variations of
of basal
basal heave
heave deformation
deformation values
values under
under freely
freely developed
developed sliding
sliding basal
basal
heave
heave failure.
failure.

When the
When the excavation
excavation width
width is B == 2H~10H,
is wide, i.e., B 2H~10H,the thebasal
basalheave
heave deformation
deformation at at
the central part is larger, and that at the side parts is smaller. It can be observed that as BB
the central part is larger, and that at the side parts is smaller. It can be observed that as
increases, the
increases, the maximum
maximum basal basal heave
heave deformation
deformation value
valuerises
risessignificantly.
significantly. ItIt is
isrecognized
recognized
that the basal heave deformation is mainly attributed to two sources: (1)
that the basal heave deformation is mainly attributed to two sources: (1) soil rebound soil rebound duedueto
thethe
to unloading
unloadingfactor from
factor fromexcavation
excavation (Source I) and
(Source (2) the
I) and deformation
(2) the deformationduedue to the
to inward
the in-
movement
ward of the of
movement support structure
the support towards
structure the excavation
towards surface surface
the excavation (Source (Source
II) [46]. When
II) [46].B
is narrow (smaller than 2H), the basal heave deformation from Source I is relatively small,
and that from Source II is comparatively large, leading to Trend I. When B is wide (greater
than 2H), the value in Source I is relatively large and contributes a greater amount to the
basal heave deformation, yielding Trend II. As B increases, the deformation from Source
I increases considerably. The rationality of using B = 2H in reference [47,48] as a critical
condition for wide and narrow pits is validated through the results of this sensitivity study,
displayed in Figure 3.

5.3.2. Influence of T on the Basal Heave Deformation

In the previous section, the freely developed sliding basal heave failure is assumed.
However, T has significant effect on basal heave deformation based on the research in
reference [49], where it is seen that the factor of safety against basal heave varies significantly
with T. Intuitively speaking, a small T tends to prevent the sliding basal heave from freely
developing, as depicted in the previous section. Define T as:

T = aTc + D (8)
√
where Tc = B/ 2 is the critical depth, a is the coefficient, and a = 0.1, 0.2, 0.3, . . . , 1.5 with
equal increment.
Since the basal heave deformation is symmetrically distributed with the center of the
pit bottom, only left half of the FEM model shown in Figure 2 is considered. Four cases
of B = H, 3H, 6H and 10H are considered. Under a specific B, 15 T values are adopted to
 T = aTc + D (8)

where Tc = B/√2 is the critical depth, a is the coefficient, and a = 0.1, 0.2, 0.3…1.5 with equal
increment.
Sustainability 2023, 15, 8985 Since the basal heave deformation is symmetrically distributed with the center of19the
9 of
pit bottom, only left half of the FEM model shown in Figure 2 is considered. Four cases of
B = H, 3H, 6H and 10H are considered. Under a specific B, 15 T values are adopted to
investigate the effect of T on basal heave deformation. For a given combination of B and
investigate the effect of T on basal heave deformation. For a given combination of B and T,
T, the basal heave deformation values at characteristic points are calculated.
the basal heave deformation values at characteristic points are calculated.
Similar to the x and y axes in Figure 3, Figure 4 plots the variation in basal heave
Similar to the x and y axes in Figure 3, Figure 4 plots the variation in basal heave
deformation values at characteristic nodes under different B and T values. In the case of B
deformation values at characteristic nodes under different B and T values. In the case of
= H, at a specific a = 0.1, as the characteristic points approach to the center point (x = 0), the
B = H, at a specific a = 0.1, as the characteristic points approach to the center point (x = 0),
maximum deformation value appears at x = −3.5. As x increases from −5 to −3.5, the defor-
the maximum deformation value appears at x = −3.5. As x increases from −5 to −3.5,
themation value increases
deformation to the maximum
value increases value, value,
to the maximum and then
anditthen
decreases when when
it decreases x is greater
x is
than −3.5.
greater than This
−3.5.variation trend has
This variation trendbeen
hasobserved for other
been observed 14 as. 14
for other Foras.the specific
For charac-
the specific
teristic pointpoint
characteristic (e.g.,(e.g.,
center point),
center as a increases
point), fromfrom
as a increases 0.1 to0.1
1.5,
toi.e.,
1.5,Ti.e.,
increases fromfrom
T increases 0.71 to
10.61 m, and the deformation value increases in a non-linear rate to an approximate
0.71 to 10.61 m, and the deformation value increases in a non-linear rate to an approximate value
of 18 mm. This variation trend has been noted for other characteristic
value of 18 mm. This variation trend has been noted for other characteristic points. points.

Sustainability 2023, 15, x FOR PEER REVIEW 10 of 20

(a) (b)

(c) (d)
Figure
Figure 4.4. Influence
Influence of
of T
T on
on basal
basal heave
heave deformation
deformation under
under different
different BB and
and T.
T. (a)
(a) B
B== H;
H; (b)
(b) B
B == 3H;
3H;
(c) B = 6H; (d) B = 10H.
(c) B = 6H; (d) B = 10H.

In the case of B == 3H,

3H, under
under aa ==0.1,
0.1,aa similar
similar variation
variation trend
trend is
is found
found to to that
that at
at BB ==H.
H.
For the
the as
asgreater
greaterthan
than0.2,
0.2,the
themaximum
maximumdeformation
deformationvaluevalue appears
appears at at
thethe
center point
center (x
point
=(x0) and
= 0) andthis variation
this variationtrend
trendis is
observed
observed totobebedifferent
differentthan
thanininthe
thecase
caseof ofBB == H.
H. For
For a
given characteristic
characteristic point,
point, the
the deformation
deformation variation
variation trend
trend with
with kk is observed
observed to be same
with that in the case of B = H. H. In
In the
the cases
cases of B ==6H 6Hand
and10H,
10H, similar
similar variation
variation trends
trends have
been noticed as those in the case of B = 3H.
The points
points with
withthe
themaximum
maximum deformation
deformation values are are
values connected in dashed
connected pink lines,
in dashed pink
as seen
lines, asinseen
Figure 4, to clearly
in Figure demonstrate
4, to clearly the variations
demonstrate of the points.
the variations of theAs a ranges
points. As abetween
ranges
between 0.1 and 1.5, the points with maximum deformation values fluctuate within a zone
about 0.25B away from the center point of pit bottom in the case of B = H. At B = 3H, the
points with maximum deformation values appear at the center point of the pit bottom.
This is true except for a = 0.1, which is located about 0.25B away from the center point. For
Sustainability 2023, 15, 8985 10 of 19

0.1 and 1.5, the points with maximum deformation values fluctuate within a zone about
0.25B away from the center point of pit bottom in the case of B = H. At B = 3H, the points
with maximum deformation values appear at the center point of the pit bottom. This is true
except for a = 0.1, which is located about 0.25B away from the center point. For the cases of
6H and 10H, a similar variation trend has been found to that at B = 3H. However, it has a
slight difference in critical a value, beyond which the points with maximum deformation
values are all at the center point of the pit bottom. To fully address the issue of the variation
trend of points with maximum deformation values, further numerical simulations are
performed regarding B = 0.5, 0.8, 2.0, 4.0, 5.0, 7.0, 8.0, 9.0H and a = 0.1, 0.2, 0.3, . . . , 1.5,
with 0.1 equal increments. The distances of characteristic points from the center point of pit
bottom are divided by 0.5B to obtain the normalized distance ranging between 0 to 1.0.
Figure 5 plots the variation of points with maximum deformation values at different
as and Bs. It can be seen that when B/H < 2, the points with maximum deformation values
Sustainability 2023, 15, x FOR PEER REVIEW
are not located at the center point of the pit bottom. When B/H ≥ 2, the points 11 with
of 20
maximum deformation values are all located at the center point of pit bottom once a ≥ 0.4.

Figure 5.
Figure 5. Variation
Variation of
of points
points with
with maximum
maximum deformation
deformation value.
value.

Based on the the aforementioned

aforementioned comparative results, it can be concluded concluded that for the the
narrow pitpit (B/H
(B/H << 2)
2) the
the influence
influence of T does not change the deformation variation trend
with
with different characteristic points.
points. However,
However,for forthe
thewide
widepitpit(B/H
(B/H ≥as2),
≥ 2), a is greater
asgreater
a is than
than 0.4, where
0.4, where the freely
the freely developed
developed basal basal
heaveheave
slidingsliding can occur,
can occur, the maximum
the maximum defor-
deformation
mation value appears
value appears at the point.
at the center center Otherwise,
point. Otherwise, the T prevents
the T prevents the fullthe full development
development of the
of
failure surface of the soil and this results in a deformation trend similar tosimilar
the failure surface of the soil and this results in a deformation trend to that
that observed
observed
from the from thepit.
narrow narrow
It is pit.
worthIt ispointing
worth pointing
out thatout thethat the critical
critical thickness
thickness of softofsoil
softlayer
soil
layer
betweenbetween the
the end of
√ end of supporting
supporting structures
structures and and
the the
top oftop
hardof hard stratum
stratum is aboutis about
0.3B 0.3B
(0.3B
≈ 0.4 ×≈ B/
(0.3B 0.4 × B/ 2).
√2).
5.4. Reliability Analysis
5.4. Reliability Analysis
Following the preliminary sensitivity analysis, E2,ur ref , ϕ0 , Eref and ϕ0 are selected as
Following the preliminary sensitivity analysis, Eref ref
2 ’ 3,50 3 ’
2,ur , φ2 , E3,50 and φ3 are selected
lognormal distributed random variables, where the single subscript or that before the coma
asthe
is lognormal
number distributed random
of the soil layer. Thevariables, wheremeans
statistics (e.g., the single subscriptdeviation)
and standard or that before the
of each
coma is the
random number
variable are of the soil
shown layer.4.
in Table The statistics
A series (e.g., means
of points and standard
with equal deviation)
horizontal of
distance at
each random variable are shown in Table 4. A series of points with equal horizontal dis-
tance at the pit bottom are selected to be the characteristic points on which the defor-
mation values are focused in order to investigate basal heave reliability. Since H is 10 m,
by referring to the documents such as Chinese Technical Code for Monitoring of Building
Foundation Excavation Engineering (GB50497-2009), the threshold value dmax can be ten-
Sustainability 2023, 15, 8985 11 of 19

the pit bottom are selected to be the characteristic points on which the deformation values
are focused in order to investigate basal heave reliability. Since H is 10 m, by referring to
the documents such as Chinese Technical Code for Monitoring of Building Foundation
Excavation Engineering (GB50497-2009), the threshold value dmax can be tentatively stated
to be in the range between 25 and 35 mm.

Table 4. Statistical parameters of random variables.

Eref
2,ur /kN/m
2 ϕ02 /(◦ ) Eref
3,50 /kN/m
2 ϕ03 /(◦ )
Mean value 73,040 20 11,300 31
The coefficients of variation 0.1 0.2 0.085 0.15
Standard deviation 7304 4 960.5 4.65
Distribution type Lognormal distribution

5.4.1. Response Surface Calibrations at Characteristic Points

Following the descriptions in Section 2, a second-order polynomial response surface
function without cross items is calibrated for each characteristic point. The calibrated
response surface function at a specific characteristic point is:

drc = a + b1 E2ur
ref
+b2 ϕ20 +b3 E3,50
ref
+b4 ϕ30 +d1 E2,ur
ref 2
+d2 ϕ20 2 +d3 E3,50
ref 2
+d4 ϕ30 2 (9)

where a, bi , di , i = 1, 2, 3, 4 are coefficients. Table 5 lists the coefficients for response surface
function calibrated at the center point and side point, as shown in Figure 1 for the cases of
B = 3H and B = 6H.

Table 5. The coefficients for calibrated response surface functions.

Center Point Side Point Center Point Side Point

Coefficients
B = 3H B = 3H B = 6H B = 6H
a 26.5248 27.02299 19.88041 21.70444
b1 −7.39 × 10−5 −7.80 × 10−5 7.53 × 10−6 −6.85 × 10−6
b2 2.75 × 10−2 −2.125 × 10−2 4.375 × 10−2 3.25 × 10−2
b3 −2.6 × 10−4 −1.4 × 10−4 1.17 × 10−4 4.04 × 10−16
b4 −2.724 × 10−1 −1.4 × 10−4 −2.3978 × 10−1 −2.7061 × 10−1
d1 3.75 × 10−10 3.75 × 10−10 −9.37 × 10−11 −2.74 × 10−22
d2 −6.3 × 10−4 −3.1 × 10−4 −9.4 × 10−4 −6.3 × 10−4
d3 1.08 × 10−8 5.43 × 10−9 −5.42 × 10−9 −1.79 × 10−20
d4 1.0637 × 10−2 1.0406 × 10−2 1.0406 × 10−2 1.0868 × 10−2

Before the calibrated response surface function can be used to predict the basal heave
deformation values at non-experimental points, a validation process must be conducted. A
total of 200 random samples were generated in order to verify the accuracy of response
surface function. Taking B = 3H and 6H as examples, Figure 6 plots the scatters of FEM
software calculation results with response surface predictions. It is noticed that almost all
the scatter points fall onto the 45◦ line, demonstrating that the predictions from response
surface function agree fairly well with those from FEM calculations. The fitted coefficient
of determination R2 = 0.9987, 0. 9838, 0.9989, 0.9992. Therefore, the calibrated response
surface functions at the points can be adopted to circumvent the time-consuming FEM
calculations.
 surface function. Taking B = 3H and 6H as examples, Figure 6 plots the scatters of FEM
software calculation results with response surface predictions. It is noticed that almost all
the scatter points fall onto the 45° line, demonstrating that the predictions from response
surface function agree fairly well with those from FEM calculations. The fitted coefficient
of determination R2 = 0.9987, 0.9838, 0.9989, 0.9992 . Therefore, the calibrated response
Sustainability 2023, 15, 8985 12 of 19
surface functions at the points can be adopted to circumvent the time-consuming FEM
calculations.

Sustainability 2023, 15, x FOR PEER REVIEW 13 of 20

(a) (b)

(c) (d)
Figure
Figure6.6.Comparison
Comparisonof of the
the FEM results for
FEM results for the
thebasal
basalheave
heavedeformation
deformationvalues
values and
and thethe
RSMRSM re-
results.
sults. (a) B = 3H; (b) B = 3H; (c) B = 6H; (d) B = 6H.
(a) B = 3H; (b) B = 3H; (c) B = 6H; (d) B = 6H.

The
Theresponse
responsesurface
surfacefunctions
functionscalibrated
calibratedusing
usingthetheremaining
remainingcharacteristic
characteristicpoints
points
are
areverified
verified thorough
thorough the similar procedure,
the similar procedure,and andthe
thefitted
fittedcoefficients
coefficients
of of determination
determination are
are greater
greater thanthan
0.9, 0.9, implying
implying that that the response
the response surface
surface functions
functions calibrated
calibrated on the on the char-
characteristic
acteristic
points canpoints can be substituted
be substituted for the
for the FEM FEM calculations
calculations in thestudy.
in the current current study.

5.4.2.Monte
5.4.2. MonteCarloCarloMethod
MethodSensitivity
SensitivityAnalysis
Analysis
AAsensitivity analysisofofNNisisgenerally
sensitivityanalysis generallyrequired
requiredtotodetermine reasonableNNvalue
determineaareasonable valueinin
MCS.When
MCS. WhenNNisistoo toosmall,
small,the obtainedPPf fhas
theobtained hasaalarge
largevariability
variabilityand
andcannot
cannotaccurately
accurately
expressthe
express thereliability
reliabilityofof the
the basal
basal heave
heave stability
stability according
according to Equation
to Equation (5). (5).
When When
N is tooN is
large, the required computational effort is prohibitively large. A reasonable N value isis
too large, the required computational effort is prohibitively large. A reasonable N value
determinedby
determined bybalancing
balancingthe theaccuracy
accuracyofofthe thecalculation
calculationresults
resultsand
andthethecalculation
calculationcost. cost.
Thesensitivity
The sensitivity analysis
analysis of N is
of N is illustrated
illustrated through
throughthe theFEMFEMmodel withB B= =HHand
modelwith TT
and = =1313 m,
1414
m, m,m,andand 1515
m,m,respectively.
respectively. TheThethreshold
threshold dmax
dmaxisistentatively
tentativelybetween
between25 25mm. mm.Figure
Figure7
7shows
showsthe thevariations
variations inin PPff with
with N.N. It
It is
is noticed
noticed that N ranges
as N
that as ranges between
between 10 1011andand101033, ,the
the
calculated P fluctuates significantly. Additionally, as N ranges between 10 3 and 10 4 (for
calculated Pf ffluctuates significantly. Additionally, as N ranges between 10 and 10 (for 3 4

the casesofofTT==1313and
thecases and1414m),m),the calculatedPP
thecalculated f fluctuatesslightly.
f fluctuates AfterNNbecomes
slightly.After becomesgreatergreater
than1010 4 , the obtained Pf tends to stable.
be stable. When N is smaller than 3
10 casein the
than 4, the obtained Pf tends to be When N is smaller than 10 3 in the of case
T = 15of
T =the
m, 15calculated
m, the calculated Pf fluctuates
Pf fluctuates considerably.
considerably. Additionally, Additionally,
when it iswhen it isthan
greater greater
103, than
the
3
obtained Pf tends to
10 , the obtained Pf be stable.
tends N stable.
to be = 10 is N
6 selected 6
= 10 istoselected
accounttofor the variable
account for thetarget Pf intarget
variable the
following reliabilityreliability
Pf in the following analysis. analysis.
 7 shows the variations in Pf with N. It is noticed that as N ranges between 101 and 103, the
calculated Pf fluctuates significantly. Additionally, as N ranges between 103 and 104 (for
the cases of T = 13 and 14 m), the calculated Pf fluctuates slightly. After N becomes greater
than 104, the obtained Pf tends to be stable. When N is smaller than 103 in the case of T = 15
m, the calculated Pf fluctuates considerably. Additionally, when it is greater than 103, the
Sustainability 2023, 15, 8985 13 of 19
obtained Pf tends to be stable. N = 106 is selected to account for the variable target Pf in the
following reliability analysis.

Sustainability 2023, 15, x FOR PEER REVIEW 14 of 20

Figure 7. Sensitivity analysis of N (at dmax = 25 mm).

Figure 7. Sensitivity analysis of N (at dmax = 25 mm).
5.4.3.
5.4.3.Influence
Influenceof ofTTand
andBBon onBasal
BasalHeave
HeaveReliability
Reliability
The
Theprevious
previousdeterministic
deterministicnumerical
numericalsimulations
simulationsdemonstrate
demonstratethatthatthethebasal
basalheave
heave
deformation is dependent on B and T. The influence of B and T on the
deformation is dependent on B and T. The influence of B and T on the basal heave reliability basal heave relia-
bility is investigated in this section using Monte Carlo method and
is investigated in this section using Monte Carlo method and with the help of the validated with the help of the
validated response surface functions shown in Section 5.4.1. The
response surface functions shown in Section 5.4.1. The influencing parameters are B = kH,influencing parameters
are
k =B1.0,
= kH,
2.0,k .=. .1.0, 2.0,10.0
, and ……, and andT= 10.0
13,and T = 16,
14, 15, 13, and
14, 15,
17 16,
m. and 17 m. Threshold
Threshold dmax is 30 dmm.max is 30
The
mm.
centerThe center
point of pitpoint of pitis bottom
bottom focusedison, focused
i.e., theon, i.e., the deformation
deformation value at thevalue
centerat point
the center
of pit
point
bottom of is
pitused
bottom is used to characterize
to characterize the failure
the failure sample and tosample and the
calculate to calculate
basal heave thefailure
basal
heave failureLet
probability. probability.
the y axis Let the y axis
represent represent the
the calculated calculated
Pf and Pf and the
x axis signify x axis signify the
parameterized
parameterized
B. Figure 8 plots B. Figure 8 plotsin
the variation thePf variation in Pf with parameterized
with parameterized B for a specific B T.
forThe
a specific
variationT.
The variation
curves under curves
differentunder
Ts aredifferent
shownTs in are showncolors
different in different colors and
and symbols. For symbols.
example,For the
example, the black
black square, square,
red circle, bluered circle,
upper blue upper
triangular, triangular,
green green lower
lower triangular, triangular,
and purple and
diamond
purple
representdiamond
the resultsrepresent
for T = the13,results
14, 15,for
16,Tand
= 13,1714,
m,15, 16, and 17 m, respectively.
respectively.

Figure
Figure8.8.Influence
Influenceof
ofTTand
andBBon
onPPf f(at
(atcenter
centerpoint
pointofofpit
pitbottom).
bottom).

ForTT=13
For thePfPvalues
= 13m,m,the f values underdifferent
under differentparameterized
parameterizedBsBsare
are0.00%,
0.00%,0.52%,
0.52%,1.12%,
1.12%,
1.46%,1.37%,
1.46%, 1.37%,1.16%,
1.16%,1.09%,
1.09%,1.00%,
1.00%,0.93%,
0.93%,and
and0.83%,
0.83%,respectively.
respectively. ItItcan
canbebeobserved
observedthat
that
asBBincreases
as increasesfrom
fromH Htoto4H, thePPf fincreases
4H,the increasessignificantly.
significantly.Then,
Then,ititvaries
variesslightly
slightlywhen
whenitit
reaches 1.37%. A similar variation trend has been found for T = 14, 15, 16 and 17 m. The
difference lies in the B at which the maximum Pf appears. This observation is attributed to
the weights in the contributions of B and T. It can be expected that when T is so large that
freely developed basal heave failure occurs for any of Bs, the Pf will increase as B increases,
Sustainability 2023, 15, 8985 14 of 19

reaches 1.37%. A similar variation trend has been found for T = 14, 15, 16 and 17 m. The
difference lies in the B at which the maximum Pf appears. This observation is attributed to
the weights in the contributions of B and T. It can be expected that when T is so large that
freely developed basal heave failure occurs for any of Bs, the Pf will increase as B increases,
which is an intuitive variation trend. Additionally, the Pf increases as T increases for the
specific B.

5.4.4. Monitoring Point Selection Based on Probabilistic Analysis

From the previous analysis, it can be concluded that the basal heave deformation
Sustainability 2023, 15, x FOR PEER REVIEW and the points with the maximum deformation value hinge on B and T, and 15
value that of the
20
maximum value is not entirely located in the center point of the pit bottom. Therefore, only
selecting the point at the center of the pit bottom or evenly placed monitoring points as
monitoring points will cause researchers to perform inaccurate assessments of the safety of
safety
the pitofbottom.
the pitIn bottom. In thisthe
this section, section, the monitoring
monitoring pointscheme
point selection selection scheme is deter-
is determined based
mined based on probabilistic basal
on probabilistic basal heave stability analysis. heave stability analysis.
The
Thebasalbasalheave
heavestability
stabilityisistreated
treatedas asaa serial
serial system
system problem,
problem,and and eacheach monitoring
monitoring
point
point at the pit bottom is treated as its subsystem. If the deformation valueat
at the pit bottom is treated as its subsystem. If the deformation value atany
anyof ofthethe
monitoring
monitoring points exceeds d , the system fails. If the center point of the pit bottomcan
points exceeds d max
max , the system fails. If the center point of the pit bottom can
represent
represent the the rest
restofofthethemonitoring
monitoring points
points to determine
to determine the basal
the basal heaveheave
failures, failures,
the evenly the
evenly
distributed distributed
monitoring monitoring
points points at the
at the pit pit bottom
bottom wouldwould be a waste
be a waste of financial
of financial resources re-
sources according to the specification. Hence, the basal heave
according to the specification. Hence, the basal heave failure probability at the center pointfailure probability at the
center
of pit point
bottom of ispitcompared
bottom isagainst
compared against
system system
failure failure probability
probability obtained by obtained
analyzing by an-the
alyzing the failure probabilities at all the monitoring points under
failure probabilities at all the monitoring points under serial system assumption. The basal serial system assump-
tion.
heave The basal at
failure heave failure point
the center at the of center point ofispit
pit bottom bottom
named is named
failure mode failure
I andmodethat at I andthe
that at the remaining
remaining points is namedpointsfailure
is named mode failure mode
II. The II. The respective
respective failure probabilities
failure probabilities are Pf1 and
are
Pf2 .PPf1f sys
andis Pthe
f2. P fsys is the
system system
failure failure probability,
probability, which is thewhich union is
ofthe union
failure mode of failure
I and II.modeUnder I
and
T = 13II. m and dTmax
Under = 13 m and
= 25 mm,dFigure
max = 259 mm, Figure 9 demonstrates
demonstrates the variations theof Pvariations
f1 , P f2 and of
P P
f sysf1, P
withf2

and Pfsys with parameterized

parameterized Bs in orange bar, Bs inblue orange bar, red
bar and bluereal
barcircles,
and red real circles, respectively.
respectively.

Figure9.9.The
Figure variationsPPf1f1, P
Thevariations ,Pf2.f2P Pf sys
. fsys with
with parameterized Bs Bs
parameterized T =T13
under
under = 13
m. m.

WhenBB==2,2,3,3,4,4,5H,
When thePPf1f1isis9.61%,
5H,the 9.61%,24.16%,
24.16%,27.98%,
27.98%,and
and26.10%,
26.10%,respectively
respectivelywhile
while
the respective
the respective PfsysP is 9.9%, 24.47%, 28.29%, and 26.37%. The negligible difference
fsysis 9.9%, 24.47%, 28.29%, and 26.37%. The negligible difference between between
PPf1f1and
andPPfsys impliesthat
fsysimplies thatthe
thecenter
centerpointpointofofthe
thepit
pitbottom
bottomcan
canrepresent
representthe
therest
restof
of the
the
characteristic points and that the optimum monitoring point is the center
characteristic points and that the optimum monitoring point is the center point. However, point. However,
for the rest of cases, i.e., B = 6, 7, 8, 9, 10H, the respective Pf1 is 23.33%, 21.43%, 19.76%,
18.57%, and 17.41%, the respective Pf2 is 23.89%, 25.40%, 22.49%, 22.86%, and 22.56%,
while the Pfsys is 23.89%, 25.40%, 22.49%, 22.86%, and 22.86%, respectively. Considerable
difference can be noticed between Pf1 and Pfsys, just as there are negligible differences be-
Sustainability 2023, 15, 8985 15 of 19

for the rest of cases, i.e., B = 6, 7, 8, 9, 10H, the respective Pf1 is 23.33%, 21.43%, 19.76%,
18.57%, and 17.41%, the respective Pf2 is 23.89%, 25.40%, 22.49%, 22.86%, and 22.56%,
while the Pfsys is 23.89%, 25.40%, 22.49%, 22.86%, and 22.86%, respectively. Considerable
difference can be noticed between Pf1 and Pfsys , just as there are negligible differences
between Pf2 and Pfsys , indicating that the center point is unable to represent the rest of
the characteristic points and that it is therefore not the optimum monitoring point. The
optimum monitoring point is favorably selected from the rest of the characteristic points.
To furtherly verify the observation in Figure 9, a flowchart in Figure 10 is designed to
calculate the conditional events E1 and E2. Where E1 = the basal heave fails at the rest of
Sustainability 2023, 15, x FOR PEER REVIEW
the characteristic points on the condition that the basal heave failure does not occur16atofthe 20

center point of the pit bottom; E2 = the basal heave fails at the center point of pit bottom
on the condition that the basal heave failure does not occur at the rest of the characteristic
points.
points. TheTherespective
respectiveconditional
conditionalfailure probabilities
failure Pf4 Pand
areare
probabilities Pf5 .PLet
f4 and ∆ f 4∆
f5. Let =f4 =P P ⁄Pfsys
f 4f4/P fsys
and∆∆f f55 =
and = PPf5f⁄5P/P representthe
fsysrepresent
fsys
therelative
relativecontributions
contributionsofof E1
E1 and
and E2 to the system
system failurefailure
probability.
probability.

Figure10.
Figure 10.The
Theflowchart
flowchartof
ofverifying
verifyingfeature
featurepoint
pointselection.
selection.

Figure
Figure 1111 plots the variations
variations of of ∆∆f4f 4and
and∆∆f5f with5 withB. B.
WhenWhen B within
B is is within thethe range
range be-
between 2H and 5H, the respective ∆ is 2.93%, 1.68%, 1.10%, and
tween 2H and 5H, the respective ∆f4f 4is 2.93%, 1.68%, 1.10%, and 0.99%. This shows that0.99%. This shows that
the
the failure
failureprobability
probability of
of the
the conditional
conditional eventevent E1 E1 isis negligibly
negligibly low.
low. That
That is, is, only
only aasmall
small
number
number of basal heave failures occur at the rest of the characteristic points forthe
of basal heave failures occur at the rest of the characteristic points for thesamples
samples
where
wherethe thebasal
basalheave
heavedoes
doesnot notfail
failatatthe
thecenter
centerpoint pointofofthe
thepit
pitbottom.
bottom. Therefore,
Therefore, thethe
center
center point of pit bottom is selected to be the optimum monitoring point inthis
point of pit bottom is selected to be the optimum monitoring point in thiscase.
case. In
In
another
anothercase whereBBisisgreater
casewhere greaterthan
than5H,5H,the the∆∆f 5f5 isis00for
forthe casesofofBB==6,6,7,7,8,
thecases 8,9,
9,10H.
10H.This
This
positively demonstrates that the basal heave does not fail at the center point of pit bottom
positively demonstrates that the basal heave does not fail at the center point of pit bottom
for the samples where the basal heave failures do not occur at the rest of the characteristic
for the samples where the basal heave failures do not occur at the rest of the characteristic
points. The optimum monitoring points are to be selected from the rest of the characteristic
points. The optimum monitoring points are to be selected from the rest of the characteris-
points.
tic points.
 where the basal heave does not fail at the center point of the pit bottom. Therefore, the
center point of pit bottom is selected to be the optimum monitoring point in this case. In
another case where B is greater than 5H, the ∆f5 is 0 for the cases of B = 6, 7, 8, 9, 10H. This
positively demonstrates that the basal heave does not fail at the center point of pit bottom
for the samples where the basal heave failures do not occur at the rest of the characteristic
Sustainability 2023, 15, 8985 16 of 19
points. The optimum monitoring points are to be selected from the rest of the characteris-
tic points.

Sustainability 2023, 15, x FOR PEER REVIEW 17 of 20

Figure
Figure11.
11.Variations of ∆𝑓4
Variations of f 4 and ∆
∆f𝑓5 withB.
5 with B.

For
Foraaspecific
specificBB==10H, 10H,different
differentTs Tsareareadopted
adoptedin inorder
orderto toinvestigate
investigatethe thevariations
variations
of
of Pf1 , Pf2 and Pf sys . Figure 12 shows the variations of Pf1 , Pf2 and Pf sys with TB for
P f1, P f2 and P fsys . Figure 12 shows the variations of Pf1 , P f2 and P fsys with T for = 10H.
B = The
10H.
respective
The respective Pf1 is 0.83%, 4.72%,4.72%,
Pf1 is 0.83%, 18.25%, 48.19%,
18.25%, 83.83%,83.83%,
48.19%, and 99.27% at T = 13,
and 99.27% at 14,
T =15, 13,16,
14,17,
15,
18
16,m.
17,The respective
18 m. Pf2 is 1.16%,
The respective Pf2 is5.75%,
1.16%,20.78%, 51.68%, 51.68%,
5.75%, 20.78%, 88.11%, and 99.60%.
88.11%, and The respec-
99.60%. The
respective
tive Pf sys is 5.75%,
Pfsys is 1.16%, 1.16%, 5.75%,
20.78%,20.78%,
51.68%,51.68%,
88.11%,88.11%, and 99.60%.
and 99.60%. It can Itbecan be noticed
noticed that
that con-
considerable
siderable difference
difference between
between Pf1 Pand
f1 and
P P
fsys is
f sys is observed
observed except
except for for
in in
thethe case
case ofofT T== 18
18 m,
m,
where a negligible
where a negligible difference difference between P and P
f2 and Pfsys
f2 is found for all the T cases. The
f sysis found for all the T cases. The variation variation
trendimplies
trend impliesthat thatthe theoptimum
optimummonitoring
monitoringpoints pointsoughtoughtto tobebeselected
selectedformformthe therest
restofof
thecharacteristic
the characteristicpoints. points.ForForthethespecial
specialcasecaseTT==18 18m,m,thethecenter
centerpointpointof ofpit
pitbottom
bottomcan can
beselected
be selectedto tobebe thethe optimum
optimum monitoring
monitoringpoint. point.Therefore,
Therefore,the theoptimum
optimumdesigndesignfor forthe
the
monitoringpoints
monitoring pointsfor forbasal
basalheave
heavestability
stabilityisisrecommended
recommendedto toaccount
accountfor forthe
theweighs
weighsin in
the contribution to the basal heave deformations of the
the contribution to the basal heave deformations of the influencing factors as excavationinfluencing factors as excavation
widthand
width andthickness
thicknessof ofsoft
soft soil
soil layer.
layer.

Figure12.
Figure Theprobabilities
12.The probabilitiesof
offailure
failurein
indifferent
differentthickness
thicknessof
ofthe
the soft
soft soil for BB ==10H.
soil for 10H.

6. Summary and Conclusions

6. Summary and Conclusions
The basal heave stability is one of the important issues in geotechnical practices. In
The basal heave stability is one of the important issues in geotechnical practices. In
addition to exploring a factor in safety-based evaluations, this paper aims to develop a
addition to exploring a factor in safety-based evaluations, this paper aims to develop a
deformation-based methodology combining non-intrusive FEM simulations and response
surface functions for assessing the basal heave stability in a probabilistic manner. The
hardening soil model embedded in commercial software package PLAXIS is adopted to
obtain the basal heave deformations at characteristic points, which are the candidate mon-
Sustainability 2023, 15, 8985 17 of 19

deformation-based methodology combining non-intrusive FEM simulations and response

surface functions for assessing the basal heave stability in a probabilistic manner. The
hardening soil model embedded in commercial software package PLAXIS is adopted
to obtain the basal heave deformations at characteristic points, which are the candidate
monitoring points. The deterministic basal heave deformation analysis is conducted in a
way that considers the influencing factors, such as excavation width (B) and thickness of
soft soil layer. The effects of B and thickness of soft soil layer on the basal heave reliability
are investigated. Finally, the recommendations for selecting optimum monitoring points
are given based on the comparisons of the failure probabilities at characteristic points
and the system failure probability. It must be noted that the proposed framework can be
applicable to excavation scenarios where deformation values are focused. The use of this
framework facilitates the design of excavation width and depth if target reliability index is
available. Most important, the optimum selection of monitoring points is helpful in cost
savings and thus favorable to sustainable underground space development. The following
may be noted:
(1) Basal heave deformations at the characteristic points are sensitive to the type of basal
heave failures, namely, the freely developed and prohibited basal heave failures. One
finds that the critical thickness of soft soil layer between the end of supporting struc-
tures and the top of a hard stratum is about 0.3B (B = excavation width). Beyond this,
the freely developed basal heave failure arises, otherwise the prohibited basal heave
failure occurs. In the case of freely developed basal heave failure, the deformations at
the center point of the pit bottom tend to be larger than those characteristic points,
but this is not the case if prohibited basal heave failure occurs.
(2) The failure probability of basal heave failure at the center point increases significantly
as B ranges within a limited value and then it begins to decrease or to vary slightly at
a certain value under a given thickness of soft soil layer. If the thickness of soft soil
layer is so sufficiently large that freely developed basal heave failure occurs for any
value of B, the failure probability of basal heave failure at the center point increases as
B increases. For a specific B, the failure probability of basal heave failure at the center
point increases as T increases.
(3) If negligible discrepancy between the failure probability of basal heave failure at the
center point and system failure probability is found for the specific excavation model,
the center point is recommended to be optimum monitoring point. However, there
exist many cases where the center point cannot fairly represent the system response.
Therefore, other candidate characteristic points are expected to be monitoring points.
In a word, we recommend that researchers use the optimum design for the monitoring
points of basal heave stability in order to account for their respective contributions to
the basal heave deformations of the influencing factors such as excavation width and
thickness of soft soil layer.

Author Contributions: Conceptualization, G.Y.; Methodology, L.L.; Software, L.X.; Investigation,

Z.X.; Resources, J.L.; Writing—original draft, K.M.; Supervision, H.Y.; Funding acquisition, Y.F. All
authors have read and agreed to the published version of the manuscript.
Funding: The work described in the paper was supported by National Natural Science Foundation of
China (Grant No. 52171264) and by Research&Development Program from China State Construction
(Grant No. CSCEC-2020-Z-49), the financial supports are gratefully acknowledged.
Data Availability Statement: Date are available upon request from corresponding author at lil-
[email protected].
Acknowledgments: Thanks to all the anonymous reviewers for their value comments on this paper,
which improved quality of our paper.
Conflicts of Interest: We declare no conflict of interest.
Sustainability 2023, 15, 8985 18 of 19

References
1. Zhang, K.; Li, J. Accident analysis for “08.11.15” foundation pit collapse of Xianghu Station of Hangzhou metro. Chin. J. Geotech.
Eng. 2010, 32, 338–342. (In Chinese)
2. Xiao, X.; Yuan, J.; Zhu, Y. Recovery work of Singapore MRT Circle Line C824 Nicoll Highway Collapse (Part II)—Implementation
of recovery option. Mod. Tunn. Technol. 2010, 47, 71–78. (In Chinese)
3. Sun, H. Research on one pit collapse in Shanghai. Chin. J. Undergr. Space Eng. 2012, 8, 1743–1746. (In Chinese)
4. Wang, L.; Liu, Y.; Long, F.; Hong, Y. Collapse of deep excavations for metro lines in soft clay. Chin. J. Geotech. Eng. 2020, 42,
1603–1611. (In Chinese)
5. Goh, A.; Kulhawy, F.H.; Wong, K.S. Reliability assessment of basal-heave stability for braced excavations in clay. J. Geotech.
Geoenviron. Eng. 2008, 134, 145–153. [CrossRef]
6. Wang, H. Improvement on safety coefficient of heave-resistant stability of excavations. Rock Soil Mech. 2014, 35, 30–36. (In Chinese)
7. Bjerrum, L. Stability of strutted excavations in clay. Geotechnique 1956, 6, 32–47. [CrossRef]
8. Wang, B.; Xia, M. Embedment depth and internal force of diaphragm wall. Chin. J. Geotech. Eng. 1983, 5, 103–114. (In Chinese)
9. Wang, H. Size effect of foundation pits and calculation method of safety factor of heave-resistant stability considering excavation
width. Rock Soil Mech. 2016, 37, 433–441. (In Chinese)
10. Wang, H.; Shen, X. Heave-resistant stability analysis method of foundation pit considering support. Rock Soil Mech. 2020, 41,
1680–1689. (In Chinese)
11. Do, T.N.; Ou, C.Y.; Lim, A. Evaluation of factors of safety against basal heave for deep excavations in soft clay using the finite
element method. J. Geotech. Geoenviron. Eng. 2013, 139, 2125–2135. [CrossRef]
12. Yu, J.; Long, Y.; Xia, X.; Gong, X. Basal stability for narrow foundation pit. J. Zhejiang Univ. 2017, 51, 2165–2174. (In Chinese)
13. Wang, C.; Zhang, T. The failure patterns and stability analysis of basal heave of long narrow foundation pits. Sci. Technol. Eng.
2017, 17, 276–281. (In Chinese)
14. Xiao, W.; Lu, Y.; Yang, Y.; Xiao, T.; Hao, Y. Stability analysis of narrow foundation pit considering size effect. Sci. Technol. Eng.
2022, 22, 4105–4112. (In Chinese)
15. Abdi, A.S.; Ou, C.Y. Numerical study of the effect of ground improvement on basal heave stability for deep excavations in
normally consolidated clays. J. Geotech. Geoenviron. Eng. 2023, 149, 04023042. [CrossRef]
16. Yang, P.; Zhang, J.; Yang, P. Study on anti-uplift stability of silt foundation pit bottom based on strength reduction method. J.
North China Inst. Sci. Technol. 2022, 19, 68–74. (In Chinese)
17. Zheng, G.; Dou, H. Stability and deformation problems during deep foundation excavation in soft soil area and their control
measures-stability problems during the whole process of excavation construction. Constr. Technol. 2011, 40, 1–6, 10. (In Chinese)
18. Wang, H.; Jia, M.; Yang, X. Influence of size effect on deformation of foundation pit and control measures. Sichuan Build. Sci. 2016,
42, 69–74. (In Chinese)
19. Zhang, T.; Baroth, J.; Dias, D. Deterministic and probabilistic basal heave stability analysis of circular shafts against hydraulic
uplift. Comput. Geotech. 2022, 150, 104922. [CrossRef]
20. Zhang, R.; Goh, A.T.C.; Li, Y.; Liu, H. A simple estimation model for basal heave stability of braced excavations in anisotropic
clay. Acta Geotech. 2022, 17, 5789–5800. [CrossRef]
21. Li, S.; NI, D.; Yan, Z.; Chen, J.; Hu, M. Strength reduction method for safety coefficient of heave-resistant stability of
asymmertrically-loaded excavations in soft soil areas. Chin. J. Geotech. Eng. 2019, 41, 13–16. (In Chinese)
22. Zhang, W.; Zhang, R.; Wu, C.; Goh, A.T.C.; Wang, L. Assessment of basal heave stability for braced excavations in anisotropic clay
using extreme gradient boosting and random forest regression. Undergr. Space 2020, 7, 233–241. [CrossRef]
23. Yodsomjai, W.; Lai, V.Q.; Banyong, R.; Chauhan, V.B. A machine learning regression approach for predicting basal heave stability
of braced excavation in non-homogeneous clay. Arab. J. Geosci. 2022, 15, 873. [CrossRef]
24. He, Y. Simplified prediction of springback deformation in soft soil foundation pit considering plastic development coefficient. J.
Yangtze River Sci. Res. Inst. 2015, 32, 94–98, 103. (In Chinese)
25. Zheng, G. Method and application of deformation control of excavations in soft ground. Chin. J. Geotech. Eng. 2022, 44, 1–36, 201.
(In Chinese)
26. Cao, L. Finite element method analysis of bottom upheaval of deep foundation pits in soft-clay ground due to excavation. Chin. J.
Geotech. Eng. 2013, 35, 819–824. (In Chinese)
27. Schanz, T.; Vermeer, P.A.; Bonnier, P.G. Formulation and verification of the Hardening-Soil Model. In Beyond 2000 in Computational
Geotechnics; Routledge: Abingdon, UK, 1999.
28. Wang, H.; Xu, M.; Song, E. A small strain constitutive model based on hardening soil model. Rock Soil Mech. 2011, 32, 39–43, 136.
(In Chinese)
29. Xue, Q. Analysis of the deep excavation deformation close to the existing operating railway based on Harding-Soil Model. J.
Railw. Eng. Soc. 2015, 32, 39–42, 92. (In Chinese)
30. Wang, W.; Wang, H.; Xu, Z. Experimental study of parameters of harding soil model for numerical analysis of excavation of
foundation pits. Rock Soil Mech. 2012, 33, 2283–2290. (In Chinese)
31. Gu, X.; Wu, R.; Liang, F.; Gao, G. On HSS model parameters for Shanghai soils with engineering verification. Rock Soil Mech. 2021,
42, 833–845. (In Chinese)
Sustainability 2023, 15, 8985 19 of 19

32. Liu, W.; Zhu, H.; He, S.; Yan, J.; Xu, C. Experimental study on parameters of hardening soil model for soils and its application in
foundation pit engineering in Nanchang. J. Civ. Environ. Eng. 2021, 43, 38–47. (In Chinese)
33. GB-50202-2002; Code for Acceptance of Construction Quality of Building Foundation. Ministry of Construction of the People’s
Republic of China: Beijing, China, 2002.
34. GB50497-2009; Technical Code for Monitoring of Building Excavation Engineering. Ministry of Construction of the People’s
Republic of China: Beijing, China, 2009.
35. Cao, J.; Gao, Y.; Liu, H. Application of Monte Carlo Method improved in reliability analysis of foundation pit system. Chin. J.
Undergr. Space Eng. 2019, 15, 1565–1572. (In Chinese)
36. Yang, L.; Xu, C. Monte Carlo simulation and reliability analysis on the deformation of foundation pit. Rock Soil Mech. 1999, 1,
16–19. (In Chinese)
37. Xu, Y.; Jiang, H.; Ma, G. Reliability analysis on stability of slopes based on response surface method. Port Waterw. Eng. 2009,
55–58. (In Chinese)
38. Jiang, K.; Meng, L. Application of slope optimization through enumerative algorithm based on response surface. Constr. Technol.
2015, 44, 77–80. (In Chinese)
39. Tan, J.; Shao, G.; Yu, J.; Li, G.; Lou, W.; Sun, Y. Reliability analysis of diaphragm wall in foundation pit based on stochastic
response surface method. Henan Sci. 2020, 38, 1831–1836. (In Chinese)
40. Li, D.Q.; Jiang, S.H.; Cao, Z.J.; Zhou, W. A multiple response-surface method for slope reliability analysis considering spatial
variability of soil properties. Eng. Geol. 2014, 187, 60–72. [CrossRef]
41. Chen, P.; Zhu, Y.; Xie, Y. Reliability analysis of response surface method based on checking point fitting. Mech. Eng. 2023, 79,
81–83, 86. (In Chinese)
42. Zhang, J.; Zhang, L.M.; Tang, W.H. New methods for system reliability analysis of soil slopes. Can. Geotech. J. 2011, 48, 1138–1148.
[CrossRef]
43. Xu, C.; Zhao, Z. Experimental study and engineering application of soil hardening model parameters of Beijing. In Proceedings of
the 2017 National Engineering Geology Academic Annual Conference, Guilin, China, 26–31 October 2017; Science Press: Beijing,
China, 2017; Volume 5. (In Chinese)
44. Hamdy, F.; Cai, F.; Ugai, K.; Hagiwara, T. Two-dimensional base stability of excavations in soft soils using FEM. Comput. Geotech.
2003, 30, 141–163.
45. Zhang, T.; Baroth, J.; Dias, D. Probabilistic basal heave stability analyses of supported circular shafts in non-homogeneous clayey
soils. Comput. Geotech. 2021, 140, 104457. [CrossRef]
46. Liu, G.; Hou, X. Residual stress analysis method for soft soil pit upheaval deformation. Tunn. Rail Transit 1996, 37, 107–115.
(In Chinese)
47. Du, T. Research on Deformation Law and Stability of Subway Long and Narrow Deep Foundation Pit. Master’s Thesis,
Beijing Jiaotong University, Beijing, China, 2012. (In Chinese)
48. Hu, M.; Gou, C.; Yan, Y.; Yuan, Z.; Chen, J.; Yan, Z. Finite element analysis of the effect of foundation pit width on its stability.
Chin. J. Ground Improv. 2020, 2, 1–8. (In Chinese)
49. Goh, A.; Zhang, W.G.; Wong, K.S. Deterministic and reliability analysis of basal heave stability for excavation in spatial variable
soils. Comput. Geotech. 2019, 108, 152–160. [CrossRef]

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual
author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to
people or property resulting from any ideas, methods, instructions or products referred to in the content.