Studia Geotechnica et Mechanica, 2025; 47(1): 17–35

Original Study Open Access

Hubert Szabowicz*, Marek Kawa, Wojciech Puła

Efficient and conservative estimation reliability

analysis of strip footing on spatially variable c - ϕ
soil using random finite element limit analysis
https://doi.org/10.2478/sgem-2025-0002
received July 24, 2024; accepted December 12, 2024.
1 Introduction
Abstract: The study deals with reliability analysis of strip Geotechnical structures, particularly the foundations,
foundation on spatially variable c - ϕ soil. The spatial are an important part of the civil engineering discipline.
variability of soil strength parameters, namely cohesion Embedded in soil or built at the surface level, they form
c and friction angle ϕ is modelled using anisotropic the base for the rest of the structure, which means that
uncorrelated random fields, generated with the Fourier their potential failure usually leads to the failure of other
series method. Random finite element limit analysis structural elements (Kawa, 2023; Kawa & Puła, 2020).
(RFELA) providing a rigorous lower and upper bound for Meanwhile, the soil surrounding the foundation is often
bearing capacity for individual Monte-Carlo simulations is highly spatially variable (Ching et al., 2018). Due to
employed. Additional use of adaptive meshing refinement limited budgets, its parameters can only be investigated at
algorithm leads to a significant reduction of the relative selected locations (Jerez et al., 2024). As a result, the design
difference between statistical moments of obtained lower of geotechnical structures is often associated with greater
and upper bound results. The influence of the horizontal uncertainty and their failure with greater consequences
and vertical scales of fluctuation and foundation depths than other branches of civil engineering. It is therefore
on the mean and standard deviation of the obtained bound not surprising that the field of structural reliability, which
moments is investigated. Additionally, the rigorousness has developed strongly in recent decades, is particularly
of the mean and standard deviation of both considered concerned with geotechnical structures, indicating that
bounds estimation is checked. As a result of the analysis, the principal source of uncertainty for them is the spatial
a novel approach based on a mixed distribution that variability of the soil parameters.
combines lower and upper bound moments is introduced. The reliability of the foundations can be analysed.
As shown, this approach offers significant benefits by In such analysis, it is possible to take into account the
providing conservative and relatively precise measures spatial variability of soil parameters, e.g., using stationary
of reliability which can be obtained in reasonable random fields (J. L. Doob, 1990), which according to many
computation time. The proposed method seems to be works (e.g., Hicks & Samy, 2002; Griffiths et al., 2009;
adequate for practical engineering reliability analysis Pieczyńska-Kozłowska et al., 2015) should be considered
of foundation bearing capacity and other limits states as anisotropic. Identification of such field consists of
problems. two parts: identification of probability distribution of
parameters and identification of so-called autocorrelation
Keywords: Strip Fundation; Random field; Random function or more precisely its parameters, scales of
Finite Element Limit Analysis; Reliability Analysis. fluctuation (SOFs), separately in the horizontal and
vertical direction. The latter quantities can be simply
understood as average dimensions of the random cluster
in the field within which the values of the given parameter
*Corresponding author: Hubert Szabowicz, Faculty of Civil Engineering, are strongly correlated. Application of random fields
Wroclaw University of Science and Technology, Wyb. Wyspianskiego 27, for the description of soil parameters is indicated as an
50-370 Wroclaw, Poland, E-mail: [email protected] appropriate approach by, i.e. annexe D of the newest
Marek Kawa, Wojciech Puła, Faculty of Civil Engineering, Wroclaw edition of the ISO 2394:2015 code dedicated to reliability
University of Science and Technology, Wyb. Wyspianskiego 27,
analysis.
50-370 Wroclaw, Poland

Open Access. © 2025 Hubert Szabowicz, Marek Kawa, Wojciech Puła, published by Sciendo. This work is licensed under the
Creative Commons Attribution alone 4.0 License.
18 Hubert Szabowicz, Marek Kawa, Wojciech Puła

A popular and universal method that employs random iterations, can quite precisely estimate the exact solution
fields to describe crucial soil parameters for numerical (their relative difference can be less than 2%). Adaptive
calculations is the random finite element method (RFEM) and nonadaptive algorithms of RFELA have been used
(Griffiths & Fenton, 1993; Fenton & Griffiths, 2008). It mostly for probabilistic assessment of foundations
consists of sequential mapping of random field values bearing capacity (e.g., Simões et al., 2014; Ali et al., 2016)
for many realisations to the part of the finite element and factors of stability for slopes until now (Ali et al.,
(FE) mesh, representing the soil medium in a numerical 2017), however, the list of problems analysed is constantly
problem model. The probabilistic results are based on expanding. Recently, the method (implemented in
a statistical analysis of considered structure responses commercial code OPTUM 2G) was, e.g., used to estimate
taken from all the mapped realisations, each computed the average diameter of randomly irregular lunar lava
as an individual Monte-Carlo simulation (MCS). The tubes (Chwała et al., 2024). Although the cited and similar
method enables probabilistic modelling of any virtually works take advantage of the high speed of the RFELA, in
geotechnical structure which can be modelled with FE by the authors’ opinion, the rigorousness of the LBs and UBs
either linear or nonlinear models. Until now, the method so far has not been fully utilised.
was used to analyse different types of foundations (e.g., This paper attempts to exploit both advantages of
Fenton & Griffiths, 2008; Kawa & Puła, 2020; Dobrzanski the RFELA method, i.e., high speed and rigorousness
& Kawa, 2021; Chen et al., 2022; Kumar et al., 2023), but of bounds estimation. Both approaches can be utilised
also other geotechnical structures such as slopes (e.g., to estimate the mean and standard deviation of the UB
Griffiths & Fenton, 2004; Huang et al., 2021), diaphragm and LB, respectively. The efficiency of the algorithm
walls (Sert et al., 2016; Kawa et al., 2021), piles (Cheng et enabled quick and accurate estimation of these moments
al., 2023), etc. The greatest disadvantage of the RFEM lies using a limited number of adaptive mesh refinements.
in its significantly high computation cost being a direct The estimated moments were used to define the mixed
result of the computation cost of a single MCS. Given the distribution, i.e., one obtained by assuming mean based
typical number of N=1000 MCSs (which in the case of the on a static approach and standard deviation based
ultimate limit state allows only for a rough estimation of on a kinematic approach. As shown, the results of the
failure probability, cf. Kawa et al., 2021) and the typical novel approach can be used for quick, conservative and
range of time needed to solve a problem realisation relatively precise reliability analysis.
(dependent on the complexity of the problem), the total The subject of this study is an analysis of the bearing
computation time on a modern PC is usually between few capacity of strip foundations on c - φ soil. Both friction
hours and several days. Such a long time significantly and cohesion are modelled by anisotropic random fields
limits the possibility of using the method in engineering generated with the Fourier series method (FSM, Jha &
practice. For this reason, more efficient alternatives, which Ching, 2013). Three foundation levels and several different
are based on the quasi-analytical kinematic method (e.g., values of the horizontal and vertical scales of fluctuation
Chwała, 2019) or limit equilibrium method (e.g., Liu et al., were analysed. The LBs and UBs of bearing capacity
2019) are still being sought. in individual realisations are determined using FELA
One of the alternatives to RFEM is random finite formulations proposed by Lyamin and Sloan (2002a) and
element limit analysis (RFELA), a more efficient Makrodimopoulos and Martin (2008), respectively. These
approach, which is also based on mapping random field formulations were implemented by the authors in Matlab
realisations into a discretised domain. It is a modification code. In each realisation, the solutions for the LBs and UBs
of RFEM, which uses finite element limit analysis (FELA, were obtained using the commercial solver MOSEK. The
Lyamin & Sloan, 2002a, 2002b) instead of the typical FE influence of both horizontal and vertical fluctuation scales
method (FEM) to solve individual simulations. The class and foundation levels on the probabilistic characteristics
of the problem is limited here to limit states. For such of bearing capacity assessments was investigated. As
problems, FELA enables the rigorous determination shown, the applied numerical algorithm allowed for quite
of both the lower (based on the static approach) and precise and rigorous estimation of mean and variance in
the upper (based on the kinematic approach) bounds both approaches. Additional use of the mentioned mixed
(LBs, UBs) of the solution (e.g., bearing capacity value) distribution enabled conservative and precise estimation
by solving constrained convex optimisation problems, of allowable load value satisfying given failure probability
which is a much faster approach than traditional elastic- condition.
perfectly plastic FEM analyses. With the additional use of
adaptive meshing refinement, these bounds, in successive
 Efficient and conservative estimation reliability analysis of strip footing on spatially variable c - ϕ soil ... 19

