Effectiveness of Random Field Approach in Serviceability Limit State Analysis of Strip Foundation

This work conducts a probabilistic inquiry on how the variability of the parameter defining soil deformability affects the settlement of the foundation located on the soil. The analysis addresses the random foundation model to relevantly estimate the probability of allowable deflection exceedance. The constitutive model parameter is based either on a single random variable or a random field. The computations incorporate direct Monte Carlo sampling and the variance reduction techniques, i.e. Stratified Sampling and Latin Hypercube Sampling. The work focuses on soil parameter modelling by random fields defined by various correlation functions. One of the analytical means is the application of a non-homogeneous function capable of reflecting the soil strata. The probabilistic methods proposed in the paper are tested on a standard example of a foundation strip featuring load eccentricity. It is proved that the modelling mode of the soil—a single variable or a random field—substantially affects the results, i.e. foundation settlements. It was detected that incorporating random fields in soil analysis allows for a valid reliability assessment of a foundation in respect to Serviceability Limit State. The relevantly adjusted correlation functions of random fields allow for a realistic subsoil analysis even in the case of a limited in situ measurement database.


Introduction
While the literature concerning standard foundations is very broad, the search is active for optimal analytical methods to precisely model the real conditions of the soil environment. Due to rising computational capability, the models are applied to consider the limitations and uncertainty of the in situ measurements. The probabilistic methods are applied in the analysis of standard cases of foundations and slope stability, complemented by various case studies. The structural assessments concerns either the foundation load-carrying capacity (ultimate limit state, ULS) or the allowable settlement exceedance (serviceability limit state, SLS). While the prior research rather regarded analytical methods, the current analytical trends are FEM-directed. The computations involve various methods, material models, and detailed solutions adjusted to specific tasks, see (Dey et al. 2019) and their references. Aldosary et al. (Aldosary et al. Abstract This work conducts a probabilistic inquiry on how the variability of the parameter defining soil deformability affects the settlement of the foundation located on the soil. The analysis addresses the random foundation model to relevantly estimate the probability of allowable deflection exceedance. The constitutive model parameter is based either on a single random variable or a random field. The computations incorporate direct Monte Carlo sampling and the variance reduction techniques, i.e. Stratified Sampling and Latin Hypercube Sampling. The work focuses on soil parameter modelling by random fields defined by various correlation functions. One of the analytical means is the application of a non-homogeneous function capable of reflecting the soil strata. The probabilistic methods proposed in the paper are tested on a standard example of a foundation strip featuring load eccentricity. It is proved that the modelling mode of the soil-a single variable or a random field-substantially affects the results, i.e. foundation settlements. It was detected that incorporating 2018) provide a state-of-the-art review of structural reliability analysis (SRA) and stochastic finite element method (SFEM), covering both technology and application aspects. It is easy to observe that various computational methods to analyse the same problem may trigger a relatively large result scatter, even on a deterministic scale, see (Phoon and Tang 2019) and their references. It justifies the use of probabilistic methods to estimate the mean values and dispersion of the results. Hence, the uncertainty of the assumed computational model is considered, the scope of the analysis is not restricted to material parameters only.
Various methods of probabilistic analysis are applied, the classical Monte Carlo (MC) simulation method is a significant one due to its simplicity and the engineering sound description. The advanced algorithms are created to accelerate the crude MC routines, e.g. variance-reduction techniques, Stratified Sampling, SS and Latin Hypercube Sampling, LHS (Dilip and Sivakumar Babu 2014).
Reliability-based robust geotechnical design (RGD) using MC simulation is elaborated (Peng et al. 2017). Response Surface Method (RSM) and Targeted Random Sampling (TRS) became popular as well. The FORM and SORM methods are also permanently developed (Low 2014). The Point Estimate Method (PEM) is also currently applied (Ahmadabadi and Poisel 2015). Their results are always uncertain, especially in the cases of several random variables due to the soil, multi-strata cases, and complex boundary conditions. The computational results are greatly distorted due to limited access to measurement data to regard the assumed, simplified model (Beer et al. 2013;Fenton et al. 2018).
The random field approach to soil analysis broadens its application presently. The works by Fenton and Griffiths (Fenton and Griffiths 2005) apply random fields in both two-dimensional and three-dimensional scales. Specific approaches are incorporated too, e.g. generation based on the Karhunen-Loeve expansion or subset simulation methodology (Ahmed and Soubra 2012). Computations involving advanced correlation functions and result verification with the use of the Weibull variable are included in Puła and Zaskórski (2015). Al-Bittar and Soubra (Al-Bittar and Soubra 2013) applied the polynomial chaos expansion in bearing capacity assessment of strip footing and conducted a sensitivity analysis of material parameters employing the Sobol index. Similar solutions were applied in Drakos and Pande (2016). Al-Bittar et al. (Al-Bittar et al. 2018) incorporated the so-called combining kriging and MC simulation (AK-MCS). Ali et al. (Ali et al. 2017) combined the adaptive FEM techniques with random fields. Ching et al. (Ching et al. 2018) proposed Young's modulus simulation using a homogenization procedure. The in situ experiments to define the soil random field are relatively rare, especially referring to correlation function. Such a complex approach may be found in Goldsworthy et al. (2005); Suchomel and Mašín 2011).
Eurocode 7 (2004) requires the employment of probabilistic methods in structural reliability assessment and safety issues. Although these computations involve advanced, dedicated software operated by users acknowledged with probabilistic methodology, their direct introduction to standard engineering computations becomes a necessity (Low and Phoon 2015;Phoon et al. 2016;Forrest and Orr 2010).
The work provides an insight on random soil parameter modelling, involving either random variables or random fields. The single random variable variant is crude, however straightforward and accessible. An advanced approach incorporates random fields. This procedure requires a precise definition of variables, their mean values, and standard deviations, but first of all their correlation function. The work proposes some correlation functions to simulate the subsoil stratification. Hence a homogenization attempt has been undertaken to successfully reflect the mechanical properties of the soil, see e.g. (Griffiths et al. 2012).
The paper focuses on probabilistic computational methods, a relatively simple example concerns the settlement of an eccentrically loaded foundation. With both standard soil material data and the load parameters assumed, the settlement zone may be predicted successfully. Therefore, the probabilistic analysis is bound to partially complement the deterministic procedure. However, the work is not intended to be a case study only, but rather to display probabilistic computational methods, featuring soil parameters in the form of random fields. The presented algorithms are universal, not restricted to foundation settlement modelling only.

