# Probabilistic analysis at the serviceability limit state of two neighboring strip footings resting on a spatially random soil

Ashraf Abd El-Modby Soliman Ahmed, Abdul-Hamid Soubra

► To cite this version:

Ashraf Abd El-Modby Soliman Ahmed, Abdul-Hamid Soubra. Probabilistic analysis at the serviceability limit state of two neighboring strip footings resting on a spatially random soil. Structural Safety, 2013, 49, pp.2-9. 10.1016/j.strusafe.2013.08.001. hal-01007203

HAL Id: hal-01007203

https://hal.science/hal-01007203

Submitted on 16 Jun 2014

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Probabilistic analysis at the serviceability limit state of two neighboring strip footings resting on a spatially random soil

Ashraf Ahmed<sup>a</sup>, Abdul-Hamid Soubra<sup>b</sup>

<sup>a</sup>University of Aswan, Abul-Rish kebly, Aswan, Egypt

<sup>b</sup>University of Nantes, Bd. de l'université, BP 152, 44603 Saint-Nazaire Cedex, France

The computation of the failure probability of geotechnical structures considering the soil spatial variability is generally performed using Monte Carlo Simulation (MCS) methodology. This method is very time-consuming especially when computing a small failure probability. As an alternative, Subset Simulation (SS) approach was proposed by Au and Beck [3] to efficiently calculate the small failure probability. In the present paper, a more efficient approach called "improved Subset Simulation (iSS)" approach is employed. In this method, the first step of the SS approach is replaced by a conditional simulation in which the samples are generated outside a hypersphere of a given radius. The efficiency of the iSS approach is illustrated here through the probabilistic analysis at the serviceability limit state (SLS) of two neighboring strip footings resting on a soil with 2D spatially varying Young's modulus. The system response is the differential settlement between the two footings. The probabilistic results have shown that the probability P<sub>e</sub> of exceeding a tolerable differential settlement computed by the iSS approach is very close to that calculated by the MCS methodology applied on the original deterministic model. The results have also shown that the use of the iSS approach has reduced the number of calls of the deterministic model by about 50% with respect to the SS approach.

Keywords:

Subset simulation

Conditional simulation

Strip footings

Differential settlement

## 1. Introduction

The classical Monte Carlo Simulation (MCS) methodology is generally used to calculate the failure probability of geotechnical problems involving random fields (e.g. [8,14,5] at ULS and [7,9] at SLS analysis). In these studies, only the mean value and the standard deviation of the system response were extensively investigated. This is because MCS requires a large number of calls of the deterministic model for the computation of the small failure probabilities. As alternative to MCS methodology, the Subset Simulation (SS) approach was proposed by Au and Beck [3] to calculate the small failure probability. The first step of the SS method is to generate a given number of realizations of the uncertain parameters using the classical MCS technique. The second step is to use the Metropolis-Hastings (M-H) algorithm to generate realizations in the direction of the limit state surface. This step is repeated until reaching the limit state surface. It should be emphasized here that in case of a small failure probability, SS requires the repetition of the second step several times to reach the limit state surface. This leads to a high number of calls of the deterministic model and consequently a high computational time. To reduce the computation time of the SS approach, Defaux et al. [6] proposed a more efficient

approach called "improved Subset Simulation (iSS)". In this method, the first step of the SS approach was replaced by a conditional simulation. In other words, instead of generating realizations directly around the origin by the classical MCS, the realizations are generated outside a hypersphere of a given radius R<sub>h</sub>. This reduces the number of realizations which are not located in the failure zone. Consequently, the number of realizations required to reach the limit state surface is significantly reduced. Notice that Defaux et al. [6] have employed the iSS to calculate the failure probability in the case where the uncertain parameters are modeled by random variables. In the present paper, the iSS is employed in the case where the uncertain parameters are modeled by random fields. This method is illustrated herein through the computation of the probability (P<sub>e</sub>) of exceeding a tolerable differential settlement between two neighboring strip footings resting on a soil with a 2D spatially varying Young's modulus. The Young's modulus is modeled by a random field. The footings are subjected to central vertical loads with equal magnitude. The random field is discretized into a finite number of random variables using the Karhunen-Loeve (K-L) expansion. The differential settlement between the two footings was used to represent the system response. The deterministic model used to compute the system response is based on numerical simulations using the commercial software FLAC.

This paper is organized as follows: a brief review of the SS approach and the Karhunen-Loeve expansion method is first
presented. Then, the iSS approach and its implementation in the case of random field problems are presented. This is followed by the illustration of the efficiency of the iSS approach through the probabilistic analysis of two neighboring strip footings resting on a spatially varying soil. The paper ends with a conclusion of the main findings.

## 2. Review of subset simulation approach

Subset simulation was proposed by Au and Beck [3] to compute the small failure probabilities. The basic idea of the subset simulation approach is that the small failure probability can be expressed as a product of larger conditional failure probabilities. Consider a failure region F defined by the condition $G < 0$ where G is the performance function and let $(s_1, ..., s_k, ..., s_{N_t})$ be a sample of $N_t$ realizations of a vector 's' composed of M random variables. It is possible to define a sequence of nested failure regions $F_1, ..., F_{j_1}, ..., F_m$ of decreasing size where $F_1 \supset ... \supset F_j \supset ... \supset F_m = F$ (Fig. 1). An intermediate failure region $F_j$ can be defined by $G < C_j$ where $C_j$ is an intermediate failure threshold whose value is larger than zero. Thus, there is a decreasing sequence of positive failure thresholds $C_1, ..., C_j, ..., C_m$ corresponding respectively to $F_1, ..., F_{j_1}, ..., F_m$ where $C_1 > ... > C_j > ... > C_m = 0$. In the SS approach, the space of uncertain parameters is divided into a number m of levels with equal number $N_s$ of realizations $(s_1, ..., s_k, ..., s_{N_s})$ where $N_t = N_s \times m$. An intermediate level j contains a safe region and a failure region defined with respect to a given failure threshold $C_j$. The conditional failure probability corresponding to this intermediate level j is calculated as follows:

$$P(F_j|F_{j-1}) = \frac{1}{N_s} \sum_{k=1}^{N_s} I_{F_j}(s_k) \quad (1)$$

where $I_{F_j}(s_k) = 1$ if $s_k \in F_j$ and $I_{F_j}(s_k) = 0$ otherwise. Notice that in the SS approach, the first $N_s$ realizations are generated using MCS methodology according to a target PDF $P_t$. The next $N_s$ realizations of each subsequent level are obtained using Markov chain method based on Metropolis-Hastings (M-H) algorithm according to a proposal PDF $P_p$. Notice that a modified M-H algorithm was proposed by [15]. This modified algorithm was used in this paper to generate the realizations of level j ($j = 1, 2, ..., m$).

The failure probability $P(F) = P(F_m)$ of the failure region F can be calculated from the sequence of conditional failure probabilities as follows:

$$P(F) = P(F_1) \prod_{j=2}^{m} P(F_j|F_{j-1}) \quad (2)$$

<figure><img src="image_2.png" /><figcaption>Fig. 1. Nested failure domain.</figcaption></figure>

where $P(F_1)$ is the failure probability corresponding to the first level of the SS approach, m is the number of levels required to reach the limit state surface and $P(F_j|F_{j-1})$ is an intermediate conditional failure probability. This equation can be regarded as a system consisting of m components (related to the m failure regions $F_1, ..., F_{j_1}, ..., F_m$) connected in parallel. Consequently, the failure probability of the failure region F is the intersection of all conditional failure probabilities of the failure regions $F_1, ..., F_j, ..., F_m$. Thus, the failure probability $P(F)$ is:

$$P(F) = P(\cap_{j=1}^{m} F_j) \quad (3)$$

where

$$\begin{aligned} P(\cap_{j=1}^{m} F_j) &= P(F_m | \cap_{j=1}^{m-1} F_j) \times P(\cap_{j=1}^{m-1} F_j) = P(F_m | F_{m-1}) \times P(\cap_{j=1}^{m-1} F_j) \\ &= \dots = P(F_1) \prod_{j=2}^{m} P(F_j | F_{j-1}) \end{aligned} \quad (4)$$

It should be noticed here that the computation of the failure probability $P(F)$ may be determined using alternatively one of the two following procedures. The first procedure consists in prescribing a sequence of $C_1, ..., C_j, ..., C_m$ so that $C_1 > ... > C_j > ... > C_m = 0$ and then, calculating the different values of $P(F_j/F_{j-1})$ at the different levels using Eq. (1). The second procedure consists in prescribing a constant conditional failure probability $P(F_j/F_{j-1})$ for the different levels and then, calculating the different $C_j$ values corresponding to these levels. The value of $C_j$ of level j is the one for which the ratio between the number of realizations for which $G < C_j$ and the number of realizations $N_s$ of this level (which is identical for the different levels), is equal to the prescribed value $P(F_j/F_{j-1})$. In this paper, the second procedure is used. Notice that, for simplicity in notations, the constant conditional failure probability $P(F_j/F_{j-1})$ will be referred to as $p_0$ later on. The algorithm of the SS approach can be described by the following steps:

(1) Generate a realization of the vector 's' of M random variables by MCS according to the target PDF $P_t$.

(2) Using the deterministic model, calculate the system response corresponding to this realization.

(3) Repeat steps 1 and 2 until obtaining a prescribed number $N_s$ of realizations of the vector 's' and the corresponding system response values. Then, evaluate the corresponding values of the performance function to obtain the vector $G_0 = \{G_0^1, ..., G_0^k, ..., G_0^{N_s}\}$. Notice that the values of the performance function of the different realizations are arranged in an increasing order in the vector $G_0$. Notice also that the subscripts '0' refer to the first level (level 0) of the subset simulation.

(4) Prescribe a constant conditional failure probability $p_0$ for all the failure regions $F_j$ ($j = 1, ..., m-1$) and evaluate the first failure threshold $C_1$ which corresponds to the failure region $F_1$ where $C_1$ is equal to the $[(N_s \times p_0) + 1]^{th}$ value in the increasing list of elements of the vector $G_0$. This ensures that the value of $P(F_1)$ will be equal to the prescribed $p_0$ value.

(5) Among the $N_s$ realizations, there are $[N_s \times p_0]$ ones whose values of the performance function are less than $C_1$ (i.e. they are located in the failure region $F_1$). These realizations are used as 'mother realizations' to generate additional $[(1-p_0)-N_s]$ realizations of the vector 's' using Markov chain method based on Metropolis-Hastings algorithm. These new realizations are located in the second level (level 1 in Fig. 1).

(6) The values of the performance function corresponding to the realizations obtained from the preceding step are listed in an increasing order and are gathered in the vector of performance function values $G_1 = \{G_1^1, ..., G_1^k, ..., G_1^{N_s}\}$.
(7) Evaluate the second failure threshold $C_2$ as the $[(N_s - \times p_0) + 1]^{th}$ value in the increasing list of elements of the vector $G_1$.

(8) Repeat steps 5-7 to evaluate the failure thresholds $C_3, C_4, \dots, C_m$ corresponding to the failure regions $F_3, F_4, \dots, F_m$. Notice that contrary to all other thresholds, the last failure threshold $C_m$ is negative. Thus, $C_m$ is set to zero and the conditional failure probability of the last level $[P(F_m/F_{m-1})]$ is calculated as follows:

$$P(F_m|F_{m-1}) = \frac{1}{N_s} \sum_{k=1}^{N_s} I_{F_m}(s_k) \quad (5)$$

where $I_{F_m}(s_k) = 1$ if the performance function $G(s_k)$ is negative and $I_{F_m}(s_k) = 0$ otherwise.

(9) The failure probability $P(F)$ is evaluated according to Eq. (2).

### 3. Discretization of random field by K-L expansion

For a Gaussian random field $E(X, \theta)$, where $X$ denotes the spatial coordinates and $\theta$ indicates the random nature of this random field, this random field can be approximated by the K-L expansion as follows [16]:

$$E(X, \theta) \approx \mu + \sum_{i=1}^{M} \sqrt{\lambda_i} \phi_i(X) \xi_i(\theta) \quad (6)$$

where $\mu$ is the mean of the random field, $M$ is the size of the series expansion, $\lambda_i$ and $\phi_i$ are the eigenvalues and eigenfunctions of the covariance function, and $\xi_i(\theta)$ is a vector of standard uncorrelated random variables.

In the present paper, the random field $E$ was assumed to follow a log-normal probability density function so that $\ln(E)$ is a normal random field with mean value $\mu_{ln}$ and standard deviation $\sigma_{ln}$. For a log-normal random field, the K-L expansion given in Eq. (6) becomes [5]:

$$E(X, \theta) \approx \exp[\mu_{\ln} + \sum_{i=1}^{M} \sqrt{\lambda_i} \phi_i(X) \xi_i(\theta)] \quad (7)$$

On the other hand, the random field $E$ was assumed to follow an exponential covariance function. This function is given as follows:

$$C[(x_1, y_1), (x_2, y_2)] = \sigma_{\ln}^2 \exp \left( - \frac{|x_1 - x_2|}{l_{\ln x}} - \frac{|y_1 - y_2|}{l_{\ln y}} \right) \quad (8)$$

where $(x_1, y_1)$ and $(x_2, y_2)$ are the coordinates of two arbitrary points in the domain over which the random field is defined and $l_{\ln x}$ and $l_{\ln y}$ are respectively the horizontal and vertical lengths over which the values of $\ln(E)$ are highly correlated. Notice that in the case of an exponential covariance function, the eigenvalues and eigenfunctions are given analytically. Their solutions are presented in Ghanem and Spanos [10]. Notice finally that the choice of the number of terms $M$ of the K-L expansion depends on the desired accuracy of the problem being treated. In the case of a Gaussian random field, the error estimate of the K-L expansion with $M$ terms can be calculated as follows ([18]):

$$err(X) = 1 - (1/\sigma) \sum_{i=1}^{M} \lambda_i \phi_i^2(X) \quad (9)$$

in which $\sigma$ is the standard deviation of the Gaussian random field. In case of a log-normal random field, one should use $\sigma_{\ln}$ instead of $\sigma$ in Eq. (9).

### 4. Improved subset simulation (iSS) approach and its implementation in the case of random fields

The basic idea of the iSS approach is to replace the first step of the SS methodology (i.e. generating realizations directly around

the origin by the classical MCS as shown in Fig. 1) by a conditional simulation [Harbitz [11] and Yonezawa et al. [19]] in which the realizations are generated outside a hypersphere of a given radius $R_h$ as shown in Fig. 2. Based on this conditional simulation, the failure probability $P(F_1)$ corresponding to the first level is calculated as follows [Harbitz [11] and Yonezawa et al. [19]]:

$$P(F_1) = (1 - \chi_M(R_h^2)) \frac{1}{N_s} \sum_{k=1}^{N_s} I_{F_1}(s_k) \quad (10)$$

where $\chi_M$ is the chi-square distribution with $M$ degrees of freedom ($M$ being the number of random variables) and $I_{F_1}(s_k) = 1$ if $s_k \in F_1$ and $I_{F_1}(s_k) = 0$ otherwise.

The advantage of using the conditional simulation is to generate realizations in the proximity of the limit state surface leading to a reduction in the number of realizations required to reach this surface. Notice finally that similar to the classical SS approach, the realizations of the remaining levels of the iSS approach are generated using the Markov chain method based on Metropolis-Hastings algorithm.

As mentioned before, this paper aims at employing the iSS approach for the computation of the failure probability in the case of a spatially varying soil property. To achieve this purpose, a link between the iSS approach and the K-L expansion through the standard normal random variables is performed. This link is somewhat similar to that presented in Ahmed and Soubra [1] and Ahmed and Soubra [2] except the fact that the generation of samples of the first step of the SS approach was performed outside a hypersphere of a given radius. It is presented herein for completeness. Based on this link, the algorithm of the iSS approach for the case of a spatially varying soil property can be described as follows:

(1) Generate a vector of $M$ standard normal random variables $\{\xi_1, \dots, \xi_M, \dots, \xi_M\}$ by MCS methodology. This vector must satisfy the condition that its norm is larger than a prescribed radius $R_h$ of a hypersphere centered at the origin of the standard space.

(2) Substitute the vector $\{\xi_1, \dots, \xi_M, \dots, \xi_M\}$ in the K-L expansion to obtain the first realization of the random field. Then, use the deterministic model to calculate the corresponding system response.

(3) Repeat steps 1 and 2 until obtaining a prescribed number $N_s$ of realizations and their corresponding system response values. Then, evaluate the corresponding values of the performance function to obtain the vector $G_0 = \{G_0^1, \dots, G_0^k, \dots, G_0^{N_s}\}$. Notice that the values of the performance function of the

<figure><img src="image_3.png" /><figcaption>Fig. 2. Samples generation outside a hypersphere of radius $R_h$.</figcaption></figure>
different realizations are arranged in an increasing order in the vector G₀. Notice also that the subscript '0' refers to the first level (i.e. level 0).

(4) Evaluate the first failure threshold C₁ of the failure region F₁ as the $[(N_s \times \bar{p}_0) + 1]^{th}$ value in the increasing list of elements of the vector G₀ where $\bar{p}_0$ is a prescribed value that represents the ratio between the number of realizations for which $G < C_1$ and the number of realizations $N_s$ (i.e. the term $\frac{1}{N_s} \sum_{k=1}^{N_s} I_{F_1}(s_k)$ in Eq. 10). Thus, among the $N_s$ realizations, there are $[N_s \times \bar{p}_0]$ ones whose values of the performance function are less than C₁ (i.e. they are located in the failure region F₁).

(5) Evaluate the conditional failure probability of the first level $P(F_1)$ using Eq. (10).