2 Methodology 𝜎𝜎𝜎𝜎𝑔𝑔𝑔𝑔2 = ln �1 +
𝜎𝜎𝜎𝜎 2 1
�, 𝜇𝜇𝜇𝜇𝑔𝑔𝑔𝑔 = ln(𝜇𝜇𝜇𝜇) − 𝜎𝜎𝜎𝜎𝑔𝑔𝑔𝑔2 (2)
𝜇𝜇𝜇𝜇2 2
2.1 Random field generation
After generation, the values of the underlying normal
The assumed soil is of c - φ type. Both soil strength
field were further transformed to lognormal field RFX(X)
parameters were described with random fields. It is
according to formula (2).
assumed that these fields are stationary and anisotropic.
Since c and φ are usually reported as negatively correlated
RFX(X)=exp(RFY(X)), (3)
and such correlation typically results in a lower standard
deviation of bearing capacity (cf. Cho & Park 2010) and
where X={x,z} denotes the coordinate of the point that
consequently lower probability of failure, the fields
belongs to the random field.
were also assumed uncorrelated for conservatism. Both
Both random fields were assumed to be characterised
probability distributions were assumed to be lognormal.
by a Gaussian autocorrelation function according to the
To generate the random fields, the FSM was utilised. This
following formula:
approach approximates the autocorrelation function
by applying a finite Fourier series. The value of the field 𝜏𝜏𝜏𝜏𝑥𝑥𝑥𝑥 2 𝜏𝜏𝜏𝜏𝑧𝑧𝑧𝑧 2
𝜌𝜌𝜌𝜌(𝜏𝜏𝜏𝜏𝑥𝑥𝑥𝑥 , 𝜏𝜏𝜏𝜏𝑧𝑧𝑧𝑧 ) = exp �−𝜋𝜋𝜋𝜋 𝜋𝜋 � + � � �� (4)
at a given point in 2D space with coordinates x and z is 𝜃𝜃𝜃𝜃𝑥𝑥𝑥𝑥 𝜃𝜃𝜃𝜃𝑧𝑧𝑧𝑧
calculated as the following double sum:

where τx, τz are, respectively, the distance between two

𝑚𝑚𝑚𝑚 𝑚𝑚𝑚𝑚
𝑚𝑚𝑚𝑚 ∙ 𝑥𝑥𝑥𝑥 𝑛𝑛𝑛𝑛 ∙ 𝑧𝑧𝑧𝑧 points on the horizontal (x) and vertical (z) axes, and θx,
𝑊𝑊𝑊𝑊(𝐗𝐗𝐗𝐗) = 𝑊𝑊𝑊𝑊(𝑥𝑥𝑥𝑥, 𝑧𝑧𝑧𝑧) = � �(𝑎𝑎𝑎𝑎𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚 + 𝑖𝑖𝑖𝑖𝑏𝑏𝑏𝑏𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚 ) exp �𝑖𝑖𝑖𝑖2𝜋𝜋𝜋𝜋 � + �� (1)
−𝑚𝑚𝑚𝑚 −𝑚𝑚𝑚𝑚
𝐿𝐿𝐿𝐿𝑥𝑥𝑥𝑥 𝐿𝐿𝐿𝐿𝑧𝑧𝑧𝑧 θz are the horizontal and vertical scales of fluctuation. The
Gaussian autocorrelation function was chosen because of
the low length of the associated Fourier series (compared
where amn, bmn are independent random variables with a
with, e.g., exponential function) which affects the time of
normal distribution N(0,σmn) with σmn dependent on the
random filed values calculation.
chosen autocorrelation function and Lx and Lz are the
dimensions of the domain over which the simulations are
performed. The length of the series is determined by the 2.2 Finite element limit analysis
numbers m and n which depend on the assumed precision
and the assumed autocorrelation function. Please note that As mentioned, the RFELA method utilises FELA to solve
the description of the random field by a continuous function, individual realisations of the given problem. The latter is
that allows the value of the field realisation at any point to one of the possible numerical methods capable of solving
be determined independently of the element mesh, is a the limit state problem, e.g., of the load-bearing capacity
mandatory condition for RFELA analysis which implements of a shallow foundation. The method requires several
the mesh adaptation procedure. The results obtained by the assumptions to be made which significantly simplify the
FELA method are strongly affected by the FE mesh fitting model while omitting the knowledge of the deformation
to the failure mechanism obtained in a given realisation. state. On the other hand, it enables a considerable
The implementation of an iterative procedure that allows reduction in calculation time compared with, e.g., the
the mesh shape to be adapted to the mechanism makes FE or the finite difference methods. Most importantly,
it possible to increase its accuracy without an excessive FELA provides rigorous lower and upper limits of bearing
increase in the number of elements. A description of the capacity.
mesh adaptation algorithm is provided in the later sections. The LB and UB solutions are determined, respectively,
The FSM algorithm generates random fields with by searching for the maximum load corresponding to the
normal distribution. Since the distribution of soil statically admissible stress field and the minimum load
parameters was assumed as lognormal, the fields with the for which the kinematically admissible displacement
parameters for the underlying normal distribution RFY(X) velocity field is developed. Both issues represent a
were first generated. The mean value (μg) and standard mathematical optimisation problem. This article uses
deviation (σg) of the underlying normal distribution can the approaches described by Lyamin and Sloan (2002a)
be calculated using the analytical formulas that link these in the context of LB estimation and Makrodimopoulos
moments to moments of the target lognormal distribution and Martin (2008) regarding UB estimation. A perfectly
{μ,σ}, namely, plastic material model with an associated flow rule is
20 Hubert Szabowicz, Marek Kawa, Wojciech Puła