The Model of the Foundation
An example of a continuous foot in plane strain conditions is investigated (Fig. 1). To minimize the impact of boundary conditions, the region is assumed with the following dimensions: 18.0 × 6.0 m. The foundation dimensions are 2.0 × 0.5 m. The subsoil is modelled beneath the footing level, the strata above this level are represented by a uniformly distributed load q = 18.0 kPa at the upper edge of the model (Fig. 1). The structural wall loading acts eccentrically. Its resultant P = 1.0 MN acts vertically, shifted by e = 0.16 m from the foundation axis (the eccentricity e = 1∕12 B, while B denotes the foundation width). The eccentric load model is assumed to produce an asymmetric solution even in a homogeneous subsoil case. The deterministic analysis of such cases is found in e.g. (Benmebarek et al. 2019) and probabilistic in e.g. (Kasama et al. 2019). Only the SLS is analysed. It is worth pointing out that the limit state function is not known due to non-linear computational characteristics, thus employing FE software is necessary here. The analysed model is not a classical case study, but rather an example to illustrate various probabilistic computational approaches proposed in the work.
Probabilistic analysis requires multiple repetitions of numerical calculations thus it is necessary to assume an optimal FE mesh density. The discretization should help to determine accurate results in a relatively short computational time for a single computational case. The sub-foundation domain is divided into two subdomains of various FE mesh densities (Fig. 1). The material of subdomain I (10 × 3 m) is described with a probabilistic model, while subdomain II is assumed deterministic. The subdomain I is described as a regular mesh of two-dimensional continuum four-node elements of the following dimensions: 0.25 × 0.1875 m (40 × 16 elements). Relatively small dimensions of subdomain I and its crude FE discretization originate from a prior parametric analysis incorporating thousands of FE models aimed at finding an optimal compromise between result accuracy and computational time. The numerical details of the FE model are collected and presented in Table 1.
In the probabilistic subdomain I the following mean value of Young's modulus E I is assumed  E I = 50, 000 kPa, while cohesion c = 10 kPa, and internal friction angle = 30 • are considered deterministic. The deterministic subdomain II exhibits constant values: E II = 30, 000 kPa, c = 10 kPa, and = 30 • . The homogeneous subdomain II exhibits reduced medium stiffness compared to subdomain I, to additionally minimize the boundary condition impact on the response, it simulates an infinite halfspace surrounding the analysed region of subsoil I, and reproduces the deterministic support of the random subdomain. Additional computations showed that the assumption of E II = 50,000 kPa for the subdomain II does not substantially affect the obtained mechanical response.
The computations assume the Mohr-Coulomb model to capture soil plastification, focusing on the foundation edge regions. Concrete is considered elastic, its parameters correspond to C30/37 concrete (PN-EN 2006+A1:2016-12 2016) of minimum compressive strength 37 MPa and Young's modulus E c = 32.8 GPa. The computations are conducted in the ZSoil commercial software (Commend et al. 2018;Truty 2018). ZSoil enables the implementation of the advanced Hardening Soil (HS) constitutive model (Cudny and Truty 2020). While the work covers only the parametric studies of random subsoil impact, the standard M-C model was chosen instead, as it is widely used in numerous other commercial packages. The advanced research on real-life structural settlement conditions makes it possible to improve the computations through the HS model application, see e.g. (Winkelmann et al. 2021).
The magnitude and eccentricity of the load were assumed deterministic, the work focuses solely on the subsoil impact on the foundation settlement. The probabilistic approach to loads requires additional assumptions and should be analyzed separately.
First of all deterministic computations are performed incorporating mean material parameters of subdomain I E I = 50, 000 kPa (Fig. 1). The following settlements of the foundation points are obtained: u l = −0.064 m(left edge), u r = −0.036 m, (right edge), and u deter = −0.050 m (midpoint). The foundation rotation, in a rigid body pattern, equals deter = 0.0142 rad. In the case of axial load action the deflection is u axial = −0.049 m. Note that the example is theoretical only, enhanced foundation settlement and rotation allow to comprehensively test the proposed algorithms. Only the plane strain case is considered, the impact of subsoil variability along the foundation is therefore neglected, see e.g. (Przewłócki 1999;Przewłócki and Górski 2001).

