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

HAL
open science

# Reliability-Based Analysis of Strip Footings Using
Response Surface Methodology

Dalia S. Youssef Abdel Massih, Abdul-Hamid Soubra

► To cite this version:

Dalia S. Youssef Abdel Massih, Abdul-Hamid Soubra. Reliability-Based Analysis of Strip Footings Using Response Surface Methodology. International Journal of Geomechanics, 2008, 8 (2), pp.134-143. 10.1061/(ASCE)1532-3641(2008)8:2(134) . hal-01007246

HAL Id: hal-01007246

https://hal.science/hal-01007246

Submitted on 17 Feb 2017

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Reliability-Based Analysis of Strip Footings Using Response Surface Methodology

Dalia S. Youssef Abdel Massih<fnref n="1" /> and Abdul-Hamid Soubra<fnref n="2" />

**Abstract** A reliability-based analysis of a strip foundation subjected to a central vertical load is presented. Both the ultimate and the serviceability limit states are considered. Two deterministic models based on numerical simulations are used. The first one computes the ultimate bearing capacity of the foundation and the second one calculates the footing displacement due to an applied load. The response surface methodology is utilized for the assessment of the Hasofer-Lind reliability indexes. Only the soil shear strength parameters are considered as random variables while studying the ultimate limit state. Also, the randomness of only the soil elastic properties is taken into account in the serviceability limit state. The assumption of uncorrelated variables was found to be conservative in comparison to the one of negatively correlated variables. The failure probability of the ultimate limit state was highly influenced by the variability of the angle of internal friction. However, for the serviceability limit state, the accurate determination of the uncertainties of the Young's modulus was found to be very important in obtaining reliable probabilistic results. Finally, the computation of the system failure probability involving both ultimate and serviceability limit states was presented and discussed.

**keywords** Shallow foundations; Bearing capacity; Foundation settlement; Serviceability; Simulation models; System reliability; Footings.

## Introduction

The commonly used approaches in the analysis and design of foundations are deterministic. The average values of the input parameters are usually considered and the uncertainties of the different parameters are taken into account via a global factor of safety which is essentially a “factor of ignorance.” A reliability-based approach for the analysis of foundations is more rational since it enables one to consider the inherent uncertainty of each input parameter. Nowadays, this is possible because of the improvement in our knowledge of the statistical properties of soil (Phoon and Kulhawy 1999).

In this paper, a reliability-based analysis of a strip foundation resting on a $c-\varphi$ soil and subjected to a central vertical load is presented. Previous investigations on the reliability analysis of foundations focused on either the ultimate or the serviceability limit state (Bauer and Pula 2000; Cherubini 2000; Griffiths and Fenton 2001; Griffiths et al. 2002; Low and Phoon 2002; Fenton and Griffiths 2002, 2003, 2005; Popescu et al. 2005; Przewlocki 2005; Youssef Abdel Massih et al. 2007). This paper considers both limit states in the analysis of foundations. Two deterministic

models based on the Lagrangian explicit finite difference code FLAC³D are used. The first one computes the ultimate bearing capacity of the foundation and the second one calculates the footing displacement due to an applied service load. The response surface methodology is utilized to find an approximation of the analytically unknown performance functions and the corresponding reliability indexes. The random variables considered in the analysis are the soil shear strength parameters $c$ and $\varphi$ for the ultimate limit state, and the soil elastic properties $E$ and $v$ for the serviceability limit state. After a brief description of the basic concepts of the theory of reliability, the two deterministic models based on numerical FLAC³D simulations are presented. Then, the probabilistic analysis and the corresponding numerical results are presented and discussed.

## Basic Reliability Concepts

Two different measures are commonly used in literature to describe the reliability of a structure: The reliability index and the failure probability.

The reliability index of a geotechnical structure is a measure of the safety that takes into account the inherent uncertainties of the input variables. A widely used reliability index is the Hasofer and Lind (1974) index defined as the shortest distance from the mean value point of the random variables to the limit state surface in units of directional standard deviations, namely $\beta = \min[R(\theta)/r(\theta)]$ (Fig. 1). Its matrix formulation is (Ditlevsen 1981)

$$ \beta_{HL} = \min_{x \in F} \sqrt{(x - \mu)^T C^{-1} (x - \mu)} \quad (1) $$

in which $x$=vector representing the $n$ random variables; $\mu$=vector of their mean values; and $C$=their covariance matrix. The minimization of Eq. (1) is performed subject to the constraint

<footnote n="1">¹Ph.D. Student, Univ. of Nantes & Lebanese Univ., BP 11-5147, Beirut, Lebanon. E-mail: dalia.youssef@univ-nantes.fr</footnote>

<footnote n="2">²Professor, Institut de Recherche en Génie Civil et Mécanique, Univ. of Nantes, UMR CNRS 6183, Bd. de l'Université BP 152, 44603 Saint-Nazaire Cedex, France (corresponding author). E-mail: abed.soubra@univ-nantes.fr</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 1. Design point and equivalent normal dispersion ellipses in space of two random variables</figcaption></figure>

$G(x) \le 0$ where the limit state surface $G(x)=0$ separates the n-dimensional domain of random variables into two regions: a failure region $F$ represented by $G(x) \le 0$ and a safe region given by $G(x) > 0$.

The classical approach for computing $\beta_{HL}$ by Eq. (1) is based on the transformation of the limit state surface into the space of standard normal uncorrelated variates. The shortest distance from the transformed failure surface to the origin of the reduced variates is the reliability index $\beta_{HL}$.

An intuitive interpretation of the reliability index was suggested in Low and Tang (1997a, 2004) where the concept of an expanding ellipse (Fig. 1) led to a simple method of computing the Hasofer-Lind reliability index in the original space of the random variables. When there are only two uncorrelated nonnormal random variables $x_1$ and $x_2$, these variables span a two-dimensional random space, with an equivalent one-sigma dispersion ellipse [corresponding to $\beta_{HL}=1$ in Eq. (1) without the min] centered at the equivalent normal mean values $(\mu_1^N, \mu_2^N)$ and whose axes are parallel to the coordinate axes of the original space. For correlated variables, a tilted ellipse is obtained. Low and Tang (1997a, 2004) reported that the Hasofer-Lind reliability index $\beta_{HL}$ may be regarded as the codirectional axis ratio of the smallest ellipse (which is either an expansion or a contraction of the $1-\sigma$ ellipse) that just touches the limit state surface to the $1-\sigma$ dispersion ellipse. They also stated that finding the smallest ellipsoid that is tangent to the limit state surface is equivalent to finding the most probable failure point.