assumed. The FE mesh consists of triangular elements for software currently available for solving many different
which a linear shape function is used for the stresses and types of optimisation problems (including SOC) (Podlich
a quadratic shape function is used for the displacement et al., 2014). The software uses the primal–dual method,
velocities when estimating the LB and UB, respectively. which allows the automatic solution of a primal and dual
As previously mentioned, the Coulomb–Mohr material problem, providing a full range of information on stresses,
model with strength parameters c and ϕ is assumed. plastic multipliers and displacement velocities for both LB
This model can easily be written in the form of a conic and UB estimates.
quadratic constraint, which enables finding the solution
of both UB and LB using second-order cone programming.
The letter appears to be a very efficient method for solving 2.3 Mesh adaptive procedure
optimisation problems.
The LB and UB problems are weakly dual, which The degree to which the mesh is adapted to the shape
means that their solutions differ by a certain gap due of the failure mechanism has a significant impact on
to discretisation error. On the other hand, strongly dual the results obtained from the optimisation problems.
problems with the same solutions can be formulated to Obviously, increasing the mesh density over the entire
each of the task (Ciria Suárez, 2004). The dual problem domain will greatly reduce the difference between the LB
to the kinematic approach involves searching for an and UB estimates, however, the number of nodes resulting
admissible stress distribution (slightly different from from such mesh is computationally prohibited. For more
the one in static approach), while the dual problem to efficient computations, elements should be densified in a
the static approach involves searching for an admissible specific manner only in the zone of interest.
failure mechanism (slightly different from the one in For mesh generation, a tool called MESH2D (Engwirda,
kinematic approach). In this article, both problems are 2014) for the Matlab environment, which allows the
formulated by searching for the stress distribution (i.e., generation of triangular meshes suitable for use in the
the original static problem and the dual to kinematic FELA method, is utilised. This generator is based on
problem are solved), because, as noted by Krabbenhoft the Delaunay method. It includes options to modify the
et al. (2005), a dual problem to the kinematic approach meshes both by subdividing the indicated elements and
could be easier to solve. requesting a density function in the form of a dimension
Both mentioned problems can be described by the expectations map. The mesh adaptation procedure
following equations provided by Podlich (2018): utilised in the calculations described below assumes
that the yielding elements should be subdivided. These
maximise: elements are indicated based on the plasticity multipliers
α obtained from both the UB and LB analyses.
subject to: (5) In the first step, the same low-density mesh is generated
BT σ=αp+p0 regardless of the realisation of the random field. Smaller
f(σ ) ≤K0 ∀ e=1,2,…,ne
e
elements occur only at the corners of the foundation due
to the use of so-called ‘element fans’, which are necessary
where α is the searched load value multiplier and σ is to reduce the influence of the singularity located at the
the global stress vector different for LB and UB analysis. specified points (Lyamin & Sloan, 2003). Performing a
Additionally, BT contains the coefficients of the static FELA analysis with such a generated mesh is extremely
equilibrium equations significantly different for the LB and time-efficient and gives a rough indication of the potential
UB limit analysis. Vectors p and p0 contain the optimised plastic zones. In the second step, a mesh density function,
and permanent load. In both cases, the Coulomb–Mohr which assumes minimum dimensions for the elements
plasticity condition is defined for all elements in the form that have yielded and maximum dimensions for all the
of SOCP (Second Order Cone Programming). others, is applied. This operation enables the refinement
On the other hand, in both approaches, the mesh of the areas where the plasticity multiplier is present
adaptation algorithm utilises parameters of plastic without an excessive increase in the number of elements.
multipliers, which are obtained by solving duals to In the next steps, the elements with the highest plastic
the mentioned problems. To solve the problems, the multiplier values are selected and subdivided using the
commercial software MOSEK run from the Matlab MESH2D tool. The limit values of the plastic multiplier are
environment was used. According to the knowledge chosen depending on the individual task. The algorithm
of the authors of this study, this is the most effective must try to match the shape of the mesh with the results
 Efficient and conservative estimation reliability analysis of strip footing on spatially variable c - ϕ soil ... 21

10 10

8 8

6 6

6 8 10 12 14 6 8 10 12 14

10 10

8 8

6 6

6 8 10 12 14 6 8 10 12 14

Figure 1: A graphical representation of the mesh adaptive procedure.

of both estimations. As mentioned, often the static and considered. The surface model is a frequently adopted
kinematic mechanisms may differ. The algorithm applied theoretical model. On the other hand, the typical direct
in the current work consists of four steps. A graphical foundation depth oscillates around 1.0 m. A value of half
representation of the discussed algorithm for a zoomed a metre was assumed as an intermediate value to present
area around the foundation is illustrated in Figure 1. a potentially nonlinear influence of embedded depth. It
was assumed that the foundation established at a depth
greater than zero is a rigid body that extended above the
2.4 Numerical model analysed with FELA ground surface (no backfill soil overlying the foundation
was assumed). Please note that the foundation material
The study investigates the effect of the foundation depth was assumed as infinitely strong (conditions of the
and the value of vertical and horizontal scales of fluctuation Mohr–Coulomb criterion were not checked within the
on the bearing capacity of a shallow foundation placed foundation).
on spatially variable soil. The first analysis results for a For each foundation level, a set of five horizontal
slightly different formulation of the FELA algorithm and a fluctuation scales θx={1.0,2.0,4.0,8.0,16.0 m} and a
different mesh adaptive procedure were presented by Puła set of three vertical scales θz={0.25,0.5,1.0 m} were
et al. (2022). The deliberation was limited to the shallow considered. The analyses included pairs formed from
foundation with D=0.0 m only. A full analysis was carried all possible combinations of these values. The values of
out in the present work. cohesion (μc,σc) and internal friction angle (μφ,σφ) were
The foundation was assumed as a rigid body with a described by the random fields realisation generated
load uniformly distributed over the top face. It was allowed for these values as described earlier. The parameters of
to rotate if such movement resulted from an optimal the lognormal distribution assumed for generation for
failure mechanism. Boundary conditions were adopted as both fields were assumed as μc=29.0kPa, σc=7.0 kPa for
indicated in Figure 2. The width and height of the model cohesion and μφ=12.41°, σφ=1.15° for friction angle. The
were assumed to be 20.0 m and 10.0 m, respectively, which soil weight was assumed to be 20 kN/m3. The specified
allowed to ensure that no failure mechanism reached the parameters were adapted after Zaskórski and Puła (2016),
area boundary. The width of the foundation was 1.0 m. who utilised them in their investigations of clay soil
Three foundation depth scenarios with D=0.0 (foundation layers. Probabilistic analysis was carried out using the
placed on the ground surface), 0.5 and 1.0 m were Monte-Carlo method. For each problem (each assumed
22 Hubert Szabowicz, Marek Kawa, Wojciech Puła

Figure 2: Model geometry with boundary condition.

pair of scales of fluctuation), 10,000 MCS was solved. The fixed vertical scale and changing horizontal fluctuation
LB and UB estimates of bearing capacity were obtained scale are presented rather, than otherwise.
for each problem realisation. A summary of the results First, the improvement in convergence of the
obtained for all considered vertical and horizontal scales obtained LB and UB results for bearing capacity with each
of fluctuations and foundation levels is reported in the subsequent adaptive mesh refinement step was analysed.
following section. All results were computed for a fixed vertical fluctuation
scale of θz=1 m with various values of the horizontal
fluctuation scale, as shown in Figures 3, 4 and 5. The
3 Results results for foundation depths of D = 0.0, 0.5 and 1.0 m are
presented in Figures 3, 4, and 5, respectively. In each case,
As a result of the performed analyses, the values of the the data are given in the form of histograms and graphs
LBs and UBs of bearing capacity from 10,000 MCS for each of the estimated probability density functions, which are
of the three considered foundation levels and 15 pairs of assumed to be lognormal. The blue lines represent the
fluctuation scales, were obtained. Only selected results results of the LB limit analysis, while the red lines represent
are presented in the following subsections. It is generally the results of the UB limit analysis. In individual columns,
known that the vertical scale of fluctuation is considerably the outcomes for the different values of the horizontal
simpler to determine due to the greater number of data scales of fluctuation, namely θx={1.0,2.0,4.0,8.0,16.0 m}.
recorded in the vertical direction during, e.g., CPTu The results obtained at each step of the mesh adaptive
tests. Typically, the horizontal scale of fluctuations is procedure are provided in separate rows, with the final
quite challenging to identify because of the considerable results being shown in the last row.
distance between sounding points. (cf. Ching et al., It can be observed that, for each case analysed,
2018). For more distanced CPTu tests, this fact can lead the mesh adaptive procedure fulfils its purpose by
to a situation where the vertical scale of fluctuation is significantly reducing the gap between the distributions of
known and the horizontal scale of fluctuation is unknown the LB and UB estimations. After four steps of adaptation,
and is, e.g., a subject of parameter study. Analysis the distributions agree almost ideally. In addition, the
with such assumption was performed in some of the character of the changes caused by varying the value of
authors’ previous work (Kawa et al. 2019, 2021). A similar fluctuation scales is very similar regardless of the assumed
assumption for presenting the result was assumed here. foundation level, namely, the histograms obtained for
For this reason on the single graph, results obtained for the subsequent horizontal scales of fluctuation oscillate
 Efficient and conservative estimation reliability analysis of strip footing on spatially variable c - ϕ soil ... 23