Single Random Variable Description
The probabilistic analysis requires the assumption of random variables, their distributions, and corresponding parameters. It is assumed that Young's modulus E is the only parameter to represent subsoil variability in subdomain I (Fig. 1). Such a simplification is justified because of the analytical domain of the SLS (as opposed to the ULS). Young's modulus is represented with a normal variable (N) with the mean E N = 50, 000 kPa, and standard deviation E N = 15, 000 kPa (the coefficient of variation equals E N = 0.3). The data on the material variability are provided in e.g. (Phoon et al. 2006;Lumb 1966).
While Gaussian variables take negative values the symmetrically bounded variable was applied here, the bounding parameter was adopted t = 3.0 . This way, the lowest Young's modulus value is 5000 kPa (features alluvia and peats), whereas the highest equals 95,000 kPa. The values beyond this interval seem unrealistic from an engineering standpoint. All material data concerning the subdomains are summarized in Table 2. A population of 10,000 samples was generated according to appropriately modified truncated Gaussian variable (N), the sample histogram is presented in Fig. 2. The estimation of mean value, standard deviation, and skewness read: Ê N = 49, 983 kPa, E N = 14, 760 kPa and ̂3 N = 0.0146. Note that the generated distribution differs from the standard type only to a slight extent. Given the truncating parameter t = 3.0 the probability of sample generation from beyond the ( The geotechnical analysis often assumes lognormal (LN) variables to represent the soil variability. Negative Young's moduli are discarded this way. Figure 2 compares the histograms generated with the use of two types of variables. To adjust the N and LN forms the following transformations are applied: To directly compare the LN and N outcomes the mean value, standard deviation, and skewness are estimated too: Ê LN = 49, 584 kPa,̂E LN = 14, 880 kPa and ̂3 LN = 0.974. Figure 2 specifies the maximum and minimum generated Young's moduli ( E lim(N) and E lim(LN) ). Note that irrationally high Young's moduli are generated. The bounded Gaussian probability density function allows for arbitrary two-side random variate limitation according to laboratory tests, in situ measurements, and engineering experience. It is worth pointing out that alternative generation methods may be employed in the computations, e.g. beta variable or hyperbolic tangent transformation (Puła and Zaskórski 2015;Fenton and Griffiths 2003). The selection of variable types to represent Young's modulus and other soil parameters is essential in geotechnical computations (Jimenez and Sitar 2009).
The first step checks the impact of Young's modulus change on the foundation settlement and slope. These parameters have been computed given various Young's moduli in the range E The results of 11 analysed cases are shown in Fig. 3. They point out a non-linear mechanical response of the subsoil. The maximum foundation settlements exceed 0.25 m, this value is unreal from an engineering viewpoint. Based on preliminary computations the truncating parameter t of Young's modulus E N variable may be corrected. However, it should be emphasized that the minimum Young's modulus E min N = E N − 3 E N = 5000 kPa does not lead to the detachment of the foundation from the subsoil, thus making the model compliant with engineering practice.