From the first-order reliability method FORM and the Hasofer-Lind reliability index $\beta_{HL}$, one can approximate the failure probability as follows

$$P_f \approx \Phi(-\beta_{HL}) \quad (2)$$

where $\Phi(\cdot)$=cumulative distribution function of a standard normal variable. In this method, the limit state function is approximated by a hyperplane tangent to the limit state surface at the design point.

## Ellipsoid Approach via Matlab

Low and Tang (1997a, 2004) showed that the minimization of the Hasofer-Lind reliability index can be efficiently carried out in the original space of the random variables. When the random variables are non-normal and correlated, the optimization approach

uses the Rackwitz-Fiessler equations to compute the equivalent normal mean $\mu^N$ and the equivalent normal standard deviation $\sigma^N$ without the need to diagonalize the correlation matrix, as shown in Low and Tang (2004) and Low (2005). Furthermore, the iterative computations of the equivalent normal mean $\mu^N$ and equivalent normal standard deviation $\sigma^N$ for each trial design point are automatic during the constrained optimization search.

In the present paper, by the Low and Tang method, one literally sets up a tilted ellipsoid in Matlab software and uses the “fmincon” command, built in the optimization tool of this software, to minimize the dispersion ellipsoid subject to the constraint that it be tangent to the limit state surface. Eq. (1) may be rewritten as (Low and Tang 1997b, 2004)

$$\beta = \min_{x \in F} \sqrt{\left[ \frac{x - \mu_x^N}{\sigma_x^N} \right]^T [R]^{-1} \left[ \frac{x - \mu_x^N}{\sigma_x^N} \right]} \quad (3)$$

in which $[R]^{-1}$=inverse of the correlation matrix. This equation will be used to set up the ellipsoid in Matlab since the correlation matrix $[R]$ displays the correlation structure more explicitly than the covariance matrix $[C]$.

## Deterministic Numerical Modeling of Bearing Capacity and Displacement of Strip Footings Using FLAC³D

FLAC³D (Fast Lagrangian Analysis of Continua 1993) is a commercially available three-dimensional finite difference code in which an explicit Lagrangian calculation scheme and a mixed discretization zoning technique are used. This code includes an internal programming option (FISH) which enables the user to add his own subroutines.

In this software, although a static (i.e. nondynamic) mechanical analysis is required, the equations of motion are used. The solution to a static problem is obtained through the damping of a dynamic process by including damping terms that gradually remove the kinetic energy from the system.

The calculation scheme invokes the equations of motion in their discretized forms to derive new velocities and displacements from stresses and forces. Then, strain rates are derived from velocities, and new stresses from strain rates. The stresses and deformations are calculated at several small timesteps (called hereafter cycles) until a steady state of static equilibrium or plastic flow is achieved. The convergence to this state may be controlled by a maximal prescribed value of the unbalanced force for all elements of the model. It should be mentioned that the application of displacements or stresses on a system creates unbalanced forces in this system. Damping is introduced in order to remove these forces or to reduce them to very small values compared to the initial ones.

## Numerical Simulations

This section focuses on the computation of the ultimate bearing capacity of the soil [ultimate limit state (ULS)] and the footing vertical displacement [serviceability limit state (SLS)] due to a central vertical footing load. Although a random soil is studied in this paper, a symmetrical velocity field is considered in both the ULS and the SLS. This is because the soil properties are modeled as random variables. Thus, each FLAC³D simulation considers a homogeneous soil. The randomness of the soil is taken into ac-
<figure><img src="image_4.png" /><figcaption>Fig. 2. Soil domain and mesh used in FLAC³D</figcaption></figure>

count from one simulation to another. A nonsymmetrical velocity field is necessary only for the computation of the reliability of a foundation resting on a spatially variable soil (i.e., where $c$ or $\varphi$ are considered as random processes).

## Ultimate Limit State—Bearing Capacity

This section focuses on the determination of the ultimate bearing capacity of a rough rigid strip footing, of breadth $B=2$ m, resting on a $c-\varphi$ soil.

Because of symmetry, only half of the entire soil domain of width $20B$ and depth $5B$ is considered. The bottom and right vertical boundaries are placed far enough from the footing and they do not disturb the soil mass in motion (i.e., velocity field) for all the soil configurations studied in this paper. A nonuniform mesh composed of 904 zones is used (Fig. 2). The region under the right half of the footing was divided horizontally into 15 zones, whose size gradually decreases from the center to the edge of the footing where very high stress gradients are developed. Beyond the edge of the footing, the domain was divided into 30 zones whose size increases gradually from the foundation edge to the right vertical boundary. Vertically, the domain was divided into 20 zones whose size decreases gradually from the bottom of the domain to the ground surface.

Since this is a two-dimensional (2D) case, all displacements in the direction parallel to the footing are fixed. For the displacement boundary conditions, the bottom boundary was assumed to be fixed and the vertical boundaries were constrained in motion in the horizontal direction.

A conventional elastic-perfectly plastic model based on the Mohr–Coulomb failure criterion is adopted to represent the soil. The soil elastic properties employed are the shear modulus $G=23$ MPa and the bulk modulus $K=50$ MPa (for which the equivalent Young's modulus and Poisson's ratio are, respectively, $E=60$ MPa and $\nu=0.3$). The values of the soil shear strength parameters used in the analysis are: $\varphi=30^\circ$, $\psi=20^\circ$, and $c=20$ kPa, where $\psi$=soil dilation angle. The soil unit weight was taken equal to $18$ kN/m$^3$. Notice that the soil elastic properties have a negligible effect on the failure load. A strip footing of half width equal to $1$ m and depth $0.5$ m is used in the analysis. It is divided horizontally into four zones. The footing is simulated by a weightless elastic material. Its elastic properties are the Young's modulus $E=25$ GPa and the Poisson's ratio $\nu=0.4$. Compared to the soil elastic properties, these values are well in excess of those of the soil and ensure a rigid behavior of the footing. The footing is connected to the soil via interface elements that follow Coulomb's law. The interface is assumed to have a friction angle equal to the soil angle of internal friction, dilation equal to that of

the soil, and cohesion equal to the soil cohesion in order to simulate a perfectly rough soil-footing interface. Normal stiffness $K_n=10^9$ Pa/m and shear stiffness $K_s=10^9$ Pa/m are assigned to this interface. These parameters do not have a major influence on the failure load.

For the computation of the ultimate bearing capacity of a rigid rough strip footing subjected to a central vertical load using FLAC³D, a displacement control method is adopted in this paper. The following procedure is performed before any simulation of the foundation loading: Geostatic stresses are first applied to the soil, then several cycles are run in order to arrive at a steady state of static equilibrium, and finally the obtained displacements are set to zero in order to obtain the footing displacement due to only the footing load.

