Undrained stability of braced excavations in clay considering the nonstationary random field of undrained shear strength

The basal heave stability of supported excavations is an essential problem in geotechnical engineering. This paper considers the probabilistic analysis of basal heave stability of supported excavations with spatially random soils by employing the random adaptive finite element limit analysis and Monte Carlo simulations to simulate all possible outcomes under parametric uncertainty. The effect of soil strength variability is investigated for various parameters, including the width and depth of the excavation ratio, strength gradient factor, and vertical correlation length. Probabilistic basal stability results have also been employed to determine the probability of design failure for a practical range of deterministic factors of safety. Considering probabilistic failure analysis, the more complete failure patterns caused by the various vertical correlation length would decrease the probability of design failure. There are different tendencies between the probability of design failure at the same safety factor with various vertical correlation lengths. These results can be of great interest to engineering practitioners in the design process of excavation problems.


Problem definition
In this study, a braced excavation in spatially random clay with a linear increase of strength with depth is defined for deterministic and stochastic analyses and shown in Fig. 1.The excavation is under plane strain conditions and has a width (B) and depth (H) of excavation and the depth of embedment (D).The clay is considered a rigid-perfectly plastic Tresca material with a mean value of undrained shear strength (μ su0 ) at a depth z = 0 and linear strength gradient ρ.As a result, the mean value of undrained shear strength at any depth can be written as μ su = μ su0 + ρz, which is a function of depth.The linear function of undrained shear strength with depth was originally proposed by Bishop 53 .By defining the linear function that relates the undrained shear strength to the random field, the function can take the random value as input and return a corresponding undrained shear strength value.As a result, a non-stationary random field can be generated.More information regarding the non-stationary random field can be found in Yi et al. 54 and Liu et al. 55 .
The constant unit weight (γ) is the objective variable (or output) for the braced excavation in spatially random clay.For deterministic analysis, the stability number (N) of this excavation problem in plane strain conditions is the normalised term of the unit weight 10 , which can be expressed as follows: This stability number results in AFELA, which can be known as a function of the normalised width of excavation (B/H), the normalised depth of embedment (D/H), the strength gradient ratio (ρH/μ su0 ), and some dimensionless coefficients describing the inherent spatial variability of soil concisely.Note that the definition of the stability number (N) and other input non-dimensional parameters are based on the previous studies by Ukritchon et al. 10 ; Keawsawasvong and Ukritchon 11 ; and Kounlavong et al. 14 .The latter will be described extensively in Section "Random field theory".

Random field theory
This study considers clay's undrained shear strength (s u ) as the spatial variability for probabilistic investigations.The spatial variability of soil properties is assigned to any location in a domain known as random fields.According to the random field theory [56][57][58][59][60][61][62] , parameters CL x and CL y are carried out to define dimensional spatial where τ xij and τ yij are the absolute horizontal and vertical distances between two discrete points, respectively; CL x and CL y are the horizontal and vertical correlation lengths, respectively.The dimensionless correlation lengths Θ X and Θ Y are defined as follows: where H is the depth of excavation.
In this paper, Karhunen-Loeve (K-L) expansion method is adopted for modelling random fields because of the advantages of exponential covariance.According to K-L expansion, the analytical solution of the eigenvalue problem for an exponential function (Eq.2) is used to generate a continuous random field.The details of the procedures can be found in Cho 63 .
The log-normal distributions commonly produce positive variables without a negative random value of s u .Note that the natural log (e = 2.718) of random variables from a normal distribution curve is carried out in the log-normal distribution.By adopting the probability density function (PDF) 39,43 , the log-normal distribution of the undrained shear strength of clay can be expressed in Eqs.(3) to (6).www.nature.com/scientificreports/ The coefficient of variation (COV) is a practical parameter for describing the inherent variation of soil properties in the field [56][57][58][59][60][61][62] .To obtain the cumulative distribution function (CDF) for a continuous random variable, the probability density function (PDF) in Eq. ( 2) is integrated, where the expression of the CDF is shown in Eq. (7).
where erfc is the complementary error function 47,48 .
The probability of design failure (P f ) for a practical design application with a factor of safety is also an interesting issue.In this work, the failure is defined to occur when the calculated result of the stability number from the probabilistic analysis with random field theory (N ran ) is less than the deterministic one (N det ) with an appropriate factor of safety (FS) as follows: Therefore, the probability of design failure (P f ) can be defined as the probability that the number of calculated random stability numbers (N ran ) is less than the "factored" deterministic value (i.e., N det /FS) over the number of realisations or the number of Monte-Carlo simulations (N c ). N c can be set as 1,000 times to ensure reliable results of stochastic analysis.Note that this number of 1,000 times the MCs is selected based on several previous studies (e.g., Kasama and Zen 35 ; Kasama and Whittle 39 ; Huang et al. 36 ; Huang et al. 40 ; Li et al. 41 ; Ali et al. 43 ; Ali et al. 45 ; Ali et al. 46 ; Wu et al. 47 ; Wu et al. 48).The probability of design failure can be approximated as follows: where P f is the probability of design failure for a given value of FS.

