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.
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.
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.
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.
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,
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%.
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.
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 .
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.
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 .
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.



















