Stress-weighted spatial averaging of random fields in geotechnical risk assessment

Reliability-based design methods have recently been widely applied to geoengineering, in which the spatial variability of soil parameters (random fields) plays an important role. Early studies from the late 60s focused on random stability of slopes. The first models considered soil strength parameters as spatial random fields; however, simple random variables are more useful in practical applications (EN 1997-1; Low and Phoon 2015, Tietje et al. 2011). Recent risk analyses widely implement random data as spatial random fields of auto- and crosscorrelated parameters or even use complete probabilistic definition with joint probability distributions (Ching et al. 2016, Javankhoshdel and Bathurst 2015, Griffith and Fenton 2001, Jiang et al. 2015, Shen et al. 2021). Data acquisition and operations on such spatial data require sophisticated numerical support, but – alternatively – some techniques of data pre-processed simplifications are also effective; the former direction can be found in some selected papers attached as references (Cho 2007, Deng et al. 2017, Farah et al. 2015, Ji et al. 2018, Jiang et al. 2014, Kim and Sitar 2013, Li et al. 2017, Liu et al. 2017, Liu et al. 2018, Shen 2021, Tietje et al. 2011), while the latter direction is presented hereafter.

Another group of boundary problems analyses random bearing capacity of footings in a random field context focusing on the most important regions, where shear resistance is generated along slip lines or slip surfaces. The undrained conditions are often analysed (frictionless materials, the Tresca model), which do not depend on stresses (Chwała 2019, Griffith and Fenton 2001, Huang et al. 2013, Li et al. 2015, Puła and Chwała 2015, Shen et al. 2021). For frictional or frictional–cohesive materials, the situation is more complex because stress fields cannot be ignored (Fenton and Griffith 2003, Liu et al. 2017).

Modelling of input data as spatial random fields complicates risk analyses; therefore, advanced analytical or numerical methods are widely in use. The fundamental difficulty is that such sophisticated models, often parameter sensitive, should be in a balance with high quality of geoengineering input data and the data look never complete; that is why, for practical design situations, the random models are sometimes reduced from spatial random fields to simple random variables (Vanmarcke 2010, Low and Phoon 2015) by making use of some averaging techniques. Note that the spatial averaging does not necessarily mean that any significant part of information is lost. In a context of interactions between adjacent regions in a subsoil mass, significant up and down fluctuations of random parameters can reduce each other; so, the volume averaging of spatially distributed soil parameters is rationally justified and economical – it reduces point variances of considered random fields (Tietje et al. 2011). In addition, the method of spatial averaging has a strong mathematical background (Vanmarcke 2010) as well as a variety of well-documented applications.

The paper focuses mainly on the subsoil shearing resistance – like τ = q_n×tan(φ) for the frictional Coulomb material – and analyses the variance reduction coefficients γ_a for the homogeneous random field tan(φ(x)), which is reduced by averaging to a corresponding random variable [tan(φ)]_a. Fields of stresses in the subsoil mass implement a new element to the averaging procedures because non-uniform (deterministic) normal stresses q_n ≠ const(x) generate a kind of subsoil inhomogeneities, in terms of the random shearing resistance τ(x). Indeed, if a normal stress q_n along a slip line is locally negligible, then τ = q_n×tan(φ) is also locally negligible and such a region does not contribute to total bearing capacity. If locally q_n ∼ 0 in shearing resistance, then the friction coefficient tan(φ) is locally not very important and especially its random fluctuations Δtan(φ) can be ignored; also, it is vice versa for a locally concentrated load q_n ≫ 0. Vanmarcke's spatial averaging applied to soil strength parameters (or to soil deformation parameters) ignores this mechanism and this is controversial. Clearly, the Vanmarcke averaging of spatial random fields is effective for many ‘stress-neutral’ random parameters like radioactive radiation, pollutant concentration, temperature, etc. For example, the spatial geometrical averaging is very useful in hydrological studies – hydrologists evaluate average rainfalls over a certain geographical region to estimate flood danger in an entire river basin (Vanmarcke 2010); but this model is insufficient for ultimate limit states, which are based on a random field tan(φ(x)) because the Coulomb shearing resistance τ is additionally ‘disturbed’ by spatially distributed normal loads q_n(x).

It is assumed that the spatial field of the normal stress q_n (x) is deterministic and physically independent of the random field tan(φ(x)). The stress field q_n(x) plays the role of a fourth source of inhomogeneity – apart from the inherent one, from deterministic trends and measurement/transformation errors – because the stress field q_n(x) can amplify (or suppress) effects of random local fluctuations. A new point of view on the spatial averaging of random fields in geoengineering is presented below, and the paper's objectives focus on simple numerical examples. The paper implements a more general and coherent concept of averaging of random fields and allows a conclusion that the standard (purely geometrical) variance reduction based only on volume integrals can be too optimistic; in a context of geotechnical structural safety, the traditional volume averaging can be dangerous.

Only few papers consider in a distinct way volume averaging in conjunction with fields of stresses, though there is a widespread agreement that the spatial regions being close to loaded places play a pivotal role in the averaging; being ‘close to loaded places’ means not only a short distance, but also local stress components of a maximal intensity. Recently, it was found by Ching and Hu (2017) and Ching et al. (2016) that the traditional spatial averaging model that treats all soil regions equally important cannot satisfactorily represent the effective Young modulus for the footing settlement problem – highly mobilised soil regions close to a footing are more significant than non-mobilised ones; the authors proposed some weight coefficients derived from a random finite element analysis. Similarly, for random shearing resistance of a rectangular sample analysed via finite elements (Ching and Phoon 2013), the authors found that the averaging over the whole sample overestimates the variance reduction, whereas the averaging along slip lines is much more representative; such a conclusion for localised shearing resistance should not be surprising. In the next paper, the authors confirmed this conclusion for a wider class of boundary problems (Ching et al. 2016).