Random adaptive finite element limit analysis
A numerical model of a braced excavation in the spatial variability of soil properties in OptumG2 is shown in Fig. 2. The boundary dimensions were defined to be at least 10 times of the width of excavation (B) in both vertical and horizontal as illustrated in Fig. 2. In the x-direction, the far vertical sides of the model are constrained.The bottom boundary domain is fixed in the x-and y-directions and the top ground surface and excavation area are free-moving surfaces.The wall is modelled by rigid plate elements, where the top of the plate boundary is activated by setting the wall to no horizontal moments and rotations to simulate the fully braced wall.All numerical models are rigorously built to ensure that the domain is adequately broad to prevent the boundary impact to achieve appropriate answers.The objective function in this study is the maximum unit weight of clay (6)   www.nature.com/scientificreports/ that results in the basal heave and the occurrence of the failure.The parameters arrangement for stochastic analysis are listed in Table 1.The adaptivity approach proposed by Lyamin et al. 44 is applied to increase the precision of upper and lowerbound solutions by merging adaptive mesh refinement with random finite element limit analysis.More information regarding RAFELA can be found in Ali et al. 43 The internal dissipation estimated from deviatoric stresses and strain rates (also known as shear power) is employed as the covariate in this adaptivity system.Three repetitions of adaptive meshing were utilised in all numerical simulations of the study, having the initial amount of 1,000 elements to the final amount of 3,000 elements.The Karhunen-Loeve (KL) expansion method is used to build a trustworthy random field for RAFELA.
Examples of the adaptive meshes are shown in Fig. 3 for deterministic analyses and typical realisation of stochastic analyses.The selected values of all parameters in Fig. 3   ρH/μ su0 = 1; COV = 60%; Θ X = 50.0;and Θ Y = 1.00.Figure 3 shows that the distribution of the undrained shear strength s u (z) for the deterministic analyses is a linearly increasing profile whereas that of the typical realisation shows a slightly ragged field of s u (z) occurring mainly in the vertical direction due to the moderate value of Θ Y = 1.0.The final adaptive meshes after three steps are also presented in Fig. 3, where the number of meshes increases in the zones that have high plastic shear strains which can reveal the failure patterns of the problem.www.nature.com/scientificreports/these values decrease with increasing normalised correlation length from Θ Y = 0.125 to 4.00.The distribution of s u (z) with a smaller value of Θ Y is more ragged in the vertical direction of random field profiles.The current study considers the six design parameters, including B/H, D/H, ρH/μ su0 , COV su , Θ X and Θ Y .Using the random stability number (N ran ) and the associated factors of safety (FS), a series of probabilistic design charts are presented for practical design uses and as decision-making considering the uncertainty of soil property.The selected range of all dimensionless parameters is shown in Table 1.Phoon and Kulhawy 51 indicated that the values of COV are about 25% to 60% for characterising inherent variation of soil properties at the site.The value of Θ X is also set to be 50 according to the previous works by Wu et al. 47 and Wu et al. 48.The values of B/H, D/H, and ρH/μ su0 are set to follow the previous works by Ukritchon et al. 10 and Keawsawasvong and Ukritchon 11 .

