Bearing capacity of spatially random rock masses obeying Hoek–Brown failure criterion

ABSTRACT This paper presents a probabilistic analysis to compute the probability density function of the bearing capacity of a strip footing resting on a spatially varying rock mass. The rock is assumed to follow the generalised Hoek–Brown failure criterion. The uniaxial compressive strength of the intact rock (σc) was considered as a random field and the geological strength index was modelled as a random variable. The uncertainty propagation methodology employed in the analysis is the sparse polynomial chaos expansion. A global sensitivity analysis based on Sobol indices was performed. Some numerical results were presented and discussed.


Introduction
The analysis and design of a strip footing resting on a rock mass obeying Hoek-Brown (HB) failure criterion are generally based on deterministic approaches (Maghous, de Buhan, and Bekaert 1998;Yang and Yin 2005;Merifield, Lyamin, and Sloan 2006;Saada, Maghous, and Garnier 2008). In this paper, the behaviour of a strip footing resting on a HB rock mass is studied using a probabilistic approach. The probabilistic approaches allow one to consider the propagation of the uncertainties from the input parameters to the system responses. Most existing probabilistic analyses on foundations consider the case of a foundation resting on a soil mass (Griffiths and Fenton 2001;Griffiths, Fenton, and Manoharan 2002;Fenton and Griffiths 2003;Popescu, Deodatis, and Nobahar 2005;Youssef Abdel Massih, Soubra, and Low 2008;Cho and Park 2010;Soubra and Youssef Abdel Massih 2010;Li, Tian, and Cassidy 2015). Only few studies may be found in literature concerning the probabilistic analysis of a foundation resting on a rock mass obeying HB failure criterion (Ching et al. 2011;Mao, Al-Bittar, and Soubra 2012). This paper fills this gap. It aims at determining the ultimate bearing capacity of a strip footing resting on a spatially varying rock mass and subjected to a vertical load. The rock mass follows the generalised HB failure criterion Brown 1980, 1997;Hoek, Carranze-Torres, and Corkum 2002;Hoek and Marinos 2007;Brown 2008). In this criterion, only intact rocks or heavily jointed rocks masses (i.e. with sufficiently dense and randomly distributed joints) can be considered. The HB failure criterion is characterised by four parameters: (i) the geological strength index (GSI), (ii) the uniaxial compressive strength of the intact rock (σ c ), (iii) the intact rock material constant (m i ) and (iv) the disturbance coefficient (D). Mao, Al-Bittar, and Soubra (2012) have modelled these four parameters as random variables and have performed a probabilistic analysis of the ultimate bearing capacity of foundations. These authors have shown that the variability of the ultimate bearing capacity is mainly due to the uniaxial compressive strength of the intact rock (σ c ) and the geological strength index (GSI). Based on this study, only these two parameters were considered herein as uncertain. The rock uniaxial compressive strength of the intact rock (σ c ) was considered as a non-Gaussian (NG) random field characterised by a square exponential autocorrelation function. The expansion optimal linear estimation (EOLE) method proposed by Li and Der Kiureghian (1993) was used to dicretise this random field. As for GSI, Ching et al. (2011) have stated that this parameter is based on engineering judgement. It characterises the overall rock mass condition and it does not represent a precise physical parameter varying in space. Thus, this parameter cannot be modelled as a random field and will be treated herein as a random variable.
As for the probabilistic method of analysis, the classical Monte Carlo simulation (MCS) methodology is generally used when dealing with random fields (Griffiths and Fenton 2001;Griffiths, Fenton, and Manoharan 2002;Fenton and Griffiths 2003;Popescu, Deodatis, and Nobahar 2005;Cho and Park 2010;Ching et al. 2011;Li, Tian, and Cassidy 2015). In spite of being a rigorous method, MCS requires a large number of calls to the deterministic model. This is not convenient in the case where a computationally expensive (finite element or finite difference) deterministic model is used. This paper makes use of an efficient probabilistic approach which significantly reduces the number of calls of the deterministic model. The sparse polynomial chaos expansion (SPCE) methodology was proposed in this regard. Notice that the SPCE is an extension of the polynomial chaos expansion (PCE). A PCE or a SPCE methodology aims at replacing the complex deterministic model by a meta-model (i.e. a simple analytical equation) (Al-Bittar and Soubra 2013Soubra , 2014Jiang et al. 2014Jiang et al. , 2015. The probability density function (PDF) of the system response can thus be easily obtained. This is because MCS is no longer applied on the original computationally expensive deterministic model, but on the meta-model. The other significant advantage of the present SPCE methodology with respect to the classical MCS method is that it allows one to easily perform a global sensitivity analysis (GSA) based on Sobol indices using the metamodel. These indices give the contribution of each uncertain parameter to the variability of the system response.
The paper is organised as follows: the next section aims at presenting the deterministic model used for the computation of the ultimate bearing capacity of a vertically loaded strip footing resting on a rock mass obeying the HB failure criterion. It is followed by a presentation of the probabilistic method used for the computation of the PDF of the ultimate bearing capacity of a strip footing resting on a HB spatially varying rock mass. Finally, the probabilistic numerical results are presented and discussed. The paper ends with a conclusion.

Deterministic model
In this section, one first presents a brief description of the generalised HB failure criterion. This is followed by a presentation of the numerical model used to compute the ultimate bearing capacity of a strip footing resting on a HB rock mass and subjected to a vertical load.

Generalised HB failure criterion
The generalised HB failure criterion only deals with intact rocks or heavily jointed rock masses Brown 1980, 1997;Hoek, Carranze-Torres, and Corkum 2002;Hoek and Marinos 2007;Brown 2008). A heavily jointed rock mass involves sufficiently dense and randomly distributed joints so that in the scale of the problem, it can be regarded as an isotropic assembly of interlocking particles. Consequently, rocks with few discontinuities cannot be considered in this framework. The generalised HB failure criterion can be described by the following equation: where σ 1 and σ 3 are, respectively, the major and minor principal stresses at failure and σ c is the uniaxial compressive strength of the intact rock material. The parameters m b , s and n are given by the following equations: In these equations, the geological strength index (GSI) characterises the quality of the rock mass and depends on its structure and its joints surface conditions (Hoek and Brown 1997). On the other hand, the parameter m i is the value of parameter m for intact rock and can be obtained from experimental tests. The parameter m i varies from 4 for very fine weak rock-like claystone to 33 for coarse igneous light-coloured rock such as granite. Finally, D is the disturbance coefficient. It varies from 0.0 for undisturbed in situ rock masses to 1.0 for very disturbed rock masses. Figure 1 presents the nonlinear generalised HB failure criterion in the (τ, σ) plan. As for the modulus of deformation of the HB rock mass, Hoek, Carranze-Torres, and Corkum (2002) have proposed the following relationship between this parameter and the generalised HB failure criterion parameters: where E m in this equation is given in GPa.