## Displacement Control Method

In this method, a controlled downward vertical velocity (i.e., displacement per timestep) is applied to the nodes of the footing. Damping of the system is introduced by running several cycles until a steady state of plastic flow is developed in the soil underneath the footing. This state is achieved when both conditions: (1) a constant footing load and (2) small values of unbalanced forces, were satisfied as the number of cycles increases. The number of cycles required to reach this state depends on the value of the applied velocity. At each cycle, the vertical footing load is obtained by using a FISH function that calculates the integral of the normal stress components for all elements in contact with the footing. The value of the vertical footing load at the plastic steady state is the ultimate footing load. The ultimate bearing capacity is then obtained by dividing this load by the footing area.

Two control parameters (the intensity of the vertical velocity and the mesh size) may greatly affect the value of the ultimate footing load. They are examined in the following sections:

An optimal vertical velocity must be chosen in order to reach a value of the ultimate bearing capacity close to the smallest most critical one (corresponding to very small velocity) with a reasonable computation time. A velocity of $2.5 \times 10^{-6}$ m/timestep downward was suggested by Yin et al. (2001) as a result of a number of verification runs. This value was tested in the present paper, and an ultimate load of $2,393.1$ kN/m was obtained at the plastic steady state after 215,000 cycles. This load corresponds to a continuous increase of the footing displacement. A smaller velocity of $10^{-6}$ m/timestep and a higher velocity of $5 \times 10^{-6}$ m/timestep were also tested. The value of the ultimate load corresponding to the smaller velocity was found equal to $2,392.7$ kN/m which is slightly smaller (i.e., more critical) than the one obtained by applying the $2.5 \times 10^{-6}$ m/timestep velocity. However, 380,205 cycles were required to achieve this value (i.e., an increase in the calculation time by 76%). For the higher velocity of $5 \times 10^{-6}$ m/timestep, a slightly greater value of $2,394.48$ kN/m was obtained (Fig. 3). The difference is smaller than 0.1% from the value obtained using the $10^{-6}$ m/timestep velocity. The necessary number of cycles to reach this value was about 107,743 which is significantly smaller than the 215,000 and the 380,205 cycles required by the two smaller velocities. Thus, the use of a vertical velocity of $5 \times 10^{-6}$ m/timestep highly reduces the computation time with a negligible deterioration in the accuracy of the solution. In this paper, this velocity is adopted for all subsequent calculations.

The effect of the mesh size on the solution was also checked. It was found that a more refined mesh under the footing does not improve the value of the ultimate footing load and may cause
<figure><img src="image_5.png" /><figcaption>Fig. 3. Load-displacement curve from displacement control method</figcaption></figure>

numerical instability. Also, a more refined mesh beyond the edge of the footing (40 zones instead of 30 horizontally and 30 zones instead of 20 in the vertical direction) improves the result (i.e., reduces the ultimate load) by only 0.24% with an increase in the calculation time by 33%. Thus, the mesh presented above will be used in all subsequent calculations.

In order to confirm the accuracy of the ultimate bearing capacity obtained by the displacement control method, incremental vertical stresses are applied to the nodes situated at the base of the footing until failure is reached. For each stress increment, damping is introduced until a steady state is obtained. This method called the load control method is found to give a closely similar result of 2,394.44 kN/m (Fig. 4) in comparison to the value of 2,394.48 kN/m obtained by the displacement control method. However, this approach is less efficient regarding the computation time. In order to compare the number of cycles required by the two methods for a given displacement of the footing, the total number of cycles was found to be about 648,000 for a vertical

<figure><img src="image_6.png" /><figcaption>Fig. 4. Load-displacement curve from load control method</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 5. (a) Normal (σ); (b) shear (τ) stresses at right half of footing base as obtained from two methods</figcaption></figure>

footing displacement of 45 cm in the load control method. The corresponding number of cycles was about 107,000 in the displacement control method.

The contact normal and shear stress distributions along the soil- footing interface as obtained by the two methods at failure are presented in Figs. 5(a and b). They show nearly identical results. Except at the footing edge which is a singular point, a quasi-uniform normal stress distribution was observed [Fig. 5(a)]. For the shear stress distribution, gradually increasing stress from the center to the edge of the footing was noticed [Fig. 5(b)]. As for the normal stress distribution, high stresses were observed at the footing edge due to the singularity at this point.

For computation of the ultimate bearing capacity of a rough rigid footing, the displacement control method was found to be the most simple and efficient one regarding the computation time. It will be used in this paper for all subsequent calculations.

## Serviceability Limit State—Vertical Displacement

For the computation of the vertical displacement of a rigid footing under an applied vertical load, it would not be interesting to apply uniform stresses directly to the surface nodes of the soil since this approach corresponds to the simulation of a flexible footing. Thus, as in the ultimate limit state, modeling the foundation by a weightless elastic material is also adopted here. An elastic-perfectly plastic model is used for the soil since it enables the development of plastic zones that may occur near the footing edges even at small service loads and it leads to more accurate solutions than a purely elastic model. The same procedure described before concerning the geostatic stresses is used here. A
<figure><img src="image_8.png" /><figcaption>Fig. 6. Vertical footing displacement versus number of cycles due to applied load $P_S=750$ kN/m</figcaption></figure>

uniform service stress is applied at the base of the footing. Damping of the system is introduced by running several cycles until a steady state of static equilibrium is reached in the soil. This state is achieved when both conditions: (1) a constant vertical displacement of the footing (Fig. 6) and (2) small values of unbalanced forces, were satisfied as the number of cycles increases.

## Reliability Analysis of Strip Footings

The aim of this paper is to perform a reliability analysis of a strip footing resting on a $c-\phi$ soil and subjected to a central vertical load. Two failure or unsatisfactory performance modes are considered in the analysis: The first one involves the ultimate limit state and emphasizes the ultimate bearing capacity of the footing and the second one considers the serviceability limit state and focuses on the maximal footing displacement. The two deterministic models presented in the previous section are used. The response surface methodology is employed to find an approximation of the analytically unknown performance functions. The cohesion $c$, the angle of internal friction $\phi$, the Young's modulus $E$, and the Poisson's ratio $\nu$ of the soil are considered as random variables. Due to the relatively low effect of the elastic modulus $E$ and the Poisson's ratio $\nu$ on the ultimate bearing capacity, only $c$ and $\phi$ will be considered as random variables while studying the ultimate limit state. Similarly, only the randomness of $E$ and $\nu$ will be taken into consideration in the analysis of the serviceability limit state. After a brief description of the performance functions used in the present analysis, the response surface methodology and its numerical implementation are presented. Then, the probabilistic numerical results based on this method are presented and discussed.