Crude Monte Carlo Method
The first step of the analysis includes the direct MC computations. The predefined series of computational procedures were performed in Python (Python manual 2019) and ZSoil packages (Commend et al. 2018). The computational time of a single realization is relatively short, allowing a multitude of automatic runs to be conducted, making the results reliable in the case of a 10,000 sample population. While Young's modulus E is assumed a bounded Gaussian variable (Fig. 2), the mean value, standard deviation, and skewness of the foundation center nodal deflection are achieved as follows: û N = −0.0530 m,  Figure 4 shows the histogram of the results, accompanied by extreme settlements based on N and LN computational variants (u lim(N) and u lim(LN) ). The minimum structural rotational slope is min = −0.0091 rad, the maximum value is max = −0.1082 rad. Comparing the distributions of Young's modulus E (Fig. 2) and the foundation settlement u (Fig. 4), a transformation is shown of diversified input variables N and LN into strongly asymmetric output variables.
The compared results N and LN are almost identical due to mean values, standard deviations, and the histogram shape (Fig. 4). The accordance of input histogram patterns is justified by their skewness coefficients (0.14% dispersion). It is worth pointing out that the positive assessment of N and LN convergence cannot be generalized onto arbitrary engineering models, as each engineering challenge should be analysed adequately. Further analysis was conducted with the use of a truncated Gaussian variable to represent soil Young's modulus E N .

Stratified Sampling
The next step attempts to reduce the number of realizations necessary to estimate the mean value and standard deviation. Stratified Sampling (SS) (Hurtado and Barbat 1998) is applied here. A population of 10,000 generated Young's moduli E N is divided into 3, 9, 27, and 81 equal strata (ES). Each interval is represented by the values close to its center. Thus the computations conducted in the previous steps are directly implemented into the next steps, hence 81 computational courses are performed. An example shows a histogram of Young's modulus E N realization domain divided into 27 equal sections of different probabilities, presented in Fig. 5, accompanied by 27 computational points. The results are presented in Table 3 and Fig. 6.
The convergence of mean value and standard deviation in the case of 3, 9, and 27 element division is sufficient (Table 3). While rising the number of intervals the solution stabilizes and the results correspond to the case of 10,000 crude MC popula- Table 3 includes the results of 3, 9, 27, and 81 samples achieved by the direct MC method. The latter compared with the SS results prove the high advantages of the variance reduction techniques.

Random Field Implementation
The software to generate random fields has been recently developed either in the form of distinct algorithms (Ahmed and Soubra 2012;Grigoriu 2003;Papaioannou and Straub 2017) or dedicated to subsequent problems Sudret and Kiureghian 2002). Nevertheless, the methodology of random fields is barely incorporated in engineering applications due to the need to apply specific software, and the difficulties in both assuming appropriate input data and interpreting results. Such an analysis usually requires an advanced background in probabilistic methods and structural The work also incorporates this computational routine. The software is employed to generate Gaussian random fields, addressed in Bielewicz and Górski (2002), combined with ZSoil (Commend et al. 2018) and Python (Python manual 2019). Note that the proposed computational algorithm makes it possible to apply any other software generating random fields and any FE program to analyse the task.
The software applied in the work conducts random variable generation (Bielewicz and Górski 2002) employing the conditional acceptance and rejection concept (Devroye 1986). The software makes it possible to simulate discrete Gaussian random fields of an arbitrary given homogeneous or non-homogeneous, two-or three-dimensional correlation function.
The conditional function f s ( u ∕ k ) is incorporated to allow for direct generation of an unknown vector u assumed the relevant random field region already specified k where c and c are a conditional correlation matrix and vector of mean values respectively, and s can be defined as: It is easy to notice that high values of the truncation parameter t (e.g. t ≈ 5) make the parameter s tend to zero.
The essential feature of the software is the generation of unknown field point values in a sequentially shifted region. Thus it is possible to consider the fields of arbitrarily large dimensions. The software has been thoroughly checked and verified to find application in many engineering tasks, e.g. silos (Górski et al. 2015) or soils (Tejchman and Górski 2011).