400 400 400 400 400

200 200 200 200 200

0 0 0 0 0

200 300 400 500 200 300 400 500 200 300 400 500 200 300 400 500 200 300 400 500

400 400 400 400 400

200 200 200 200 200

0 0 0 0 0

200 300 400 500 200 300 400 500 200 300 400 500 200 300 400 500 200 300 400 500

400 400 400 400 400

200 200 200 200 200

0 0 0 0 0

200 300 400 500 200 300 400 500 200 300 400 500 200 300 400 500 200 300 400 500

400 400 400 400 400

200 200 200 200 200

0 0 0 0 0

200 300 400 500 200 300 400 500 200 300 400 500 200 300 400 500 200 300 400 500

Figure 3: Upper (red lines) and lower (blue lines) bound of bearing capacities of a shallow foundation D=0.0 m.

around the same mean values, while simultaneously being In most of the instances, the values of the coefficient of
characterised by greater values of standard deviations. As variation may unexpectedly change their relative order
the results obtained for other analysed vertical scales of during the convergence process.
fluctuation are qualitatively very similar, they were not Figure 9 illustrates the estimation relative difference
presented in this form. between the upper and lower approach in the mean,
In Figures 6–8, the convergence of mean values, standard deviation and coefficient of variation as a
standard deviations and coefficients of variation of function of the step of the mesh adaptive procedure.
LBs μLB, σLB, νLB and LBs μUB, σUB, νUB of bearing capacity The results corresponding to different foundation levels
calculated for obtained series of results are presented as are presented in different columns. The advantage of
a function of the mesh adaptive procedure step. In the applying an algorithm that adapts the mesh topography
individual drawing, the results obtained for different to the failure mechanism is clearly visible in this figure.
values of vertical fluctuation scale namely θz={1.0,0.5,0.25 For the mean value, the differences between the lower and
m} are presented. Similar to the previous illustrations, upper estimates for the initial calculations with sparse
results for different horizontal scales of fluctuation mesh are in the range of 15.0%–23.5% and then, in the
are presented in individual columns. Again, the blue last step of mesh adaption, fall to the range 2.0%–6.1% for
colour indicates the LB and the red colour the UB of the foundation level D=0.0 m, 2.7%–7.2% for D=0.5 m and
bearing capacity. The results corresponding to different 3.2%–8.6% for D=1.0 m. The analogous behaviour can be
foundation levels, D=0.0, 0.5 and 1.0 m, are marked with observed for the relative differences in standard deviation
continuous, dashed and dotted lines, respectively. As estimation. The algorithm reduces the error from 19.7% to
observed, although the value and rate of convergence 37.2% for the initial meshing to, respectively, 3.0%–12.3%
depend on the specific case of foundation depth or for D=0.0 m, 3.2%–7.5% for D=0.5 m and 2.1%–5.0%
fluctuation scale, in each scenario, both moments for D=1.0 m in the final calculation step. Once again,
exhibit lower values for the LB estimation and higher differing rates of convergence for the mean and standard
values for the UB estimation. While this behaviour is deviation lead to unexpected behaviour in the coefficient
intuitive for the mean values, it is less so for the standard of variation. In some cases, the relative difference in
deviation. Additionally, since the rate of convergence can this statistic can be higher in the final step than in the
differ between the mean and standard deviation, this preceding one. Moreover, for greater depths, negative
behaviour does not extend to the coefficient of variation. values of the relative difference may occur indicating that
24 Hubert Szabowicz, Marek Kawa, Wojciech Puła

400 400 400 400 400

200 200 200 200 200

0 0 0 0 0

300 400 500 600 300 400 500 600 300 400 500 600 300 400 500 600 300 400 500 600

400 400 400 400 400

200 200 200 200 200

0 0 0 0 0

300 400 500 600 300 400 500 600 300 400 500 600 300 400 500 600 300 400 500 600

400 400 400 400 400

200 200 200 200 200

0 0 0 0 0

300 400 500 600 300 400 500 600 300 400 500 600 300 400 500 600 300 400 500 600

400 400 400 400 400

200 200 200 200 200

0 0 0 0 0

300 400 500 600 300 400 500 600 300 400 500 600 300 400 500 600 300 400 500 600

Figure 4: Upper (red lines) and lower (blue lines) bound of bearing capacities of a shallow foundation D=0.5 m.

400 400 400 400 400

200 200 200 200 200

0 0 0 0 0

300 400 500 600 700 300 400 500 600 700 300 400 500 600 700 300 400 500 600 700 300 400 500 600 700

400 400 400 400 400

200 200 200 200 200

0 0 0 0 0

300 400 500 600 700 300 400 500 600 700 300 400 500 600 700 300 400 500 600 700 300 400 500 600 700

400 400 400 400 400

200 200 200 200 200

0 0 0 0 0

300 400 500 600 700 300 400 500 600 700 300 400 500 600 700 300 400 500 600 700 300 400 500 600 700

400 400 400 400 400

200 200 200 200 200

0 0 0 0 0

300 400 500 600 700 300 400 500 600 700 300 400 500 600 700 300 400 500 600 700 300 400 500 600 700

Figure 5: Upper (red lines) and lower (blue lines) bound of bearing capacities of a shallow foundation D=1.0 m.

the coefficient of variation is higher for the LB results than again, the LB estimations are indicated by blue and the
for the UB results. UB estimations by red lines. In individual columns, the
The mean values, standard deviations and results for each foundation depth are presented.
coefficients of variation for the last step of the mesh The results presented in Figure 10 confirm the previous
adaptive procedure as a function of the horizontal scale observations. The mean values seem to be practically
of fluctuation have been presented in Figure 10. Once independent of the value of the horizontal fluctuation
 Efficient and conservative estimation reliability analysis of strip footing on spatially variable c - ϕ soil ... 25

Figure 6: Diagrams of mean value, standard deviation and coefficient of variation of bearing capacities of a shallow foundation as a function
of the step of adaptive mesh refinement procedure for vertical scale of fluctuation θz=0.25 m.

Figure 7: Diagrams of mean value, standard deviation and coefficient of variation of bearing capacities of a shallow foundation as a function
of the step of adaptive mesh refinement procedure for vertical scale of fluctuation θz=0.5 m.
26 Hubert Szabowicz, Marek Kawa, Wojciech Puła

Figure 8: Diagrams of mean value, standard deviation and coefficient of variation of bearing capacities of a shallow foundation as a function
of the step of adaptive mesh refinement procedure for vertical scale of fluctuation θz=1.0 m.

Figure 9: Diagrams of the relative difference in mean, standard deviation and coefficient of variation between the upper and lower bound
estimations as a function of following mesh adaptation procedure steps for all considered foundation levels.
 Efficient and conservative estimation reliability analysis of strip footing on spatially variable c - ϕ soil ... 27

Figure 10: Diagrams of mean value, standard deviation and coefficient of variation as a function of the horizontal scales of fluctuations θx for
all considered foundation levels and vertical scales of fluctuations θz.