## Performance Functions

Two performance functions are used in this reliability analysis. The first one is defined with respect to the ultimate bearing capacity of the soil. It is given as follows

$$G = P_u/P_S - 1 \quad (4)$$

where $P_u$=ultimate foundation load calculated by $FLAC^{3D}$ and $P_S$=applied footing load. The second performance function, defined with respect to a prescribed admissible footing displacement, is given as follows

$$G = u_{\max} - u \quad (5)$$

where $u$=vertical displacement of the footing calculated by $FLAC^{3D}$ due to a service load $P_S$; and $u_{\max}$=maximal admissible vertical displacement.

## Response Surface Method

If the performance function is an explicit function of the random variables, the reliability index can be calculated easily. In the $FLAC^{3D}$ model, the closed form solution of the performance function is not available. Thus, the determination of the reliability index is not straightforward. An algorithm based on the response surface methodology proposed by Tandjiria et al. (2000) is used in this paper with the aim to calculate the reliability index and the corresponding design point. The basic idea of this method is to approximate the performance function by an explicit function of the random variables, and to improve the approximation via iterations. The approximate performance function used in this study has a quadratic form. It uses a second-order polynomial with squared terms but no cross terms. The expression of this approximation is given by

$$G(x) = a_0 + \sum_{i=1}^{n} a_i \cdot x_i + \sum_{i=1}^{n} b_i \cdot x_i^2 \quad (6)$$

where $x_i$=random variables; $n$=number of the random variables; and $(a_i, b_i)$=coefficients to be determined. In this paper, two random variables are considered for each limit state (i.e., $n=2$). They are characterized by their mean values $\mu_i$ and their coefficients of variation $\sigma_i$. A brief explanation of the algorithm used is as follows:

1. Evaluate the performance function $G(x)$ at the mean value point $\mu$ and the $2n$ points each at $\mu \pm k\sigma$ where $k=1$ in this paper;

2. The above $2n+1$ values of $G(x)$ can be used to solve Eq. (6) for the coefficients $(a_i, b_i)$. This obtains a tentative response surface function;

3. Solve Eq. (1) to obtain a tentative design point and a tentative $\beta_{HL}$ subject to the constraint that the tentative response surface function of step 2 be equal to zero; and

4. Repeat steps 1–3 until convergence. Each time step 1 is repeated, the $2n+1$ sampled points are centered at the new tentative design point of step 3.

## Numerical Implementation of Response Surface Method

As described in the previous section, the determination of the Hasofer-Lind reliability index requires: (1) the determination of the coefficients $(a_i, b_i)$ of the tentative response surface via the resolution of Eq. (6) for the $2n+1$ sampled points; and (2) the minimization of the Hasofer-Lind reliability index subject to the constraint that the tentative response surface function be equal to zero. These two operations which constitute a single iteration
were done using the optimization toolbox available in Matlab 7.0 software. Several iterations were performed until convergence of the Hasofer-Lind reliability index.

Notice that the determination of the performance function at the 2n+1 sampled points was performed using deterministic FLAC³D calculations. The results of these computations constitute input parameters for the determination of the coefficients ($a_i, b_i$) of the tentative response surface using Matlab 7.0. Also, the value of the design point determined using the minimization procedure in Matlab 7.0 is an input parameter for the determination of the performance function at the 2n+1 sampled points in FLAC³D. Therefore, an exchange of data between FLAC³D and Matlab 7.0 in both directions was necessary to enable an automatic resolution of the iterative algorithm for the determination of the Hasofer-Lind reliability index. The link between FLAC³D and Matlab 7.0 was performed using text files and FISH commands.

## Numerical Results

For the ultimate limit state, different values of the coefficients of variation of the angle of internal friction and cohesion are presented in the literature. For most soils, the mean value of the effective angle of internal friction is typically between 20 and 40°. Within this range, the corresponding coefficient of variation as proposed by Phoon and Kulhawy (1999) is essentially between 5 and 15%. For the effective cohesion, the coefficient of variation (COV) varies between 10 and 70% (Cherubini 2000). For the coefficient of correlation, Harr (1987) has shown that a correlation exists between the effective cohesion $c$ and the effective angle of internal friction $\phi$. The results of Wolff (1985) ($\rho_{c,\phi} = -0.47$), Yuceman et al. (1973) ($-0.49 \le \rho_{c,\phi} \le -0.24$), Lumb (1970) ($-0.7 \le \rho_{c,\phi} \le -0.37$), and Cherubini (2000) ($\rho_{c,\phi} = -0.61$) are among the ones cited in the literature. In this paper, the illustrative values used for the statistical moments of the shear strength parameters and their coefficient of correlation $\rho_{c,\phi}$ are given as follows: $\mu_c = 20$ kPa, $\mu_\phi = 30°$, $COV_c = 20\%$, $COV_\phi = 10\%$, and $\rho_{c,\phi} = -0.5$. These values are within the range of values cited above. For the probability distribution of the random variables, $c$ is assumed to be lognormally distributed while $\phi$ is assumed to be bounded and a beta distribution is used (Fenton and Griffiths 2003). The parameters of the beta distribution are determined from the mean value and standard deviation of $\phi$. It should be mentioned that the soil elastic properties (i.e., $K$ and $G$ or $E$ and $\nu$) considered as deterministic in the present ultimate limit state have no effect on the value of the ultimate bearing capacity. Higher values of these properties, $G=100$ MPa and $K=133$ MPa (for which $E=240$ MPa and $\nu=0.2$), were checked. No change was observed in the value of the ultimate bearing capacity. Furthermore, a reduction by 50% in the number of cycles necessary to reach failure was noticed (i.e., a reduction in the computation time by half). Consequently, these values will be used in all subsequent calculations when studying the ultimate limit state. The CPU time required for each simulation was found to be about 15 min on a Centrino 2.0 GHz PC.

For the serviceability limit state, soils with small values of Young's modulus are used in this paper. In such soils, the variability of the compressibility characteristics is very large (Bauer and Pula 2000). A lognormal distribution is used for $E$ with a mean value of 60 MPa (Nour et al. 2002). For the coefficient of variation, some values proposed and used by several authors are listed in Table 1. A value of 15% is used in this paper. Regarding the Poisson's ratio, there is no available information about its