Numerical model
The deterministic model used to calculate the ultimate bearing capacity of a strip footing resting on a HB rock mass and subjected to a vertical load was based on the commercial numerical code FLAC 3D . A footing of breadth B = 1 m was considered in the analysis. For this calculation, a rock mass of 20 m wide by 6 m deep was considered ( Figure 2). In contrast to the homogeneous rock case where only one-half of the rock mass domain may be considered in the analysis, the entire rock domain shown in Figure 2 was considered herein. This is because the random field creates a non-symmetrical failure mechanism. A control displacement approach was used in this paper. The rigid strip footing was modelled by a prescribed uniform velocity for all the rock nodes in contact with the footing. A value of 10 −6 m/time step was chosen for this velocity since a smaller value was proved to negligibly decrease the value of the ultimate bearing capacity. The mesh was refined near the footing edges where high stress gradient may occur. 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. The rock behaviour was modelled by an elastic perfectly plastic model obeying the generalised HB failure criterion. It should be emphasised here that an associated flow rule was considered in the study in order to be able to compare the obtained results with those given in literature (Merifield, Lyamin, and Sloan 2006;Mao, Al-Bittar, and Soubra 2012). For this purpose, the confining stress at constant volume s cv 3 must be properly selected. In fact, beyond the value of s cv 3 , no volume changes are expected to appear. This means that when s cv 3 /s c is very small, the case of a deformation at constant volume is rapidly reached and the model can be considered to follow a non-associated flow rule with a zero dilation angle. On the contrary, the case of a large value of s cv 3 /s c means that the deformation at the constant volume cannot be reached easily and thus the model can be considered to follow an associated flow rule. In the present paper, a large value of s cv 3 /s c = 2 was selected. This value was chosen since greater values have led to the same value of the ultimate bearing capacity. The present deterministic model was validated by comparison of its results with those provided by Merifield, Lyamin, and Sloan (2006) and Mao, Al-Bittar, and Soubra (2012) for different configurations of the rock parameters. The results are presented for the case of a weightless material. The value of the Poisson ratio adopted in this paper is 0.3. Table 1 presents a comparison between the results obtained from the present deterministic model and those given by Merifield, Lyamin, and Sloan (2006) and Mao, Al-Bittar, and Soubra (2012). It should be mentioned here that the results given by Merifield, Lyamin, and Sloan (2006) present the average values between the upper and lower bound solutions of the limit analysis theory. On the other hand, Mao, Al-Bittar, and Soubra (2012) present only an upper bound solution of the ultimate bearing capacity. Table 1 shows that the present numerical model provides slightly more critical values of the ultimate bearing capacity. This model will be used to perform the probabilistic analysis.   Merifield, Lyamin, and Sloan (2006) and by Mao, Al-Bittar, and Soubra (2012)