scale (however, slight increase for greater scale value studies (Pieczyńska-Kozłowska et al., 2015;, Kawa & Puła,
can be observed), although significantly dependent on 2020).
the foundation level. However, the standard deviation is Tables 1, 2 and 3 contain the calculated values of the
strongly related to the value of the horizontal fluctuation mean, standard deviation and coefficient of variation of
scale and the foundation level. This is caused by the the estimated limit-bearing capacities for D=0.0, 0.5 and
magnitude of the failure mechanism that occurs when the 1.0 m, respectively.
yield criterion is reached. The greater the foundation level,
the wider and deeper the failure mechanism will extend,
which means that the soil parameters are averaged over a 3.1 Conservative estimation of the allowable
much wider area. However, the standard deviation for the load value
only isotropic case considered, namely for the horizontal
and vertical scales of fluctuation θx=θz=1.0 m, for each The results described in the previous chapter are based
foundation depth, are very close to each other regardless on 10,000 MCS realisations. By using RFELA for a single
of the foundation depth. As the value of the horizontal fluctuation scale, these amount of calculations were
fluctuation scale increases, the values of the standard carried out in about 8–10 hours. It is worth noting that,
deviations increase with different rates. The greater the although the lower estimation results always estimate
foundation depth, the greater the increase in the standard the load conservatively, even such a large number of
deviation for successive scales of fluctuation values. Still, realisations do not allow a sufficiently accurate estimation
since the increase in the mean value of the bearing capacity of the probability of failure for the ultimate limit state of
for the deeper foundation (regardless of fluctuations scale) the foundation. According to EN-1990, for the standard
is even greater, the resulting coefficients of variations 50-year reference period, the allowable reliability index
for the deeper foundation present smaller values. The should be greater than β=3.8 which corresponds to the
obtained trends are similar to the values obtained in other failure probability lower than pf=7.23∙10-5. According
28 Hubert Szabowicz, Marek Kawa, Wojciech Puła

Table 1: The mean values, standard deviations and coefficients of variation for the lower and upper bound of bearing capacity estimated
based on results from the last (fourth) step of the adaptive mesh procedure – D=0.0 m.

Scales of fluctuations Lower-bound estimation Upper-bound estimation

Mean value Standard Coefficient of Mean value Standard Coefficient of
(kPa) deviation (kPa) variation (%) (kPa) deviation (kPa) variation (%)

θx =1.0 m θz =0.25m 252.11 16.25 6.44 267.60 18.24 6.82

θz =0.5 m 254.02 22.25 8.76 264.53 24.05 9.09

θz =1.0 m 258.28 28.77 11.14 264.19 29.80 11.28

θx =2.0 m θz =0.25 m 253.51 20.82 8.21 268.67 23.19 8.63

θz =0.5 m 255.50 28.19 11.03 265.44 30.29 11.41

θz =1.0 m 259.00 36.39 14.05 264.32 37.49 14.18

θx =4.0 m θz =0.25 m 255.01 24.58 9.64 270.63 27.30 10.09

θz =0.5 m 258.76 33.90 13.10 268.74 36.19 13.47

θz =1.0 m 263.80 43.09 16.34 269.03 44.39 16.50

θx =8.0 m θz =0.25 m 258.78 27.83 10.75 275.21 30.71 11.16

θz =0.5 m 263.71 36.93 14.01 274.16 39.31 14.34

θz =1.0 m 269.69 49.62 18.40 275.12 51.20 18.61

θx =16.0 m θz =0.25 m 260.23 28.48 10.94 277.44 31.81 11.46

θz =0.5 m 254.02 39.49 14.81 278.15 42.09 15.13

θz =1.0 m 274.18 50.94 18.58 279.80 52.57 18.79

Table 2: The mean values, standard deviations and coefficients of variation for the lower and upper bound of bearing capacity estimated
based on results from the last (fourth) step of the adaptive mesh procedure – D=0.5m.

Scales of fluctuations Lower-bound estimation Upper-bound estimation

Mean value Standard Coefficient of Mean value Standard Coefficient of
(kPa) deviation (kPa) variation (%) (kPa) deviation (kPa) variation (%)

θx =1.0 m θz =0.25 m 374.15 15.73 4.20 401.09 16.79 4.19

θz =0.5 m 379.86 21.81 5.74 398.32 22.95 5.76

θz =1.0 m 383.55 28.25 7.36 397.40 29.39 7.39

θx =2.0 m θz =0.25 m 375.96 19.79 5.26 402.26 21.20 5.27

θz =0.5 m 381.39 27.42 7.19 398.58 28.83 7.23

θz =1.0 m 384.83 36.53 9.49 397.09 37.92 9.55

θx =4.0 m θz =0.25 m 376.99 24.40 6.47 403.00 26.23 6.51

θz =0.5 m 384.16 34.17 8.90 400.92 35.96 8.97

θz =1.0 m 387.88 45.13 11.64 399.40 46.74 11.70

θx =8.0 m θz =0.25 m 380.13 28.46 7.49 405.97 30.59 7.53

θz =0.5 m 388.19 39.09 10.07 404.65 41.01 10.14

θz =1.0 m 393.19 54.23 13.79 404.47 56.00 13.85

θx =16.0 m θz =0.25 m 380.82 30.02 7.88 406.25 32.29 7.95

θz =0.5 m 390.91 42.87 10.97 407.12 44.85 11.02

θz =1.0 m 397.74 57.36 14.42 408.76 59.19 14.48

 Efficient and conservative estimation reliability analysis of strip footing on spatially variable c - ϕ soil ... 29

Table 3: The mean values, standard deviations and coefficients of variation for the lower and upper bound of bearing capacity are estimated
based on results from the last (fourth) step of the adaptive mesh procedure – D=1.0 m.

Scales of fluctuations Lower-bound estimation Upper-bound estimation

Mean value Standard Coefficient of Mean value Standard Coefficient of
(kPa) deviation (kPa) variation (%) (kPa) deviation (kPa) variation (%)

θx =1.0 m θz =0.25 m 476.12 16.68 3.50 517.23 17.42 3.37

θz =0.5 m 485.51 22.77 4.69 514.61 23.37 4.54

θz =1.0 m 491.12 30.03 6.11 512.43 30.66 5.98

θx =2.0 m θz =0.25 m 478.43 20.74 4.34 518.28 21.61 4.17

θz =0.5 m 487.44 28.43 5.83 514.32 29.29 5.69

θz =1.0 m 493.39 38.40 7.78 511.90 39.23 7.66

θx =4.0 m θz =0.25 m 478.85 25.47 5.32 518.15 26.62 5.14

θz =0.5 m 490.31 35.75 7.29 516.22 36.59 7.09

θz =1.0 m 496.35 47.37 9.54 513.48 48.38 9.42

θx =8.0 m θz =0.25 m 481.79 30.08 6.24 520.81 31.47 6.04

θz =0.5 m 494.37 41.32 8.36 519.78 42.54 8.18

θz =1.0 m 501.01 56.89 11.35 517.61 58.12 11.23

θx =16.0 m θz =0.25 m 482.99 32.12 6.65 521.34 33.73 6.47

θz =0.5 m 497.53 46.35 9.32 522.34 47.79 9.15

θz =1.0 m 507.52 62.26 12.27 523.68 63.66 12.16