<table><caption>**Table 1.** Values of Coefficient of Variation (COV) of Young's Modulus Proposed by Several Authors</caption><thead><tr><th>Authors</th><th>COV of Young's modulus (%)</th></tr></thead><tbody><tr><td>Phoon and Kulhawy (1999)</td><td>30</td></tr><tr><td>Bauer and Pula (2000)</td><td>15</td></tr><tr><td>Nour et al. (2002)</td><td>40-50</td></tr><tr><td>Baecher and Christian (2003)</td><td>2-42</td></tr></tbody></table>

random variation. Some authors have suggested that the randomness can be neglected in an analysis of settlement taking place in the case of elastic soil. Others have stated that $\nu$ changes with a relatively narrow interval. In this paper, $\nu$ is considered as a lognormally distributed variable with a coefficient of variation of 5%. Its mean value is taken equal to 0.3. For the correlation coefficient of these two parameters, there is no information available. The results reported by some researchers (Bauer and Pula 2000) lead to the conclusion that this correlation is negative. In this paper, the cases of uncorrelated and correlated soil elastic properties with $\rho_{E,\nu} = -0.5$ are considered. The CPU time required for each serviceability limit state simulation was found to be about 10 min on a Centrino 2.0 GHz PC.

## Ultimate Limit State

### Graphical Representation of Successive Tentative Response Surfaces

Fig. 7 shows the evolution of the tentative response surface in the standard space ($u_1, u_2$) for a footing applied load equal to 775 kN/m (i.e., a safety factor of 3.1). The two equations used for the transformation of each ($c, \phi$) of the limit state surface from the physical space to the standardized normal uncorrelated space ($u_1, u_2$) are (e.g., Lemaire 2005)

<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Evolution of tentative response surface</figcaption></figure>
<table><caption>**Table 2. Reliability Index and Design Point for Uncorrelated and Correlated Shear Strength Parameters**</caption><thead><tr><th rowspan="2">P<sub>S</sub><br>(kN/m)</th><th rowspan="2">F</th><th colspan="5">ρ<sub>c,φ</sub>=0.0</th><th colspan="5">ρ<sub>c,φ</sub>=-0.5</th></tr><tr><th>c<sup>*</sup><br>(kPa)</th><th>φ<sup>*</sup><br>(°)</th><th>β<sub>HL</sub></th><th>F<sub>c</sub></th><th>F<sub>φ</sub></th><th>c<sup>*</sup><br>(kPa)</th><th>φ<sup>*</sup><br>(°)</th><th>β<sub>HL</sub></th><th>F<sub>c</sub></th><th>F<sub>φ</sub></th></tr></thead><tbody><tr><td>750</td><td>3.19</td><td>14.12</td><td>20.86</td><td>3.49</td><td>1.42</td><td>1.52</td><td>17.08</td><td>19.34</td><td>4.62</td><td>1.17</td><td>1.64</td></tr><tr><td>1,150</td><td>2.08</td><td>16.30</td><td>24.27</td><td>2.12</td><td>1.23</td><td>1.28</td><td>18.26</td><td>23.08</td><td>2.71</td><td>1.10</td><td>1.35</td></tr><tr><td>1,550</td><td>1.54</td><td>17.90</td><td>26.61</td><td>1.21</td><td>1.12</td><td>1.15</td><td>18.64</td><td>26.43</td><td>1.53</td><td>1.07</td><td>1.16</td></tr><tr><td>1,780</td><td>1.35</td><td>18.52</td><td>27.71</td><td>0.81</td><td>1.08</td><td>1.10</td><td>20.04</td><td>27.22</td><td>1.00</td><td>0.99</td><td>1.12</td></tr><tr><td>1,950</td><td>1.23</td><td>18.89</td><td>28.45</td><td>0.55</td><td>1.06</td><td>1.07</td><td>19.92</td><td>28.14</td><td>0.67</td><td>1.00</td><td>1.08</td></tr><tr><td>2,395</td><td>1.00</td><td>19.61</td><td>30.00</td><td>0.00</td><td>1.01</td><td>1.00</td><td>19.61</td><td>30.00</td><td>0.00</td><td>1.01</td><td>1.00</td></tr></tbody></table>

$$u_1 = \left( \frac{c - \mu_c^N}{\sigma_c^N} \right) \qquad (7)$$

$$u_2 = \frac{1}{\sqrt{1 - \rho^2}} \left[ \left( \frac{\varphi - \mu_{\varphi}^N}{\sigma_{\varphi}^N} \right) - \rho \left( \frac{c - \mu_c^N}{\sigma_c^N} \right) \right] \quad (8)$$

where $ρ$=coefficient of correlation of $c$ and $φ$; and $μ_c^N$, $μ_{φ}^N$, $σ_c^N$, and $σ_{φ}^N$=equivalent normal means and standard deviations of the random variables $c$ and $φ$. They are determined from the translation approach using the following equations

$$\frac{c - \mu_c^N}{\sigma_c^N} = \Phi^{-1}[F_c(c)] \qquad (9)$$

$$\frac{\varphi - \mu_{\varphi}^N}{\sigma_{\varphi}^N} = \Phi^{-1}[F_{\varphi}(\varphi)] \qquad (10)$$

where $F_c$ and $F_φ$=non-Gaussian cumulative distribution functions of $c$ and $φ$; and $Φ^{-1}(·)=inverse$ of the standard normal cumulative distribution.

A convergence criterion on the reliability index was adopted. It considers that convergence is reached when a difference smaller than $10^{-2}$ between two successive reliability indexes is achieved. One can notice that this criterion is reached after only four iterations. Thus, only 20 numerical simulations by $FLAC^{3D}$ were necessary. The corresponding CPU time required is about $20 \times 15 = 300$ min (i.e., 5 h). A value of 4.35 was found for the reliability index. This value corresponds to a failure probability of $6.84 \times 10^{-6}$ calculated by the $FORM$ approximation.

## Reliability Index and Design Point

Table 2 presents the Hasofer-Lind reliability index and the corresponding design point for different values of the vertical applied load $P_S$ (i.e., safety factor $F=P_u/P_S$) varying from small values up to the deterministic ultimate load. The cases of correlated ($ρ_{c,φ}=-0.5$) and uncorrelated ($ρ_{c,φ}=0$) shear strength parameters are considered.

The reliability index decreases with the increase of the applied load $P_S$ (i.e., the decrease of the safety factor $F=P_u/P_S$) until it vanishes for an applied load equal to the deterministic ultimate load. This case corresponds to a deterministic state of failure for which $F=1$ using the mean values of the random variables and the failure probability is equal to 50%. The comparison of the results of correlated variables with those of uncorrelated variables shows that the reliability index corresponding to uncorrelated variables is smaller than the one of negatively correlated variables. One can conclude that the hypothesis of uncorrelated shear strength parameters is conservative in comparison to the one of negatively correlated parameters. For instance, when the safety