Probabilistic analysis
First, the discretisation of the NG (log-normal in the present paper) random field of σ c is briefly presented. It is followed by a brief presentation of the SPCE methodology used for the probabilistic analysis.

Discretisation of the random field
Consider a 2D NG random field Z NG described by (i) a NG (log-normal in the present paper) marginal cumulative distribution function F G for which μ and σ are, respectively, its mean and standard deviation values and (ii) a square exponential isotropic autocor- , which gives the values of the correlation between two arbitrary points (x, y) and (x ′ , y ′ ). Notice that this function is given as follows: where a is the autocorrelation distance. The EOLE method proposed by Li and Der Kiureghian (1993) to discretise a random field is used herein. In this method, one should first define a stochastic grid composed of s grid-points (or nodes) and determine the NG autocorrelation matrix S NG , which gives the correlation between each grid point of the stochastic mesh and the other grid-points of this mesh using Equation (6). The NG autocorrelation matrix S NG should then be transformed into the Gaussian space using the Nataf transformation (1962). As a result, one obtains a Gaussian autocorrelation matrix S G that can be used to discretise the Gaussian random field Z as follows: where μ ln Z and σ ln Z are the mean and standard deviation values of the underlying Gaussian random field (i.e. ln(Z )); (l j , f j ) are the eigenvalues and eigenvectors of the Gaussian autocorrelation matrix S G ; V is the correlation vector between the value of the field at an arbitrary point (x, y) and its values at the different grid-points and j j is a vector of N standard normal random variables, where N is the number of terms (expansion order) retained in the EOLE method. Notice that both f j and V are vectors of dimension s, where s is the number of the grid-points of the stochastic mesh. Notice also that the number N in Equation (7) is obtained (i) by sorting the eigenvalues l j ( j = 1, … , s) in a descending order and (ii) by choosing the number N of eigenmodes that leads to a variance of the error which is smaller than a prescribed value (10% in this paper). Notice finally that the variance of the error for EOLE is given as follows: where Z(x, y) andZ(x, y)are, respectively, the exact and the approximate values of the random field at a given point (x, y) and (f j ) T is the transpose of the eigenvector f j . Once the Gaussian random field is obtained, it should be transformed into the NG (log-normal in the present paper) space by exponentiating the approximated Gaussian random fieldZ(x, y) given by Equation (7).