Soils display inherent inhomogeneities caused by long-term geological processes, unknown stress history, soil moisture, chemical reactions, etc.; therefore, soil parameters have a random nature, which can be described in terms of random variables or spatial random fields (time effects are not considered hereafter). Such local random fluctuations that happen in situ can be amplified during soil parameter evaluation, that is, by soil testing (usually invasive), measurement errors, calibration procedures, transformation models and others (Phoon and Kulhawy 1999a, 1999b); therefore, in geotechnical risk assessment, the term uncertainty of geotechnical parameters is more adequate than the term inherent inhomogeneity. Reported fluctuations of geotechnical parameters are relatively great, much greater than in cases of man-made construction materials. Considering a subsoil layer as a macro-homogeneous mass, the uncertainty can be expressed using coefficients of variations (c.o.v.) called ν = σ/μ, so a standard deviation σ over an expected value μ > 0. This distribution-free parameter is the minimal information about subsoil uncertainty; it can be sufficient only in simple practical applications.

Numerous authors present the results, literature reviews and discussions on ‘typical’ values of c.o.v. (%), like the values reviewed by Phoon and Kulhawy (1999a). Clearly, the data from many countries, regions, geological units, climate zones and even from different geological companies differ substantially and may be incomparable. This explains inevitable discrepancies in published data. Reported c.o.v. are sometimes based on estimated inherent inhomogeneity, but sometimes they are based on a total uncertainty that includes measurement errors, etc. Moreover, algorithms of de-trending can be different, non-standard testing procedures and individual data transformation methods are sometimes implemented (Phoon and Kulhawy 1999a). For an illustrative character of this paragraph, generalised ‘most typical’ values can be derived based on the review papers by Phoon and Kulhawy (1999a, 1999b): for the internal friction coefficient ν_tan(φ) ∼ 5%–15% (non-cohesive soils), ν_tan(φ) ∼ 25% (cohesive soils, even up to 50% is possible; slightly less for the angle φ itself), for the cohesion ν_c ∼ 20%–30% (for undrained shear strength, even up to 70% is possible) and for the deformation modulus ν_D ∼ 30%–50%. The presented values refer to the mean values of inherent variability; the intervals for global uncertainty can be still wider, at least 10%–20% due to measurement errors.

The following negative crosscorrelation coefficient is usually reported: −0.75 < r_tan(φ);c < 0 and the value r_tan(φ);c ∼ −0.3 can be recommended as a conservative estimation. Note that deterministic physical arguments do not always confirm the assumed negative crosscorrelation between φ or tan(φ) and c. Indeed, local fluctuations of soil moisture cause simultaneous decrease (or increase) of both tan(φ) and c; this fact can support suggestion about a positive correlation. On the other hand, local fluctuations of fine particle content result in random changes of tan(φ), which are simultaneously opposite to changes of c; this fact can support the conclusion about a negative correlation. Tests on drained and undrained samples also suggest a negative correlation (c_u > c′, φ_u ∼ 0^o < φ′). A review of previous papers and study on the role of negative crosscorrelation can be found in Javankhoshdel and Bathurst (2015).

Reliability-based design methods that use only the random variable format can cause an overdesign of foundations or other geoengineering structures. There are several possible reasons for this: ultimate limit state conditions may be oversimplified (too restrictive), some improvements of soil parameters may happen during construction, the database about geotechnical failures may be still poor – or the values of the c.o.v. may be ‘too large’, even though they are correct as point estimations. Indeed, it is a well-established fact that point variances Var{P(x)} = σ² = const(x) of a homogeneous random field of soil parameters P(x) are usually not representative (overestimated) in geotechnical safety analyses because they ignore the interactions between adjacent regions in the considered soil mass. Bearing in mind dimensions of such adjacent regions – expressed in metres or dozens of metres – local random fluctuations can compensate each other in a considered soil mass as a whole, if total (summarised) effects are of interest. Formally, these are autocorrelation and crosscorrelation functions (between different parameters), which cause such reduction effects. This observation is a pivot of the monographies by Vanmarcke (2010, 1983) and numerous next papers.

Simple reliability models based on random variables that use only the expected values μ and standard deviations σ – instead of spatial random fields of subsoil parameters – are the most popular tool in safety evaluation, including the Eurocodes (EN 1997-1). In this context, a special attention is paid to averaging techniques called as spatial random homogenisation or volume integrals, which can transform spatially fluctuating random fields to some equivalent random variables. Since for a homogeneous random field, the expected value E{P(x)} = μ = const, any spatial averaging will not change this value, but the averaging yields an ‘effective’ variance Var_a = γ_a×σ² with a reduction coefficient γ_a ≤ 1; equivalently, but in terms of standard deviations σ_a = Γ_a×σ, therefore, the symbol Γ_a = √γ_a is also in use (Vanmarcke 2010, Tietje et al. 2011).

Consider a homogeneous random field of a soil parameter (or several soil parameters in a vectorial version) denoted as P(x), where a point x = (x₁, …, x_n) belongs to a domain D ⊂ Eⁿ in n-dimensional Euclidean space. The assumed homogeneity means that for every two points x, y, the expected values and the variances of two random variables P(x), P(y) are the same, E{P(x)} = E{P(y)} = μ = const, Var{P(x)} = Var{P(y)} = σ² = const < + ∞; the covariance function equals Cov{P(x);P(y)} = E{[P(x) − μ]×[P(y) − μ]} = σ²×ρ(x − y), so the corresponding correlation coefficient ρ depends only on the difference Δx = x − y or simply ρ = ρ(Δx). Note that the approach is also effective for some quasi-homogeneous random fields when homogeneous random fluctuations happen around a functional deterministic trend μ + f(x), like for a subsoil stiffness and strength which often increase with depth (Li et al. 2015). In principle, the homogeneity requirement concerns the homogeneous random fluctuations and not so much deterministic trends; so, during the pre-processing, the model should be de-trended (Bagińska et al. 2016, Li et al. 2015). Another technique is used by Shen et al. (2021); the authors recommend a modulating amplitude A(x), which is a deterministic function, but modifies the considered random field of subsoil parameters P(x) = A(x)×H(x). The randomness is still generated by a homogeneous random field H(x). In particular, subsoil resistance derived from vertical cone penetration tests (CPT) can have the form of a stochastic process q_c(z) with increasing both the expected value and the standard deviation (proportionally), but the model is not more complicated.