(6) Prescribe an intermediate constant conditional failure probability $P(F_j|F_{j-1})$ for all the remaining failure regions $F_j$ where $j = 2, 3, ..., m-1$. Although the $P(F_j|F_{j-1})$ value can be arbitrary chosen, it is recommended to be chosen equal to the value of $\bar{p}_0$ used in step 4 to facilitate the implementation of the iSS approach. Notice that, for simplicity in notation, $P(F_j|F_{j-1})$ for $j = 2, 3, ..., m-1$ will be referred to as $p_0$ in the remaining sections of this paper.

(7) The different vectors of random variables {$\bar{\xi}_1, ..., \bar{\xi}_i, ..., \bar{\xi}_M$} corresponding to the realizations that are located in the failure region F₁ (from step 4) are used as 'mother vectors' to generate additional $[(1-p_0) \times N_s]$ vectors of random variables {$\bar{\xi}_1, ..., \bar{\xi}_i, ..., \bar{\xi}_M$} using the Markov chain method based on Metropolis-Hastings algorithm. These new vectors are substituted in the K-L expansion to obtain the corresponding random field realizations. Thus, one obtains the $N_s$ realizations of level 1.

(8) The values of the performance function corresponding to the realizations of level 1 are listed in an increasing order and are gathered in the vector of performance function values $G_1 = \{G_1^1, ..., G_1^k, ..., G_1^{N_s}\}$.

(9) Evaluate the second failure threshold $C_2$ as the $[(N_s \times p_0) + 1]^{th}$ value in the increasing list of elements of the vector $G_1$.

(10) Repeat steps 7 and 8 to evaluate the failure thresholds $C_3, C_4, ..., C_m$ corresponding to the failure regions $F_3, F_4, ..., F_m$ by using each time the vectors of random variables {$\bar{\xi}_1, ..., \bar{\xi}_i, ..., \bar{\xi}_M$} corresponding to the realizations that are located in the failure region $F_j$ as mother vectors to generate the additional vectors in this region. Notice that contrary to all other thresholds, the last threshold $C_m$ is negative. Thus, $C_m$ is set to zero and the conditional failure probability of the last level $[P(F_m/F_m-1)]$ is calculated as:

$$P(F_m|F_{m-1}) = \frac{1}{N_s} \sum_{k=1}^{N_s} I_{F_m}(s_k) \qquad (11)$$

where $I_{F_m} = 1$ if the performance function $G(s_k)$ is negative and $I_{F_m} = 0$ otherwise.

(11) Finally, the failure probability $P(F)$ is evaluated according to Eq. (2) in which $P(F_1)$ is calculated using Eq. (10) and the failure probability of the last level is calculated using Eq. (11).

Notice that, a normal PDF was used as a target probability density function $P_r$. However, a uniform PDF was chosen as a proposal probability density function $P_p$. The intermediate failure probability $p_0$ was chosen equal to 0.1.

## 5. Probabilistic analysis of two neighboring strip footings

The efficiency of the iSS approach is illustrated in this paper through the probabilistic analysis of two neighboring strip footings

resting on a soil with a spatially varying Young's modulus and subjected to equal vertical loads. Indeed, due to the soil spatial variability, the two footings exhibit a differential settlement $\delta$. The differential settlement $\delta$ was used to represent the system response. It is calculated as follows: $\delta = |\delta_1 - \delta_2|$ where $\delta_1$ and $\delta_2$ are the settlements (computed at the footing centers) of the two footings. The Young's modulus was modeled by a random field and it was assumed to follow a log-normal probability density function. Its mean value and coefficient of variation are respectively $\mu_E = 60$ MPa and $COV_E = 15\%$. The random field was discretized using K-L expansion. It was assumed to follow an exponential covariance function. Although an isotropic random field is often assumed in literature, a two-dimensional (2D) anisotropic random field with horizontal and vertical autocorrelation lengths (denoted by $l_{\text{ln x}}$ and $l_{\text{ln y}}$, respectively) was used herein. This is because the horizontal autocorrelation length tends to be larger than the vertical autocorrelation length. A ratio of $l_{\text{ln x}}$ to $l_{\text{ln y}}$ of 1 to 10 for these autocorrelation lengths is usually found in practice [4]. Notice however that a wide range of values of the autocorrelation lengths was considered herein in order to explore some interesting features related to the autocorrelation lengths. The performance function used to calculate the probability $P_e$ of exceeding a tolerable differential settlement is defined as follows:

$$G = \delta_{\max} - \delta \qquad (12)$$

where $\delta_{\max}$ is a prescribed tolerable differential settlement and $\delta$ is the computed differential settlement due to the soil spatial variability.

In the following subsections, the deterministic model used to calculate the differential settlement will be presented. Then, the validation of the iSS approach in the case of random fields will be performed by comparison of its results with those obtained by MCS methodology. Finally, the effect of the autocorrelation lengths on the $P_e$ value in both cases of isotropic and anisotropic random fields will be presented and discussed.

### 5.1. Deterministic model

The deterministic model used to calculate the differential settlement $\delta$ is based on numerical simulations using FLAC. For this computation, two footings (each of width $b = 2$ m) were considered in the analysis (Fig. 3). Each footing is subjected to a central vertical load $P_a = 1000$ kN/m (i.e. a uniform vertical applied pressure $q_a = 500$ kN/m²). The two footings' centers are separated by a distance $D = 4$ m. This small distance was chosen in order to obtain a small soil domain that requires relatively small computation time. The small computation time helps to validate the results obtained by the iSS approach by comparison with those obtained by MCS methodology using a large number of calls of the deterministic model. An optimal non-uniform but symmetrical mesh composed of 1290 zones was employed. For the displacement boundary conditions, the bottom boundary was assumed to be

<figure><img src="image_4.png" /><figcaption>Fig. 3. Soil domain and mesh used in the numerical simulations.</figcaption></figure>
fixed and the vertical boundaries were constrained in motion in the horizontal direction. The soil behavior was modeled by a conventional elastic-perfectly plastic model based on Mohr-Coulomb failure criterion in order to take into account the possible plasticification that may take place near the footing edges even under the service loads. The strip footings were modeled by a linear elastic model. They are connected to the soil via interface elements. The values of the different parameters of the soil, footings and interfaces are given in Table 1 where $\psi$ is the dilation angle, E and v are the Young's modulus and Poisson ratio and $K_n$ and $K_s$ are the normal and shear stiffness of the interface.