Homogeneous Random Field
The most important issue in the random field modelling is to assume a relevant correlation function of the problem (Sudret 2008). Parameters governing the .  The random field description of standard engineering materials (subsoil, concrete, composites, etc.) usually employs an exponential correlation function related to the first-order autoregression function or Markov process. It takes the following form in the univariate case, according to Vanmarcke (1983) where |x| is the distance between the field points, and d x is the so-called damping parameter (Knabe et al. 1998). The correlation function (3) Figure 7 presents the impact of the assumed damping parameter value d x (Eq. 5) on the correlation between the points.
Upon comparison of the diagrams in Fig. 7, the initial adjustment of the field parameters relative to the problem becomes possible (Sudret 2008). The information on correlation function is usually supplemented by an additional measure-the scale of fluctuation-proposed by Vanmarcke (1983) In the case of correlation function (Eq. 5) the fluctuation scale is = 2∕d x . Different methods of field classification exist, e.g. based on incorporating information entropy (Walukiewicz et al. 1995). Relevant classification methods allow for the appropriate selection of correlation functions according to the problem. However, linking the function with the insitu measurements is very difficult, and requires further in-depth analysis.
While defining a correlation function referring to a specified linear section Δx (e.g. FE element dimension) it is necessary to assess local averages integrating Eq. (5). Uniform division yields a formula (Knabe et al. 1998) where L x is the distance between centers of the averaging intervals.
Note that a sufficiently small Δx yields small variations in computed results according to Eqs. (5) and (7). The application of Eq. (7) makes the random field values and the FE mesh assumption independent. That was the way to minimize the impact of a sparse sub-foundation region division (Fig. 1). Equation (7) makes it possible to capture the differences in both horizontal and vertical dimensions of finite elements.
The work employs two correlation function types. The first generalizes Eq. (5) for a two-dimensional case, i.e. a combination of two coordinates where d x and d y are damping parameters capturing the correlation decay, |x| and |y| are the distances between field points along these axes.
The field defined by Eq. (8) may be considered homogeneous and anisotropic in a general case. The averaged version makes it possible to apply Eq. (8) in FE computations takes its form from the product of relevant Eq. (7). It seems that only moderately and strongly correlated fields capture the soil conditions related to the natural homogeneity of soils within respective subdomains. Therefore, a mean correlation range i.e. d x = d y = 0.3 m −1 (Fig. 7) is assumed between Young's moduli in the points of the subdomain I (Fig. 1). The assumption d x = d y simplifies the issue, while the real situation suggests a variable correlation in vertical and horizontal scales. Note that damping parameters distinction in vertical and horizontal directions simulates the soil stratification, thus this action is not limited to non-homogeneity display The field is modelled by the Gaussian bounded function of mean value E = 50, 000 kPa, standard deviation E = 15, 000 kPa, and bounding parameter t = 3.0 (Eq. 4). Thus the basic random field data are identical to the case of a one-variable field (Chap. 3), to compare these two sets of results.
A population of 10,000 field realizations (vectors) i (i = 1, 2, ..., 10000) was generated according to the subdomain I division into 16 × 40 finite elements (Fig. 1). It should be pointed out that appropriate modelling of a soil medium by Eq. (8) or other correlated fields requires a broader FE domain than the one in the model (Fig. 1), and a denser FE discretization. However, the parametric study to compare various field patterns makes it relevant to use the assumed FE model. This does not hold in the case-study type analysis. An example of generated field realizations is shown in Fig. 8. The impact of the assumed decay parameters d x and d y on the generated fields is available in Zyliński et al. (2020), including strongly correlated (d x = d y = 0.1 m −1 ) and uncorrelated (white noise) cases.
The next important analytical stage performs the correctness check of a generated field. The estimator of the covariance matrix ̂ E is where i is a vector to represent a single generated field and ̂ is a vector of mean values of the generated fields. The norm of the covariance matrix A global error of generation, concerning the corresponding theoretical value (Eq. 8) is estimated: The error points out the possibility to apply the generated fields in the assessment of foundation settlements (Walukiewicz et al. 1995). Additionally, Fig. 9 compares theoretical (Eq. 8) and generated (Eq. 9) correlation coefficients between the first and subsequent 40 upper row elements belonging to the subdomain I (Fig. 1). The estimation of the generated matrix ̂ E shown in Fig. 9 is almost perfect. First of all, a comparative analysis is conducted employing the MC method regarding 10,000 random field realizations. Every operation is performed automatically by the author-made FORTRAN-coded software, cooperating with commercial platforms (ZSoil and Python). The histogram of the obtained foundation settlements is shown in Fig. 10, the estimators are included in Table 4 with other analytical variants. The difference is substantial between the diagrams in Fig. 10 (Young's modulus represented by a random field RF) and Fig. 4 (Young's modulus represented by a single random variable E N ). While the mean values of the results are similar, their standard deviations and the deflection ranges differ significantly.
The computations concerning a single random variable prove that the selected MC variance reduction method (i.e. SS) leads to the results of engineering sound accuracy. It is possible to conduct the related analysis in the random field domain (Hurtado and Barbat 1998;Tejchman and Górski 2011). Two reduction techniques have been used: SS (two variants) and LHS. These methods may be applied while the generated random fields i (i = 1, 2, ..., 10000) are appropriately classified. Two field-classifying parameters have been specified, i.e. mean values of each field E i and maximum difference between Young's moduli of a single field The parameter E extr i denotes the field value gap (a lower parameter E extr i yields a "flat" field shape). Moreover, the parameter E extr i is independent of the mean values of the fields E i , adding substantial information on the generated realization domain. The joint probability density functions are presented in Fig. 11, a dot represents a single random vector i defined by its mean value E i and the field value gaps E extr i . This is the so-called anthill-type representation, it yields that the classification does not trigger parameter correlation, hence the basic requirement to apply the Stratified and LHS methods is fulfilled.
The determined intervals of classification mean values are 21, 906 ≤ E i ≤ 81, 279 kPa and 35, 538 ≤ E extr i ≤ 89, 801 kPa referring to mean values and the field value gaps, respectively. These two variables are represented by marginal distributions, each one is divided into nine sections of ES (Fig. 11a) and equal probabilities (Fig. 11b).
First of all classical SS is incorporated, assuming ES division SS-ES (Fig. 11a). Hence the Fig. 10 Histogram of foundation settlements u RF determined by the MC simulation method, homogeneous random field description, the population of 10,000 random fields)  . 11 The anthill-type distribution of the generated 10,000 homogeneous random fields, the cases of equal strata (a) and equal probability (b) computations consider both parameters classifying the random field. The entire range is split into 9 × 9 = 81 sections (Fig. 11a), however, due to the sample distribution in 11 sections at the edge regions no samples are present. Thus the computations cover 70 samples only, chosen from the central parts of each section. The probabilities assumed for selected sections differ substantially, affected by the number of samples. The results are included in Table 4 (the SS-ES method). Next, the computations are performed by the SS method, dividing the domain into equal probability sections SS-EP (Fig. 11b). This approach directs the computations into samples located near the mean values, reducing the contribution of the anthill region margins. The results are presented in Table 4 (the SS-EP method). Note that taking a higher impact of anthill central region samples (Fig. 11b) into account results in a slight computational variation only.
The next step checks out the possibility of further sample space reduction. The LHS is applied here. Nine out of the 81 subregions are randomly sampled from the interval limited by vertical and horizontal lines, see Fig. 11b. The samples selected for the computation are chosen from the central sections of the zones marked by boxes (Fig. 11b). The computational results are included in Table 4 (the LHS method). It may be stated that the LHS method allowing for multiple reductions of the number of analysed samples (in this case 81 to 9) is proved successful.
To assess the impact of the assumed field classification the standard SS computations are provided, regarding a single classification parameter-mean values of selected random fields. Nine samples are chosen from the centers of the intervals, separated by vertical strands, see Fig. 11b. Here, the second variable classifying the field (field value gap) is not taken into consideration. This computational variant resembles the analysis of the soil defined by a single random variable (Chap. 3.2). The results included in Table 4 (the SS method) slightly differ from the LHS method variant.
The computational results of foundation settlement resting on a subsoil and defined by a single random variable (Table 3) and random field (Table 4) draw to the conclusion that both cases yield similar mean values ( û N = −0.0530 m, and û RF = −0.0529 m, respectively). However the standard deviations differ substantially: ̂u N = −0.0130 m and ̂u RF = −0.0074 m.
A low random field result scatter is confirmed by the histograms in Fig. 10, compared to the ones in Figs. 4 and 6 (a single random variable model for a field). The lower bound for the single variable solution is u lim(N) = −0.27 m (Fig. 4) while the random field yields u lim(RF) ≈ −0.12 m (Fig. 10). It proves that the random field results are more realistic in the considered case here. The results in Table 4 do not differ substantially. The most accurate solution to the direct MC case (Table 4, MC method) was achieved in a classical SS algorithm, splitting the generated region into nine sections of equal probabilities (Table 4, SS method). The application of LHS with the sample number left unchanged reduces the mean value and standard deviation to a slight extent only (Table 4, LHS method). However, this conclusion is not general as a relatively simple example is analysed.