to Kolmogorov’s inequality, 4,194,304 realisations are allowable load estimation. For this reason, adopting a
required to obtain a credible (with 95% confidence and distribution with a lower mean based on the results of
a relative error of less than 10%) Monte-Carlo estimator the LB (static) approach and a higher standard deviation
of such a low probability (Kawa, 2023). Such a number based on the UB (kinematic) approach appears to be a
of simulations is not feasible in a reasonable time. This safe and conservative approximation of the probabilistic
number could be probably significantly reduced by using bearing capacity behaviour. Such a distribution with
one of the variance reduction techniques, e.g., Subset the following moments will be referred to as the mixed
Simulation (Au & Beck, 2001), but using such a method is solution.
out of the scope of the current work.
As mentioned above, as the consequence of the μ=min{μLB,μUB}=μLB
(6)
rigorousness of FELA, the results described above are σ=max{σLB,σUB}=σUB
consistently characterised by a lower mean value and lower
standard deviation for LB bearing capacity estimations Based on the distributions with values of mean and
and greater mean and greater standard deviation for UB standard deviation presented above including upper-,
bearing capacity estimations. While this relationship lower- and mixed-approach distributions, the values of
is less intuitive in the case of standard deviation, it was the reliability-based design loads for the analysed shallow
clearly observed in all considered cases. It is important foundations were estimated. By following the mentioned
to note that the lower mean value of bearing capacity EN-1990 guidelines, pf=7.23∙10-5 (β=3.8) was assumed as an
distribution with fixed standard deviation always lead to allowable probability condition. For a given distribution,
greater probability of ultimate limit state failure and thus the allowable reliability-based design load value can be
conservative estimation of allowable load. On the other obtained as:
hand, an increase in the standard deviation of structural
response results in an increased probability of failure, qd=Φ-1(pf) (7)
which in the consequence also lead to conservative
30 Hubert Szabowicz, Marek Kawa, Wojciech Puła

Figure 11: Diagrams of allowable loads and safety factors as a function of the horizontal scales of fluctuations θx with a constant vertical
scale of fluctuation for all considered foundation levels.

where Φ-1 denotes the inverse cumulative distribution The results from Figure 11 are also summarised in Tables
function of the considered distribution 4, 5 and 6.
Figure 11 illustrates the results in the form of diagrams Please note that the approach utilising the distribution
of the allowable load qd and the safety factor, understood estimated based on the LB usually provides a lower (more
as the ratio of the mean value to the value of the allowable conservative) estimate of the allowable load than the
load. Both variables were plotted as a function of the analogous approach based on the UB. However, since the
horizontal scale of fluctuation. Subsequent vertical scales rate of convergence differs between the LB, UB, mean and
of fluctuation have been marked with solid, dash-dotted standard deviation, it may be possible that the LB results
and dashed lines, respectively. Once again, results for will estimate the mean value of the real solution better
different foundation levels are presented in different than the standard deviation. If such a standard deviation
columns. The red colour represents the upper estimate, is too low, the allowable load distribution based on the
the blue the lower estimate and the green the mixed LB may not be a conservative one. It should be noted that
approach. As the horizontal scale of fluctuation increases, for the different scenarios, the mixed approach always
a decrease in the allowable load value of the foundation provides a conservative and quite precise (as indicated
and an increase in the safety factor can be observed. The by the small differences between the values provided
safety factor determines the reduction in the ultimate by the different approaches) estimation. Moreover,
bearing capacity that must occur for the assumption this precision could be further improved (if the current
of an acceptable probability of structural failure to be precision level in the given problem is not satisfactory)
satisfied. For the same values of the fluctuation scales, by increasing the number of steps in the adaptive mesh
increasing the foundation level from D=0.0 to 1.0 m means refinement procedure. Since the moments of the upper
a decrease in the factor by 10%–25%. Again, the general and lower approaches can be accurately estimated even
trends are similar to the values obtained in other studies for lower numbers of realisations, the mixed approach
(Pieczyńska-Kozłowska et al., 2015;, Kawa & Puła, 2020). provides a rigorous, conservative and effective approach
 Efficient and conservative estimation reliability analysis of strip footing on spatially variable c - ϕ soil ... 31

Table 4: Allowable loads for β=3.8 obtained for lower LB, upper UB and mixed approach estimated distributions – D=0.0 m.

Scales of fluctuations Design value Safety Factor

Horizontal Vertical LB UB MIXED LB UB MIXED
[kPa] [kPa] [kPa] [-] [-] [-]
θx =1.0 m θz =0.25 m 196.99 206.12 191.09 1.28 1.30 1.32

θz =0.5 m 181.53 186.61 176.61 1.40 1.42 1.44

θz =1.0 m 168.33 171.25 165.75 1.53 1.54 1.56

θx =2.0 m θz =0.25 m 185.03 192.94 178.46 1.37 1.39 1.42

θz =0.5 m 167.21 171.17 161.95 1.53 1.55 1.58

θz =1.0 m 150.77 153.06 148.29 1.72 1.73 1.75

θx =4.0 m θz =0.25 m 176.14 183.70 169.00 1.45 1.47 1.51

θz =0.5 m 156.28 160.01 151.00 1.66 1.68 1.71

θz =1.0 m 140.51 142.40 137.87 1.88 1.89 1.91

θx =8.0m θz =0.25 m 171.18 179.21 163.95 1.51 1.54 1.58

θz =0.5 m 153.78 157.81 148.49 1.71 1.74 1.78

θz =1.0 m 132.59 134.16 129.61 2.03 2.05 2.08

θx =16.0 m θz =0.25 m 170.89 178.54 162.62 1.52 1.55 1.60

θz =0.5 m 150.67 155.24 145.06 1.77 1.79 1.84

θz =1.0 m 133.86 135.50 130.80 2.05 2.06 2.10

for reliability analysis. With these features, the proposed the structure response. This means that 2D analysis of
approach seems to be applicable in engineering practice. strip foundation based on random field lead to crude but
conservatism estimation of bearing capacity probability
distribution. A full 3D analysis of spatially variable soil
3.2 Research significance
was beyond the scope of the current work. However, we
believe that the proposed method for estimating failure
As has been mentioned in the introduction, the application
probability or allowable load values based on the mixed
of the FELA approach to the estimation of the bearing
approach remains applicable to 3D analyses.
capacity of a shallow foundation has a recognisable history.
However, until now, the rigorousness of this technique has
not been fully exploited. The main motivation for writing
this paper is to show a relatively simple application 4 Conclusions
of this rigorousness for the conservative estimation of
probability-based design values. As has been shown in In the present work, probabilistic analysis of bearing
the above section, such an application is possible. capacity of strip foundation on spatially variable c - ϕ soil
It should be noted that the current study is subject to was performed. Variability of both strength parameters,
certain limitations. One of the most obvious is considering i.e., cohesion c and angle of internal friction φ was
the structure in 2D as a plain strain problem. As a described with lognormal stationary anisotropic random
consequence, the value of the out-of-plane fluctuation fields generated with FSM. The analysis included a few
scale was assumed to be infinite. The issue of modelling different values of horizontal and vertical fluctuation
a problem of spatially variable soil, which is, in general, a scales and different levels of foundation. The bearing
3D problem, in 2D space is well described in the literature capacity value was estimated using the FELA method
for different issues (e.g., Liu et al. 2022). As shown in the with an adaptive mesh refinement procedure within the
paper (Chwała & Puła, 2020), omitting the fluctuation in Monte-Carlo framework. For each problem (one pair of
out-of-plane direction can strongly influence the bearing fluctuation scales and one foundation level), 10,000
capacity result by increasing the standard deviation of simulations were computed. The results of the analysis
32 Hubert Szabowicz, Marek Kawa, Wojciech Puła

Table 5: Allowable loads for β=3.8 obtained for lower LB, upper UB and mixed approach estimated distributions – .

Scales of fluctuations Design value Safety factor

Horizontal Vertical LB UB MIXED LB UB MIXED
[kPa] [kPa] [kPa] [-] [-] [-]
θx =1.0 m θz =0.25 m 318.66 341.82 315.19 1.17 1.17 1.19

θz =0.5 m 304.96 319.53 301.46 1.25 1.25 1.26

θz =1.0 m 289.25 299.34 285.95 1.33 1.33 1.34

θx =2.0 m θz =0.25 m 307.40 328.85 303.02 1.22 1.22 1.24

θz =0.5 m 289.56 302.10 285.45 1.32 1.32 1.34