In order to calculate the differential settlement for a given random field realization, (i) the coordinates of the center of each element of the mesh were calculated; then, the K-L was used to calculate the value of the Young's modulus at each element center, (ii) geostatic stresses were applied to the soil, (iii) the obtained displacements were set to zero in order to obtain the footings displacements due to only the footings applied loads and finally, (iv) the service loads were applied to the footings and the vertical displacements at the footings centers ($\delta_1$ and $\delta_2$) due to these loads are calculated. The differential settlement is calculated as the absolute difference between $\delta_1$ and $\delta_2$.

## 5.2. Validation of the iSS approach

This section presents a validation of the proposed iSS approach. Notice that for all the probabilistic analyses performed in this paper, the tolerable differential settlement $\delta_{max}$ was assumed equal to $3.5 \times 10^{-3}$ m. Notice also that the horizontal and vertical autocorrelation lengths $l_{ln x}$ and $l_{ln y}$ were normalized with respect to the distance D between the centers of the two footings (i.e. $L_{ln x} = l_{ln x}/D$ and $L_{ln y} = l_{ln y}/D$). The numerical results have shown that this assumption is valid when the ratio $D/b$ is constant where b is the footing width. Notice finally that, for the small values of $L_{ln x}$ and $L_{ln y}$ used in this paper [(($L_{ln x} = L_{ln y} = 1$ in case of isotropic random field) and ($L_{ln x} = 0.5, L_{ln y} = 0.25$ in case of anisotropic random field)], 100 terms of K-L expansion were used. This number was necessary to achieve an error estimate less than 15% (Fig. 4a). For only one case where $L_{ln x} = L_{ln y} = 0.25$, a number of 500 terms was required to achieve such an error estimate (Fig. 4b).

### 5.2.1. Selection of the optimal number $N_s$ of realizations per level of iSS approach

In order to determine the optimal number of realizations $N_s$ to be used per level, different values of $N_s$ were considered. For each $N_s$ value, the failure thresholds $C_1$, $C_2$, etc. were calculated and presented in Table 2 when the radius $R_h$ of the hypersphere is equal to zero (i.e. for the classical SS approach). This table shows that the failure threshold value decreases with the successive levels until reaching a negative value at the last level. Table 3 presents the $P_e$ values and the corresponding values of the coefficient of variation for the different number of realizations $N_s$. As expected, the coefficient of variation of $P_e$ decreases with the increase in the number of realizations $N_s$. It should be noted that for each $N_s$ value presented

<caption>**Table 1**</caption>

Shear strength and elastic properties of soil, footing, and interface.

<table><thead><tr><th>Variable</th><th>Soil</th><th>Footing</th><th>Interface</th></tr></thead><tbody><tr><td>c</td><td>20 kPa</td><td>N/A</td><td>20 kPa</td></tr><tr><td>φ</td><td>30°</td><td>N/A</td><td>30°</td></tr><tr><td>ψ = 2/3 φ</td><td>20°</td><td>N/A</td><td>20°</td></tr><tr><td>E</td><td>60 MPa</td><td>25 GPa</td><td>N/A</td></tr><tr><td>v</td><td>0.3</td><td>0.4</td><td>N/A</td></tr><tr><td>K<sub>h</sub></td><td>N/A</td><td>N/A</td><td>1 GPa</td></tr><tr><td>K<sub>s</sub></td><td>N/A</td><td>N/A</td><td>1 GPa</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Error estimate versus the number of eigenmodes for different values of $L_{ln x}$ and $L_{ln y}$.</figcaption></figure>

in Table 2, $P_e$ corresponding to each level j was calculated by the iSS approach as follows:

$$P(F_j) = P(F_1) \times P(F_2|F_1) \times \dots \times P(F_j|F_{j-1}) \quad (13)$$

These $P_e$ values were compared to those computed by the crude MCS methodology using a number $N = 30,000$ of realizations (Fig. 5). Notice that at a given level j, the $P_e$ value is calculated by MCS methodology as follows:

$$P(F_j) = \frac{1}{N} \sum_{k=1}^{N} I_F(G_k) \quad (14)$$

in which, $G_k$ is the value of the performance function at the kth realization and $I_F = 1$ if $G_k < C_j$ and $I_F = 0$ otherwise. The comparison has shown that for $N_s \ge 1000$ realizations, the $P_e$ value computed by the iSS approach at the different levels is very close to that computed by the crude MCS methodology (Fig. 5a and b). Thus, $N_s = 1000$ realizations will be used in all the probabilistic analyses performed in this paper. Notice that when $N_s = 1000$ realizations, the coefficient of variation of $P_e$ by the iSS approach is $COV_{P_e} = 31.5\%$. A quasi similar value of COV (COV = 31.3%) was obtained by MCS methodology (but when using 30,000 realizations).

As was explained before, the optimal $N_s$ value was determined by comparing the values of $P_e$ obtained by both iSS and MCS. It should be mentioned here that although the computation time of the 30,000 realizations by MCS is significant (about 145 days), this number of realizations remains insufficient to assure an accurate $P_e$ value with a small value of $COV_{P_e}$. As an alternative approach to MCS, one may determine the optimal $N_s$ value by successively increasing $N_s$ and comparing the $P_e$ values given by the iSS approach for each $N_s$ value. The $N_s$ value beyond which $P_e$ converges (i.e. slightly varies with the increase of $N_s$) is the optimal one. In the present analysis, it was found that $P_e$ converges when $N_s = 1000$ realizations. This is because the final $P_e$ values (corresponding to $C=0$) are equal to $3.65 \times 10^{-4}$ and $3.67 \times 10^{-4}$ for $N_s = 1000$ and 1200 realizations, respectively. As a conclusion, this alternative procedure seems to work well for the determination of the optimal $N_s$ value.

### 5.2.2. Selection of the optimal radius $R_h$ of the hypersphere

When $R_h = 0$, four levels were required to reach the limit state surface $G=0$. This means that a total number of $N_t = 1000 + (900 \times 3) = 3700$ realizations were required to calculate $P_e$ with the iSS approach. Thus, for the same accuracy, the number of realizations (and consequently, the computation time) is reduced by 87.7% with respect to MCS when $R_h = 0$ (i.e. when the classical SS is used). This number can again be reduced by increasing $R_h$ (i.e. by using iSS). Table 4 shows that, when $R_h$ increases, the total number of realizations decreases. When $R_h = 11.5$, only two levels are required. Thus, the total number of realizations is $N_t = 1000 + 900 = 1900$ realizations. As a conclusion, the number of realizations (and consequently,
<caption>**Table 2**</caption>

Evolution of the failure threshold with the different levels of the iSS approach and with the number of realizations $N_s$ ($R_h = 0$, $L_{ln x} = 2.5$ and $L_{ln y} = 0.25$).

<table><caption>**Table 3**</caption><thead><tr><th rowspan="2">Failure threshold C<sub>j</sub> for each level j</th><th colspan="6">Number of realizations N<sub>s</sub> per level</th></tr><tr><th>200</th><th>400</th><th>600</th><th>800</th><th>1000</th><th>1200</th></tr></thead><tbody><tr><td>C<sub>1</sub></td><td>0.00191</td><td>0.00199</td><td>0.00189</td><td>0.00204</td><td>0.00191</td><td>0.00195</td></tr><tr><td>C<sub>2</sub></td><td>0.00103</td><td>0.00099</td><td>0.00096</td><td>0.00110</td><td>0.00102</td><td>0.00103</td></tr><tr><td>C<sub>3</sub></td><td>0.00041</td><td>0.00032</td><td>0.00021</td><td>0.00037</td><td>0.00036</td><td>0.00034</td></tr><tr><td>C<sub>4</sub></td><td>-0.00009</td><td>-0.00051</td><td>-0.00036</td><td>-0.00037</td><td>-0.00039</td><td>-0.00038</td></tr></tbody></table>

Values of $P_e$ and $COV_{Pe}$ versus the number $N_s$ of realizations per level.

<table><thead><tr><th rowspan="2"></th><th colspan="6">Number of realizations N<sub>s</sub> per level</th></tr><tr><th>200</th><th>400</th><th>600</th><th>800</th><th>1000</th><th>1200</th></tr></thead><tbody><tr><th>P<sub>e</sub> × (10<sup>-4</sup>)</th><td>1.85</td><td>3.48</td><td>4.63</td><td>2.36</td><td>3.65</td><td>3.67</td></tr><tr><th>COV<sub>Pe</sub></th><td>0.669</td><td>0.505</td><td>0.385</td><td>0.348</td><td>0.315</td><td>0.285</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>(a) $N_s = 1000$ realizations per level</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(b) $N_s = 1200$ realizations per level</figcaption></figure>

**Fig. 5.** Comparison between $P_e$ computed by iSS and that computed by MCS methodology at each level of iSS ($R_h = 0$, $L_{ln x} = 2.5$ and $L_{ln y} = 0.25$).

the computation time) required by the SS approach could be reduced by 48.6% by employing the iSS approach.

It is to be mentioned here that the radius $R_h$ should be carefully chosen. If $R_h$ is very small, the number of levels of the iSS approach will be equal to the number of levels of the classical SS approach and consequently the time cost will remain constant. On the other hand, if $R_h$ is very large, the hypersphere might overlap with the failure region F leading to unsampled area in the failure region which leads to inaccurate value of the failure probability. This issue can be overcome (i) by calculating an approximate value of the failure probability using a simple and fast approach and then (ii) by computing the corresponding approximate value of the radius $R_h$ to be used in the iSS approach.

In this paper, the approximate value $P_{eapp}$ of the probability of exceeding a tolerable differential settlement was calculated by the Collocation-based Stochastic Response Surface Method (CSRSM) using a small number of random variables and a small Polynomial Chaos Expansion (PCE) order [Huang et al. [13] and Huang and Kou [12]]. It should be mentioned here that a high order of the PCE is not necessary herein since an approximate $P_e$ value is sought. A second order PCE was thus used to approximate the system response. Concerning the number of standard normal random variables (number of terms in the K–L expansion), a small number

was firstly selected and the corresponding $P_{eapp}$ value was calculated. This number was then successively increased until $P_{eapp}$ converges to a constant value as shown in Table 5. This table indicates that $P_{eapp}$ converges to a value of $6.46 \times 10^{-4}$ when the number of standard normal random variable is equal to 6. In this case, the number of collocation points is equal to 28 according to the concept of matrix invertibility by Sudret [17].

After the determination of $P_{eapp}$, the corresponding approximate radius $R_{happ}$ of the hypershere can be determined. Notice that $R_{happ}$ represents the distance between the origin and the last failure region corresponding to $C_m$. This means that only the first level of the iSS approach will likely be required to reach the last failure threshold. Thus, $P_{eapp}$ can be supposed equal to $P(F_1)$ and $C_1 = C_m = 0$. Consequently, one obtains:

$$P_{eapp} = P(F_1) = (1 - \chi_M(R_h^2)) \frac{1}{N_s} \sum_{k=1}^{N_s} I_{F_1}(s_k) \quad (15)$$

where $I_{F_1}(s_k) = 1$ if $G(s_k) \le 0$ and $I_{F_1}(s_k) = 0$ otherwise. Since the term $\frac{1}{N_s} \sum_{k=1}^{N_s} I_{F_1}(s_k)$ in Eq. (15) is equal to $\bar{p}_0$ as mentioned before in Section 4 (see step 4), Eq. (15) can be rewritten as follows:

$$\chi_M(R_{happ}^2) = 1 - \frac{P_{eapp}}{\bar{p}_0} \quad (16)$$

from which:

$$R_{happ} = \sqrt{\chi_M^{-1}\left(1 - \frac{P_{eapp}}{\bar{p}_0}\right)} \quad (17)$$

in which $\chi_M^{-1}(\cdot)$ is the inverse of the chi-square CDF. By using Eq. (17), for the case studied herein where $\bar{p}_0 = 0.1$ and $M = 100$, the approximate radius $R_{happ}$ corresponding to the approximate value of $P_{eapp} = 6.46 \times 10^{-4}$ is equal to 11.77. After the determination of $R_{happ}$, the iSS approach can be used with $R_h$ slightly smaller than $R_{happ}$ and $M = 100$ terms to rigorously discretize the random field.

### 5.3. Parametric study

This section aims at presenting a parametric study showing the effect of the autocorrelation lengths on the $P_e$ value in both cases of isotropic and anisotropic random fields.

#### 5.3.1. Effect of the autocorrelation length on $P_e$ in the case of isotropic random field

Fig. 6 shows the effect of the autocorrelation length on the $P_e$ value in the case of an isotropic random field. This figure indicates

<caption>**Table 4**</caption>

Effect of the radius of the hypersphere on the number of realizations required to calculate $P_e$ ($L_{ln x} = 2.5$ and $L_{ln y} = 0.25$).

<table><thead><tr><th rowspan="2"></th><th rowspan="2">MCS</th><th rowspan="2">R<sub>h</sub> = 0 (Classical SS)</th><th colspan="3">iSS</th></tr><tr><th>R<sub>h</sub> = 10</th><th>R<sub>h</sub> = 11</th><th>R<sub>h</sub> = 11.5</th></tr></thead><tbody><tr><th>P<sub>e</sub> (×10<sup>-4</sup>)</th><td>3.40</td><td>3.65</td><td>3.58</td><td>3.36</td><td>3.45</td></tr><tr><th>Number of levels</th><td>-</td><td>4</td><td>3</td><td>3</td><td>2</td></tr><tr><th>Number of realizations</th><td>30000</td><td>3700</td><td>2800</td><td>2800</td><td>1900</td></tr><tr><th>Computation time (minutes)</th><td>210000</td><td>25900</td><td>19600</td><td>19600</td><td>13300</td></tr></tbody></table>
<table><caption>**Table 5**
Effect of the number of standard normal random variables on $P_{e,opp}$ value for the case where $L_{\ln x} = 2.5$ and $L_{\ln y} = 0.25$.</caption><thead><tr><th>Number of standard normal random variables <i>M</i></th><th><i>P</i><sub>e,opp</sub></th><th>Number of collocation points according to Sudret [17]</th></tr></thead><tbody><tr><td>3</td><td>0.00</td><td>10</td></tr><tr><td>4</td><td>8.63 × 10<sup>-5</sup></td><td>15</td></tr><tr><td>5</td><td>8.73 × 10<sup>-9</sup></td><td>21</td></tr><tr><td>6</td><td>6.46 × 10<sup>-4</sup></td><td>28</td></tr><tr><td>7</td><td>5.98 × 10<sup>-4</sup></td><td>36</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Effect of the autocorrelation length on $P_e$ (isotropic random field).</figcaption></figure>

that $P_e$ presents a maximum value when $L_{\ln x} = L_{\ln y} = 1$ (i.e. when the autocorrelation length is equal to the distance between the centers of the two footings). This can be explained by the fact that when the autocorrelation length is very small, one obtains a highly heterogeneous soil in both the vertical and the horizontal directions with a great variety of high and small values of the Young's modulus beneath the footings (Fig. 7a). In this case, the soil under the footings contains a mixture of stiff and soft soil zones. Due to the high rigidity of the footings, their movements are resisted by the numerous stiff soil zones in the soil mass; the numerous soft soil zones being of negligible effect on the footings displacements. This leads to small values of the footings displacements (i.e. to a small differential settlement) and thus, to a small value of $P_e$. On the other hand, when the autocorrelation length is large, the soil tends to be homogenous (Fig. 7b). This means that the differential settlement tends to be very small (close to zero) which leads to a very small value of $P_e$. For the intermediate values of the autocorrelation length, there is a high probability that one footing rests on

<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Grey-scale representation of the random field to show the effect of $L_{\ln x} = L_{\ln y}$ in case of isotropic random field.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Effect of $L_{\ln x}$ on $P_e$ when $L_{\ln y} = 0.25$.</figcaption></figure>

a stiff soil zone and the other on a relatively soft soil zone (Fig. 7c). This leads to a high differential settlement and thus to a high $P_e$ value. In this case, $P_e$ presents a maximum.

### 5.3.2. Effect of the autocorrelation lengths on $P_e$ in the case of anisotropic random field

Fig. 8 shows the effect of $L_{\ln x}$ on $P_e$ when $L_{\ln y} = 0.25$. This figure shows that $P_e$ presents a maximum value when the autocorrelation length is equal to the distance between the two footings centers (i.e. when $L_{\ln x} = 1$). This observation can be explained as follows:

For the very small values of $L_{\ln x}$ compared to $L_{\ln y}$, one obtains a vertical multilayer composed of thin sub-layers where each sub-layer may have a high or a small value of the Young's modulus (Fig. 9a). The sub-layers with high values of the Young's modulus prevent the movements of both footings and thus lead to a small value of $P_e$. On the other hand, when $L_{\ln x}$ is very large compared to $L_{\ln y}$, one obtains a horizontal multilayer (i.e. the soil tends to the case of a one-dimensional vertical random field) for which each sub-layer may have a high or a small value of the Young's modulus (Fig. 9b). This leads to the same displacement for both footings and thus to a very small value of $P_e$. Finally, for intermediate values of $L_{\ln x}$ (cf. Fig. 9c), the horizontally extended stiff layers of Fig. 9b become less extended leading to a high probability that the footings rest on soil zones with different values of the Young's modulus. This leads to a greater differential settlement and consequently a greater value of $P_e$.

The effect of $L_{\ln y}$ is presented in Fig. 10 when $L_{\ln x} = 2.5$. This figure also presents the $P_e$ value corresponding to the case of a one-dimensional horizontal random field with $L_{\ln x} = 2.5$. In this case, the soil was considered to be spatially varying only in the horizontal

<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Grey-scale representation of the random field to show the effect of $L_{\ln x}$ in case of anisotropic random field.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Effect of $L_{in y}$ on $P_e$ when $L_{in x} = 2.5$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Grey-scale representation of the random field to show the effect of $L_{in y}$ in case of anisotropic random field.</figcaption></figure>

direction while it was considered to be homogeneous in the vertical direction. Fig. 10 shows that the $P_e$ value increases with the increase in $L_{in y}$. This can be explained as follows: when $L_{in y}$ is very small, the two footings rest on a horizontal multilayer composed of thin sub-layers where each sub-layer may have a high or a small value of the Young's modulus (Fig. 11a). This means that $\delta_1$ and $\delta_2$ are almost equal. Thus, the differential settlement $\delta$ is very small which results in a small value of $P_e$. On the other hand, when $L_{in y}$ is very large, the soil tends to the case of a one-dimensional horizontal random field. In this case, one obtains vertically extended stiff sub-layers adjacent to vertically extended soft sub-layers (Fig. 11b). For the chosen value of $L_{in x}$, there is a high probability that one footing rests on a vertical stiff layer and the other one rests on a vertical soft layer which leads to a high differential settlement and thus to a great value of $P_e$.

## 6. Conclusion

This paper presents an efficient method, called improved subset simulation (iSS), to perform a probabilistic analysis of geotechnical structures that involve spatial variability of the soil properties. This method is an improvement of the subset simulation approach presented by Ahmed and Soubra [2] in the case of a spatially varying soil. It allows one to calculate the small failure probabilities using a reduced number of calls of the deterministic model. This was made possible by generating the samples of the first step of the SS approach outside a hypersphere whose radius was determined by a simple fast and approximate approach. The iSS approach was illustrated through the probabilistic analysis at SLS of two neighboring strip footings resting on a soil with spatially varying Young's modulus. The differential settlement between the two footings was used to represent the system response. The probability $P_e$ (i.e. the probability of exceeding a tolerable differential settlement) calculated by

the improved subset simulation approach was found very close to that computed by Monte Carlo Simulation methodology or the subset simulation approach with a significant reduction in the number of calls of the deterministic model. The use of the iSS approach has reduced the number of calls of the deterministic model by about 50% with respect to the SS approach by Ahmed and Soubra [2].

A parametric study to investigate the effect of the autocorrelation lengths on $P_e$ in both cases of isotropic and anisotropic random fields has shown that:

(1) In case of an isotropic random field, the probability $P_e$ of exceeding a tolerable differential settlement presents a maximum value when the autocorrelation length is equal to the distance $D$ that separates the two footings centers.

(2) In case of an anisotropic random field, $P_e$ significantly increases with the increase of the vertical autocorrelation length (for a given value of the horizontal autocorrelation length) and then, it attains an asymptote which corresponds to the case of a horizontal one-dimensional random field. On the other hand, for a given value of the vertical autocorrelation length, $P_e$ presents a maximum when the horizontal autocorrelation length is equal to the distance $D$ that separates the two footings centers.

## References

[1] Ahmed A, Soubra A-H. Subset simulation and its application to a spatially random soil. Geotechnical risk assessment and management, GeoRisk 2011. Atlanta, Georgia: ASCE; 2011. p. 209-216.

[2] Ahmed A, Soubra A-H. Probabilistic analysis of strip footings resting on a spatially random soil using subset simulation approach. Georisk Assessm Manage Risk Eng Syst Geohazard 2012;6(3):188-201.

[3] Au SK, Beck JL. Estimation of small failure probabilities in high dimensions by subset simulation. Probab Eng Mech 2001;16:263-77.

[4] Baecher GB, Christian JT. Reliability and statistics in geotechnical engineering. New York: Wiley; 2003.

[5] Cho SE, Park HC. Effect of spatial variability of cross-correlated soil properties on bearing capacity of strip footing. Int J Numer Anal Meth Geomech 2010;34:1-26.

[6] Defaux G, Meister E, Pendola M. (2010). Mécanique probabiliste et intégrité des cuves des réacteurs nucléaires. JFMS'10, Toulouse, France, p. 23.

[7] Fenton GA, Griffiths DV. Probabilistic foundation settlement on a spatially random soil. J Geotech Geoenviron Eng ASCE 2002;128(5):381-90.

[8] Fenton GA, Griffiths DV. Bearing capacity prediction of spatially random c-φ soils. Can Geotech J 2003;40:54-65.

[9] Fenton GA, Griffiths DV. Three-dimensional probabilistic foundation settlement. J Geotech Geoenviron Eng ASCE 2005;131(2):232-9.

[10] Ghanem R, Spanos PD. Stochastic finite elements – a spectral approach. New York: Springer; 1991.

[11] Harbitz A. An efficient sampling method for probability of failure calculation. Struct Saf 1986;3(2):109-15.

[12] Huang S, Kou X. An extended stochastic response surface method for random field problems. Acta Mech Sin 2007;23:445-50.

[13] Huang S, Mahadervan S, Rebba R. Collocation-based stochastic finite element analysis for random field problems. Probab Eng Mech 2007;22:194-205.

[14] Popescu R, Deodatis G, Nobahar A. Effect of random heterogeneity of soil properties on bearing capacity. Probab Eng Mech 2005;20:324-41.

[15] Santoso AM, Phoon KK, Quek ST. Modified Metropolis-Hastings algorithm with reduced chain-correlation for efficient subset simulation. Probabilistic Engineering Mechanics 2011;26(2):331-41.

[16] Spanos PD, Ghanem R. Stochastic finite element expansion for random media. J Eng Mech ASCE 1989;115(5):1035-53.

[17] Sudret B. Global sensitivity analysis using polynomial chaos expansion. Reliabil Eng Syst Saf 2008;93:964-79.

[18] Sudret B, Berveiller M. Stochastic finite element methods in geotechnical engineering. In: Phoon KK, editor. Reliability-based design in geotechnical engineering: computations and applications. Taylor & Francis; 2008. p. 260-96.

[19] Yonezawa M, Okuda S, Park YT. Structural reliability estimation based on simulation within limited sampling region. Int J Prod Econ 1999;60(61):607-12.