SPCE for the system response
In this section, one first presents the PCE and then its extension, the SPCE. The PCE methodology allows one to replace a complex deterministic model whose input uncertain parameters are random variables by a metamodel. Thus, the random system response may be easily calculated (when performing the probabilistic analysis by MCS) using a simple analytical equation. Within the PCE methodology, the system response Γ of a given model with M random variables can be expressed by a PCE as follows: where P is the number of terms retained in the truncation scheme, j = {j i } i=1,...,M is a vector of M independent standard random variables that represent the M random variables, a β are unknown coefficients to be computed and C b are multivariate Hermite polynomials. These multivariate Hermite polynomials can be obtained from the product of one-dimensional Hermite polynomials as follows: where α i (i = 1, … , M ) are a sequence of M non-negative integers and H a i (.) is the a i th one-dimensional Hermite polynomial. The expressions of the onedimensional Hermite polynomials are given in Appendix 1. The coefficients a β of the PCE are computed in this paper using a non-intrusive technique where the deterministic calculations are done using the finite difference software FLAC 3D treated as a black box.
The most used non-intrusive method is the regression approach (Isukapalli, Roy, and Georgopoulos 1998;Huang, Liang, and Phoon 2009;Mollon, Dias, and Soubra 2009;Blatman and Sudret 2010;Li et al. 2011;Mao, Al-Bittar, and Soubra 2012;Soubra 2013, 2014;Jiang et al. 2014Jiang et al. , 2015. This method is used in the present work. In fact, for a PCE of order p, only the multivariate polynomials C b of degree less than or equal to p should be retained. This leads to a number P of the unknown PCE coefficients (see Equation (9)) equal to ((M + p)!)/(M!p!). It should be noticed that the number of the PCE coefficients to be computed grows dramatically with the size M of the input random vector and the PCE order p. When dealing with a random field as is the case in the present paper (and especially when considering small values of the autocorrelation distance), the discretisation of the random field by EOLE may lead to a significant number of random variables, which makes the determination of the PCE coefficients unfeasible because of the significant increase in the number of calls of the deterministic model. To address such problems, the SPCE developed by Blatman and Sudret (2010) is used herein. Indeed, Blatman and Sudret (2010) have shown that the number of significant terms in a PCE is relatively small since the multivariate polynomials C b corresponding to high-order interaction (i.e. those resulting from the multiplication of the H a i with increasing α i values) are associated with very small values for the coefficients a β . Based on this observation, these authors have proposed a so-called hyperbolic truncation scheme to determine the significant C b terms. This scheme suggests that the q-norm a q of the retained C b term should be less than or equal to the order p of the PCE. The q-norm is given by ( 1 1 ) where q is a coefficient (0 < q < 1). In this formula, q can be chosen arbitrarily. Blatman and Sudret (2010) have shown that sufficient accuracy is obtained for q ≥ 0.5. The proposed SPCE methodology leads to a SPCE that contains a small number of unknown coefficients, which can be calculated from a reduced number of calls of the deterministic model with respect to the classical PCE methodology. Notice that the SPCE methodology as proposed by Blatman and Sudret (2010) is based on an iterative procedure to arrive to a minimal number for the SPCE coefficients. This procedure is presented in a flowchart in Appendix 2. It is used in this paper to build up a SPCE for the system response. For more details on this procedure, the reader may refer to Al-Bittar and Soubra (2013Soubra ( , 2014 and Blatman and Sudret (2010). Once the coefficients a β are computed, the PDF of the system response and the corresponding statistical moments (mean μ, standard deviation σ, skewness δ u , and kurtosis κ u ) can be easily calculated with no additional cost using the meta-model. The two next subsections are, respectively, devoted to (i) the method used for the computation of the coefficients a β of the SPCE using the regression approach and (ii) the GSA based on Sobol indices.

Computation of the SPCE coefficients by the regression approach
. . , j M )} of the standard normal random vector ξ. These realisations are called experimental design (ED) and can be obtained from MCSs. We note G = {G(j (1) ), . . . , G(j (K) )} the corresponding values of the response determined by deterministic calculations. The computation of the SPCE coefficients using the regression approach is performed using the following equation: in which the matrix η is defined by where J is the number of the retained SPCE coefficients. Notice that in order to ensure the numerical stability of the treated problem in Equation (12), the size K of the ED must be selected in such a way that the matrix (h T h) −1 is well-conditioned. This implies that the rank of this matrix should be larger than or equal to the number of unknown coefficients. This test was systematically performed while solving the linear system of equations of the regression approach to arrive to the numerical stability. Notice finally that the quality of the output approximation via a SPCE closely depends on the SPCE order p. To ensure a good fit between the meta-model and the true deterministic model (i.e. to obtain the optimal SPCE order), one successively increases the SPCE order until a target accuracy was obtained. In this paper, the coefficient of determination Q 2 (see Blatman and Sudret 2010) is used. This coefficient is more efficient than the classical coefficient of determination R 2 since it allows one to check the capability of the meta-model of correctly predicting the model response at any point which does not belong to the ED.

Global sensitivity analysis
Once the SPCE coefficients are determined, a GSA based on Sobol indices can be easily performed. Notice that the first-order Sobol index of a given random variable ξ i (i = 1, … , M) gives the contribution of this variable in the variability of the system response. The first-order Sobol index is given as follows (Saltelli, Chan, and Scott 2000;Sobol ′ 2001): where Y is the system response (i.e. the ultimate bearing capacity in the present paper), E(Y|j i ) is the expectation of Y conditional on a fixed value of j i and Var denotes the variance. In the present paper, the system response is represented by a SPCE. Thus, by replacing Y in Equation (14) with the SPCE expression, one obtains the Sobol index as a function of the different terms of the SPCE as follows (Sudret 2008): where a b are the obtained SPCE coefficients, C b are the multivariate Hermite polynomials, E[.] is the expectation operator and E[(C b ) 2 ] is given by Sudret (2008) as follows: where the α i are the same sequence of M non-negative integers {a 1 , . . . , a M } used in Equation (10). It should be mentioned here that Equation (15) is similar to Equation (14), in which Y is the system response (i.e. the ultimate bearing capacity in the present paper). Equation (15) was obtained by replacing the system response Y by the PCE expression given by Equation (9). The mathematical derivation of Equation (15) was presented by Sudret (2008). Notice also that I i , which appears in the numerator of Equation (15), denotes the set of indices β for which the corresponding C b terms are only functions of the random variable ξ i (i.e. they only contain the variable ξ i ), and I i (i = 1, … , M) in the denominator of the same equation regroup all the indices β for which the corresponding C b terms are functions of all the random variables ξ i (i = 1, … , M).
In the present paper, where both a random variable (GSI) and a random field (σ c ) are involved, the Sobol index of the random field (σ c ) is computed as the sum of the Sobol indices of the different variables that represent this field.
In order to illustrate the construction of a PCE and the derivation of the equations providing Sobol indices, an example of a PCE of order p = 3 using only M = 2 random variables (ξ 1 and ξ 2 ) is presented in Appendix 1.

Numerical results
The aim of this section is to present the probabilistic numerical results obtained from the analysis at the ultimate limit state of a shallow strip footing resting on a HB spatially varying rock mass and subjected to a vertical load. The uniaxial compressive strength of the intact rock (σ c ) was considered as a log-normal random field characterised by a square exponential autocorrelation function. Its mean value and coefficient of variation (referred to in this paper as reference values) were taken as follows: m s c = 10 MPa, COV s c = 25%. On the other hand, GSI was modelled as a log-normally distributed random variable with a mean value and a coefficient of variation given as follows: m GSI = 25, COV GSI = 10%. As for the autocorrelation distance (a) of the random field σ c , the reference value adopted is 2 m; however, a range of 0.5-100 m was considered for the parametric study. For the different values of the autocorrelation distance (a), the total number N of random variables (or eigenmodes) that should be used to discretise the random field of σ c within the prescribed value of 10% for the variance of the error is presented in Table 2. Notice that the intact rock material constant (m i ) and the disturbance coefficient (D) were assumed to be deterministic since the probabilistic ultimate bearing capacity was found not sensitive to the variability of these parameters (Mao, Al-Bittar, and Soubra 2012). Their corresponding values were respectively m i = 8 and D = 0.3. Notice finally that the probabilistic results are presented in the case of a weightless rock mass. The deterministic numerical model was presented in the previous section. It should be noted that the size of a given element in the deterministic mesh depends on the autocorrelation distance of the spatially varying rock property. Der Kiureghian and Ke (1988) have suggested that the length of the largest element of the deterministic mesh in a given direction (horizontal or vertical) should not exceed 0.5 times the autocorrelation distance in that direction. In order to respect this criterion for the different values of the autocorrelation distance, a refinement of the deterministic mesh was performed in FLAC 3D for the very small values of the autocorrelation distance (<1 m).
The following subsections are organised as follows: first, the probabilistic results of the reference case (i.e. a = 2 m) are presented and discussed. This is followed by a presentation of the parametric study. The aim of this parametric study is to show the effect of the different governing statistical parameters on both the PDF of the ultimate bearing capacity and the Sobol indices of the uncertain parameters (i.e. σ c and GSI).

Probabilistic results for the reference case
The different steps of the probabilistic computation are presented herein for the reference case where a = 2 m (similar steps were done when performing the parametric study). These steps can be summarised by two stages as follows: In a first stage, discretise the random field σ c into a finite number of random variables using EOLE and determine a realisation of the two uncertain parameters σ c and GSI using M random variables as follows: . Define the stochastic grid: Li and Der Kiureghian (1993) have shown that the variance of the error (Equation (8)) is large at the boundaries of the stochastic domain. This problem can be solved by using a stochastic domain Ω RF that extends beyond the boundaries of the physical domain Ω. In this paper, a uniform stochastic grid of dimensions Ω RF = [21 m, 7 m] was used while the size of the physical domain was Ω = [20 m, 6 m]. On the other hand, Li and Der Kiureghian (1993) have shown that the number of grid-points in the stochastic grid strongly depends on the autocorrelation distances. These authors have shown that a ratio of about l RF /a = 1/5 provides a sufficient accuracy in terms of the variance of the error, where l RF is the typical element size in the stochastic grid, and a is the autocorrelation distance. In this paper, the number of grid-points in the stochastic grid was chosen as follows: six grid-points were considered within each autocorrelation distance with a minimum of six grid-points in that direction when the autocorrelation distance is larger than the size of the stochastic domain. Thus, a fine stochastic mesh was used for a highly heterogeneous soil and a coarse stochastic mesh was used for a slightly heterogeneous soil. . Calculate the autocorrelation matrix S NG using Equation (6) (notice that the dimension of this matrix depends on the value of the autocorrelation distance a). Then, compute its corresponding eigenmodes (i.e. eigenvalues and eigenvectors). Finally, identify for the random field σ c its N largest eigenmodes l j and w j (j = 1, … , N ), for which the variance of the error is smaller than a threshold of say 10%. As may be seen from Table 2, for smaller values of the autocorrelation distance of the random field, the number of random variables (which is equal to the number N of eigenmodes) increases. The total number of random variables M to be introduced in the probabilistic analysis is equal to the number N of terms retained in EOLE method (to discretise the random field σ c ) plus one (which represents the random variable GSI), that is, M = N + 1. . For each realisation, compute the values of σ c at the centroid of each element of the deterministic mesh using Equation (7), the value of GSI being the same for all the elements of the deterministic mesh and it is obtained by generating a realisation from the PDF describing this parameter. Figure 3 presents a realisation of the random field σ c . As may be seen from this figure, the isotropy in the autocorrelation distance is well presented within this realisation.
In a second stage, an arbitrary number of realisations K = 200 (with M random variables for each realisation) was performed using MCS technique to construct the ED. The number K is an initial (arbitrary) value since the iterative algorithm by Blatman and Sudret (2010) suggests to automatically add other simulations (an arbitrary number of 100 realisations was taken in this paper) each time the regression problem is ill-posed. The algorithm stops if either the target accuracy Q 2 TARGET is achieved or if the SPCE order p reached a maximal prescribed value fixed by the user. In this paper, a target accuracy Q 2 TARGET = 0.999, a value of q = 0.7 (see Equation (11)), and a maximal SPCE order equal to 5 were used. Notice that for the reference case (i.e. a = 2 m), the algorithm has stopped when the target accuracy Q 2 TARGET was reached. The corresponding order of the SPCE was equal to 3. In this case, where 35 random variables were needed for the random field σ c (see Table 2) and a single random variable was used to represent GSI, the classical PCE leads to P = 9139 unknown coefficients and thus a minimum of 9139 calls of the deterministic model were needed to accurately represent the ultimate bearing capacity by a meta-model. By using the SPCE iterative algorithm by Blatman and Sudret (2010), only P = 127 unknown coefficients were retained and only 500 calls of the deterministic model were found to be largely sufficient to construct the meta-model. Consequently, a significant reduction in the number of calls of the deterministic model can be obtained using the SPCE. This greatly facilitates the solution of problems involving random fields.
In order to show the efficiency of the proposed SPCE methodology, a comparison with the crude MCS method was performed. A total number of simulations of 10,000 was considered in the computation. This number is largely sufficient to lead to convergence of the mean, standard deviation, skewness and kurtosis of the ultimate bearing capacity as may be seen from Figure 4. Figure 5 presents the histogram obtained from the 10,000 MCSs (performed on the original expensive finite difference code FLAC 3D ) and the PDF obtained using the SPCE metamodelling technique based on only 500 simulations. As may be seen from this figure and from Table 3, the SPCE metamodelling technique is an efficient alternative to the crude MCS technique because the results show good agreement between the two methods (in terms of the mean, standard deviation, skewness and kurtosis of the ultimate bearing capacity) with a much smaller number of calls of the deterministic model when using the SPCE. Figure 6 depicts the values of Sobol indices as given by the obtained SPCE for (i) the random variable GSI and (ii) the 35 random variables representing the random field σ c . The first random variable ξ 1 corresponds to GSI and its Sobol index was found to be equal to 0.66. However, the last 35 random variables [i.e. ξ i for i = 2, … , 36] are those corresponding to the σ c random field. The sum of their Sobol indices gives the weight of the random field σ c in the variability of the ultimate bearing capacity. This sum was found to be equal to 0.34. Figure  6 shows that only six random variables (ξ 2 , ξ 4 , ξ 6 , ξ 8 , ξ 9 , ξ 12 ) of the σ c random field are the most influential (they involve 89% of the variance of σ c ). This can be explained by the fact that the system response (i.e. the ultimate bearing capacity) is a quantity that depends on the average distribution of the spatially varying rock property over the entire domain and it is therefore quite insensitive to small-scale fluctuations of σ c . Notice that the first eigenmodes provide the average distribution of σ c over the rock domain; however, the remaining eigenmodes give the small-scale fluctuations around this average distribution.

Probabilistic parametric study
The aim of this section is to study the effect of the different statistical governing parameters (autocorrelation distance of σ c and coefficient of variation of both σ c and GSI) (i) on the PDF of the ultimate bearing capacity and (ii) on Sobol indices. Figure 7 provides the PDFs of the ultimate bearing capacity for different values of the autocorrelation distance of σ c (a = 0.5, 1, 2, 5, 10, 50, 100 m) and for the case where σ c is modelled as a random variable (case of a homogenous rock mass). Table 4 presents the four statistical moments and the number of needed simulations to construct the meta-model using the SPCE technique for the cases presented in the same figure. The PDF and the statistical moments corresponding to a great value of the autocorrelation distance (a = 100 m) are similar to those given by the case of a random variable. This is because the case of a random variable can be considered as the limiting case of a random field with an infinite value of the autocorrelation distance. Figure 7 shows that the PDF is less spread out when the autocorrelation distance decreases. For the very large values of the autocorrelation distance (i.e. a = 100 m), the coefficient of variation of the ultimate bearing capacity tends to a constant maximal value (see Table 4), which is the value corresponding to the case of a random variable. In this case, the different values of σ c (at the different locations in a given realisation) are perfectly correlated. This means that for a given simulation, a single value of σ c is affected to the entire rock mass. This value of σ c is chosen according to the prescribed PDF of σ c and it may vary (from a realisation to another one) in the range of values imposed by this PDF. This observation is also valid for GSI. This leads to a large variability of the ultimate bearing capacity. To conclude, the large value of the variability is because one obtains a large variety of homogeneous rock masses with low, intermediate and high values of σ c and GSI from simulation to another one. The decrease in the autocorrelation distance of the random field σ c from infinity to a finite value (moderate or small where a ≤ 10 m) limits the correlation between the values of this random field (in a given simulation) to a finite zone which leads to several zones with different values of σ c over the entire rock mass. This means that in a single simulation, one obtains a set of weak and strong zones for which the position may change from one simulation to another. The case of moderate-to-small values of the autocorrelation distance leads to a decrease in the variability of the   ultimate bearing capacity since the presence of the rock mass heterogeneity (weak and strong zones) will produce a somewhat close global behaviour of the footing because of the averaging phenomenon over the zone of possible failure mechanism. The decrease in the variability of the ultimate bearing capacity becomes the most significant for the case of a very small value of the autocorrelation distance because the rapid change in the values of σ c from one element to another of the mesh leads to quasi-similar ultimate bearing capacity for all the realisations. For these very small values of the autocorrelation distance, the rock mass can be considered as having a somewhat deterministic value of σ c with only GSI being a random variable. Figure 8 and Table 4 show that the probabilistic mean value of the ultimate bearing capacity presents a minimum when the autocorrelation distance is nearly equal to the footing breadth B (i.e. in our case when a = 1 m). Notice that the minimal probabilistic mean was also observed by other authors (Griffiths, Fenton, and Manoharan 2002;Fenton and Griffiths 2003) when considering the bearing capacity of foundations resting on a soil mass. For the very large values of the autocorrelation distance (a = 100 m), the probabilistic mean tends to the one of the homogenous rock mass (case of random variables) as may be seen from Table 4. On the other hand, for the very small values of the autocorrelation distance, the probabilistic mean becomes greater than the minimal value because the weakest path becomes increasingly tortuous and its length is also longer (cf. Fenton and Griffiths 2003). As a result, the failure mechanism will start to look for shorter path cutting through higher values of σ c . Table 4 also shows the impact of the autocorrelation distance on both the skewness and the kurtosis of the PDF. For small values of the autocorrelation distance, the skewness and kurtosis of the response are small, which means that the PDF of the response is not far from a Gaussian one in these cases. Notice however that these moments increase when the autocorrelation distance increases, which means that for great values of the autocorrelation distance, the shape of the PDF of the output becomes far from a Gaussian one. Finally, Figure 9 and Table 5 show the effect of the autocorrelation distance on the Sobol indices of the random field σ c and the random variable GSI. The results show that for very large values of the autocorrelation distance (i.e. a = 100 m), the variability of the ultimate bearing capacity is mainly due to σ c . Similar results were obtained by Mao, Al-Bittar, and Soubra (2012), where   Figure 8. Influence of the autocorrelation distance on the probabilistic mean of the ultimate bearing capacity.

Effect of the autocorrelation distance
the uncertain parameters were modelled by random variables. It should be emphasised here that σ c is the most weighted parameter in the variability of the ultimate bearing capacity only in the case of very large values of the autocorrelation distance and in the case of random variables. Indeed, Figure 9 shows that the decrease in the autocorrelation distance of σ c reduces its weight in the variability of the ultimate bearing capacity and increases the weight of GSI. Although this result was impossible to be detected when a simplified modelling (i.e. random variables) of the uncertain rock parameters was used, it can be explained by the fact that the small values of the autocorrelation distance increase the rock mass heterogeneity (i.e. one obtains a set of weak and strong zones), which will produce a somewhat close global behaviour of the footing from one simulation to another because of the averaging phenomenon over the zone of possible failure mechanism. The expected decrease in the variability of the ultimate bearing capacity with the decrease in the autocorrelation distance of σ c is reflected herein by a decrease in the weight of σ c in the variability of this response. For the limiting case of a very small value of the autocorrelation distance, σ c can be seen as a deterministic value, which implies that in this case the variability of the ultimate bearing capacity is only due to GSI (i.e. S(GSI) tends to one).

Effect of the coefficients of variation
The effect of the coefficients of variation (COVs) of the random field σ c and the random variable GSI is studied Figure 9. Influence of the autocorrelation distance on the Sobol indices of GSI and σ c .  Figure 10. Influence of the coefficients of variation COVs of the random variable GSI and the random field σ c . on the PDF of the ultimate bearing capacity: (a) influence of COV(GSI); (b) influence of COV(σ c ). and presented in Figure 10 and Tables 6 and 7. Notice that in this study, the adopted value of the autocorrelation distance of the random field σ c is the reference value (i.e. a = 2 m). Figure 10 and Table 6 show that the variability of the ultimate bearing capacity increases when the coefficient of variation of either the random field σ c or the random variable GSI increases, with the increase being more significant for the GSI parameter (see Table 6). This is because an increase in the coefficient of variation of σ c by 50% (with respect to the reference case where COVσ c = 25%) increases the COV of the ultimate bearing capacity by only about 13.9%, while increasing the coefficient of variation of GSI by 50% (with respect to the reference case where COV GSI = 10%) increases the COV of the ultimate bearing capacity by about 34.9%. Table 6 also shows that the probabilistic mean value of the ultimate bearing capacity slightly decreases when the coefficients of variation increase. From Table 7, one can see that an increase in the coefficient of variation of a rock parameter increases its Sobol index and thus its weight in the variability of the ultimate bearing capacity; this automatically reduces the contribution of the other uncertain parameter. This increase is more significant for σ c . This is because an increase in the coefficient of variation of GSI by 50% (with respect to the reference case where COV GSI = 10%) increases its Sobol index by only about 24.3%, while increasing the coefficient of variation of σ c by 50% (with respect to the reference case where COVσ c = 25%) increases its Sobol index by about 50%.

Conclusions
A probabilistic analysis of a vertically loaded strip footing resting on a spatially varying rock mass has been performed to compute the PDF of the ultimate bearing capacity. The rock was assumed to follow the generalised HB failure criterion. The uniaxial compressive strength of the intact rock (σ c ) was considered as a random field and the Geological Strength Index (GSI) was modelled as a random variable. The aim of this paper is to propagate the uncertainties from the input parameters (GSI and σ c ) to the system response (i.e. the ultimate bearing capacity). The uncertainty propagation methodology employed in the analysis makes use of a non-intrusive approach to build up a SPCE for the system response.
In addition to the uncertainty propagation, a GSA based on Sobol indices was performed. The deterministic model was based on numerical simulations using FLAC 3D software. The probabilistic numerical results were presented in the case of a weightless rock mass. A good agreement with the probabilistic results obtained from MCS was observed in terms of the mean, standard deviation, skewness and kurtosis of the ultimate bearing capacity. The parametric study has shown that the variability of the ultimate bearing capacity increases with the increase in the coefficients of variation of the rock parameters σ c and GSI, with the increase being more significant for the GSI parameter. Sobol indices have shown that for the very large values of the autocorrelation distance, the variability of the ultimate bearing capacity is mainly due to σ c ; however, in the case of very small values of the autocorrelation distance, GSI is the most weighed variable. It was also shown that an increase in the coefficient of variation of a rock parameter (σ c or GSI) increases its Sobol index and thus its weight in the variability of the system response and decreases the weight of the other parameter; the increase is more significant for σ c . With a decrease in the autocorrelation distance of the uniaxial compressive strength of the intact rock (σ c ), a less spread-out PDF of the ultimate bearing capacity was obtained; however, the probabilistic mean value of the ultimate bearing capacity presents a minimum. This minimum was obtained when the autocorrelation distance is equal to the footing breadth B. Small values of the autocorrelation distance lead to small values of the skewness and the kurtosis of the system response. Thus, a PDF of the system response that is not far from a Gaussian one was obtained in these cases.

Disclosure statement
No potential conflict of interest was reported by the authors.