Joint probability distributions of the random field P(x) are not specified, if the simplest distribution-free second-order model is used. Moreover, it is assumed that the probabilistic averaging (ensemble averaging) and spatial averaging (volume averaging) coincide. This is only a useful hypothesis that cannot be proved in frames of distribution-free theories. In geotechnical practice of soil testing, one cannot carry out any sequence of independent measurements at the same point – the list of fully non-invasive geotechnical tests is very short. Replacing the term ‘at the same point’ by a more liberal ‘at the same place’ can be controversial; this difficulty leads to multiscale models, ambiguous definitions of ‘representative volume elements’ (Ostoja-Starzewski 2006) and nested micro-meso inhomogeneous models (Jaksa et al. 1999, Ching and Phoon 2013), which are not discussed here.

Several types of the decay functions ρ(Δx) became popular; four of them are discussed in detail by Vanmarcke (2010), Bagińska et al. (2016) and Oguz et al. (2019), the squared exponential (Gaussian) in particular (equation 1a). Due to natural sedimentation processes over vast areas, the geological model of the subsoil can reveal a horizontally layered structure like for transversally isotropic continua; hence, the following squared exponential representation of the correlation function ρ(x;y) = ρ(x − y) is useful: (1a) ρ=exp{−(x1−y1)2+(x2−y2)2(dh)2}⋅exp{−(x3−y3)2(dv)2} \rho = \exp \left\{ { - {{{{\left( {{x_1} - {y_1}} \right)}^2} + {{\left( {{x_2} - {y_2}} \right)}^2}} \over {{{\left( {{d_{\rm{h}}}} \right)}^2}}}} \right\} \cdot \exp \left\{ { - {{{{\left( {{x_3} - {y_3}} \right)}^2}} \over {{{\left( {{d_{\rm{v}}}} \right)}^2}}}} \right\} where the symbols x₃, y₃, d_v are used for vertical direction. Alternatively, the linear exponential (Markovian) function (equation 1b) is in use: (1b) ρ=exp{−|x1−y1|+|x2−y2|dh′}⋅exp{−|x3−y3|dv′} \rho = \exp \left\{ { - {{\left| {{x_1} - {y_1}} \right| + \left| {{x_2} - {y_2}} \right|} \over {{d_{{\rm{h'}}}}}}} \right\} \cdot \exp \left\{ { - {{\left| {{x_3} - {y_3}} \right|} \over {{d_{{\rm{v'}}}}}}} \right\} Although the equations (1a) and (1b) apply to general 3D situations, the points of interest x, y are often situated on certain specific subsets like along slide lines (surfaces) in ultimate limit states or limit equilibrium studies.

In data processing, sometimes a negative autocorrelation can happen that is beyond the scope of the exponential functions (1a and 1b); if such negative values are not incidental, and so if they are confirmed by significance tests, then a wider class of cosine exponents like exp{−|Δx|/d′}×cos(|Δx|/d″) can be considered. Such a ‘more chaotic’ behaviour seems to be possible for anthropogenic soils, for example for lignite mine dumping soils (Bagińska et al. 2016) or low-quality sandy backfills (Brząkała, 1981); in the former case, averaged locally negative values up to ρ ∼ −0.2 are reported (but not used, assumed d″ = +∞) and in the latter one, the values are up to ρ ∼ −0.5 (horizontal variability of a subsoil vertical response over a meso-scale laboratory field 4.2 m × 2.0 m, several dozens of testing points). The problem of negative autocorrelation for nearshore sea bottom soils is recently discussed by Oguz et al. (2019). The values derived from CPT tests in stiff, overconsolidated clays can also be negative (up to ρ ∼ −0.2) (Jaksa et al. 1999). In the context of risk assessment and averaging, note that ignoring the negative values of ρ is usually on the safe side (Fig. 1b).

Figure 1

The decay constants d_v, d_h (m) in equation (1a) control the speed of convergence to zero for increasing distances and they are in close relation to so-called scales of fluctuations δ_i (Vanmarcke 2010, Bagińska et al. 2016, Pieczyńska-Kozłowska et al. 2017, Tietje et al. 2011), which are good intuitive measures of random fluctuations. The following definitions of δ and δ′ are in use for unidirectional random variations due to equations (1a) and (1b), respectively: (2a) δ=∫−∞+∞ρ(r)dr=∫−∞+∞exp(−r2d2)dr=πd \delta = \int_{ - \infty }^{ + \infty } {\rho (r)dr = \int_{ - \infty }^{ + \infty } {\exp \left( { - {{{r^2}} \over {{d^2}}}} \right)dr = \sqrt \pi d} } (2b) δ′=∫−∞+∞ρ(r)dr=2∫0+∞exp(−rd′)dr=2d′ \delta ' = \int_{ - \infty }^{ + \infty } {\rho (r)dr = 2\int_0^{ + \infty } {\exp \left( { - {r \over {d'}}} \right)dr = 2d'} } Empirical (discrete) correlation function ρ_e(r) can be applied to approximate a distance r_o such that ρ_e(r_o) ≈ 0. In terms of both exponential functions ρ(r), this could be assumed as ρ(r) ≅ 0.05. If so, then r_o=√3d and r_o=3d‘; thus, d ∼ √3d‘. In this way, the following working hypothesis is derived: δ ∼ 1.5δ‘. A wider statistical analysis of the distance r_o and its confidence limits were presented by Jaksa et al. (1999). For d ≫ 0, the model turns to a degenerate spatial random field – just a random variable; for d ∼ 0, the model turns to a very chaotic random field with uncorrelated succeeding values.

There exist several sophisticated methods of evaluation of the δ-value; selected techniques and numerical results were presented by Bagińska et al. (2016), Lloret-Cabot et al. (2014), Pieczyńska-Kozłowska et al. (2017), Phoon and Kulhawy (1999a) and Vanmarcke (2010). Standard CPT testing is most often used – in vertical direction, first of all (Oguz et al. 2019); a very original horizontal CPT testing has been presented by Jaksa et al. (1999). Although the CPT method provides a very large amount of data – continuous realisations of the soil resistance, as expected for a stochastic process q_c(z) – the interpretation of results encounters difficulties.

