## 2 Random Spatial FE Models of the Subsoil  
The issue of the numerical uncertainty definition of the subsoil is analysed in various FE models of a spatial subsoil block presented in Figure 1. The analysis is performed in ZSoil(r) software [31; 32; 33] coupled with Python(r) programming language [34] to aid the series of repetitive parametric computations. The representative random subsoil volume area is modelled as a block with dimensions of \(10\times 10\times 10\) m (all finite elements have dimensions equal to \(1\times 1\times 1\) m), distinguished in pink in Figure 1.  
The boundary conditions at the interface with the surrounding subsoil are reflected by an additional layer of elements whose Young's modulus is \(E_{bc}=560\) kPa and Poisson's ratio is \(\
u_{bc}=0.25\), distinguished in green in Figure 1.  
These parameters are established based on a preliminary numerical test. First, a simple compression test of the boundary element stratum is performed, assuming [PERSON]'s modulus \(E_{bc,t}=60,000\) kPa and Poisson's ratio \(\
u_{bc,t}=0\), so that the boundary element soil stratum does not contribute to lateral displacements of the main subsoil block. The lateral displacement of the boundary layer \(\delta_{bc}=0.000803\) m is determined. Subsequently, this value is substituted as the target lateral displacement of the global subsoil block FE model, with the lateral contribution (\(\
u_{bc}=0.25\)) of the boundary layer assumed. This allows determining the required effective stiffness parameter of the boundary layer as \(E_{bc}=560\) kPa, ensuring proper subsoil deformability in all directions. The total number of finite elements is equal to \(n_{el,tot}=12\times 12\times 11=1584\), whereas only \(n_{el}=10\times 10\times 10=1000\) constitute the regarded subsoil. The subsoil loading \(q=170\) kPa is situated on the soil block indirectly, represented by a flexible metal sheet (plate) with a width \(w=6.0\) m, thickness \(t=0.01\) m, and material parameters corresponding to S235 steel ([PERSON]'s modulus \(E_{s}=210\) MPa, [PERSON]'s ratio \(\
u_{s}=0.30\)). The subsoil is governed by the Coulomb-Mohr material model. In the task, the Young's modulus of the subsoil \(E\) is considered the only random variable, with a Gaussian probability density function (PDF) of the mean value \(\mu_{E}=28,000\) kPa and standard deviation \(\sigma_{E}=7000\) kPa. While the coefficient of variation is \(v=0.25\), there is a \(0.00317\%\) probability that the generated \(E_{i}\) value will be negative; therefore, a truncation parameter of \(t=4.0\) (\(E_{i}\in\mu_{E}\pm t\sigma_{E}\)) is introduced, and hence no negative Young's moduli appear in the subsequent generation. Note that a relatively large Young's modulus variability of the subsoil was purposefully assumed to introduce a significant impact on the scatter of the numerical results of the settlements. Thus, the random variable distribution is not realistic. However, in most design cases in real-life engineering, the data on subsoil variability acquired from the in situ surveys are limited. Hence, suggestions in design codes and literature guidelines vicariously reflect the random parameters, adopting their high variability. As an example,  
Figure 1: The FE model of the subsoil block, featuring the adopted indirect loading. variability \(v=0.3\) is incorporated in the reliability assessment of foundations resting on piles [35] or in stability analysis of a coastal cliff [36]. Both cohesion \(c=10\) kPa and internal friction angle \(\phi=30^{\circ}\) are assumed deterministic but, in fact, can be assumed random. The algorithms proposed in the paper do not limit the number of material parameters to be considered random variables. However, to compare the influence of the adopted FE models on the variability of the computed settlements, [PERSON]'s modulus \(E\) is assumed the only random variable. Various random FE models of the subsoil block are assembled in the ZSoil(r) software, representing the global randomness of the FE model in different ways.