Non-homogeneous Random Field
In the second case, a second-order random field is applied where x is the so-called frequency parameter (Knabe et al. 1998), its use allows for positive and negative values of Eq. (12). Figure 12 presents possible correlation function shapes assuming a constant decay parameter d x = 0.15 m −1 and various frequency parameters x = 0.5, 1.5, 3.0 [m −1 ].
In a two-dimensional field case, it is possible to freely shape the layout of quasi-strata. The effects of multi-layered soil on foundation settlement are (12) (x) = e −d x |x| cos( x |x|), Fig. 12 Graphical presentation of non-homogeneous correlation function (Eq. 12) analysed in Kuo et al. (2004);Huang et al. 2018). As an example, taking a cosine function affected by the y parameter only, the field is quasi-stratified in the vertical direction Taking various values of y it is possible to simulate the samples of quasi-stratified layout-the higher the parameter value, the higher density of quasi-strata of diverse Young's modulus E. While the parameters y show small values in subdomain I (Fig. 1) no stratification occurs. However, to achieve limited Young's moduli dispersion in particular quasi-strata, it is suggested to increase the correlation of the stratum points. Thus, with regard to the computations conducted in Chap. 4.1. the half-value decay parameter d x = d y = 0.15 m −1 and frequency y = 1.5 m −1 are assumed. It is worth pointing out that the generated fields show high variability, e.g. Figure 13a shows a stratified feature while Fig. 13b makes it hard to observe. Many other generated field types show a homogeneous pattern due to the strong correlation between points. Nevertheless, it is possible to apply Eq. (13) to map the soil parameters which depend on geological processes. Note that this stratification does not reflect e.g. the in situ experimental results, instead, the averaged values only. Therefore the MC computations including stratification simulation may be regarded as material parameter homogenization in the sub-foundation zone.
On the other hand, Fig. 14 shows examples of fields generated with the help of a full two-dimensional correlation function The generation assumes the prior parameter values: d x = d y = 0.15 m −1 and x = y = 1.5 m −1 . Appropriate variants of parameters x and y allows for shaping the quasi-strata slopes.
The computation of foundation settlements employs the field defined by Eq. (12), regarding only horizontal stratification. 10,000 Realizations are generated here. In this case, no crude MC computations have been conducted. The SS variant with 9 × 9 = 81 samples has been instead employed, in the variant of equal probability intervals. Two variables defining the field are assumed: mean values E (NH)i and maximum differences between Young's moduli E extr (NH)i . The following results are obtained: mean settlements û NH = −0.0521 m and standard deviation ̂u NH = −0.0033 m . While the mean value does not exceed the previous results for both the single variable case (û N = −0.0530 m) and the homogeneous field case (û RF = −0.0529 m) the standard deviation according to these models is far beyond ( ̂u N = −0.0130 m and ̂u RF = −0.0074 m, respectively). It may be clarified based on field mean values-Young's moduli. In the single value case, the bounded Gaussian variable, the range is 5165 kPa ≤ E N ≤ 94, 810 kPa (Fig. 2), the mean Young's modulus of homogeneous random fields (14) (x) = e −d x |x|−d x |y| cos( x |x|) cos( y |y|). 21, 906 ≤ E i ≤ 81, 279 kPa (Fig.  11), and 43, 052 kPa ≤ E (NH)i ≤ 56, 926 kPa in the case of non-homogeneous fields. The non-homogeneous field shows the lowest dispersion due to the features of a strongly correlated function. All the obtained results are summarized in Table 5.