In the common opinion δ_h ≫ δ_v, like the ratio δ_h/δ_v = 10 used by Tietje et al. (2014); such proportion δ_h ≫ δ_v can be suggested by a regular sedimentation of geological units. Again, for the same reasons as in cases of c.o.v., there is no agreed opinion about representative values of the scales of fluctuations d. Moreover, in published data, the distinction between definitions in equations (2a) and (2b) is not always respected. For example, results of CPT soundings have been presented with the best fitted Gaussian curve (equation 1a) in the paper by Jaksa et al. (1999) – the authors report measured cone tip resistance q_c for which δ_v ∼ 0.17 m (average value) with a c.o.v. on the level of 30%. Independently, another value r_o ≅ 0.28 m can be estimated from the authors’ empirical plot of ρ_e(r), used in place of the least-square fitting; in this way, the following values can be derived: δ_v = √πd_v=√πr_o/√3 ≅ r_o ≅ 0.28 m. Similar estimation δ_v ∼ 0.8×5.3/20 = 0.21 m can be found from the same data using the Vanmarcke proposal δ_v ∼ 0.8d, where d denotes the average distance between intersections of the fluctuating resistance q_c and its parabolic trend function (Vanmarcke 2010, Jaksa et al. 1999). The horizontal scale of fluctuations δ_h was also analysed by Jaksa et al. (1999) and its nested structure was documented – the estimated results depend on the scale of sampling: δ_h ∼ 0.15 m for CPT sampling at 5 mm intervals (so, δ_h ≅ δ_v), but δ_h ∼ 1–2 m for CPT sampling at 0.5 and 1.0 m. The authors’ conclusion says that ‘both the small-scale vertical and lateral correlation distances are a manifestation of the CPT and not the soil itself’ (Jaksa et al. 1999). On the other hand, the review by Phoon and Kulhawy (1999a) is in sharp contrast to the conclusions presented by Jaksa et al. (1999), but the former review is focused on a wider spectrum of testing methods, not only on the CPT sounding. For different geotechnical parameters, as well as for different testing techniques, there is generally δ_v ∼ 1.0–2.0 m (less than 1.0 m only for CPT tests), but δ_h is expressed in dozens of metres; so, the ratio δ_h/δ_v ∼ 10 is confirmed. The references presented by Shen et al. (2021) for the offshore soil in the North Sea report δ_v ∼ 0.48–7.14 m and δ_h ∼ 24.62–66.48 m.

Standard spatial averaging, called as a geometric one, has a well-documented mathematical background presented by Vanmarcke (2010); also, it has many practical applications (Vanmarcke 2010, Chwała 2019, Tietje et al. 2011, Ching et al. 2016, Puła and Chwała 2015). For a frictional material, where tan(φ) > 0 and c = 0 (kPa), consider two intervals dL extracted from a sliding line and assume the two random variables shown in Fig. 1a: a random variable tan(φ₁) on the left part dL and a random variable tan(φ₂) on the right one.

