Reliability-Based Analysis of Strip Footings Using Response Surface Methodology

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 reliabilitybased 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 − 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 3D 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 for the ultimate limit state, and the soil elastic properties E and 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 3D 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 ␤ = min͓R͑͒ / r͔͑͒ ͑Fig.1͒.Its matrix formulation is ͑Ditlevsen 1981͒ in which xϭvector representing the n random variables; ϭvector of their mean values; and Cϭtheir covariance matrix.The minimization of Eq. ͑1͒ is performed subject to the constraint G͑x͒ Յ 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͒ Յ 0 and a safe region given by G͑x͒ Ͼ 0.
The classical approach for computing ␤ 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 ␤ HL .
An intuitive interpretation of the reliability index was suggested in Low andTang ͑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 ␤ HL = 1 in Eq. ͑1͒ without the min͔ centered at the equivalent normal mean values ͑ 1 N , 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 ␤ HL may be regarded as the codirectional axis ratio of the smallest ellipse ͑which is either an expansion or a contraction of the 1 − ellipse͒ that just touches the limit state surface to the 1− 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 ␤ HL , one can approximate the failure probability as follows where ⌽͑•͒ϭ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 N and the equivalent normal standard deviation 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 N and equivalent normal standard deviation 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͒ 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 3D
FLAC 3D ͑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 3D simulation considers a homogeneous soil.The randomness of the soil is taken into ac-

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 − 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 = 0.3͒.The values of the soil shear strength parameters used in the analysis are: = 30°, = 20°, and c = 20 kPa, where ϭ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 = 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 3D , 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ϫ 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 ϫ 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ϫ 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 ϫ 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 ϫ 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 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 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 elasticperfectly 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  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 − 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 , the Young's modulus E, and the Poisson's ratio of the soil are considered as random variables.Due to the relatively low effect of the elastic modulus E and the Poisson's ratio on the ultimate bearing capacity, only c and will be considered as random variables while studying the ultimate limit state.Similarly, only the randomness of E and 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 ͑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 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 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 i and their coefficients of variation i .A brief explanation of the algorithm used is as follows: 1. Evaluate the performance function G͑x͒ at the mean value point and the 2n points each at ± k where k = 1in 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 ␤ 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 3D 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 3D .Therefore, an exchange of data between FLAC 3D 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 3D 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 .The results of Wolff ͑1985͒ ͑ c, = −0.47͒,Yuceman et al. ͑1973͒ ͑−0.49Յ c, Յ −0.24͒, Lumb ͑1970͒ ͑−0.7Յ c, Յ −0.37͒, and Cherubini ͑2000͒ ͑ c, = −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 c, are given as follows: c = 20 kPa, = 30°, COV c = 20%, COV = 10%, and c, = −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 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 .It should be mentioned that the soil elastic properties ͑i.e., K and G or E and ͒ 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 = 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 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 changes with a relatively narrow interval.In this paper, 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 E, = −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.

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 , ͒ of the limit state surface from the physical space to the standardized normal uncorrelated space ͑u 1 , u 2 ͒ are ͑e.g., Lemaire 2005͒ 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 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ϫ 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 ϫ 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 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-

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 .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 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.

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 as follows Table 3 shows that for uncorrelated soil elastic properties, the partial safety factors F E and F 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 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 .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-  trast, compared to E, the uncertainties in 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 where P f ͑U പ 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͒͒ Յ P f ͑P പ S͒ Յ P͑A͒ + P͑B͒ ͑16͒ where ␤ U and ␤ S ϭreliability indexes corresponding to the ultimate and the serviceability failure modes, respectively, and US ϭcor-relation between the two failure modes.␣ U and ␣ S for both modes are given by 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 US = 0 since the two failure modes are independent.The probability of the intersection P f ͑U പ 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 sys .
The system reliability index can be approximated using the FORM approximation as follows Table 4 presents the system failure probability P fsys and the corresponding reliability index ␤ 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 sate.In this case, the system reliability index is meaningless.

Conclusions
A reliability-based analysis of a strip footing resting on a c − 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-  cal simulations using the Lagrangian explicit finite difference code FLAC 3D 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 serviceability 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 * 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 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 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 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.

Fig. 1 .
Fig. 1.Design point and equivalent normal dispersion ellipses in space of two random variables

Fig. 5 .
Fig. 5. ͑a͒ Normal ͑͒; ͑b͒ shear ͑͒ stresses at right half of footing base as obtained from two methods

Fig. 6 .
Fig. 6.Vertical footing displacement versus number of cycles due to applied load P S = 750 kN/ m

Fig. 9 .
Fig. 9. Effect of variability of soil shear strength parameters on failure probability

Fig. 10 .
Fig. 10.Effect of variability of soil elastic properties on failure probability

Table 1 .
Values of Coefficient of Variation ͑COV͒ of Young's Modulus Proposed by Several Authors Fig. 7. Evolution of tentative response surface

Table 2 .
Reliability Index and Design Point for Uncorrelated and Correlated Shear Strength Parameters Fig. 8. General layout of dispersion ellipse for different correlation coefficients clusion is similar to that found by Youssef Abdel Massih et al. ͑2008͒.

Table 3 .
Reliability Index and Design Point for Uncorrelated and Correlated Soil Elastic Properties