θz =1.0 m 267.29 275.22 263.59 1.44 1.44 1.46

θx =4.0 m θz =0.25 m 294.26 314.12 288.81 1.28 1.28 1.31

θz =0.5 m 273.07 284.18 268.20 1.41 1.41 1.43

θz =1.0 m 247.97 254.68 244.02 1.56 1.57 1.59

θx =8.0 m θz =0.25 m 285.32 304.16 279.22 1.33 1.33 1.36

θz =0.5 m 263.68 274.16 258.67 1.47 1.48 1.50

θz =1.0 m 231.18 237.32 227.17 1.70 1.70 1.73

θx =16.0 m θz =0.25 m 281.51 299.53 275.09 1.35 1.36 1.38

θz =0.5 m 256.46 266.60 251.48 1.52 1.53 1.55

θz =1.0 m 228.21 234.00 224.17 1.74 1.75 1.77

Table 6: Allowable loads for β=3.8 obtained for lower LB, upper UB and mixed approach estimated distributions – .

Scales of fluctuations Design value Safety factor

Horizontal Vertical LB UB MIXED LB UB MIXED
[kPa] [kPa] [kPa] [-] [-] [-]
θx =1.0 m θz =0.25 m 416.52 454.87 414.07 1.14 1.14 1.15

θz =0.5 m 405.84 432.65 403.93 1.20 1.19 1.20

θz =1.0 m 388.64 407.56 386.73 1.26 1.26 1.27

θx =2.0 m θz =0.25 m 405.40 441.97 402.59 1.18 1.17 1.19

θz =0.5 m 389.96 413.64 387.31 1.25 1.24 1.26

θz =1.0 m 366.12 381.61 363.75 1.35 1.34 1.36

θx =4.0 m θz =0.25 m 390.72 425.73 387.11 1.23 1.22 1.24

θz =0.5 m 370.82 393.48 368.37 1.32 1.31 1.33

θz =1.0 m 344.08 357.64 341.39 1.44 1.44 1.45

θx =8.0 m θz =0.25 m 379.37 413.29 375.19 1.27 1.26 1.28

θz =0.5 m 358.78 379.78 355.39 1.38 1.37 1.39

θz =1.0 m 323.79 336.16 320.73 1.55 1.54 1.56

θx =16.0 m θz =0.25 m 374.40 406.94 369.61 1.29 1.28 1.31

θz =0.5 m 347.94 367.66 344.07 1.43 1.42 1.45

θz =1.0 m 316.59 328.09 313.23 1.60 1.60 1.62

 Efficient and conservative estimation reliability analysis of strip footing on spatially variable c - ϕ soil ... 33

allowed to introduce mixed distribution with mean values a conservative and precise estimate of allowable load.
based on lower and standard deviation on UBs which can Due to these features, the proposed mixed approach
be further used for efficient and conservative reliability with the moments obtained based, e.g., on a limited
analysis of the considered problem. number of realisations obtained with FELA, seems to
The following conclusions can be drawn from the be an adequate candidate for a method possible to
paper: apply in engineering practice.
– Foundations are parts of the structures associated
with greater failure consequences than the other
element. Meantime, they are the only parts of the
structure in contact with strongly spatially variable
References
soil. For this reason, this part of the structure might be [1] Ali, A., Lyamin, A. V., Huang, J., Li, J. H., Cassidy, M. J., & Sloan,
recommended for reliability analysis. While the spatial S. W. (2017). Probabilistic stability assessment using adaptive
variability of the foundation can be described with limit analysis and random fields. Acta Geotechnica, 12(4),
the random field, the complicacy of such description 937–948. https://doi.org/10.1007/s11440-016-0505-1
and, more importantly, the high computation cost [2] Ali, A., Lyamin, A. V., Huang, J., Sloan, S. W., & Cassidy, M.
(2016). Effect of Spatial Correlation Length on the Bearing
of associated calculations are the main reasons that
Capacity of an Eccentrically Loaded Strip Footing. In H. W.
such analysis is not popular in engineering practice. Huang, J. Li, J. Zhang, & Chen J.B. (Eds.), APSSRA.
– The long computation time associated with [3] Au, S.-K., & Beck, J. L. (2001). Estimation of small failure
probabilistic modelling of foundations is strongly probabilities in high dimensions by subset simulation.
associated with the assumed method for solving Probabilistic Engineering Mechanics, 16(4), 263–277. https://
doi.org/10.1016/S0266-8920(01)00019-4
individual realisations of random fields. RFELA
[4] Chen, X.-J., Fu, Y., & Liu, Y. (2022). Random finite element
provides a much more effective way of solving analysis on uplift bearing capacity and failure mechanisms of
problems with soil parameter variability modelled square plate anchors in spatially variable clay. Engineering
with random fields. While the class of tasks solved Geology, 304, 106677. https://doi.org/10.1016/j.
with this method is limited to the limit state problems, enggeo.2022.106677
the foundation’s ultimate limit state (bearing capacity) [5] Cheng, P., Guo, J., Yao, K., & Chen, X. (2023). Numerical
investigation on pullout capacity of helical piles under
belongs to this category. Since the computation time
combined loading in spatially random clay. Marine
benefits are significant the method seems to be of Georesources & Geotechnology, 41(10), 1118–1131. https://doi.
particular importance. org/10.1080/1064119X.2022.2120843
– The analysis performed in this paper focused on [6] Ching, J., Wu, T. J., Stuedlein, A. W., & Bong, T. (2018).
single-layer c - ϕ soil with both parameters defined as Estimating horizontal scale of fluctuation with limited CPT
soundings. Geoscience Frontiers, 9(6), 1597–1608. https://doi.
random fields. The influence of horizontal and vertical
org/10.1016/j.gsf.2017.11.008
scales, and depth of foundation on estimated bearing [7] Chwała, M., Komatsu, G., & Haruyama, J. (2024). Structural
capacity results were analysed. As shown, the mean stability of lunar lava tubes with consideration of variable
value slightly increases and the standard deviation cross-section geometry. Icarus, 411. https://doi.org/10.1016/j.
significantly increases with the scale of fluctuation. icarus.2023.115928
As a result, according to a simple reliability analysis [8] Chwała, M., & Puła, W. (2020). Evaluation of shallow
foundation bearing capacity in the case of a two-layered soil
performed based on estimated distributions the
and spatial variability in soil strength parameters. PLoS One,
allowable load decreases and the factor of safety 15(4), e0231992.
increases with the increasing scale of fluctuation. [9] Cho, S. E., & Park, H. C. (2010). Effect of spatial variability of
The effect is less significant for the greater depth of cross‐correlated soil properties on bearing capacity of strip
the foundation. Similar observation comes from other footing. International Journal for Numerical and Analytical
Methods in Geomechanics, 34(1), 1-26.
studies (Pieczyńska-Kozłowska et al., 2015; Kawa &
[10] Ciria Suárez, H. (2004). Computation of Upper and Lower
Puła, 2020). Bounds in Limit Analysis using Second-order Cone
– The most important finding of the current work is Programming and Mesh Adaptivity [Master of Science].
introducing a mixed distribution with mean value Massachusetts Institute of Technology.
based on LB estimation results and standard deviation [11] Dobrzanski, J., & Kawa, M. (2021). Bearing capacity of
based on UB estimation results. As shown the eccentrically loaded strip footing on spatially variable cohesive
soil. Studia Geotechnica et Mechanica, 43(4), 425–437. https://
identification of the distribution moments is possible
doi.org/10.2478/sgem-2021-0035
in reasonable time and the obtained distributions, [12] EN-1990, Basis of Structural Design. (2002).
using the rigorousness of FELA approaches, provide
34 Hubert Szabowicz, Marek Kawa, Wojciech Puła

