Deformation-Based Basal Heave Reliability Analysis and Selection on Monitoring Points for General Braced Excavations

: 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.3 B ( 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 signiﬁcantly 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 sufﬁciently 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 inﬂuencing factors such as excavation width and thickness of soft soil layer. The proposed framework is applicable to basal heave reliability analysis for braced excavations where deformation values are focused.


Introduction
Basal heave stability has seen witnessed its rise importance within the large number of braced excavations [1][2][3][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 from worldwide researchers. Traditionally, basal heave stability is evaluated via factor of safety against basal heave. The methods for calculating the factor of safety against basal heave are classified into three categories, i.e., (1) the Terzaghi method based on bearing capacity [5][6][7], (2) the slip circle method [8][9][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 reduction technique is widely adopted owing to its advantage in the consideration of complicated braces and construction stages [12][13][14][15][16]. In numerical methods, the factor of safety against basal heave is similar to that used in slope stability analysis. It is a factor by which the shear strength is reduced to bring the pit system to a limit equilibrium state. It has been found that, before the pit system reaches the limit equilibrium state, fairly substantial deformation occurs at the locations adjacent to the pit walls. This results in the disfunction of pit from the perspective of serviceability limit state. Therefore, the factor of safety by 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][21][22][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 reliability 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 development. 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 methodology 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.

Deformation-Based Basal Heave Failures
The braced excavations are mainly conducted in the urban area with populated highrise 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.

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: E oed (tangent 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). k 0 is the coefficient of lateral pressure, σ 3/k 0 is effective vertical stress. p ref is the reference stress. c and ϕ are effective cohesion and effective internal friction angle. Previous studies [28][29][30][31][32] ,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.

Limit State Function in Terms of Basal Heave Deformation
From the perspective of serviceability limit design, the deformation values at characteristic 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: 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. d max 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. d c (X) is the predicted deformation value obtained using finite element simulation or by using response surface function with sufficient accuracy.

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 d c (X i ), 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 P f as: 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 P f is equal to the 'true' value. Traditionally, the choice of N is dependent on the accuracy of simulated P f . Let δ P f denote the coefficient of variation in simulated P f , and the following equation relates the target P f , N and δ P f [36]: where δ P f is the coefficient of variation of failure probability, P f is the target failure probability, and N is the total number of samples. A sensitivity study of N on P f is recommended to select a proper N before the MCS is applied.

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 computational 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 d c . If the calibrated response surface function has a fairly good accuracy, it is adopted instead of finite element simulation to obtain d c , thereby reducing the computational cost. The second-order polynomial response surface function without cross items has been widely used in the reliability analysis of geotechnical engineering problems [37][38][39][40][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 d c where d c r (X) is the predicted basal heave deformation value by response surface function, X = (x 1 , x 2 , . . . , x n ) is the vector of n random variables, and a, b i , and d i (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., X 1 = (µ x1 , µ x2 , . . . , µ xn ). For the ith random variable x i , 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 x i , 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 X 1 , X 2 , . . . , X 2n+1 are input into the finite element model in PLAXIS 2D to obtain the corresponding basal heave deformation values (d c (X 1 ), d c (X 2 ), . . . , d c (X 2n+1 )). By equating d c (X i ) with d c r (X i ), 2n + 1 linear equations regarding coefficients (a, b i , d i ), 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.

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 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 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 firstly designed to facilitate the implementation of part I. The statistics of soil parameters involved in the excavations are determined according to the engineering geological survey reports and the relevant literature to execute the second part. Part I is conducted by PLAXIS 2D CE V20 Update 4, as shown in the left half of Figure 1. This includes the establishment of FEM model, the selection of the characteristic points, the execution of the model, and the use of outputs for the deformation values. Part II starts with the experimental points design for the calibration of the response surface functions, followed by the repeated execution of deterministic analysis to obtain the benchmark deformation values at experimental points. The response surface functions are calibrated and validated using self-developed code in Matlab. The final step of part II is to conduct MCS to obtain the number of failure samples and to calculate the failure probability.
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 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 firstly designed to facilitate the implementation of part I. The statistics of soil parameters involved in the excavations are determined according to the engineering geological survey reports and the relevant literature to execute the second part. Part I is conducted by PLAXIS 2D CE V20 Update 4, as shown in the left half of Figure 1. This includes the establishment of FEM model, the selection of the characteristic points, the execution of the model, and the use of outputs for the deformation values. Part II starts with the experimental points design for the calibration of the response surface functions, followed by the repeated execution of deterministic analysis to obtain the benchmark deformation values at experimental points. The response surface functions are calibrated and validated using self-developed code in Matlab. The final step of part II is to conduct MCS to obtain the number of failure samples and to calculate the failure probability.

Finite Element Modeling
Consider a schematic model of braced excavation shown in Figure 2. The soil layers 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: 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 by EA, is 4.8 × 10 6 kN/m. The EA and EI of diaphragm walls are 12 × 10 6 kN/m and 120 × 10 3 kN·m 2 , respectively. The groundwater table is at a depth of −3 m. Table 1 lists the five construction stages of the braced excavations. The means of soil parameters are shown in Table 2.
tively. A width coefficient k is introduced to study the effect of B on the basal heave deformations: 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 by EA, is 4.8 × 10 6 kN m ⁄ . The EA and EI of diaphragm walls are 12 × 10 6 kN m ⁄ and 120 × 10 3 kN•m 2 , respectively. The groundwater table is at a depth of −3 m. Table 1 lists the five construction stages of the braced excavations. The means of soil parameters are shown in Table 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 (R inter ) to consider the strength parameters between structure and soil. The default set R inter = 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. The mesh density has a significant effect on the calculation results and the computational 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 deformation 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 CPU@3.60 GHz and 16.0 GB RAM. It can be noticed from Table 3   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 significant 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.

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 T c proposed in the literature [44] is used in this study. T c = B/ √ 2 is the minimum depth to ensure a freely developed sliding basal heave failure with D = 0 [44]. In this section, T = T c + 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.  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). When the excavation width is wide, i.e., B = 2H~10H, the basal heave deformation at the central part is larger, and that at the side parts is smaller. It can be observed that as B increases, the maximum basal heave deformation value rises significantly. It is recognized that the basal heave deformation is mainly attributed to two sources: (1) soil rebound due to the unloading factor from excavation (Source I) and (2) the deformation due to the inward movement of the support structure towards the excavation surface (Source II) [46]. When the excavation width is wide, i.e., B = 2H~10H, the basal heave deformation at the central part is larger, and that at the side parts is smaller. It can be observed that as B increases, the maximum basal heave deformation value rises significantly. It is recognized that the basal heave deformation is mainly attributed to two sources: (1) soil rebound due to the unloading factor from excavation (Source I) and (2) the deformation due to the inward movement of the support structure towards the excavation surface (Source II) [46]. When 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.

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: where T c = 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  investigate the effect of T on basal heave deformation. For a given combination of B and T, 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 deformation values at characteristic nodes under different B and T values. In the case of B = H, at a specific a = 0.1, as the characteristic points approach to the center point (x = 0), the maximum deformation value appears at x = −3.5. As x increases from −5 to −3.5, the deformation value increases to the maximum value, and then it decreases when x is greater than −3.5. This variation trend has been observed for other 14 as. For the specific characteristic point (e.g., center point), as a increases from 0.1 to 1.5, i.e., T increases from 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 points. 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 investigate the effect of T on basal heave deformation. For a given combination of B and T, 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 deformation values at characteristic nodes under different B and T values. In the case of B = H, at a specific a = 0.1, as the characteristic points approach to the center point (x = 0), the maximum deformation value appears at x = −3.5. As x increases from −5 to −3.5, the deformation value increases to the maximum value, and then it decreases when x is greater than −3.5. This variation trend has been observed for other 14 as. For the specific characteristic point (e.g., center point), as a increases from 0.1 to 1.5, i.e., T increases from 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 points.   In the case of B = 3H, under a = 0.1, a similar variation trend is found to that at B = H. For the as greater than 0.2, the maximum deformation value appears at the center point (x = 0) and this variation trend is observed to be different than in the case of B = H. For a given characteristic point, the deformation variation trend with k is observed to be same with that in the case of B = H. In the cases of B = 6H and 10H, similar variation trends have been noticed as those in the case of B = 3H.
The points with the maximum deformation values are connected in dashed pink lines, as seen in Figure 4, to clearly demonstrate the variations of the points. As a ranges between  Based on the aforementioned comparative results, it can be concluded that for the narrow pit (B/H < 2) the influence of T does not change the deformation variation trend with different characteristic points. However, for the wide pit (B/H ≥ 2), as a is greater than 0.4, where the freely developed basal heave sliding can occur, the maximum deformation value appears at the center point. Otherwise, the T prevents the full development of the failure surface of the soil and this results in a deformation trend similar to that observed from the narrow pit. It is worth pointing out that the critical thickness of soft soil layer between the end of supporting structures and the top of hard stratum is about 0.3B (0.3B ≈ 0.4 × B/√2).

Reliability Analysis
Following the preliminary sensitivity analysis, E 2,ur ref , φ 2 ' , E 3,50 ref and φ 3 ' are selected as lognormal distributed random variables, where the single subscript or that before the coma is the number of the soil layer. The statistics (e.g., means and standard deviation) of each random variable are shown in Table 4. A series of points with equal horizontal distance at 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 d max can be ten- Based on the aforementioned comparative results, it can be concluded that for the narrow pit (B/H < 2) the influence of T does not change the deformation variation trend with different characteristic points. However, for the wide pit (B/H ≥ 2), as a is greater than 0.4, where the freely developed basal heave sliding can occur, the maximum deformation value appears at the center point. Otherwise, the T prevents the full development of the failure surface of the soil and this results in a deformation trend similar to that observed from the narrow pit. It is worth pointing out that the critical thickness of soft soil layer between the end of supporting structures and the top of hard stratum is about 0.3B (0.3B ≈ 0.4 × B/ √ 2).

Reliability Analysis
Following the preliminary sensitivity analysis, E ref 2,ur , ϕ 2 , E ref 3,50 and ϕ 3 are selected as lognormal distributed random variables, where the single subscript or that before the coma is the number of the soil layer. The statistics (e.g., means and standard deviation) of each random variable are shown in Table 4. A series of points with equal horizontal distance at 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 d max can be tentatively stated to be in the range between 25 and 35 mm. 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: where a, b i , d i , 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. 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 R 2 = 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 R 2 = 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. The response surface functions calibrated using the remaining characteristic points are verified thorough the similar procedure, and the fitted coefficients of determination are greater than 0.9, implying that the response surface functions calibrated on the characteristic points can be substituted for the FEM calculations in the current study.

Monte Carlo Method Sensitivity Analysis
A sensitivity analysis of N is generally required to determine a reasonable N value in MCS. When N is too small, the obtained Pf has a large variability and cannot accurately express the reliability of the basal heave stability according to Equation (5). When N is too large, the required computational effort is prohibitively large. A reasonable N value is determined by balancing the accuracy of the calculation results and the calculation cost. The sensitivity analysis of N is illustrated through the FEM model with B = H and T = 13 m, 14 m, and 15 m, respectively. The threshold d max is tentatively between 25 mm. Figure  7 shows the variations in Pf with N. It is noticed that as N ranges between 10 1 and 10 3 , the calculated Pf fluctuates significantly. Additionally, as N ranges between 10 3 and 10 4 (for the cases of T = 13 and 14 m), the calculated Pf fluctuates slightly. After N becomes greater than 10 4 , the obtained Pf tends to be stable. When N is smaller than 10 3 in the case of T = 15 m, the calculated Pf fluctuates considerably. Additionally, when it is greater than 10 3 , the obtained Pf tends to be stable. N = 10 6 is selected to account for the variable target Pf in the following reliability analysis. The response surface functions calibrated using the remaining characteristic points are verified thorough the similar procedure, and the fitted coefficients of determination are greater than 0.9, implying that the response surface functions calibrated on the characteristic points can be substituted for the FEM calculations in the current study.

Monte Carlo Method Sensitivity Analysis
A sensitivity analysis of N is generally required to determine a reasonable N value in MCS. When N is too small, the obtained P f has a large variability and cannot accurately express the reliability of the basal heave stability according to Equation (5). When N is too large, the required computational effort is prohibitively large. A reasonable N value is determined by balancing the accuracy of the calculation results and the calculation cost. The sensitivity analysis of N is illustrated through the FEM model with B = H and T = 13 m, 14 m, and 15 m, respectively. The threshold d max is tentatively between 25 mm. Figure 7 shows the variations in P f with N. It is noticed that as N ranges between 10 1 and 10 3 , the calculated P f fluctuates significantly. Additionally, as N ranges between 10 3 and 10 4 (for the cases of T = 13 and 14 m), the calculated P f fluctuates slightly. After N becomes greater than 10 4 , the obtained P f tends to be stable. When N is smaller than 10 3 in the case of T = 15 m, the calculated P f fluctuates considerably. Additionally, when it is greater than 10 3 , the obtained P f tends to be stable. N = 10 6 is selected to account for the variable target P f in the following reliability analysis. 7 shows the variations in Pf with N. It is noticed that as N ranges between 10 1 and 10 3 , the calculated Pf fluctuates significantly. Additionally, as N ranges between 10 3 and 10 4 (for the cases of T = 13 and 14 m), the calculated Pf fluctuates slightly. After N becomes greater than 10 4 , the obtained Pf tends to be stable. When N is smaller than 10 3 in the case of T = 15 m, the calculated Pf fluctuates considerably. Additionally, when it is greater than 10 3 , the obtained Pf tends to be stable. N = 10 6 is selected to account for the variable target Pf in the following reliability analysis.  For T =13 m, the Pf values under different parameterized Bs are 0.00%, 0.52%, 1.12%, 1.46%, 1.37%, 1.16%, 1.09%, 1.00%, 0.93%, and 0.83%, respectively. It can be observed that as B increases from H to 4H, the Pf increases significantly. Then, it varies slightly when it 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, For T = 13 m, the P f values under different parameterized Bs are 0.00%, 0.52%, 1.12%, 1.46%, 1.37%, 1.16%, 1.09%, 1.00%, 0.93%, and 0.83%, respectively. It can be observed that as B increases from H to 4H, the P f increases significantly. Then, it varies slightly when it 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 P f 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 P f will increase as B increases, which is an intuitive variation trend. Additionally, the P f increases as T increases for the specific B.

Monitoring Point Selection Based on Probabilistic Analysis
From the previous analysis, it can be concluded that the basal heave deformation value and the points with the maximum deformation value hinge on B and T, and that the 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 the pit bottom. In this section, the monitoring point selection scheme is determined based on probabilistic basal heave stability analysis.
The basal heave stability is treated as a serial system problem, and each monitoring point at the pit bottom is treated as its subsystem. If the deformation value at any of the monitoring points exceeds d max , the system fails. If the center point of the pit bottom can represent the rest of the monitoring points to determine the basal heave failures, the evenly distributed monitoring points at the pit bottom would be a waste of financial resources according to the specification. Hence, the basal heave failure probability at the center point of pit bottom is compared against system failure probability obtained by analyzing the failure probabilities at all the monitoring points under serial system assumption. The basal heave failure at the center point of pit bottom is named failure mode I and that at the remaining points is named failure mode II. The respective failure probabilities are P f1 and P f2 . P fsys is the system failure probability, which is the union of failure mode I and II. Under T = 13 m and d max = 25 mm, Figure 9 demonstrates the variations of P f1 , P f2 and P fsys with parameterized Bs in orange bar, blue bar and red real circles, respectively.    When B = 2, 3, 4, 5H, the P f1 is 9.61%, 24.16%, 27.98%, and 26.10%, respectively while the respective P fsys is 9.9%, 24.47%, 28.29%, and 26.37%. The negligible difference between P f1 and P fsys implies that the center point of the pit bottom can represent the rest of the characteristic points and that the optimum monitoring point is the center point. However, for the rest of cases, i.e., B = 6, 7, 8, 9, 10H, the respective P f1 is 23.33%, 21.43%, 19.76%, 18.57%, and 17.41%, the respective P f2 is 23.89%, 25.40%, 22.49%, 22.86%, and 22.56%, while the P fsys is 23.89%, 25.40%, 22.49%, 22.86%, and 22.86%, respectively. Considerable difference can be noticed between P f1 and P fsys , just as there are negligible differences between P f2 and P fsys , 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 the characteristic points on the condition that the basal heave failure does not occur at the 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. The respective conditional failure probabilities are P f4 and P f5 . Let ∆ f 4 = P f 4 /P fsys and ∆ f 5 = P f 5 /P fsys represent the relative contributions of E1 and E2 to the system failure probability. points. The respective conditional failure probabilities are Pf4 and Pf5. Let ∆ f4 = P f4 P fsys ⁄ and ∆ f5 = P f5 P fsys ⁄ represent the relative contributions of E1 and E2 to the system failure probability.  Figure 11 plots the variations of ∆ f4 and ∆ f5 with B. When B is within the range between 2H and 5H, the respective ∆ f4 is 2.93%, 1.68%, 1.10%, and 0.99%. This shows that the failure probability of the conditional event E1 is negligibly low. That is, only a small number of basal heave failures occur at the rest of the characteristic points for the samples 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 points. The optimum monitoring points are to be selected from the rest of the characteristic 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 points. The optimum monitoring points are to be selected from the rest of the characteristic points. Figure 11. Variations of ∆ 4 and ∆ 5 with B. Figure 11. Variations of ∆ f 4 and ∆ f 5 with B.
For a specific B = 10H, different Ts are adopted in order to investigate the variations of P f1 , P f2 and P fsys . Figure 12 shows the variations of P f1 , P f2 and P fsys with T for B = 10H. The respective P f1 is 0.83%, 4.72%, 18.25%, 48.19%, 83.83%, and 99.27% at T = 13, 14, 15, 16, 17, 18 m. The respective P f2 is 1.16%, 5.75%, 20.78%, 51.68%, 88.11%, and 99.60%. The respective P fsys is 1.16%, 5.75%, 20.78%, 51.68%, 88.11%, and 99.60%. It can be noticed that considerable difference between P f1 and P fsys is observed except for in the case of T = 18 m, where a negligible difference between P f2 and P fsys is found for all the T cases. The variation trend implies that the optimum monitoring points ought to be selected form the rest of the characteristic points. For the special case T = 18 m, the center point of pit bottom can be selected to be the optimum monitoring point. Therefore, the optimum design for the monitoring points for basal heave stability is recommended to account for the weighs in the contribution to the basal heave deformations of the influencing factors as excavation width and thickness of soft soil layer. For a specific B = 10H, different Ts are adopted in order to investigate the variations of Pf1, Pf2 and Pfsys. Figure 12 shows  . It can be noticed that considerable difference between Pf1 and Pfsys is observed except for in the case of T = 18 m, where a negligible difference between Pf2 and Pfsys is found for all the T cases. The variation trend implies that the optimum monitoring points ought to be selected form the rest of the characteristic points. For the special case T = 18 m, the center point of pit bottom can be selected to be the optimum monitoring point. Therefore, the optimum design for the monitoring points for basal heave stability is recommended to account for the weighs in the contribution to the basal heave deformations of the influencing factors as excavation width and thickness of soft soil layer.

Summary and Conclusions
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 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-

Summary and Conclusions
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 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 structures 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.