Reliability Estimation
Determining the settlement distribution for 10,000 realizations of a single random variable and random field case (the homogeneous variant, Eq. 8) allows for reliability assessment incorporating the classical MC. The direct MC approach estimates the probability of failure p f in the entire domain of , using the expression applying the indicator function I 0∕1 where LS( ) is the limit state function, and NS denotes the number of samples used in approximation.
The impact of both computational variants on the foundation reliability assessment referring to allowable midpoint deflection is presented in Fig. 15. The computations assume the limit settlements variable −0.28 m ≤ u lim ≤ −0.05 m. The highest difference between the solutions may be observed while comparing the settlement range of zero damage probability.
In the case of a single variable description it reads u lim(N) = −0.27 m while in the random homogeneous field case it equals u lim(RF) = −0.12 m (Fig. 15), i.e. less than half the value. As an example, the level u lim = −0.09 m triggers the estimated reliabilities p f (N) = 0.0183 and p f (RF) = 0.0015, respectively, the ratio is greater than 10 here. Above the level of u lim = −0.08 m the rapid growth of failure probability occurs in both analytical cases. A comparative look into two diagrams and analysing of prior results yield a conclusion that the definition incorporating random fields is more realistic.
It is worth pointing out that a reliability estimation employing a direct MC method is not possible in the case of real engineering structures due to the necessity to perform thousands of FE computations. In these cases, other well-known methods may be applied, such as e.g. RSM, TRS, PEM, and many more addressed in detail in a great many applicationoriented works. The commercial software packages are bound to support these computations.  15 Reliability estimation with the use of a single random variable and the random field in terms of a limit settlement function