Results and discussion
As the first step of the investigation, the deterministic stability number of fully braced supported excavations in homogeneous clays under plane strain conditions from the present study is compared with existing solutions from the limit equilibrium method (LEM) by Terzaghi 1 and Bjerrum and Eide 2 , FEM by Goh 3 , and FELA by Ukritchon et al. 10 Fig. 5 shows that the solutions by Terzaghi 1 and Goh 3 provide an overestimation of the stability number.The present solutions using AFELA fairly agree with the previous solutions by Ukritchon et al. 10 and Bjerrum and Eide 2 .Note that Terzaghi 1 employed the LEM in this analysis by assuming the failure lines of excavations while Goh 3 carried out the solutions from FEM models with very coarse mesh distribution.This can lead to the large difference from the present study and Terzaghi 1 and Goh 3 since this study used the FELA with the mesh adaptivity approach to obtain more accurate limit state solutions of the excavation problem.
All RAFELA results of the mean of random stability numbers μN ran are presented in Tables 2 and 3 for COV = 25% and 60%, respectively.Some results are selected to portray the effect of soil strength variability for all considered parameters including width and depth of the excavation ratio, strength gradient factor and dimensionless vertical correlation length.Figure 6 presents the mean stability numbers μN ran against the dimensionless vertical correlation length Θ Y for various values of COV su , ρH/μ su0 , and D/H while the value of B/H is set to be 1.0.The mean of the random stability number significantly increases and tends to be stable corresponding with the lower and further values of the dimensionless vertical correlation length Θ Y because a lower value of Θ Y causes a rougher spatial variability of s u .A higher value of Θ Y would lead to a smoother spatial variability of s u dealing with a maintained value of random stability number N ran .When the value of Θ Y reaches 4.00, the completely smooth random field of soils with linearly increasing strength can be obtained.Hence, the value of the mean stability numbers μN ran come closer to the deterministic stability numbers N det .As expected, different values of the mean stability numbers μN ran between ρH/s u and COV can be observed, and the higher values of ρH/s u and COV are considered.The larger and lower values of the mean stability numbers μN ran are also observed.
In order to determine accurate probability of design failure (P f ), 1,000 simulations were conducted to achieve convergence of P f , ensuring that the P f has stabilised, which was demonstrated in Fig. 7 as an arbitrary case of braced excavations.Figure 8  , the effect of COV su on the P f values is found to be significant.A large value of COV su can result in a higher value of P f for the same giving value of FS.Moreover, the different P f among FS of COV = 60% are higher than that of COV = 25%.Based on the results in Fig. 8, the designed FS can be derived for no failures.For example, the excavation would be stable during the www.nature.com/scientificreports/construction process considering the designed FS > 1.6 at any Θ Y and COV su = 25%, while the distribution of s u is very high fluctuating at the site (COV = 60%).Hence, the designed FS should be increased, FS > 3.0 at any Θ Y .From the perspective of modelling material properties and loads in engineering, which are typically positive, the log-normal distribution has the advantage of allowing only positive values.Therefore, the results of all simulations should be considered to fit with the log-normal distribution in this study.The random variable X is said to follow a lognormally distribution if ln(X) follows a normally distribution.The standard deviation, the  coefficient of variation and the mean are derived through the transformation of the parameters of the normal distribution as defined in Eqs. ( 4), ( 5  distribution is smooth and wide when the values of Θ Y is small.However, the curve of the probabilistic density function becomes narrow and left skewed when Θ Y becomes very larger (e.g., Θ Y = 4.00).The dotted lines in the histograms in Fig. 9 represent the deterministic values of the stability number.The number of random stability numbers N ran , which are lesser than the deterministic one N det (or FS < 1), is higher when the value of Θ Y is larger.It should be significantly interesting from the results in Fig. 9 that over half of the number of random stability numbers N ran could be observed, which are lower than the value of deterministic analysis even though the spatial variability of s u becomes to be completely smooth.This implies that ignoring random fields is not suitable for assessing The effect of the strength gradient factor ρH/μ su0 on the mean of random stability numbers is demonstrated next.In Fig. 11, the tendency between the mean stability numbers μN ran and ρH/μ su0 is a linearly increasing line.An increase in ρH/μ su0 causes an increase in the mean stability number, and the gaps between COV = 25 and 60% increase with increasing ρH/μ su0 because the larger values of the undrained shear strength are simulated, and greater values of random stability number are produced.Figure 12 shows the variation of μN ran with different values of the width of excavation ratio B/H.Unlike the effect of ρH/μ su0 , the tendency between the mean stability numbers μN ran and B/H is the nonlinear decreasing line.The curves of the mean stability numbers μN ran quickly go down at the lesser B/H (from B/H = 0.2 to 2.0) and tend to be flat at the further H/B.The reduced levels of the mean stability numbers μN ran are significantly different among ρH/μ su0 at B/H = 0.2 to 2.0, which confirms that an excavation with a high area has less stability than a narrow area.The effects of the depth of excavation ratio D/H on the random stability number is illustrated in Fig. 13.The relationship between μN ran versus D/H is linear, where an increase in the depth of excavation ratio D/H yields an increase in the excavation stability.Hence, the level of the excavation stability would be increased corresponding to the larger value of ρH/μ su0 , and remain unchanged at a very small value of ρH/μ su0 for any D/H.Finally, the effect of B/H on the probability of design failure (P f ) is shown in Fig. 14.The effect of ρH/μ su0 and D/H on P f is very small and is not presented here for brevity.The results in Fig. 14 have shown that for all chosen values of FS from 1.0 to 3.0, an increase in B/H causes an increase in P f .When COV is larger, the value of P f is also larger for all giving values of FS.Relating to practical applications based on the findings in Fig. 14, the failure probability of braced excavation is not affected by B/H if the designed FS is great enough value when the depth of excavation is fixed at the site.

Conclusion
This study adopted the RAFELA method to investigate the effect of a nonstationary random field of undrained shear strength on the failure probability of fully braced excavation.The quantitative findings from this study can be highlighted as follow: 1.In comparison to the results obtained from the LEM and FEM, the stability number calculated in the present study using the AFELA method was found to be underestimated.Nonetheless, these results align with the deterministic analysis conducted in a previous study utilising the FELA method.2. In general, the mean stability number of stochastic analysis μN ran is smaller than the deterministic stability number.The larger vertical correlation length leads to a further increase in the mean stability number μN ran due to the more distributions of smooth spatial variability s u were produced by random field theory.3. Considering probabilistic failure analysis, the more complete failure patterns caused by the smaller vertical correlation length would decrease the probability of design failure.The failure probability was not affected by simulating the higher values of dimensionless vertical correlation length and B/H.Besides, the increase in FS would drastically decrease the failure probability, it is suggested that the FS = 1.6 can guarantee relatively high-level safety braced excavation with COV su = 25% for most cases.The FS should be raised with the greater COV su .For example, FS = 3.0 with COV su = 60% in this study.4. Based on the statistical analysis of undrained stability in braced excavations, it can be observed that disregarding the spatial variability of undrained shear strength leads to an overestimation of both the mean stability numbers and the probability of designed failures (P f ).Therefore, it becomes essential to consider the variability of soil properties for reliable analysis and design of braced excavations, particularly in highly non-homogeneous soil conditions.Moreover, the P f is highly sensitive to variations in different vertical correlation lengths and the FS. 5.The limitation of this study is that only the 2D plane strain condition is employed in the analysis.Hence, only the stability solutions of 2D braced excavations are considered in the present study, which differs from the real-world case for 3D braced excavations.Future works may include the 3D stability analysis of braced excavations in clay considering the nonstationary random field of undrained shear strength.

Fundings
The work was supported by Faculty of Engineering Research Fund, Thammasat University and the Thailand Science Research and Innovation Fundamental Fund fiscal year 2023.
horizontal and vertical directions to capture the scale of fluctuations of random soils.These parameters represent a characterisation of spatial variability of s u in the soil domain.A large dimensional spatial correlation length yields the smoothly varying field of random soils, whereas a low one makes a ragged random field.Every point in the random field becomes independent when the dimensional spatial correlation length approximates zero.In the random field context, CL x and CL y are generally incorporated through a correlation function ρ which can be assumed by the exponential function as shown in Eq. (2a).
shows examples of the variations of the probability of design failure (P f ) versus the dimensionless vertical correlation length Θ Y , where B/H is fixed as 0.25.Note that other parameters in Fig. 8 are D/H = 1.0, ρH/μ su0 = 1.0, and COV su = 25% and 60%.The various values of FS are considered by varying from FS = 1.0 to 3.0.The results indicate the curves are different trends between FS = 1.0 and other FS.The P f tends to decrease when FS = 1.0, while the increasing values of P f are generated by FS > 1.0 according to increasing values of Θ Y for both COV = 25 and 60%.The FS = 1.0 is well known as a critical state for simulating stable analysis.In addition, by comparing the results between Figs. 8(a) and 8(b)
stability excavation.Examples of failure mechanisms of the basal heave stability problem are shown in Fig. 10 for the different values of Θ Y and COV su .The distributions of shear dissipation are employed to show the pattern of failures, where the values of B/H = 1.0,D/H = 1.0, and ρH/μ su0 = 1.0.Figure10shows that considering a smaller value of Θ Y and a larger value of COV, the failure pattern is not symmetric because of the presence of more generated variation and rough spatial variability of s u .Figure10indicates that the pattern of failures becomes more completely symmetric when the values of Θ Y are larger.The failure patterns of the stochastic case with Θ Y = 4.00 for both COV su = 25 and 60% are analogous to that of the deterministic case.The findings in Fig.10also demonstrate no differences in the failure patterns between a lower and higher values of COV without considering spatial variability of s u .Thus, the failure mechanisms of the basal heave excavation problem mainly depend on a characterisation of spatial variability of s u .