factor is equal to 3.2 (i.e., $P_S=750$ kN/m), the reliability index increases by 32% if the variables $c$ and $φ$ are considered as negatively correlated.

The values of the design points corresponding to different values of the vertical applied load can give an idea about the partial safety factors of each of the strength parameters $c$ and $tan φ$ as follows

$$F_c = \frac{\mu_c}{c^*} \qquad (11)$$

$$F_{\varphi} = \frac{\tan(\mu_{\varphi})}{\tan \varphi^*} \qquad (12)$$

Table 2 shows that for uncorrelated shear strength parameters, the values of $c^*$ and $φ^*$ at the design point are smaller than their respective mean values and increase with the increase of the applied load. Consequently, the partial safety factors $F_c$ and $F_φ$ decrease with the increase of the applied load. They tend to 1 when $P_S=P_u$. For negatively correlated shear strength parameters, $c^*$ slightly exceeds the mean for some values of the applied load. This can be explained by the counterclockwise rotation of the dispersion ellipse due to the negative correlation (Fig. 8). The position of the design point, which is the point of tangency between the ellipse and the limit state surface, changes from that found for uncorrelated soil shear strength parameters. A higher $c^*$ (respectively, a lower $φ^*$) is found. Consequently, $c^*$ can become greater than the mean value for a negative correlation. This con-

<figure><img src="image_10.png" /><figcaption>**Fig. 8.** General layout of dispersion ellipse for different correlation coefficients</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Effect of variability of soil shear strength parameters on failure probability</figcaption></figure>

clusion is similar to that found by Youssef Abdel Massih et al. (2008).

## Sensitivity of Failure Probability to Variability of Soil Shear Strength Parameters

In order to study the effect of the variability of the soil shear strength parameters on the failure probability, Fig. 9 shows the FORM failure probability versus the coefficient of variation of $c$ and $\phi$. For each curve, the coefficient of variation of a parameter is held to the same constant value given in the introduction of the section “Numerical Results” and the coefficient of variation of the second parameter is varied over the range 10–40%. The results show that the failure probability is highly influenced by the coefficient of variation of the angle of internal friction; the greater the scatter in $\phi$ the higher the failure probability of the foundation. This means that the accurate determination of the distribution of this parameter is very important in obtaining reliable probabilistic results. In contrast, the coefficient of variation of $c$ does not significantly affect the failure probability.

## Serviceability Limit State

### Reliability Index and Design Point

The threshold value of the settlement is $u_{max}=0.1$ m. Table 3 presents the Hasofer–Lind reliability index and the corresponding design point for different values of the vertical applied load $P_S$. The cases of correlated and uncorrelated soil elastic properties are considered. The reliability index decreases with the increase of the applied load $P_S$. A comparison of the results of correlated soil elastic properties with those of uncorrelated ones shows that, as in

the ultimate limit state, the hypothesis of uncorrelated soil elastic properties is conservative in comparison to the one of negatively correlated properties.

By comparing Tables 2 and 3, one can notice that for small values of the applied load, the reliability index of the ultimate limit state is significantly smaller than that of the serviceability limit state. Thus, for small values of the applied load, the ultimate failure mode is predominant and will have the highest contribution in the determination of the system failure probability. The difference between the reliability indexes of the two failure modes becomes smaller for higher values of the applied load. Consequently, when the applied load increases, the two failure modes (i.e., the ultimate and the serviceability ones) will have nearly similar contributions in the computation of the system failure probability (see another interpretation in the section “System Failure Probability”).

As for the ultimate limit state, the values of the design point allow us to calculate the partial safety factors of $E$ and $\nu$ as follows

$$F_E = \frac{\mu_E}{E^*} \qquad (13)$$

$$F_v = \frac{\mu_v}{\nu^*} \qquad (14)$$

Table 3 shows that for uncorrelated soil elastic properties, the partial safety factors $F_E$ and $F_v$ decrease with the increase of $P_S$. They become equal to 1 when $P_S$ is equal to the load that leads to the maximal prescribed foundation settlement $u_{max}$ for the mean values of the soil elastic properties. For negatively correlated soil elastic properties, $F_v$ is found smaller than 1. The same interpretation given in the ultimate limit state for negatively correlated variables remains valid in the present case.

## Sensitivity of Failure Probability to Variability of Soil Elastic Properties

As for the ultimate limit state, the effect of the variability of the soil elastic properties on the failure probability is shown in Fig. 10 where FORM failure probability is plotted versus the coefficient of variation of $E$ and $\nu$. For each curve, the coefficient of variation of a parameter is held to the same constant value given in the introduction of the section “Numerical Results” and the coefficient of variation of the second parameter is varied over the range 5–35%. The results show that the failure probability of the serviceability limit state is highly influenced by the coefficient of variation of the Young’s modulus; the greater the scatter in $E$ the higher the failure probability of the foundation. This means that an accurate determination of the distribution of this parameter is very important in obtaining reliable probabilistic results. In con-

<table><caption>**Table 3.** Reliability Index and Design Point for Uncorrelated and Correlated Soil Elastic Properties</caption><thead><tr><th rowspan="2">P<sub>S</sub><br>(kN/m)</th><th colspan="5">ρ<sub>E,v</sub>=0.0</th><th colspan="5">ρ<sub>E,v</sub>=-0.5</th></tr><tr><th>E<sup>*</sup><br>(MPa)</th><th>ν<sup>*</sup></th><th>β<sub>HL</sub></th><th>F<sub>E</sub></th><th>F<sub>v</sub></th><th>E<sup>*</sup><br>(MPa)</th><th>ν<sup>*</sup></th><th>β<sub>HL</sub></th><th>F<sub>E</sub></th><th>F<sub>v</sub></th></tr></thead><tbody><tr><td>750</td><td>19.42</td><td>0.280</td><td>7.60</td><td>3.10</td><td>1.07</td><td>16.57</td><td>0.337</td><td>8.80</td><td>3.62</td><td>0.89</td></tr><tr><td>1,150</td><td>33.71</td><td>0.288</td><td>3.87</td><td>1.78</td><td>1.04</td><td>31.06</td><td>0.319</td><td>4.46</td><td>1.93</td><td>0.94</td></tr><tr><td>1,550</td><td>49.75</td><td>0.296</td><td>1.21</td><td>1.21</td><td>1.01</td><td>48.37</td><td>0.306</td><td>1.41</td><td>1.24</td><td>0.98</td></tr><tr><td>1,780</td><td>59.33</td><td>0.300</td><td>0.00</td><td>1.01</td><td>1.00</td><td>59.33</td><td>0.300</td><td>0.00</td><td>1.01</td><td>1.00</td></tr></tbody></table>
<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Effect of variability of soil elastic properties on failure probability</figcaption></figure>