[13] Engwirda, D. (2014). Locally Optimal Delaunay-refinement and of Soil Properties in Bearing Capacity and Settlement Analysis
Optimisation-based Mesh Generation. of a Strip Footing Using Random Finite Element Method.
[14] Fenton, G. A., & Griffiths, D. V. (2008). Risk assessment in Civil Engineering Infrastructures Journal. https://doi.org/
geotechnical engineering. John Wiley & Sons. DOI:10.22059/CEIJ.2023.360871.1930
[15] Griffiths, D. V., & Fenton, G. A. (1993). Seepage beneath [30] Liu, Y., Chen, X., & Hu, M. (2022). Three-dimensional large
water retaining structures founded on spatially random soil. deformation modeling of landslides in spatially variable and
Géotechnique, 43(4), 577–587. https://doi.org/10.1680/ strain-softening soils subjected to seismic loads. Canadian
geot.1993.43.4.577 Geotechnical Journal, 60(4), 426-437.
[16] Griffiths, D. V., & Fenton, G. A. (2004). Probabilistic Slope [31] Liu, X., Wang, Y., & Li, D. Q. (2019). Investigation of slope
Stability Analysis by Finite Elements. Journal of Geotechnical failure mode evolution during large deformation in spatially
and Geoenvironmental Engineering, 130(5), 507–518. https:// variable soils by random limit equilibrium and material point
doi.org/10.1061/(ASCE)1090-0241(2004)130:5(507) methods. Computers and Geotechnics, 111, 301–312. https://
[17] Griffiths, D. V., Huang, J., & Fenton, G. A. (2009). Influence doi.org/10.1016/j.compgeo.2019.03.022
of Spatial Variability on Slope Reliability Using 2-D Random [32] Lyamin, A. V., & Sloan, S. W. (2002a). Lower bound limit
Fields. Journal of Geotechnical and Geoenvironmental analysis using non‐linear programming. International Journal
Engineering, 135(10), 1367–1378. https://doi.org/10.1061/ for Numerical Methods in Engineering, 55(5), 573–611. https://
(ASCE)GT.1943-5606.0000099 doi.org/10.1002/nme.511
[18] Hicks, M. A., & Samy, K. (2002). Influence of heterogeneity on [33] Lyamin, A. V., & Sloan, S. W. (2002b). Upper bound limit
undrained clay slope stability. Quarterly Journal of Engineering analysis using linear finite elements and non-linear
Geology and Hydrogeology, 35(1), 41–49. https://doi. programming. International Journal for Numerical and
org/10.1144/qjegh.35.1.41 Analytical Methods in Geomechanics, 26(2), 181–216. https://
[19] Huang, L., Cheng, Y. M., Li, L., & Yu, S. (2021). Reliability doi.org/10.1002/nag.198
and failure mechanism of a slope with non-stationarity and [34] Lyamin, A. V., & Sloan, S. W. (2003). Mesh generation for lower
rotated transverse anisotropy in undrained soil strength. bound limit analysis. Advances in Engineering Software, 34(6),
Computers and Geotechnics, 132. https://doi.org/10.1016/j. 321–338. https://doi.org/10.1016/S0965-9978(03)00032-2
compgeo.2020.103970 [35] Makrodimopoulos, A., & Martin, C. M. (2008). Upper bound
[20] ISO 2394:2015, General principles on reliability for structures. limit analysis using discontinuous quadratic displacement
(2015). fields. Communications in Numerical Methods in Engineering,
[21] J. L. Doob. (1990). Stochastic processes. Wiley-Interscience. 24(11), 911–927. https://doi.org/10.1002/cnm.998
[22] Jerez D. J. & Chwała M. & Jensen H. A. & Beer M. (2024). [36] Pieczyńska-Kozłowska, J. M., Puła, W., Griffiths, D. V., &
Optimal borehole placement for the design of rectangular Fenton, G. A. (2015). Influence of embedment, self-weight and
shallow foundation systems under undrained soil conditions: A anisotropy on bearing capacity reliability using the random
stochastic framework. Reliability Engineering & System Safety. finite element method. Computers and Geotechnics, 67,
doi.org/10.1016/j.ress.2023.109771. 229–238. https://doi.org/10.1016/j.compgeo.2015.02.013
[23] Jha, S. K., & Ching, J. (2013). Simulating Spatial Averages of [37] Podlich, N. C. (2018). The Development of Efficient Algorithms
Stationary Random Field Using the Fourier Series Method. for Large-Scale Finite Element Limit Analysis [Doctor of
Journal of Engineering Mechanics, 139(5), 594–605. https:// Philosophy]. University of Newcastle.
doi.org/10.1061/(ASCE)EM.1943-7889.0000517 [38] Podlich, N. C., Lyamin, A. V., & Sloan, S. W. (2014). A
[24] Kawa, M. (2023). Zastosowania pól losowych do opisu Comparison of Conic Programming Software for Finite
anizotropowych ośrodków gruntowych w wybranych Element Limit Analysis. Applied Mechanics and Materials,
zagadnieniach geoinżynierii. Oficyna Wydawnicza Politechniki 553, 439–444. https://doi.org/10.4028/www.scientific.net/
Wrocławskiej (in Polish). AMM.553.439
[25] Kawa, M., Baginska, I., & Wyjadlowski, M. (2019). Reliability [39] Puła, W., Szabowicz, H., & Kawa, M. (2022). Efficient and
analysis of sheet pile wall in spatially variable soil including conservative estimation of failure probability of strip footing
CPTu test results. Archives of civil and mechanical engineering, on spatially variable soil using random finite element
19, 598-613. limit analysis. In J. Huang, D. V. , Griffiths, S.-H. Jiang, A.
[26] Kawa, M., & Puła, W. (2020). 3D bearing capacity probabilistic Giacomini, & R. Kelly (Eds.), 8th International Symposiumon
analyses of footings on spatially variable c–φ soil. Acta Geotechnical Safety and Risk (ISGSR) (pp. 303–308). https://
Geotechnica, 15(6), 1453–1466. https://doi.org/10.1007/ doi.org/10.3850/978-981-18-5182-7_04-007-cd
s11440-019-00853-3 [40] Sert, S., Luo, Z., Xiao, J., Gong, W., & Juang, C. H. (2016).
[27] Kawa, M., Puła, W., & Truty, A. (2021). Probabilistic analysis Probabilistic analysis of responses of cantilever wall-supported
of the diaphragm wall using the hardening soil-small (HSs) excavations in sands considering vertical spatial variability.
model. Engineering Structures, 232. https://doi.org/10.1016/j. Computers and Geotechnics, 75, 182–191. https://doi.
engstruct.2021.111869 org/10.1016/j.compgeo.2016.02.004
[28] Krabbenhoft, K., Lyamin, A. V., Hjiaj, M., & Sloan, S. W. (2005). [41] Simões, J. T., Neves, L. C., Antão, A. N., & Guerra, N. M. C.
A new discontinuous upper bound limit analysis formulation. (2014). Probabilistic analysis of bearing capacity of shallow
International Journal for Numerical Methods in Engineering, foundations using three-dimensional limit analyses.
63(7), 1069–1088. https://doi.org/10.1002/nme.1314 International Journal of Computational Methods, 11(2). https://
[29] Kumar, V., Burman, A., Portelinha, F. H. M., Kumar, M., Burman, doi.org/10.1142/S0219876213420085
A., Portelinha, F. H. M., & Das, G. (2023). Influence of Variation
 Efficient and conservative estimation reliability analysis of strip footing on spatially variable c - ϕ soil ... 35

[42] Zaskórski, L., & Puła, W. (2016). Calibration of characteristic

values of soil properties using the random finite element
method. Archives of Civil and Mechanical Engineering, 16(1),
112–124. https://doi.org/10.1016/j.acme.2015.09.007