Looking for one ‘effective’ parameter, the random vector (tan(φ₁);tan(φ₂)) is replaced by one averaged random variable tan(φ)_ga acting along 2dL. Assume the homogeneity condition: E{tan(φ_i)} = μ, Var{tan(φ_i)} = σ²; moreover, Cov{tan(φ₁);tan(φ₂)} = σ²×ρ_1;2. The values dL are geometric weighting coefficients and the intuitively averaged (homogenised) friction coefficient [tan(φ)]_ga on 2dL can be calculated as a new random variable: (3) [tan(φ)]ga=[tan(φ1)⋅dL+tan(φ2)⋅dL]/(2dL)=12⋅tan(φ1)+12⋅tan(φ2) \matrix{ {{{\left[ {\tan (\varphi )} \right]}_{{\rm{ga}}}} = \left[ {\tan ({\varphi _1}) \cdot dL }\,+ \right.} \cr {\left. {\tan ({\varphi _2}) \cdot dL} \right]/(2dL) = {1 \over 2} \cdot \tan ({\varphi _1}) + {1 \over 2} \cdot \tan ({\varphi _2})} \cr } Therefore, E{[tan(φ)]_ga} = μ, Var{[tan(φ)]_ga} = σ²×(1 + ρ_1;2)/2 < Var{tan(φ_i) = σ². The variance reduction factor equals γ_ga = Var{[tan(φ)]_ga}/Var{tan(φ_i) = (1 + ρ_1;2)/2 = 1 −(1 − ρ_1;2)/2; the lower index g refers to the purely geometric sense of equation (3). Note that potentially, the randomness of [tan(φ)]_ga in equation (3) can even disappear, if always both random fluctuations of tan(φ₁) and tan(φ₂) happen as opposite numbers, thus taking the hypothetic correlation coefficient ρ_1;2 = −1.

Generalisations of the geometric averaging are straightforward, if for a longer interval AB, the sum in equation (3) is replaced by an integral over this line AB: (4) [tan(φ)]ga=∫ABtan(φ(x))dx|AB|=∫ABtan(φ(x))dx∫AB1dx {\left[ {\tan (\varphi )} \right]_{{\rm{ga}}}} = {{\int_A^B {\tan \left( {\varphi (x)} \right)dx} } \over {\left| {AB} \right|}} = {{\int_A^B {\tan \left( {\varphi (x)} \right)dx} } \over {\int_A^B {1dx} }} where the interval length equals |AB|; the definition is similar for surfaces in spatial cases (Vanmarcke 2010).

Ultimate limit state for random shearing is a start point of this approach due to the context of the random vector (tan(φ₁);tan(φ₂)). For the frictional material, the two variables tan(φ₁), tan(φ₂) should be considered in conjunction with two normal deterministic loadings in Fig. 1a, q_n1 = const > 0, q_n2 = const > 0. The random shearing resistance along AB can be expressed as: (5a) 2dT=qn1⋅tan(φ1)⋅dL+qn2⋅tan(φ2)⋅dL 2dT = {q_{n1}} \cdot \tan ({\varphi _1}) \cdot dL + {q_{n2}} \cdot \tan ({\varphi _2}) \cdot dL Looking for one ‘equivalent’ random variable [tan(φ)]_sa on the same interval AB: (5b) 2dT=qn1⋅[tan(φ)]sa⋅dL+qn2⋅[tan(φ)]sa⋅⋅ dL=[tan(φ)]sa⋅[qn1+qn2]⋅dL \matrix{ {2dT = } \hfill & {{q_{n1}} \cdot {{\left[ {\tan \left( \varphi \right)} \right]}_{{\rm{sa}}}} \cdot dL + {q_{n2}} \cdot {{\left[ {\tan \left( \varphi \right)} \right]}_{{\rm{sa}}}} \cdot } \hfill \cr {} \hfill & { \cdot \;dL = {{\left[ {\tan \left( \varphi \right)} \right]}_{{\rm{sa}}}} \cdot \left[ {{q_{n1}} + {q_{n2}}} \right] \cdot dL} \hfill \cr } Therefore, on comparing the two equations (5a) and (5b): (6) [tan(φ)]sa=qn1⋅tan(φ1)+qn2⋅tan(φ2)qn1+qn2=w1⋅tan(φ1)+w2⋅tan(φ2) \eqalign{ & {\left[ {\tan (\varphi )} \right]_{{\rm{sa}}}} = {{{q_{n1}} \cdot \tan ({\varphi _1}) + {q_{n2}} \cdot \tan ({\varphi _2})} \over {{q_{n1}} + {q_{n2}}}} = \cr & {w_1} \cdot \tan ({\varphi _1}) + {w_2} \cdot \tan ({\varphi _2}) \cr} where the dimensionless weight coefficients w_i > 0 equal w₁ = q_n1/(q_n1 + q_n2), w₂ = q_n2/(q_n1 + q_n2), w₁ + w₂ = 1.

In particular, the local contribution of tan(φ₁) is negligible if q_n1 ≪ q_n2, w₁ ≪ w₂, so the randomness is localised, almost not averaged. If a (potentially possible) physical dependence between q_n and tan(φ) is ignored, then E{[tan(φ)]_sa} = w₁×μ + w₂×μ = μ and Var{[tan(φ)]_sa} = [1 − 2×w₁×w₂×(1 − ρ_1;2)]×σ².

Therefore, the variance reduction factor equals γ_sa = Var{[tan(φ)]_sa}/σ² = 1 − 2×w₁×(1 − w₁)×(1 − ρ_1;2); the lower index s refers to the stress-weighted averaging. The variance reduction factor γ_sa defined above is sensitive to q_n gradients and depends on the ratio q_n1/q_n2, what is not true for the variance reduction factor γ_ga = 1 −(1 − ρ_1;2)/2, which is stress-independent. The following inequality γ_ga ≤ γ_sa is true, but γ_ga = γ_sa if and only if w₁ = w₂ = 1/2, like for five dots depicted in Fig. 1b. To conclude, the geometric averaging has only an intuitive character, weakly related to the shearing resistance of frictional materials; the overestimation of the beneficial variance reduction γ_ga ≤ γ_sa can be dangerous in safety evaluation. Parabolic functions γ_sa(w₁) are depicted in Fig. 1b for five different correlation coefficients ρ = ρ_1;2.

Passage to a discrete multiterm sum for n > 2 random variables is straightforward; the same applies for a second-order stochastic process tan(φ(x)) with an autocorrelation function ρ(x;y) along the interval AB, where E{tan(φ(x))} = μ = const, Var{tan(φ(x)} = σ² = const. The shearing resistance T is a random variable generated by the stationary stochastic process tan(φ(x)) in the formula: (7a) T=∫ABqn(x)⋅tan(φ(x))dx T = \int_A^B {{q_n}(x) \cdot \tan \left( {\varphi (x)} \right)dx} Looking for an equivalent simplified model, the stochastic process tan(φ(x)) is replaced by a random variable [tan(φ)]_sa = const(x). If so, the same random variable T has the following form: (7b) T=∫AB[tan(φ)]sa⋅qn(x)dx=[tan(φ)]sa⋅∫ABqn(x)dx \matrix{ {T = \int_A^B {{{\left[ {\tan (\varphi )} \right]}_{{\rm{sa}}}} \cdot {q_n}(x)dx = } } \cr {{{\left[ {\tan (\varphi )} \right]}_{{\rm{sa}}}} \cdot \int_A^B {{q_n}(x)dx} } \cr } Both alternative expressions for the same random variable T formulate the background for the definition of [tan(φ)]_sa; if equations (7a) and (7b) are compared, the following definition can be obtained: (8) [tan(φ)]sa=∫ABqn(x)⋅tan(φ(x))dx∫ABqn(x)dx {\left[ {\tan (\varphi )} \right]_{{\rm{sa}}}} = {{\int_A^B {{q_n}(x) \cdot \tan (\varphi (x))dx} } \over {\int_A^B {{q_n}(x)dx} }} Therefore, (9a) E{[tan(φ)]sa}=∫ABqn(x)⋅E{tan(φ(x))}dx∫ABqn(x)dx=∫ABqn(x)⋅μdx∫ABqn(x)dx=μ \matrix{ {E\left\{ {{{\left[ {\tan (\varphi )} \right]}_{{\rm{sa}}}}} \right\} = {{\int_A^B {{q_n}(x) \cdot E\left\{ {\tan (\varphi (x))} \right\}dx} } \over {\int_A^B {{q_n}(x)dx} }} = } \cr {{{\int_A^B {{q_n}(x) \cdot \mu dx} } \over {\int_A^B {{q_n}(x)dx} }} = \mu } \cr } (9b) Var{[tan(φ)]sa}=∫AB∫ABqn(x)⋅qn(y)⋅Cov{tan(φ(x));tan(φ(y))}dxdy[∫ABqn(x)dx]2 \matrix{ {Var\left\{ {{{\left[ {\tan (\varphi )} \right]}_{{\rm{sa}}}}} \right\} = } \cr {{{\int_A^B {\int_A^B {{q_n}(x) \cdot {q_n}(y) \cdot {Cov}\left\{ {\tan (\varphi (x));\tan (\varphi (y))} \right\}dxdy} } } \over {{{\left[ {\int_A^B {{q_n}(x)dx} } \right]}^2}}}} \cr } The dimensionless variance reduction factor equals (10) γsa=Var{[tan(φ)]sa}Var{tan(φ(x))}=Var{[tan(φ)]sa}σ2==∫AB∫ABqn(x)⋅qn(y)⋅ρ(x;y)dxdy[∫ABqn(x)dx]2≤1 \matrix{ {{\gamma _{{\rm{sa}}}} = {{{Var}\left\{ {{{\left[ {\tan (\varphi )} \right]}_{{\rm{sa}}}}} \right\}} \over {{Var}\left\{ {\tan (\varphi \left( x \right))} \right\}}} = {{{Var}\left\{ {{{\left[ {\tan (\varphi )} \right]}_{{\rm{sa}}}}} \right\}} \over {{\sigma ^2}}} = } \cr {{ = }{{\int_{A}^{B} {\int_{A}^{B} {{q_n}(x) \cdot {q_n}(y) \cdot \rho (x;y)dxdy} } } \over {{{\left[ {\int_A^B {{q_n}(x)dx} } \right]}^2}}} \le 1} \cr } where ρ(x;y) = ρ(x − y) because of the assumed homogeneity; the expression (10) has only loose ties with the approximate variance reduction factor proposed for practical applications by Vanmarcke (2010): Γ² = min{1; δ/|AB|}.

Only if q_n = const along AB, then the geometric averaging coincides with stress-weighted averaging, but generally the latter one is more representative.

Simple numerical results are presented in Fig. 2; a fixed interval AB of length|AB|=7 m is loaded in a localised way on [0;7]: q_n(x) = 0 on [0;S], 0 ≤ S < 7, and q_n(x) = q = const > 0 on (S;7], so in the Heaviside form, q_n(x) = q×[H(x − S) − H(x − 7)]; the Gaussian autocorrelation ρ(x;y) = exp{−[(x−y)/d]²}, 0 ≤ x, y ≤ 7 is assumed. Both methods coincide if and only if S = 0 (four dots on the vertical axis), but γ_ga < γ_sa elsewhere; the geometric spatial averaging is insensitive to the length S, which describes very different stresses along AB.

Figure 2

As a practical example, consider the standard GEO-sliding stability condition (EN 1997-1) for a massive shallow foundation BxLxh, where −B/2 ≤ x₁ ≤ +B/2 and −L/2 ≤ x₂ ≤ +L/2. Define normal (vertical) contact stresses under the footing as q_n(x₁,x₂) = q_n(x) ≥ 0. The foundation base BxL is horizontal and ‘rough’, that is, the stabilising horizontal force can be calculated from the simple frictional model as T = V×tan(φ), if the subsoil is cohesionless; note that this deterministic model does not depend on the eccentricities of the deterministic vertical load V called as e_B and e_L, respectively.

In the random variable format, the friction coefficient tan(φ) in the expression T = V×tan(φ) has a mean value μ and a variance σ², so E{T} = Vμ, Var{T} = V₂σ². For a second-order homogeneous isotropic random field tan(φ(x₁,x₂)) = tan(φ(x)) and the random shearing resistance, T can be expressed following equations (7a) and (7b), that is: (11) T=∫BxLqn(x_)⋅tan(φ(x_))dx_ and,on the other hand, T=[tan(φ)]sa⋅∫BxLqn(x_)dx_ \matrix{ {T = \int_{BxL} {{q_n}(\underline x) \cdot \tan \left( {\varphi (\underline x)} \right)d{\underline x}\,\,{\rm{and}},} } \cr {{\rm{on}}\,{\rm{the}}\,{\rm{other}}\,{\rm{hand}},\,\,T = {{\left[ {\tan (\varphi )} \right]}_{{\rm{sa}}}} \cdot \int_{BxL} {{q_n}(\underline x)d{\underline x}} } \cr } Therefore, (12) [tan(φ)]sa=∫BxLqn(x_)⋅tan(φ(x_))dx_∫BxLqn(x_)dx_and eventually,γsa=∫BxL∫BxLqn(x_)⋅qn(y_)⋅ρ(x_;y_)⋅dx_dy_[∫BxLqn(x_)dx_]2 \matrix{ {{{\left[ {\tan (\varphi )} \right]}_{{\rm{sa}}}} = {{\int_{BxL} {{q_n}({\underline x} ) \cdot \tan \left( {\varphi ({\underline x} )} \right)d{\underline x} } } \over {\int_{BxL} {{q_n}({\underline x} )d{\underline x} } }}{\rm{and}}\,{\rm{eventually,}}} \cr {{\gamma _{{\rm{sa}}}} = {{\int_{BxL} {\int_{BxL} {{q_n}({\underline x} ) \cdot {q_n}({\underline y}) \cdot \rho ({\underline x} ;{\underline y}) \cdot d{\underline x} d{\underline y}} } } \over {{{\left[ {\int_{BxL} {{q_n}({\underline x} )d{\underline x} } } \right]}^2}}}} \cr } Note that the vertical coordinate x₃ (or simply z) is not used in this model.

The presented numerical example employs the horizontally isotropic Gaussian autocorrelation function ρ(x;y) = exp{−[(x₁ − y₁)/d_h]² − [(x₂ − y₂)/d_h]²}; following the recommendations of the Eurocode 7 (EN 1997-1), the vertical contact stress q_n under the footing can be assumed as a double linear function like q_n(x) = a_o + a₁x₁ + a₂x₂ (kPa). Consider a symmetrical case B = L = 10 m, V =10 MN, q = V/(B×L) = 100 kPa, e_B = e_L = e (m) ≥ 0.

If both eccentricities 0 ≤ e < 10/12 m, then under the corners: q₁ = 100 + 120×e, q₂ = q₃ = 100, q₄ = 100 − 120×e > 0, so q_n(x) = 100 + 12×e×(x₁ + x₂) (kPa).

If both eccentricities e > 10/12, then a gap appears under foundation (zero contact stress, locally); for a maximal value e = 30/12 m, which is the maximal acceptable by the Eurocode 7, there is q_n(x) = max{0;60×(x₁ + x₂)} (kPa).

The values of the variance reduction factor γ_sa are depicted in Fig. 3 for e ≤ 30/12 m.

Figure 3

The geometric spatial averaging and the variance reduction factor γ_ga do not depend on the eccentricities e and they overestimate the variance reduction effects (four dots in Fig. 3). Special attention should be paid to the curves for d ∼ 20 m and d ∼ 10 m, which correspond to the most realistic values of the horizontal correlation parameter d_h: –

for d_h = 10 m, the variance reduction γ_sa is on the level of 74%–87%, whereas γ_ga = 74%;

The integral equations (7)–(10) apply not only to shearing resistance, but also to settlement analysis and to any linear operator, in fact. In a simplified linear form suggested by the Eurocode 7 (EN 1997-1), the foundation settlement w (m) can be estimated as follows: (13a) w=∫0zmaxqn(z)⋅D(z)dz w = \int_0^{{z_{\max }}} {{q_n}(z) \cdot D(z)dz} The symbol D(z) > 0 (1/kPa) denotes a random deformation modulus – stationary stochastic process as a function of the depth x₃ = z > 0 under foundation.

The following equivalent [D]_sa = const(z) as a corresponding random variable can be introduced: (13b) w=∫0zmax[D]sa⋅qn(z)dz=[D]sa⋅∫0zmaxqn(z)dz \matrix{ {w = \int_0^{{z_{\max }}} {{{\left[ D \right]}_{{\rm{sa}}}} \cdot {q_n}(z)dz} = } \cr {{{\left[ D \right]}_{{\rm{sa}}}} \cdot \int_0^{{z_{\max }}} {{q_n}(z)dz} } \cr } Equations (13a) and (13b) define the random variable [D]_sa.

Equations (13a) and (13b) are focused on a finite interval 0 ≤ z ≤ z_max under a shallow footing; taking into account a founding depth h_f under the ground surface, the value of z_max can be estimated from the equation γ×(h_f/B + z_max/B)×B×0.2 = q_n(z_max) (EN 1997-1), where g×(h_f + z) denotes vertical stresses in situ for a dry soil and γ (kN/m³) is the subsoil unit weight. The average vertical stress q_n under a rectangular footing BxL vanishes with depth z; usually, the integrated Steinbrenner formula is in use (a double integral of the Boussinesq solution) or just the following simplified estimation: q_n(z) = V/[(B + z)×(L + z)] = q/[(1 + z/B) × (1 + β×z/B)] for β = B/L ≤ 1, q = V/(B×L) = const > 0. The deterministic depth z_max depends on the β coefficient as in Fig. 4. It is assumed that random fluctuations of D(z) do not change the deterministic load function q_n(z) and vice versa.

Figure 4

The stochastic parameter D(z) (1/kPa) needs a comment. Clearly, the random deformation modulus D(z) corresponds to 1/E(z) or 1/M(z) by making use of the Young modulus E or the oedometric modulus M of the subsoil; both are assumed as stochastic processes. As far as distribution-free models are considered, the first two moments E{D(z)} = μ = const(z), Var{D(z)} = σ² = const(z), Cov{D(z₁);D(z₂)} = σ²×ρ(z₁ − z₂) cannot be derived from the corresponding moments E{M(z)}, Cov{M(z₁);M(z₂)} – the model should be completed by some joint probability distributions of the stochastic process M(z); in particular, correlated lognormal random variables M(z₁), M(z₂) can be useful. Alternatively, the first-order perturbation is acceptable for a ‘small’ variance of M(z); also, an adequate data pre-processing can be recommended to apply the sample statistics directly to 1/M, not to M, because the stationarity of both 1/M and M has different conditions.

Equations (13a) and (13b) yield the variance of the averaged deformation modulus Var{[D]_sa} and then the variance reduction factor equals: (14) γsa=Var{[D]sa}Var{D(z)}=∫0zmax∫0zmaxqn(z1)⋅qn(z2)⋅ρ(z1;z2)dz1dz2[∫0zmaxqn(z)dz]2≤1 {\gamma _{{\rm{sa}}}} = {{{Var}\left\{ {{{\left[ D \right]}_{{\rm{sa}}}}} \right\}} \over {{Var}\left\{ {D(z)} \right\}}} = {{\int_0^{{z_{\max }}} {\int_0^{{z_{\max }}} {{q_n}({z_1}) \cdot {q_n}({z_2}) \cdot \rho ({z_1};{z_2})d{z_1}d{z_2}} } } \over {{{\left[ {\int_0^{{z_{\max }}} {{q_n}(z)dz} } \right]}^2}}} \le 1 Spatial horizontal random fluctuations are much less significant for such rectangular footing BxL; they are not considered here (δ_h ≫ d_v).

The numerical example in Fig. 4 uses the Gaussian autocorrelation function ρ(z₁;z₂) = exp{−[(z₁ − z₂)/d_v]²}, 0 ≤ z₁, z₂ ≤ z_max; the subsoil unit weight γ = 18 kN/m³, q = 250 kPa, h_f/B = 2 and the foundation width B = 1 m. For such fixed data, the effective dimensionless z_max/B depends on the dimensionless shape coefficient β = B/L ≤ 1 as in Fig. 4; the shape ratio β = 0 corresponds to an infinite beam and β = 1 corresponds to a square footing.

The geometric variance reduction factor γ_ga in Fig. 5 is also derived from equation (14) for the same z_max, but for q_n = const(z). Differences between the factors γ_sa and γ_ga in Fig. 5 are significant.

Figure 5

For perfectly cohesive materials, the shearing resistance in Fig. 1a equals 2dT = c₁×dL + c₂×dL and is replaced by 2dT = [c]_a×dL + [c]_a×dL; hence, [c]_a = (c₁ + c₂)/2. For such a case, the geometric spatial averaging (Vanmarcke 2010) is representative because the model is independent of the normal stress q_n > 0. There is only one averaged cohesion: (15) [c]a=[c]sa=[c]ga=∬Dc(x_)dx_∬D1dx_ {\left[ c \right]_a} = {\left[ c \right]_{{\rm{sa}}}} = {\left[ c \right]_{{\rm{ga}}}} = {{\int\!\!\!\int_D {c({\underline x} )d{\underline x} } } \over {\int\!\!\!\int_D {1d{\underline x} } }} and only one dimensionless reduction coefficient γ_a = γ_sa = γ_ga.

To be precise, there exist some rare exceptions because this is not the soil cohesion (or adhesion) itself which generates the shearing resistance, but the soil cohesion (or adhesion) on a real contact with the foundation. As in section 6, a gap can appear under a rigid foundation for ‘great’ eccentricities e, hence the averaging is limited (localised) only to a soil-foundation direct contact (q_n > 0); to a certain extent, for a rigid foundation, the model is not completely independent of q_n, so γ_ga can be less than γ_sa also for perfectly cohesive soils. This aspect is not analysed below.

Although the Tresca model with tan(φ) ≡ 0, c > 0 found certain practical applications to undrained bearing capacity, this is, however, a peculiar case. Generally, the shearing ultimate limit states are stress dependent, not only for the Coulomb material. In a general case of τ(x) = q_n(x)×tan(φ(x)) + c(x), two homogeneous random fields are considered for a frictional–cohesive soil. Both autocovariance functions are completed by a crosscovariance function Cov{tan(φ(x));c(y)} = σ_tan(φ)×σ_c×r_tan(φ);c×ρ_tan(φ);c(x − y). More or less significant distinctions between the decay functions ρ_tan(φ)(x − y), ρ_c(x − y) and ρ_tan(φ);c(x − y) are possible. If q_n(x) > 0, then the usual definition of a local shear coefficient tan(ψ(x)) = tan(φ(x)) + c(x)/q_n(x) > tan(φ(x)) reduces the situation to an equivalent cohesionless material. The new random field tan(ψ(x)) is not homogeneous if q_n(x) ≠ const(x); in particular, E{tan(ψ(x))} = μ_tan(φ) + μ_c/q_n(x), Var{tan(ψ(x))} ≠ const(x). By making use of the procedure defined by equations (7) and (8) for the frictional soil, (16a) T(tan(φ(x_));c(x_))=∬Dqn(x_)⋅tan(ψ(x_))dx_ T\left( {\tan \left( {\varphi ({\underline x} )} \right);c({\underline x} )} \right) = \int\!\!\!\int_D {{q_n}({\underline x} ) \cdot \tan \left( {\psi ({\underline x} )} \right)d{\underline x} } (16b) T([tan(ψ)]sa)=∬D[tan(ψ)]sa⋅qn(x_)dx_=[tan(ψ)]sa⋅∬Dqn(x_)dx_ \matrix{ {T\left( {{{\left[ {\tan (\psi )} \right]}_{{\rm{sa}}}}} \right) = \int\!\!\!\int_D {{{\left[ {\tan (\psi )} \right]}_{{\rm{sa}}}} \cdot {q_n}({\underline x} )d{\underline x} = } } \cr {{{\left[ {\tan (\psi )} \right]}_{{\rm{sa}}}} \cdot \int\!\!\!\int_D {{q_n}({\underline x} )d{\underline x} } } \cr } Therefore, (17) [tan(ψ)]sa=∬Dqn(x_)⋅tan(φ(x_))dx_+∬Dc(x_)dx_∬Dqn(x_)dx_=[tan(φ)]sa+[c]a⋅∬D1dx_∬Dqn(x_)dx_ \matrix{ {{{\left[ {\tan (\psi )} \right]}_{{\rm{sa}}}} = {{\int\!\!\!\int_D {{q_n}({\underline x} ) \cdot \tan \left( {\varphi ({\underline x} )} \right)d{\underline x} + \int\!\!\!\int_D c ({\underline x} )d{\underline x} } } \over {\int\!\!\!\int_D {{q_n}({\underline x} )d{\underline x} } }} = } \cr {{{\left[ {\tan (\varphi )} \right]}_{{\rm{sa}}}} + {{\left[ c \right]}_a} \cdot {{\int\!\!\!\int_D {1d{\underline x} } } \over {\int\!\!\!\int_D {{q_n}({\underline x} )d{\underline x} } }}} \cr } Here, the variance reduction factor (equation 10) is less useful because it is a position-dependent function (in the denominator): (18) γsa=Var{[tan(ψ)]sa}Var{tan(ψ(x_))}′x_∈D {\gamma _{{\rm{sa}}}} = {{{Var}\left\{ {{{\left[ {\tan (\psi )} \right]}_{{\rm{sa}}}}} \right\}} \over {{Var}{{\left\{ {\tan (\psi ({\underline x} ))} \right\}}^\prime }}}{\underline x} \in D

Alternatively, the vectorial random field (tan(φ(x));c(x)) can be replaced by the equivalent random vector ([tan(φ)]_sa;[c]_a) by making use of the previously defined averages. Depending on a stress level of q_n(x) – not only on gradients – either the value [tan(φ)]_sa or [c]_a dominates in the averaging; thus, the interpretation becomes more complex. To overcome this difficulty, the effect of the variance reduction can be studied in terms of the random shearing resistance T, not in terms of its strength parameters. Therefore, (19) γsa=Var{T([tan(φ)]sa;[c]a)}Var{T(tan(φ);c)}=Var{[tan(φ)]sa⋅∬Dqn(x_)dx_+[c]a⋅∬D1dx_}Var{tan(φ)⋅∬Dqn(x_)dx_+c⋅∬D1dx_} \matrix{ {{\gamma _{{\rm{sa}}}} = {{{Var}\left\{ {T\left( {{{\left[ {\tan (\varphi )} \right]}_{sa}};{{\left[ c \right]}_a}} \right)} \right\}} \over {{Var}\left\{ {T\left( {\tan (\varphi );c} \right)} \right\}}} = } \cr {{{{Var}\left\{ {{{\left[ {\tan (\varphi )} \right]}_{{\rm{sa}}}} \cdot \int\!\!\!\int_D {{q_n}({\underline x} )d{\underline x} + {{\left[ c \right]}_a} \cdot \int\!\!\!\int_D {1d{\underline x} } } } \right\}} \over {{Var}\left\{ {\tan (\varphi ) \cdot \int\!\!\!\int_D {{q_n}({\underline x} )d{\underline x} + c \cdot \int\!\!\!\int_D {1d{\underline x} } } } \right\}}}} \cr } The random solution T(tan(φ);c) in the denominator operates simply on the (correlated) random variables tan(φ), c, yielding from the random fields tan(φ(x)), c(x) when ignoring their spatial effects, that is, for ρ(x;y) ≡ 1.