trast, compared to *E*, the uncertainties in *v* have a minor effect on the failure probability.

## System Failure Probability

The system failure probability under the two failure modes involving the ultimate and the serviceability limit states of the footing is given by

$$P_{f\text{sys}} = P_f(U \cup S) = P_f(U) + P_f(S) - P_f(U \cap S) \quad (15)$$

where $P_f(U \cap S)$=failure probability under the ultimate and the serviceability failure modes; $P_f(U)$=failure probability under only the ultimate failure mode; and $P_f(S)$=failure probability under only the serviceability failure mode. The failure probability of the intersection is given as follows (Lemaire 2005)

$$\max(P(A), P(B)) \le P_f(P \cap S) \le P(A) + P(B) \quad (16)$$

where

$$P(A) = \Phi(-\beta_U)\Phi\left(-\frac{\beta_S - \rho_{US}\beta_U}{\sqrt{1 - \rho_{US}^2}}\right) \quad (17)$$

$$P(B) = \Phi(-\beta_S)\Phi\left(-\frac{\beta_U - \rho_{US}\beta_S}{\sqrt{1 - \rho_{US}^2}}\right) \quad (18)$$

$$\rho_{US} = \langle \alpha_U \rangle \{ \alpha_S \} \quad (19)$$

$\beta_U$ and $\beta_S$=reliability indexes corresponding to the ultimate and the serviceability failure modes, respectively, and $\rho_{US}$=cor-

relation between the two failure modes. $\alpha_U$ and $\alpha_S$ for both modes are given by

$$\alpha_{U_i} = - \left. \frac{\partial \beta_U}{\partial u_{U_i}} \right|_{\{u_{U_i}^*\}} = - \frac{u_{U_i}^*}{\beta_U} \quad (20)$$

$$\alpha_{S_i} = - \left. \frac{\partial \beta_S}{\partial u_{S_i}} \right|_{\{u_{S_i}^*\}} = - \frac{u_{S_i}^*}{\beta_S} \quad (21)$$

where $u_{U_i}^*$ and $u_{S_i}^*$=standard uncorrelated normal variables at the design points [see Eqs. (7) and (8)]. Here, it was found that $\rho_{US}=0$ since the two failure modes are independent.

The probability of the intersection $P_f(U \cap S)$ is set equal to its lower limit [i.e., $\max(P(A), P(B))$] in order to obtain the higher limit of the system failure probability $P_{f\text{sys}}$.

The system reliability index can be approximated using the FORM approximation as follows

$$\beta_{\text{sys}} = - \Phi^{-1}(P_{f\text{sys}}) \quad (22)$$

Table 4 presents the system failure probability $P_{f\text{sys}}$ and the corresponding reliability index $\beta_{\text{sys}}$ for different values of the applied load. Four cases are considered: They are the combinations of correlated and uncorrelated shear strength parameters with correlated and uncorrelated soil elastic properties. From this table, it can be seen that even for the system reliability, the assumption of uncorrelated parameters is conservative in comparison to the one of negatively correlated variables. For small values of the applied load, where the ultimate failure mode is predominant, one can notice that the system reliability index is equal to the reliability index of the ultimate failure mode. When the applied load increases, the system reliability depends on both failure modes and the system reliability index is smaller than the ones corresponding to a single failure mode. Finally, one can notice that a negative system reliability index is found for an applied load of 1,780 kN/m corresponding to the deterministic failure state of the serviceability mode. This corresponds to a system failure probability higher than 50% (i.e., higher than the one corresponding to a zero reliability index). This value is due to the combination of the failure probability of the serviceability limit state (i.e., 50%) and the one of the ultimate limit state. In this case, the system reliability index is meaningless.

## Conclusions

A reliability-based analysis of a strip footing resting on a $c-\varphi$ soil and subjected to a central vertical load is presented. Both the ultimate and the serviceability limit states are used to characterize the footing behavior. Two deterministic models based on numeri-

<table><caption>**Table 4.** System Failure Probability and Reliability Index</caption><thead><tr><th rowspan="2">P<sub>S</sub><br>(kN/m)</th><th colspan="2">ρ<sub>c,φ</sub>=0.0<br>ρ<sub>E,v</sub>=0.0</th><th colspan="2">ρ<sub>c,φ</sub>=-0.5<br>ρ<sub>E,v</sub>=-0.5</th><th colspan="2">ρ<sub>c,φ</sub>=-0.5<br>ρ<sub>E,v</sub>=0.0</th><th colspan="2">ρ<sub>c,φ</sub>=0.0<br>ρ<sub>E,v</sub>=-0.5</th></tr><tr><th>P<sub>fsys</sub><br>(%)</th><th>β<sub>sys</sub></th><th>P<sub>fsys</sub><br>(%)</th><th>β<sub>sys</sub></th><th>P<sub>fsys</sub><br>(%)</th><th>β<sub>sys</sub></th><th>P<sub>fsys</sub><br>(%)</th><th>β<sub>sys</sub></th></tr></thead><tbody><tr><td>750</td><td>0.02</td><td>3.49</td><td>2.00 × 10<sup>-4</sup></td><td>4.62</td><td>2.00 × 10<sup>-4</sup></td><td>4.62</td><td>0.02</td><td>3.49</td></tr><tr><td>1,150</td><td>1.70</td><td>2.12</td><td>0.34</td><td>2.71</td><td>0.34</td><td>2.70</td><td>1.70</td><td>2.12</td></tr><tr><td>1,550</td><td>21.30</td><td>0.79</td><td>13.70</td><td>1.09</td><td>16.90</td><td>0.96</td><td>18.30</td><td>0.90</td></tr><tr><td>1,780</td><td>60.45</td><td>-0.26</td><td>57.93</td><td>-0.20</td><td>57.93</td><td>-0.20</td><td>60.45</td><td>-0.26</td></tr></tbody></table>
cal simulations using the Lagrangian explicit finite difference
code FLAC³D are employed. The first one computes the ultimate
bearing capacity of the foundation and the second one calculates
the footing displacement due to an applied service load. The
Hasofer-Lind reliability index is adopted here for the assessment
of the foundation reliability. The response surface methodology is
used to find an approximation of the analytically unknown limit
state surfaces and the corresponding reliability indexes. Only the
soil shear strength parameters are considered as random variables
while studying the ultimate limit state. Also, the randomness of
only the soil elastic properties is taken into account in the service-
ability limit state. The main conclusions of this paper can be
summarized as follows:

1. The hypothesis of uncorrelated parameters was found to be conservative in comparison to the one of negatively correlated variables;

2. For uncorrelated shear strength parameters, the values of $c^*$ and $\phi^*$ at the design point are found smaller than their respective mean values and increase with the increase of the applied load $P_S$. Consequently, the partial safety factors $F_c$ and $F_\phi$ decrease with the increase of the applied load. They tend to 1 when $P_S=P_u$. For negatively correlated shear strength parameters, $c^*$ slightly exceeds the mean for some values of the applied load;

3. For uncorrelated soil elastic properties, the partial safety factors $F_E$ and $F_v$ decrease with the increase of $P_S$. They tend to 1 when $P_S$ is equal to the load that leads to the maximal prescribed foundation settlement $u_{max}$ for the mean values of the soil elastic properties. For negatively correlated soil elastic properties, $F_v$ is found smaller than 1;

4. The failure probability is found to be highly influenced by the uncertainties of the angle of internal friction for the ultimate limit state and by those of the Young's modulus for the serviceability limit state; and

5. For small values of the applied load, the ultimate limit state is predominant in the computation of the system failure probability. Consequently, the system reliability index is found equal to that of the ultimate limit state. For higher values of the applied load, the system reliability index depends on both limit states. It is smaller than the ones corresponding to a single failure mode. Thus, both failure modes have to be considered in the reliability analysis of foundations for high values of the applied load.

## Acknowledgments

The writers would like to thank the Lebanese National Council for Scientific Research (CNRSL) and the French organization EGIDE for providing the financial support for this research.

## References

Baecher, G., and Christian, J. (2003). *Reliability and statistics in geotechnical engineering*, Wiley, U.K.

Bauer, J., and Pula, W. (2000). "Reliability with respect to settlement limit-states of shallow foundations on linearly-deformable subsoil." *Comput. Geotech.*, 26, 281–308.

Cherubini, C. (2000). "Reliability evaluation of shallow foundation bearing capacity on c', φ' soils." *Can. Geotech. J.*, 37, 264–269.

Ditlevsen, O. (1981). *Uncertainty modelling: With applications to multi-dimensional civil engineering systems*, McGraw-Hill, New York.

Fenton, G. A., and Griffiths, D. V. (2002). "Probabilistic foundation settlement on spatially random soil." *J. Geotech. Geoenviron. Eng.*, 128(5), 381-390.

Fenton, G. A., and Griffiths, D. V. (2003). "Bearing capacity prediction of spatially random C-φ soils." *Can. Geotech. J.*, 40, 54-65.

Fenton, G. A., and Griffiths, D. V. (2005). "Three-dimensional probabilistic foundation settlement." *J. Geotech. Geoenviron. Eng.*, 131(2), 232-239.

FLAC³D. (1993). *Fast Lagrangian analysis of continua*, ITASCA Consulting Group, Inc., Minneapolis.

Griffiths, D. V., and Fenton, G. A. (2001). "Bearing capacity of spatially random soil: the undrained clay Prandtl problem revisited." *Geotechnique*, 51(4), 351-359.

Griffiths, D. V., Fenton, G. A., and Manoharan, N. (2002). "Bearing capacity of rough rigid strip footing on cohesive soil: Probabilistic study." *J. Geotech. Geoenviron. Eng.*, 128(9), 743-755.

Harr, M. E. (1987). *Reliability-based design in civil engineering*, McGraw-Hill, New York.

Hasofer, A. M., and Lind, N. C. (1974). "Exact and invariant second-moment code format." *J. Engrg. Mech. Div.*, 100(1), 111-121.

Lemaire, M. (2005). *Fiabilité des structures*, Hermès, Lavoisier, Paris (in French).

Low, B. K. (2005). "Reliability-based design applied to retaining walls." *Geotechnique*, 55(1), 63-75.

Low, B. K., and Phoon, K. K. (2002). "Practical first-order reliability computations using spreadsheet." *Proc., Probabilistics in Geotechnics: Technical and Economic Risk Estimation*, Verlag Gluckauf GmbH. Essen, Graz, Austria, 39-46.

Low, B. K., and Tang, W. H. (1997a). "Efficient reliability evaluation using spreadsheet." *J. Eng. Mech.*, 123(7), 749-752.

Low, B. K., and Tang, W. H. (1997b). "Reliability analysis of reinforced embankments on soft ground." *Can. Geotech. J.*, 34, 672-685.

Low, B. K., and Tang, W. H. (2004). "Reliability analysis using object-oriented constrained optimization." *Struct. Safety*, 26, 68-89.

Lumb, P. (1970). "Safety factors and the probability distribution of soil strength." *Can. Geotech. J.*, 7, 225-242.

Nour, A., Slimani, A., and Laouami, N. (2002). "Foundation settlement statistics via finite element analysis." *Comput. Geotech.*, 29, 641-672.

Phoon, K.-K., and Kulhawy, F. H. (1999). "Evaluation of geotechnical property variability." *Can. Geotech. J.*, 36, 625-639.

Popescu, R., Deodatis, G., and Nobahar, A. (2005). "Effect of random heterogeneity of soil properties on bearing capacity." *Probab. Eng. Mech.*, 20, 324-341.

Przewlocki, J. (2005). "A stochastic approach to the problem of bearing capacity by the method of characteristics." *Comput. Geotech.*, 32, 370-376.

Tandjiria, V., Teh, C. I., and Low, B. K. (2000). "Reliability analysis of laterally loaded piles using response surface methods." *Struct. Safety*, 22, 335-355.

Wolff, T. H. (1985). "Analysis and design of embankment dam slopes: A probabilistic approach." Ph.D. thesis, Purdue Univ., Lafayette, Ind.

Yin, J.-H., Wang, Y.-J., and Selvadurai, P. S. (2001). "Influence of non-associativity on the bearing capacity of a strip footing." *J. Geotech. Geoenviron. Eng.*, 127(11), 985-989.

Youssef Abdel Massih, D. S., Soubra, A.-H., and Low, B. K. (2008). "Reliability-based analysis and design of strip footings against bearing capacity failure." *J. Geotech. Geoenviron. Eng.*, in press.

Yuceman, M. S., Tang, W. H., and Ang, A. H. S. (1973). "A probabilistic study of safety and design of earth slopes." *Civil engineering studies, Structural Research Series* 402, University of Illinois Press, Urbana, Ill.