Conclusions
The work addresses the analysis of an eccentrically loaded foundation strip. The subsoil under the foundation is defined in random terms through a single random parameter-Young's modulus. The computations are conducted by a joint effort of the author's procedures and commercial software, e.g. ZSoil and Python. Appropriately prepared, automated algorithms made it possible to generate and compute thousands of foundation models. The probabilistic analysis makes it possible to estimate the mean values of foundation settlements, their standard deviations, and reliability in the function of the assumed limit settlement threshold. First, Young's modulus is defined by a single variable. In both cases of variable types, Gaussian and log-normal their impact on final results is analysed. It is demonstrated that the bounded Gaussian pattern allows to limit the generation domain and to discard excessively small and large Young's moduli E.
The major goal of the work is to compare the soil definitions through a single random variable and random fields. The homogeneous field simulates the possible subsoil parameter variability while the nonhomogeneous field allows for introducing an additional factor, i.e. soil stratification. The selection of correlation function parameters makes it possible to capture strata thickness and their orientation. The random fields are more realistic than random variables in the light of soil modelling.
It should be pointed out that the random field generation is bound to average the detected soil condition parameters, not mapping them in their entirety. The probabilistic analysis is directed onto multiple computations of models with various random fields. Hence medium homogenization is conducted. The computations prove that the Gaussian-modeled Young's modulus, complemented by various correlation functions relevantly reflects the subsoil performance at various stiffness ranges. This approach shows key importance in the case of limited and partial in situ data concerning the subsoil under the designed foundation.
The work applies various MC versions. It is proved that for random fields in a single-variable modelling variant (Young's modulus) a relevant classification of generated fields allows for using SS and LHS methods. This action further reduces the computational effort reducing the population to several samples only.
The computations are performed due to theoretical examples to provide a detailed parametric analysis. However, the proposed algorithms may be applied directly e.g. in foundation analysis of dedicated, special structures of high sensitivity to non-uniform settlements.
Author Contributions All authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by KŻ and JG. The first draft of the manuscript was written by KŻ, KW, and JG and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.
Funding The authors declare that no funds, grants, or other support were received during the preparation of this manuscript.

Data Availability
The datasets generated and analysed during the current study are not publicly available, but are available from the corresponding author on reasonable request.

Conflict of interest
The authors have no relevant financial or non-financial interests to disclose.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.