This paper presents a comprehensive statistical analysis to evaluate the uncertainties inherited in the lower–upper matrix (LU) decomposition-based random field theory, which is commonly used in geotechnical engineering to model the random spatial variability of soil properties. To facilitate this study, several realisations of Monte Carlo simulation coupled with LU decomposition-based random field theory are made for varying values of input statistical parameters. Thereafter, the statistical parameters are estimated for each simulated random field using the method of maximum likelihood. The results obtained in this study indicate that the LU decomposition-based random field theory does not yield the random field as it is anticipated to be for given input parameters. Ideally, the statistical parameters estimated from the simulated random fields should be nearly equal to the values of input parameters; however, there is significant randomness in the statistical parameters of the simulated random fields. Furthermore, this paper also presents the analysis to evaluate the factors affecting the uncertainties in the random fields simulated by LU decomposition-based theory by varying the input statistical and geometrical parameters.
Notation
intercept on axis
slope of trend line of variation of along the depth
coefficient of variation of X
average of vertical distances between intersection of soil property curve and fitted trend line
- L
soil domain depth
cone tip resistance in MPa
- X
soil property
percentage relative error
standard deviation of log-normal distribution of X
input scale of fluctuation
mean of log-normal distribution of X
mean value of X
correlation matrix
lag distance
Introduction
The soil is formed through several intricate natural processes that unfold over extended geological time periods. Due to this, the geotechnical properties vary over short distances of just a few metres, which results in an inherent spatial variability of the ground (Baker, 1984; Fenton and Vanmarcke, 1990; Phoon and Kulhawy, 1999a, 1999b). Exact characterisation of soil spatial variability requires extensive geotechnical test data to obtain the complete information of the soil properties of the ground of interest, which is practically near to impossible. The presence of considerable inherent variability and limited geotechnical data gives rise to randomness in the quantification and modelling of soil spatial variability; therefore, treating the ground as a spatially varying random field is the most appropriate approach. Several studies have been done to assess the influence of spatially varying random ground on the behaviour of geotechnical engineering systems, including shallow foundations, slope stability, and pile foundation (Fenton and Griffiths, 2002, 2003, 2004a, 2004b, 2005; Haldar and Babu, 2008). These studies revealed that the soil spatial variability significantly affects the geotechnical system response; hence, it cannot be ignored. Therefore, the characterisation and modelling of spatial variability become of utmost importance in the risk assessment of geotechnical engineering systems. Over the period of time, several methods for probabilistic characterisation and modelling of the random spatial variable ground are developed. The method of moments and the method of maximum likelihood are commonly used methods to characterise the spatial variability of the ground (Cami et al., 2020; DeGroot and Baecher, 1993; Fenton, 1999a, 1999b). There are several mathematical models available, such as LU decomposition, fast Fourier transformation, and local average subdivision methods to model the random spatial variable ground (Fenton, 1994; Fenton and Vanmarcke, 1990; Jha and Ching, 2013). These mathematical models are generally known as random field theory.
Among these random field theories, LU decomposition-based random field theory is mostly used in practice (Cheng et al., 2019; Haldar and Babu, 2008; Li et al., 2015; Singhai and Vikash, 2022; Singhal et al., 2023). The LU decomposition-based random field theory follows Cholesky’s decomposition method coupled with Monte Carlo simulation, and it requires three parameters, namely, mean, variance, and scale of fluctuation, to model the stationary spatially varying random soil domain. Although the LU decomposition-based theory is commonly used to assess the risk associated with geotechnical engineering systems due to random spatially varying ground, a critical evaluation of the performance of this theory is yet to be done. Considering this, the LU decomposition-based random field theory is analysed in this paper, wherein the following points are critically evaluated.
Does the theory yield the simulated random fields as it is expected to for a given set of input parameters? Are the statistical parameters characterising the simulated spatially varying ground the same as those given as the input parameters to the random field theory, or are there some uncertainties involved in these?
Are the random fields simulated by the LU decomposition-based theory stationary?
What are the factors that enhance the randomness in the simulated random fields significantly?
The evaluation of these points is important, as the uncertainties in the random fields generated using LU decomposition-based random field theory will lead to an error in the risk assessment of geotechnical engineering systems. To achieve this, a comprehensive stochastic analysis of the uncertainties involved in modelling the soil spatial variability is performed, which is described in subsequent sections.
Problem statement
The inherent soil spatial variability is concisely described by mean, coefficient of variation, and scale of fluctuation. Based on an extensive literature review, Phoon and Kulhawy (1999a) summarised the values of these three statistical parameters for various mechanical properties of soils. They reported that the coefficient of variation of soils varies from 10% to 50%, and the scale of fluctuation lies between 0.2 and 6 m. Details of these can be found elsewhere (Phoon and Kulhawy, 1999a). In geotechnical engineering, the spatial variability of ground is commonly characterised by cone penetration test data, as it provides a continuous profile of cone tip resistance, , and sleeve friction, (Onyejekwe et al., 2016). Phoon and Kulhawy (1999a) reported that the mean value of of the soils ranges from 0.5 to 30 MPa, depending on the soil type. Considering these, this study considers a soil domain of depth (L) 20 m with as a soil property with a varying degree of coefficient of variation and scale of fluctuation. The cone tip resistance () is denoted hereafter by the term X for brevity. Numerous sets of 2000 one-dimensional (1-D) random fields are realised for X with a given set of input parameters. The soil domain is discretised into 100 equidistant nodes for the simulation of random fields. The realisations are made using the Monte Carlo simulation technique with the input parameters: mean, = 25 MPa; coefficient of variation, , varying from 10% to 30%; and the scale of fluctuation, , varies from 0.5 to 5 m. The soil property, X, is assumed to follow a log-normal distribution as it is always non-negative. The statistical analysis performed in this paper is categorised into two major sections, viz, the estimation of statistical parameters characterising the generated random fields, and the evaluation of uncertainties in the statistical parameters of the random fields. The former section is primarily focused on finding out the suitable method for the estimation of parameters of soil spatial variability, followed by statistical hypothesis testing. To accomplish this, a comparative analysis is made among the mostly used statistical estimation methods, namely, the method of moments and the method of maximum likelihood. The comparison between these two methods is done to determine the most appropriate method to characterise the soil spatial variability. Besides this, the goodness-of-fit test, the K–S test (Kolmogorov–Smirnov test) is performed to verify whether the generated random fields follow the given log-normal distribution or not. The latter section of the paper evaluates the randomness in statistical parameters of 1-D random fields to assess the potential factors influencing the uncertainties in the random fields generated by LU decomposition-based random field theory. Several rigorous statistical analyses are made on various sets of random field realisations by varying some of their input parameters (, , and L) and observing its impact on their output statistical parameters. Numerical codes for Monte Carlo simulation coupled with random field theory, method of moments, and method of maximum likelihood are developed in this study using the Python programming environment to perform these statistical analyses.
LU decomposition-based random field theory
The LU decomposition-based random field theory considers that the spatial correlation of the soil property, X, follows the Markovian exponentially decaying function, which is defined as
where is the correlation coefficient defining the correlation of the soil property at and nodes, is the absolute distance between these two nodes, and is the scale of fluctuation. For n number of nodes, the correlation matrix, [], which consists of all the correlation coefficients, is written as,
As the correlation matrix is a positive symmetric matrix, it can be decomposed into a product of lower and upper triangular matrices using the Cholesky decomposition technique, thus it is written as
where [L] is the lower triangular matrix.
The soil property, X, is assumed to follow a log-normal distribution, whose probability density function (pdf) is defined as
where and are the parameters of the log-normal distribution, which are related to the parameters of the normal distribution and as shown.
To generate a random field of the soil property, the LU decomposition-based random field theory introduces a vector G of size 1, which is defined as
where z is the standard random normal variate with zero mean and unit variance. The spatially varying random field of soil property, X, is then determined by utilising the vector, G, and log-normal distribution parameters (, ) as
Thereafter, a Monte Carlo simulation is performed to generate distinct spatially varying soil properties in each realisation.
Estimation of statistical parameters defining soil spatial variability
Three parameters are required to estimate the randomness in soil spatial variability: mean, variance, and scale of fluctuation. These parameters are generally estimated using the method of moments and the method of maximum likelihood. This section describes both the methods in brief, estimates the statistical parameters of the simulated random fields, followed by a quantitative comparison between these two methods to select the suitable method out of these two for further analysis.
Method of moments
In the method of moments, the parameters defining the distribution are estimated using the information on its moments. The parameters of the log-normal distribution, and , are determined as the first moment, , and the second moment, Var(), of logarithmic of the observation (, , , , ) of the random variable, X, respectively. The observation (, , , , ) is the value of X at different nodes along the depth obtained from each realisation of the random field. In the method of moments, the scale of fluctuation, , is generally determined using three methods, namely, the approximate method, the autocorrelation function method, and the semivariogram method. These methods are described below.
Approximate method
The approximate method is also known as the rule of thumb method for estimation of the scale of fluctuation, which was given by Fenton and Vanmarcke (1990). It assumes that the correlation function follows the squared exponential decay model. The method involves the fitting of the soil property profile into a trend line using a statistical method such as linear regression analysis. Furthermore, the average of the vertical distances, , between the intersection of the soil property curve and the fitted trend line is determined. The scale of fluctuation is then estimated as
Figure 1 illustrates the rule of thumb method to estimate the scale of fluctuation.
The diagram shows soil property xi along the horizontal axis and depth z increasing downward. A straight line indicates the trend, while a wavy curve shows deviations from this trend. At a depth z i, the horizontal distance between the curve and trend is labelled w z i. The profile is divided into five vertical segments labelled d 1, d 2, d 3, d 4, and d 5 using dashed horizontal lines. At the bottom, the mean spacing is defined as d bar equals 1 over 5 times the sum of d i.Illustration of the rule of thumb for estimation of scale of fluctuation
The diagram shows soil property xi along the horizontal axis and depth z increasing downward. A straight line indicates the trend, while a wavy curve shows deviations from this trend. At a depth z i, the horizontal distance between the curve and trend is labelled w z i. The profile is divided into five vertical segments labelled d 1, d 2, d 3, d 4, and d 5 using dashed horizontal lines. At the bottom, the mean spacing is defined as d bar equals 1 over 5 times the sum of d i.Illustration of the rule of thumb for estimation of scale of fluctuation
Autocovariance function fitting
In the autocovariance function fitting, the covariance between the random vectors and , which are subsets of the realisations of the random variable X, is determined as,
where, and are the vectors of equal size separated by the lag distance and and are the mean of and , respectively.
For a given lag distance , where , and is the distance between two successive nodes, and can be written as,
Thus, for a given lag distance can be determined as,
In the above formula, random fields are assumed to be stationary and therefore and are considered to be the same as , which is the mean of X. Furthermore, is normalised with to obtain the sample correlation coefficient . is the covariance between the random vectors with zero lag distance, which is defined as,
Therefore,
The normalised correlation coefficient is plotted against the lag distance, , to determine the appropriate correlation function and its parameter. Generally, the Morkov model is considered a correlation function for spatially varying ground. The Markov model, which is an exponentially decaying correlation function, is defined as,
where is called scale of fluctuation. The scale of fluctuation is determined by performing least square regression analysis to fit the variation of with for the random field X obtained from the theoretical function as given in Equation 14.
Semivariogram fitting
The semivariogram fitting is an alternative method to determine the correlation between the random vectors and . It is defined as,
where is semivariogram, which can be written as below for a given lag distance ,
The difference between the autocovariance and semivariogram fitting methods is that the former depends on the mean , whereas the latter does not depend on it. However, for the stationary process and V() are related as follows
The scale of fluctuation in the semivariogram fitting is determined using the regression analysis as similar to the autocovariance function fitting. For the least square regression analysis, the error for autocovariance and semivariogram fitting is defined as,
where and ) are the autocorrelation coefficients at given obtained from the theoretical function (Equation 14) and that from the simulated random field, respectively; and are the semivariogram coefficients obtained from the theoretical function (Equation 17) and that from the simulated random field, respectively. The error is minimised by performing an iteration over . The value of for which the error is minimum is considered as the scale of fluctuation for a given random field. Figures 2(a) and 2(b) show the variation of error with for autocovariance and semivariogram function fitting, respectively. Figures 3(a) and 3(b) depict the variation of with obtained from the simulated random field and from the theoretical function with the estimated for autocovariance and semivariogram function fitting, respectively. From Figure 3, it can be clearly observed that the autocovariance method fits the against relatively better than the semivariogram function fitting. The estimated value of for semivariogram function fitting is also higher than the input value of ; whereas, it is close to the input value of for the autocovariance method.
The left chart plots error versus theta in metres. Error decreases from about 1.2 to a minimum near 0.95 at theta about 1.54, then increases again. The right chart shows a similar trend, with minimum error near 0.96 at theta about 2.91. Both curves are U-shaped, and the minimum points are marked with dots and vertical dashed lines.Variation of error with respect to estimated using (a) autocovariance function fitting () and (b) semivariogram method ()
The left chart plots error versus theta in metres. Error decreases from about 1.2 to a minimum near 0.95 at theta about 1.54, then increases again. The right chart shows a similar trend, with minimum error near 0.96 at theta about 2.91. Both curves are U-shaped, and the minimum points are marked with dots and vertical dashed lines.Variation of error with respect to estimated using (a) autocovariance function fitting () and (b) semivariogram method ()
The left chart shows the correlation rho of tau versus lag tau in metres. The theoretical curve decreases smoothly to near 0, while the observed curve fluctuates around zero with small oscillations. The right chart shows the variance V of tau versus lag tau. The theoretical curve rises quickly and stabilises near 25, while the observed curve varies widely between about 10 and 45, with peaks around 16 to 18 before dropping near 0 at high lag.Comparison of curves between simulated data and theoretical function along estimated of using (a) autocovariance function fitting and (b) semivariogram method
The left chart shows the correlation rho of tau versus lag tau in metres. The theoretical curve decreases smoothly to near 0, while the observed curve fluctuates around zero with small oscillations. The right chart shows the variance V of tau versus lag tau. The theoretical curve rises quickly and stabilises near 25, while the observed curve varies widely between about 10 and 45, with peaks around 16 to 18 before dropping near 0 at high lag.Comparison of curves between simulated data and theoretical function along estimated of using (a) autocovariance function fitting and (b) semivariogram method
Method of maximum likelihood
In the method of maximum likelihood, the statistical parameters defining uncertainties are estimated by maximising the likelihood of observations under the assumed distribution. The likelihood of producing the observation ( is the vector and superscript T is the matrix transpose) given the parameters is defined as (Fenton, 1999a),
where is the covariance matrix, is the determinant of , and is the vector of means defined as ; . The expression for likelihood, as shown in Equation 19, considers that the observations follow the normal distribution. However, the simulation of the spatially varying random field is done considering that the realisations of X follow a log-normal distribution. Therefore, the data obtained from the simulation is transformed by taking the natural logarithmic to use it for the maximum likelihood estimation. Equation 19 is further simplified to estimate the parameters as shown below. Details of the mathematical derivations of the maximum likelihood method can be found in Fenton (1999a).
where is logarithm of , is the estimator of , and is the determinant of correlation matrix. is determined by setting the partial derivative of with respect to equal to zero, which yields,
where, is a vector defined as , and is estimator of , which is obtained by setting the partial derivative of with respect to equal to zero,
The estimation of requires partial derivative of (as given in Equation 20) with respect to . However, the closed-form solution does not exist for the same as differentiation of the determinant of is not possible. Therefore, an iterative approach has been used to estimate the wherein the iteration over has been done to maximise the . The value of and have been updated in each iteration, and this process has been continued till the global maximum of is obtained.
Table 1 shows the estimated using the method of moments and the method of maximum likelihood, respectively. The is the average value of estimated for all the simulated 2000 random fields with = 25 MPa, = 10, and the input = 2 m. Two cases of the soil domain with the length of 20 and 50 m are considered for estimation of . Table 1 depicts that the estimated from the approximate method for both the L equals to 20 and 50 m are much lesser than the input . The value of in this case is between 0.75 and 0.81; however, given as input to generate the random fields is equal to 2. The semivariogram yields the value of relatively higher than the input . The autocorrelation function fitting gives the value of closer to the input . From Table 1, it can be clearly observed that the method of maximum likelihood estimates the value of much closer to the input than any of the three methods of the method of moment. This indicates that the method of maximum likelihood performs better than the method of moment in estimation of the scale of fluctuation, . Therefore, the method of maximum likelihood is used for further analysis.
Estimated values of using the method of moments and method of maximum likelihood
| L: m | Estimated values of : m | |||
|---|---|---|---|---|
| Method of moments | Method of maximum likelihood | |||
| Approximate method | Autocovariance function fitting | Semivariogram fitting | ||
| 20 | 0.756 | 1.501 | 3.303 | 1.843 |
| 50 | 0.805 | 1.815 | 2.947 | 1.944 |
| L: m | Estimated values of | |||
|---|---|---|---|---|
| Method of moments | Method of maximum likelihood | |||
| Approximate method | Autocovariance function fitting | Semivariogram fitting | ||
| 20 | 0.756 | 1.501 | 3.303 | 1.843 |
| 50 | 0.805 | 1.815 | 2.947 | 1.944 |
Results and discussion
The analysis of uncertainties in the LU decomposition-based random field theory begins with performing the hypothesis testing for the following hypotheses.
The simulated random field of the soil property, X, follows the log-normal distribution.
The simulated random field of the soil property, X, follows the same distribution that is given as input to the random field theory for the realisation of the random fields.
The K–S test is performed to test the hypotheses mentioned above. The one-sample test is done for the first hypothesis, whereas the two-sample test is done for the second hypothesis. These tests are performed on the total 2000 random fields simulated with the input parameters, = 25 MPa, = 30, and = 2 m. Table 2 shows the percentage of acceptance of the hypothesis at different significance levels. From this table, it can be observed that the percentage of acceptance of the first hypothesis is in the range of 95%–99% for the given significance levels. This indicates that most of the simulated random fields follow the log-normal distribution. The result of the two-sample test shows that the percentage of acceptance of the second hypothesis is in the range of 43%–58% for the given values of significance level. This result indicates that approximately half of the total number of simulated random fields do not follow the distribution given as an input to the random field theory. Based on the results of one-sample and two-sample tests, it can be inferred that although simulated random fields follow the log-normal distribution, the statistical parameters of a major proportion of the simulated random fields are not the same as their input values. This indicates that the generated random fields are not the same as those expected to be; thereby, it results in the randomness in the statistical parameters of the simulated random fields.
Evaluation of randomness in the simulated random fields
To evaluate the randomness in the simulated random fields, a set of 2000 realisations of spatial random field are made with the input parameters: equals to 25 MPa, = 20, and = 2 m. Figures 4–6 show the histogram and best-fitted pdf of , and , respectively. Where, , , and are the mean, coefficient of variation, and scale of fluctuation estimated from the realised random fields using maximum likelihood estimation, respectively. These figures indicate that the statistical parameters (, , and ) characterising the simulated random fields are not constant as these are expected to be, instead these behave as the random variables following skewed normal distribution with their distribution parameters shown in the respective figures. The distribution parameters shown in Figures 4–6 are , , and , where is the mean, is the standard deviation, and is the skewness of the input statistical parameter, respectively. To further investigate the factors governing the randomness in output parameters, an elaborate stochastic analysis was conducted using 12 distinct sets of 2000 realisations done with varying input parameters. The scale of fluctuation was varied as 0.5, 1, 2, and 5 m, and the coefficient of variation was varied as 10%, 20%, and 30% with 25 MPa.
The chart plots f of mu sub x versus mu sub x. A histogram spans about 21 to 32, with the highest frequency near 26. A smooth curve overlays the bars. The peak occurs at about 25.89 with a spread sigma of 1.64 and a skewness of about 0.32. The distribution is approximately symmetric with slight right skew.Histogram and best-fitted probability density function (pdf) of
The chart plots f of mu sub x versus mu sub x. A histogram spans about 21 to 32, with the highest frequency near 26. A smooth curve overlays the bars. The peak occurs at about 25.89 with a spread sigma of 1.64 and a skewness of about 0.32. The distribution is approximately symmetric with slight right skew.Histogram and best-fitted probability density function (pdf) of
The chart shows f of C O V sub x percent versus C O V sub x percent. Values range from about 27 to 43, with a peak near 33.5. A smooth curve overlays the histogram. The mean is about 33.51 with a sigma of 2.59 and a skewness of about 0.19. The distribution is moderately spread and nearly symmetric.Histogram and best-fitted probability density function (pdf) of
The chart shows f of C O V sub x percent versus C O V sub x percent. Values range from about 27 to 43, with a peak near 33.5. A smooth curve overlays the histogram. The mean is about 33.51 with a sigma of 2.59 and a skewness of about 0.19. The distribution is moderately spread and nearly symmetric.Histogram and best-fitted probability density function (pdf) of
The chart plots f of theta versus theta in metres. Values range from about 0 to 5.5 with a peak near 1.8. A fitted curve overlays the histogram. The mean is about 1.84, the sigma is 0.64, and the skewness is about 1.02. The distribution shows noticeable right skew with a longer tail at higher values.Histogram and best-fitted probability density function (pdf) of
The chart plots f of theta versus theta in metres. Values range from about 0 to 5.5 with a peak near 1.8. A fitted curve overlays the histogram. The mean is about 1.84, the sigma is 0.64, and the skewness is about 1.02. The distribution shows noticeable right skew with a longer tail at higher values.Histogram and best-fitted probability density function (pdf) of
Effect of input on randomness in , , and
Figure 7 shows the best-fitted distribution of of the random fields generated with m and varying from 10% to 30%. This figure indicates that the mean and standard deviation of increases with increase in the value of . The percentage relative error () of the mean of varies from 0.88% to 8.04% as varies from 10% to 30%; where the percentage relative error () is determined as,
The chart plots f of mu sub x versus mu sub x for three cases of C O V sub x equal to 10 percent, 20 percent, and 30 percent. The 10 percent curve peaks sharply near 25.2 with a narrow spread. The 20 percent curve peaks near 25.9 with a moderate spread. The 30 percent curve peaks near 27.0 with the widest spread. Increasing C O V sub x results in broader and slightly right-shifted distributions.Effect of variation of on randomness in
The chart plots f of mu sub x versus mu sub x for three cases of C O V sub x equal to 10 percent, 20 percent, and 30 percent. The 10 percent curve peaks sharply near 25.2 with a narrow spread. The 20 percent curve peaks near 25.9 with a moderate spread. The 30 percent curve peaks near 27.0 with the widest spread. Increasing C O V sub x results in broader and slightly right-shifted distributions.Effect of variation of on randomness in
It can also be observed from Figure 7 that the standard deviation of varies from 0.8 to 2.57 with an increase in the value of from 10% to 30%. Similar observations can also be made from Figure 8, which shows that the percentage relative error in mean of varies from 65.4% to 70.7% and standard deviation of increases from 1.2 to 4.27 as the increases from 10% to 30%. Figures 7 and 8 indicate that the randomness in and increases with increase in the value of . However, Figure 9 shows that is minimally sensitive towards the variation in as in mean of varies from 9.5% to 7.5% and standard deviation of varies from 0.63 to 0.66 for varying from 10% to 30%.
The chart plots f of C O V sub x versus C O V sub x percent for three cases of C O V sub x equal to 10 percent, 20 percent, and 30 percent. The 10 percent curve peaks near 16.5 with a narrow spread. The 20 percent curve peaks near 33.5 with a moderate spread. The 30 percent curve peaks near 51.2 with the widest spread. As C O V sub x increases, the distribution shifts to higher values and becomes broader.Effect of variation of on randomness in
The chart plots f of C O V sub x versus C O V sub x percent for three cases of C O V sub x equal to 10 percent, 20 percent, and 30 percent. The 10 percent curve peaks near 16.5 with a narrow spread. The 20 percent curve peaks near 33.5 with a moderate spread. The 30 percent curve peaks near 51.2 with the widest spread. As C O V sub x increases, the distribution shifts to higher values and becomes broader.Effect of variation of on randomness in
The chart plots f of theta versus theta in metres for C O V sub x values of 10 percent, 20 percent, and 30 percent. All three curves peak near 1.8 and have very similar shapes and spreads. The mean values range from about 1.81 to 1.85, and sigma values are close to 0.63 and 0.66. The distributions almost overlap, indicating minimal influence of C O V sub x on theta distribution shape.Effect of variation of on randomness in
The chart plots f of theta versus theta in metres for C O V sub x values of 10 percent, 20 percent, and 30 percent. All three curves peak near 1.8 and have very similar shapes and spreads. The mean values range from about 1.81 to 1.85, and sigma values are close to 0.63 and 0.66. The distributions almost overlap, indicating minimal influence of C O V sub x on theta distribution shape.Effect of variation of on randomness in
Effect of input on randomness in , , and
Figure 10 depicts the best-fitted distribution of for different values of ranging from 0.5 to 5 m for the given . This figure indicates that the scatteredness in increases as increases. The value of of varies from 1.39 to 2.21 as varies from 0.5 to 5 m. However, the percentage relative error in mean of decreases from 15.28% to 0.52% with an increase in the value of from 0.5 to 5 m. Figure 11 shows the probability distribution of of the simulated random fields with input 0.5, 1, 2, and 5 m for the given . Analysis of these distributions indicates that the percentage relative error decreases from 207.3% to 8.45% and the standard deviation decreases from 5.51 to 1.59 as the value of increases from 0.5 to 5 m. It is to be noted that the coefficient of variation measures the deviation of the values of the random variable around its mean. Thus, it can be inferred from the observations made from Figure 11 that the simulated random fields at lower value of are more scattered around their mean and the scatteredness decreases with increase in the value of . Effect of increase in dispersion with decrease in the value of can be observed in the percentage relative error in the mean of and the standard deviation of . Although the uncertainties in , which are measured herein with the relative percentage error and standard deviation, are relatively high, but their effect on the randomness in are minimal. Figure 12 indicates that the accuracy as well as precision in increases with decrease in value of input . The relative percentage error in the mean of decreases from 25% to 2% and the value of decreases from 1.41 to 0.13 as the values of input varies from 5 to 0.5 m. This observation indicates that the simulated random fields are relatively more accurate at the lower value of as their scale of fluctuation () are close to the value of for which realisations are made.
The graph plots probability density against mean soil property mu x bar. Four curves correspond to theta values of 0.5 metres, 1 metre, 2 metres, and 5 metres. The curves shift left as theta increases, indicating decreasing mean values from about 28.82 to 25.13. Standard deviation remains similar for smaller theta and increases slightly for larger theta. Each curve is labelled with mu, sigma, and epsilon percent values.Best-fitted probability density functions (pdfs) of for varying at
The graph plots probability density against mean soil property mu x bar. Four curves correspond to theta values of 0.5 metres, 1 metre, 2 metres, and 5 metres. The curves shift left as theta increases, indicating decreasing mean values from about 28.82 to 25.13. Standard deviation remains similar for smaller theta and increases slightly for larger theta. Each curve is labelled with mu, sigma, and epsilon percent values.Best-fitted probability density functions (pdfs) of for varying at
The graph shows probability density versus the coefficient of variation percent. Four curves represent theta values of 0.5 metres, 1 metre, 2 metres, and 5 metres. The peak shifts from about 61.46 percent at theta 0.5 metre to about 21.69 percent at theta 5 metre. Spread decreases as theta increases, indicating reduced variability. Each curve includes values of mean, standard deviation, and epsilon percent.Best-fitted probability density functions (pdfs) of for varying at
The graph shows probability density versus the coefficient of variation percent. Four curves represent theta values of 0.5 metres, 1 metre, 2 metres, and 5 metres. The peak shifts from about 61.46 percent at theta 0.5 metre to about 21.69 percent at theta 5 metre. Spread decreases as theta increases, indicating reduced variability. Each curve includes values of mean, standard deviation, and epsilon percent.Best-fitted probability density functions (pdfs) of for varying at
The graph plots probability density against theta zero in metres. Four curves correspond to theta values of 0.5 metres, 1 metre, 2 metres, and 5 metres. The peaks move from about 0.49 to 3.75 as theta increases. The distributions widen significantly for larger theta, indicating increased spread. Each curve includes labelled mean, standard deviation, and epsilon percent values.Best-fitted probability density functions (pdfs) of for varying at
The graph plots probability density against theta zero in metres. Four curves correspond to theta values of 0.5 metres, 1 metre, 2 metres, and 5 metres. The peaks move from about 0.49 to 3.75 as theta increases. The distributions widen significantly for larger theta, indicating increased spread. Each curve includes labelled mean, standard deviation, and epsilon percent values.Best-fitted probability density functions (pdfs) of for varying at
Effect of input on randomness of
In the above section, effect of input value of on the randomness in is analysed. However, the normalised value of with respect to the depth of the domain (L), defined as , can be considered to be a more appropriate parameter as compared with . Therefore, effect of on the parameters defining randomness in is further investigated. The values of ranging from 0.01 to 0.2 are considered for this analysis. This range has been chosen based on the practical values of reported in the literature (Fei et al., 2022; Onyejekwe et al., 2016). Two cases are considered to obtain the varying values of . In the first case, is varied from 0.5 to 2 m, and L is kept constant; and in the second case, is kept constant, and L is varied from 10 to 50 m. Figures 13 and 14 depict the effect of on the percentage relative error in mean of and the standard deviation of for the first and the second cases, respectively. From these figures, it can be clearly observed that the effect of on the variation of parameters defining randomness is not significant for both the cases. Figures 13 and 14 indicate that the randomness in is very small as the percentage relative error is found to be around 5% and standard deviation is around 0.2 for , which are increased to around 18% and 0.80, respectively, as is increased to 0.2. It can be therefore inferred from these observations that the randomness in is less for the lower value of and it increases with increase in the values of . Figures 15(a) and 15(b) present an empirical relationship for the variation of the percentage relative error in the mean of and the standard deviation of with , respectively, which is obtained by performing least square regression analysis on the combined results of both the cases. These empirical relationships are useful in determining an optimum value of for a given domain to obtain the expected random fields using LU decomposition-based random field theory by setting a tolerable limit for and . In addition to this, these relationships also yield the amount of uncertainty in the random fields simulated for an input value of .
The left chart plots epsilon in mean of theta percent against theta over L with values 0.025, 0.05, 0.1, and 0.25. Three lines for C O V x 10 percent, 20 percent, and 30 percent all increase, reaching about 24 to 25 percent at 0.25. The right chart shows sigma of theta percent versus theta over L with the same x values. All three lines rise steadily from about 0.14 to about 1.45, with minimal difference between C O V x cases.The effect of for case 1, (, is varying and L is kept constant (L = 20)) on (a) in mean of and (b) in mean of
The left chart plots epsilon in mean of theta percent against theta over L with values 0.025, 0.05, 0.1, and 0.25. Three lines for C O V x 10 percent, 20 percent, and 30 percent all increase, reaching about 24 to 25 percent at 0.25. The right chart shows sigma of theta percent versus theta over L with the same x values. All three lines rise steadily from about 0.14 to about 1.45, with minimal difference between C O V x cases.The effect of for case 1, (, is varying and L is kept constant (L = 20)) on (a) in mean of and (b) in mean of
The left chart shows epsilon in the mean of theta percent plotted against theta over L. Points for lengths 8, 10, 20, 30, 40, and 50 metres increase from about 3 percent at 0.025 to about 20 percent at 0.25. Larger lengths correspond to smaller epsilon at low theta over L. The right chart shows the sigma of theta percent versus theta over L. Values increase from about 0.42 to about 0.92 as theta over L increases, with shorter lengths generally giving higher sigma.The effect of for case 2, (L is varying and , are kept constant (, )) on (a) in mean of and (b) in mean of
The left chart shows epsilon in the mean of theta percent plotted against theta over L. Points for lengths 8, 10, 20, 30, 40, and 50 metres increase from about 3 percent at 0.025 to about 20 percent at 0.25. Larger lengths correspond to smaller epsilon at low theta over L. The right chart shows the sigma of theta percent versus theta over L. Values increase from about 0.42 to about 0.92 as theta over L increases, with shorter lengths generally giving higher sigma.The effect of for case 2, (L is varying and , are kept constant (, )) on (a) in mean of and (b) in mean of
The left plot shows epsilon in mean of theta percent against theta over L with scattered points and a fitted line labelled y equals 84.07 x plus 0.55. Values increase linearly from about 3 to about 22 percent. The right plot shows sigma of theta percent versus theta over L with a fitted line labelled y equals 2.13 x plus 0.38. Sigma increases from about 0.42 to about 0.92, showing a clear linear relationship.Empirical relationship derived using linear regression analysis for (a) in mean of against and (b) in mean of against
The left plot shows epsilon in mean of theta percent against theta over L with scattered points and a fitted line labelled y equals 84.07 x plus 0.55. Values increase linearly from about 3 to about 22 percent. The right plot shows sigma of theta percent versus theta over L with a fitted line labelled y equals 2.13 x plus 0.38. Sigma increases from about 0.42 to about 0.92, showing a clear linear relationship.Empirical relationship derived using linear regression analysis for (a) in mean of against and (b) in mean of against
Evaluation of stationary assumption
LU decomposition-based random field theory implicitly considers that the simulated spatially varying random fields are stationary. The stationarity assumption infers that the joint probability distribution function does not change with respect to space. In a simple way, it emphasises that the mean and variance should remain constant throughout the space. To evaluate this assumption, the soil domain is divided into four equal and successive subdomains, each having a depth of 5 m. Furthermore, the mean and the variance are estimated for each subdomain, and their values are shown at mid point of the subdomains in Figures 16(a) and 16(b). These figures show the variation of the mean and the variance along the depth estimated for one of the realisations, respectively. It can be clearly observed from Figure 16 that the values of the mean and the variance are not the same for all the subdomains, thereby indicating the non-stationary behaviour of soil property X along the spatial dimension. For detailed analysis, the subdomain mean values of X are fitted into a linear trend line using regression analysis as
where is the intercept on axis and is the slope of the trend line of the variation of along the depth.
The left plot shows mean mu in kilopascal versus depth from 0 to 20 metres. Data points increase with depth from about 22 to about 26 kilopascal, with a straight trend line. The right plot shows sigma in kilopascal versus depth. Values increase from about 6 to about 8.5 kilopascals with depth. Both plots display scattered measurements aligned with upward linear trends.Variation and linear regression fitting of (a) mean values of X and (b) standard deviation values of X along each sub-domain
The left plot shows mean mu in kilopascal versus depth from 0 to 20 metres. Data points increase with depth from about 22 to about 26 kilopascal, with a straight trend line. The right plot shows sigma in kilopascal versus depth. Values increase from about 6 to about 8.5 kilopascals with depth. Both plots display scattered measurements aligned with upward linear trends.Variation and linear regression fitting of (a) mean values of X and (b) standard deviation values of X along each sub-domain
For a stationary field, should be equal to the input value of mean and should be equal to zero. However, likelihood of single point prediction of or to analyse the stationarity assumption is practically zero in such random scenario. Therefore, a reasonable range of slope of trend line ranging from to , for which the value of ranges from to 0.1, around its mean is assumed to be the stationary field condition for statistical evaluation of stationarity assumption. Figure 17 shows the histogram of and its best-fitted probability distribution for the random field simulated for = 2 m, , and MPa. A hatched portion is also shown in this figure to depict the area corresponding to the assumed range of for stationarity assumption. Thus, this area yields likelihood of the simulated random fields to be stationary, which is termed as . Figure 18 shows the variation of with . This figure shows that the probability to obtain a stationary random field is reduced on increasing the value of following the exponentially decaying trend. Hence, this observation indicates that the likelihood of the simulated random fields to be stationary is more at the lower value of .
The chart shows the probability density p of a 1 versus the value of slope a 1. A histogram spans about negative 1.2 to positive 1.2, with the highest bars near 0. A smooth curve overlays the bars and peaks at 0, indicating a symmetric normal distribution. The central bin around 0 is highlighted, and frequencies decrease evenly towards both tails.Statistical variation of the regression parameter, (assumed reasonable range of stationary condition shown by hatched region)
The chart shows the probability density p of a 1 versus the value of slope a 1. A histogram spans about negative 1.2 to positive 1.2, with the highest bars near 0. A smooth curve overlays the bars and peaks at 0, indicating a symmetric normal distribution. The central bin around 0 is highlighted, and frequencies decrease evenly towards both tails.Statistical variation of the regression parameter, (assumed reasonable range of stationary condition shown by hatched region)
The chart plots p of S against theta over L. The curve starts near 0.38 at low theta over L and decreases sharply to about 0.2 by 0.1. It then gradually declines and stabilises around 0.135 at 0.5. The trend shows a rapid initial drop followed by a slight decrease.Variation of probability of stationary random field, P(S) with
The chart plots p of S against theta over L. The curve starts near 0.38 at low theta over L and decreases sharply to about 0.2 by 0.1. It then gradually declines and stabilises around 0.135 at 0.5. The trend shows a rapid initial drop followed by a slight decrease.Variation of probability of stationary random field, P(S) with
Conclusions
A comprehensive statistical analysis to evaluate the uncertainties associated with LU decomposition-based random field theory is done in this paper. The first part of the study analyses the suitable method to characterise the spatially varying ground, and the later part evaluate the uncertainties associated in modelling the random spatial variable ground using LU decomposition-based random field theory. The results obtained in this study indicate that the method of maximum likelihood performs better than the method of moments in estimation of the scale of fluctuation, . The one-way and two-way K–S hypothesis test performed on the simulated random fields indicate that the LU decomposition-based random field theory does not yield the random fields as these are expected to be for given input parameters. The mean, variance, and the scale of fluctuation estimated from the simulated random fields are different from those given as the input parameters to generate the random fields. The results obtained in this study indicate that the statistical parameters estimated from the simulated random fields behave as a random variable. It is observed that the effect of on the randomness in and are significant; however, it is minimal in . The randomness in is less for the lower value of and it increases with increase in the value of . The likelihood to obtain a stationary random field is more at the lower value of , and it reduces following exponentially decaying trend as increases. The results obtained from this study clearly highlight the limitations of LU decomposition-based random field theory, which is commonly used in geotechnical engineering to generate spatially variable random fields for risk assessment analyses for ground spatial variability. The findings of this study can be utilised to select the optimum values of statistical parameters, thereby reducing uncertainties and enabling the simulation of more accurate random fields.
Acknowledgements
Financial support from Shiv Nadar Institution of Eminence, Delhi-NCR, India, is gratefully acknowledged. Any opinions, findings, and conclusions or recommendations expressed in this materials are these of authors and do not necessarily reflect the views of Shiv Nadar Institution of Eminence, Delhi-NCR.

