Spatial variability of medium dense sand deposit

: This paper presents a study of the spatial variability of the cone resistance in a medium dense sand deposit in Norway. Spatial variability studies have been done earlier for clay deposits, but there are few results available for sands. To bridge this gap and to compare the spatial variability of sand and clay, the in situ cone resistance measured with the piezocone at the Øysand benchmark site near Trondheim was studied. Corrected cone resistances q t derived from the measured cone resistance and pore pressure measurements were used to establish the autocorrelation structure of the Øysand deposit. The q t data were examined visually and conclusions were drawn on the data that should be analyzed. A depth interval of 7 to 12 m was selected for spatial variability further analysis. The scale of fluctuation was identified. Several autocorrelation functions were considered, and the single exponential function was found to be the one offering the best autocorrelation. The spatial variability in the vertical and horizontal directions was analyzed statistically, using three different approaches, the auto-correlation fitting, maximum likelihood estimation and simplified Vanmarcke method. The results indicate short autocorrelation distances of 3 m or less in the horizontal direction, suggesting a very variable sand at the Øysand site. These horizontal autocorrelation distances are much shorter than those obtained for clays. In the vertical direction the scale of fluctuation was less than one meter, as obtained for other soils. This paper presented a study of the spatial variability of the cone resistance in the medium dense Øysand sand in Norway. The corrected cone resistance q t derived from the cone resistance and pore pressure measurements was used to establish the autocorrelation structure. The q t data were first analyzed visually. A depth interval of 7 to 12 m was selected for further analysis. To study spatial variability numerically the scale of fluctuation was identified. Several autocorrelations functions were considered, and the single exponential function was found to be the one offering the best autocorrelation. The spatial variability in the vertical and horizontal directions was analyzed statistically, using three different approaches, i.e. auto-correlation fitting, maximum likelihood estimation and simplified Vanmarcke method. The results indicate short autocorrelation distances of 3.5 m or less in the horizontal direction, suggesting a very variable sand at the Øysand site. In the vertical direction the scale of fluctuation was less than one meter, as obtained for other soils. The CPTU data at the Øysand presented some very interesting spatial variability aspects. The study of vertical and horizontal spatial variability of the corrected cone resistance q t led to the following observations.


Introduction
The engineering properties of a soil always exhibit variation from point to point within a soil volume. It is of interest to quantify the spatial variability in order to assess its importance, help plan future soil investigations and, if possible, reduce the uncertainty in a property by spatial averaging. This paper quantifies the spatial variability of the cone resistance in a sand deposit. The spatial variability of clays has been studied earlier. However, little information exists on the spatial variability of sands and silts.
The piezocone test has become one of the most frequently used methods to obtain the characteristics of a soil in situ. Piezocone tests in sand, even when closely spaced, can show large variability. Quantifying the distance over which the soil parameters in a sand are correlated can help plan site investigations and reduce the uncertainty in the soil parameters derived from the piezocone test. Reducing the uncertainty can have a significant impact on the reliability of geotechnical structures (e.g. [1][2][3][4]). In this paper, the variability in the corrected cone resistance is analyzed through first a visual inspection of the data, and second a calculation of the spatial correlation structure.
The spatial correlation structure of a geotechnical property is characterized by random field modelling [5]. The scale of fluctuation (SOF), describing the random field, defines the distance over which there is a significant correlation of a geotechnical parameter. It is recognized that soil properties are anisotropic by nature and the horizontal and vertical SOFs will probably differ (e.g. [1,5,6]). However, 3D spatial variability characterization remains a challenging task, due to limited site-specific investigation data, especially in the horizontal direction, and the lack of suitable characterization method [4,7,8].
The Norwegian GeoTest Sites (NGTS) project established five research sites in Norway in 2016. The sites are referred to as the sand, soft clay, quick clay, silt and permafrost sites. The establishment of the sites was funded by the Research Council of Norway and the aim of the project is to establish, characterize, share digital data and manage the test sites over the next 20 years [9]. Each of the sites have been carefully and extensively characterized to offer well documented data sets in five typical geo-materials for further testing, physical model testing, constitutive model development and modelling.
The objective of this study is to do a characterization of the inherent spatial variability of the sand in Øysand 1 . The results of piezocone tests (CPTU, but also denoted CPT in some of the figures and tables in this paper) from the Øysand benchmark site, near Trondheim in Norway, were utilized to establish a correlation model. The paper describes the soil profile at the test site (Section 2) and provides an overview of the results of the cone penetration tests selected for spatial variability analysis (Sections 3 and 4). In Sections 3 and 4, no calculation of the spatial variability is done. Section 5 summarizes briefly concepts from random field theory and the analysis approaches used in the study. The result of the analyses for Øysand sand in terms of scale of fluctuations in the vertical and horizontal directions are compared and discussed in Section 6. Conclusions are drawn in Section 7.

The Øysand site
The Øysand research site is located in central Norway, approximately 15 km south of Trondheim. The locality sits on the south side of the Gaula River, at the head of the Gaulosen, an arm off the main Trondheim Fjord. Over the past thousand years, the river has mostly prograded westwards in the fjord. The ground surface at the site is at an elevation of 2.7 m above mean sea level.
The fluvial and deltaic deposit at Øysand consists of a 20-25 m fine silty sand with occasional high gravel content [10,11]. Figure 1 (top) presents a schematic longitudinal cross-section of a deltaic deposit, depicting its characteristic tripartite architecture (topset, foreset and bottomset) [12]. At Øysand, the stratigraphy features a general coarsening upward sequence as typically observed in deltaic deposits with topset, foreset and bottomset units (Figure 1 (middle)). The layers in these units can have different properties, geometry, fall and dip that can be linked to the depositional history at the mount of the river delta. Figure 1 (bottom) illustrates the layering at Øysand, based on the results of two resistivity cone tests [10]. Figure 2 presents a snapshot of the stratigraphy and index properties of the soils at the site, as obtained from in situ and laboratory tests 2 . The borehole log is for Borehole OYSB09, which is located very close to CPTU OYSC09. Figures 3 and 4 illustrate two of the piezocone tests on the soil behavior classification charts [13], as interpreted from the pore pressure ratio B q , and from the friction ratio F r [14]. Figure 3, with data from CPTU OYSC42 (called CPT8 later in the spatial variability analysis), indicates that for this borehole, there is a fairly clean sand between 4 and about 10 m. Below 10 or 11m, the sand becomes more silty. Figure 4, with data from CPTU OYSC40 (called CPT6 later in the spatial variability analysis), indicates a clean sand to silty sand between 5 and 20 m. Note that the data points are so close to each other that the colors seem darker than in reality. Many of the CPT showed a change in gradient in the CPTU measured values from depth 13 or 14 m. The chart indicates a clean sand to silty sand. As the depth increases, the sand contains more silt and clay.
Based on all the boreholes and CPTUs (NGI [14,15]), a "more uniform" layer of clean sand to lightly silty sand was found between depths of 7 and 12 m was identified. Above this unit, the fine to coarse sand contained some gravel. Below 12 m, the sand contained an increasing portion of silty sand and sandy silt with traces of organic material. For this reason, and because the study was to address specifically medium dense sands, the spatial variability analyses in this paper concentrated on the sand between depths 7 to 12 m. This depth interval is shown with the zone between two pink horizontal lines in Figure 2. 2 The symbols used in Figure 1 are defined either in the figure text above each profile or at the bottom of the table. Other symbols not directly defined are γ s , the density of solids and D 10 and D 60 , the particle diameter (in mm) on the grain size distribution curve with 10 and 60% of the particles by weight passing. Figure 1. Schematic-longitudinal-cross-section-of-a-Gilbert-type-delta (top) [12]; typical stratigraphy at the Øysand research site (middle, not to scale); results of two resistivity cone soundings (ERT01 and ERT05) at the Øysand site (bottom) [10]. Figure 5 (left) illustrates the locations of the piezocone tests at the Øysand site. There are a large quantity of tests in a cluster. These are shown in more detail in Figure 5 (right). The original testing plan was to have a minimum distance of 1.5 m to 2 m between two adjacent tests, but several tests were run at closer intervals.  The CPTU tests in or close to the cluster were used for the spatial variability analysis. Figure 5 (right) shows the 13 piezocone tests considered in this study: the piezocones from the cluster plus three CPTU tests 10 to 15 m away. For the spatial variability analyses, the tests were renumbered CPT1 to CPT13 ( Figure 5 (right) and as will be listed in the second column of Table 1). All the tests except piezocone tests OYSC50, OYSC51 and OYSC52 (CPT11, 12 and 13) were run within an area of 18 m by 15 m (270 m 2 ). Table 1 gives for each test the cone type, the zero drifts on the cone resistance q c measured, the sleeve friction f s and the pore pressure just behind the cone u 2 , the measurement interval in mm and the air temperature during testing. The air temperature was taken from meteorological records for all tests except tests with cone types 7 and 8 for which the air temperature was measured on site [16].

Overview of cone penetration tests and corrected cone resistance q t
The selection of the CPTU tests was based on proximity and tests run under as similar conditions as possible. The air temperatures were between 16 and 18 °C, except for three tests where the temperature was 10 °C. To avoid introducing a statistical uncertainty, only the CPTUs with 10 and 20 mm measurement interval were used to evaluate the scale of fluctuation. Several types of cones were used during the characterization of the sand at the Øysand site. Tests with the seismic cone were excluded. The selected CPTU tests for spatial variability analysis were run with the cone types 1, 2, 6 and 7. Tests with cone type 1 were performed approximately 1.5 m away from other boreholes. The tests with cone type 7 were performed 2 m away from other boreholes.   1) Original number of CPTU, in site reports and in Figure 5, left.
2) Numbering of 13 CPTU selected for spatial variability analysis. Numbering, in Figure 5, right, is used in rest of paper. 3) Representative air temperature used to correct measured results. 4) It is assumed that the zero shift occurred at the end of test when the probe hit the harder layer. Table 2 lists the characteristics of the piezocone at the Øysand site, including the zero drifts, the measurement intervals and some of the observations made. Tests CPT 1 to 13 are those that are used in the spatial variability analysis. As mentioned, only the sand between depths 7 and 12 m was included in the spatial variability analysis.

Visual inspection of corrected cone resistance q t
Before the spatial variability analyses were started, the cone resistance data were looked at in detail. This section examines the cone data without any calculation of the spatial structure. The corrected cone resistance q t was judged to be a more representative parameter to compare the CPTU results than the measured cone resistance, because the derived q t accounts for the pore pressure generated during cone penetration and acting behind the cone tip. Figure 6 presents the corrected cone resistance, q t , versus depth for each of the selected 13 CPTU tests and Figure 7 superposes the 13 tests in a single graph. There is very wide variability in the corrected cone resistance derived from the cone penetration measurements and the measured pore pressure u 2 , as q t varies between 1.5 and 4.5 MPa or more in the depth interval 7 to 12 m. This is not unexpected and links to the depositional history at the site. The layers found in the delta foreset ( Figure 1) can have varying geometry and varying content of clay, silt and sand.
While examining the data, it was observed that the 13 CPTU seemed slightly out of phase depth-wise because of the structure (strike and dip) in the foreset of the deltaic deposits. A simple depth adjustment was made by shifting the corrected cone resistance q t data in each CPTU sounding up or down in order to match the peaks and troughs in the corrected cone resistances. Three figures illustrate this depth-wise phase difference and adjustment: Figure 8, Cone resistance q t versus depth once they are "in phase". Figure 9, Comparison of q t -values for 13 CPTU brought in phase depth-wise. Figure 10, Polynomial surface of depth adjustments bringing the 13 CPTU's "in phase". A comparison of Figures 7 and 9 shows how much better the data from the 13 CPTU tests fit together once the depths have been adjusted. Table 3 lists the coordinates and the ground elevation at the location of each of the CPTs in the spatial variability analysis, as well the depths from the ground to the peak values of q t for each CPT. These depths were used to do the depth-wise adjustment of the cone resistance values. Figure 10 presents an interpolation of the relative depths. The three-dimensional figure shows the relative easting and northing coordinates on the horizontal plane and the depth to the peak value of the corrected cone resistance (in m). The bi-harmonic spline interpolation in Matlab was used for the interpolation of the depth to the peak q t value based on the known location of peak values in each CPT. One could also have done such interpolation analysis with kriging techniques, e.g. [17], and obtain an even better spatial representation pf the spatial representation of the variation of the depth-wise adjustment that had to be made.
The spatial variability analyses will be done with the corrected cone resistance q t from the measured cone resistance and pore pressure vs depth values (Figure 6), and for the corrected cone resistance q t adjusted depth-wise ( Figure 8).

Random field theory
The uncertainty in soil properties estimates results from many different sources associated with data scatter and systematic errors [1,18]. The spatial correlation of a soil property, g(z), can be modelled as the sum of a deterministic trend component, t(z), and a randomly varying fluctuation term w(z) [2,19]: Vanmarcke [19] proposed that the spatial variability can be modeled as a random field. Among random field models, stationary random fields are widely used due to their simplicity and because stationary random fields are possibly the only practical type of random field that can be characterized from limited data [21]. The removal of trends is the most common means for achieving stationary fields [2]. Linearly de-trended CPT data were used to characterize spatial variability. The scale of fluctuation δ is used to quantify the spatial variability. Small values of the scale of fluctuation δ indicate rapid fluctuations about the mean trend and large values of δ are indicative of long-distance correlations. The scale of fluctuation is often denoted as SOF. Details on the approach and examples of its use for cone resistance profiles in mixed materials can be found in [3] and [20].

Auto-correlation functions
Four autocorrelation functions are commonly used in geotechnical applications: the single-exponential (SNX), the cosine exponential (CSX), the second-order Markov (SOM), and the squared exponential (SQX) models [19,21,22]. Table 4 lists the formulations of the autocorrelation models and the relations to scale of fluctuation and Figure 11 illustrates the four autocorrelation functions as a function of normalized separation distance. In Figure 11, ∆h is the horizontal separation distance and δ h is the scale of fluctuation in the horizontal direction. The functions in Figure 11 illustrate that for two CPTUs with a horizontal separation distance of 2δ h , the auto-correlations are small. Table 4. Autocorrelation functions (after [21]).  The four function models were tested with the Øysand CPTU data, and the single exponential function (SNX) was found to give the most consistent results. All four functions gave good convergence with the ACF method (Section 5.3). With the MLE method (Section 5.3), the SNX model showed the best convergence. In the rest of the paper, the results with the SNX model are shown.

Earlier results of spatial variability
Tables 5 and 6 list the scales of fluctuation obtained earlier. Table 5 lists earlier spatial variability results per soil type, Table 6 list earlier results per soil property and testing method. There are far fewer evaluations of horizontal correlation distance than evaluations of vertical correlation distance, because of the lack of data horizontally. The earlier data indicate that the horizontal autocorrelation distance is much shorter for sand than for clay. The vertical autocorrelation distance is usually shorter than the horizontal value. This reflects the dependence of the shear strength on the effective stress (or thickness of the overburden). Offshore clay Horizontal 317 [40] Offshore clay Horizontal 4500 [6] Clay Vertical 0.35-0.6 [41]  * s u = undrained shear strength; q c = cone resistance; q t = corrected cone tip resistance; q c1N = dimensionless, stress−normalised cone resistance; f s = sleeve friction; F R = stress−normalised friction ratio; I c = CPT soil behaviour classification index; N = SPT blow count; w = water content; w L = liquid limit; γ´ = submerged unit weight; γ = unit weight; e = void ratio; σ´p = preconsolidation stress; K S = subgrade modulus; D R = relative density; n = porosity. ** VST = vane shear testing; CPT = cone penetration testing; CPTU = piezocone testing; SPT = standard penetration testing; PLT = plate load testing.

Analysis methods used for Øysand CPTU data
Several statistical techniques are available in the geotechnical literature for the estimation of random field parameters in spatial variability characterization, such as autocorrelation model fitting (AMF) (e.g. [1,7]) and maximum likelihood estimation (MLE) [8,42]. Three approaches were used to characterize the vertical and horizontal variability of the cone resistance q t of Øysand sand.

Auto-correlation model fitting (AMF) method for vertical SOF characterization:
The basic principle of the autocorrelation function (ACF) approach is to fit a plausible theoretical autocorrelation function and evaluate the scale of fluctuation based on the parameters in the model (Table 4). There are two key parameters required to characterize a zero-mean stationary random field model using the ACF method, namely (1) the scale of fluctuation (δ) and (2) the auto-correlation function (ACF) ( Table 4).
Maximum Likelihood Estimation (MLE) method for vertical and horizontal SOF characterization: The Maximum Likelihood Estimation (MLE) method has the objective to quantify the soil parameters by maximize the likelihood function. However, this method has limitations computationally when large number of data are to be analyzed, because a very large matrix needs to be solved. There is on-going research to improve the computation efficiency of the calculation process of MLE [8,41].
Simplified Vanmarcke method for horizontal SOF characterization: Vanmarcke [19] proposed an approximate relationship between the scale of fluctuation and the average distances between the intersections of the soil property curves and the trend line. This method is very simple to estimate the scale of fluctuation and is widely used for the reliability analysis of geo-structures [3].

Results of spatial variability analyses
The spatial variability analyses were run in the vertical and horizontal directions, and for two sets of data: the corrected cone resistance q t inferred from the CPTU measurements of cone resistance and pore pressure, and the corrected cone resistance q t after adjusting the depths to phase the data depth-wise.  Table 7). The coefficient of variation 3 (COV) was 22%. For the fitting of the ACF function, the least square method was used on the linearly detrended data.   Table 7 presents the results of the spatial variability analyses. The horizontal scale of fluctuation is between 1.9 and 2.3 m, indicating high variability of the soil over short distance. The COV of the horizontal scale of fluctuation was 42 to 50%, which is very high. The values of scale of fluctuation were very close in the northing and easting directions. It is expected that there is link between spatial variability and the direction of delta progradation which mostly was towards the west [43].

Corrected cone resistance q t inferred directly from measurements
The vertical scale of fluctuation was close to 0.5 m for the AMF method, but as much as 1.6 m for the maximum likelihood evaluation (MLE) method. The MLE method tended to give higher values than the other methods, and was not always stable, when the matrices involved in the calculations became very large.  Figure 13. Example of the horizontal SOF at 8m depth below ground surface, using simplified Vanmarcke method (circles indicate relative location of the CPT's in the Northing direction).  Table 8 presents the results of the spatial variability analyses for the phased data over an adjusted "common" relative depth of 4 meters ( Figure 9). The average vertical scale of fluctuation, δ v , from the 13 CPTU tests from the AMF method was 0.27 m with a standard deviation (SD) of 0.03 m ( Table 8). The coefficient of variation (COV) was only 11%. The reduction in the COV is a direct result of the phasing with depth of the observed cone resistance data. For the fitting of the ACF function, the least square method was used on the linearly detrended data. The horizontal scale of fluctuation is between 1.9 and 3.2 m, indicating high variability of the soil over short distance. The COV of the horizontal scale of fluctuation was 32% in the northing direction and as high as 83% in the easting direction. The value of scale of fluctuation varied by more than one meter in the northing and easting directions. For this dataset, there is a link between spatial variability and the direction of delta progradation, with less autocorrelation in the direction of delta progradation and more correlation perpendicular to the progradation.

Corrected cone resistance q t after phasing the data depth-wise
The vertical scale of fluctuation was 0.8 m for the maximum likelihood evaluation (MLE) method.

Discussion of results
The CPTU data at the Øysand presented some very interesting spatial variability aspects. The study of vertical and horizontal spatial variability of the corrected cone resistance q t leads to the following observations. Only Table 8 with the depth-adjusted CPT-values was considered as it represents a more reliable result for scale of fluctuation in the sand layer than Table 7: -The vertical scale of fluctuation in Øysand sand was between 0.3 and 0.8 m (or 1 m rounded off for practice). The vertical scale of fluctuation for Øysand is similar to values obtained earlier for both sand and clay (Tables 5 and 6). -The horizontal scale of fluctuation was very short, between 1.9 and 3.2 m (depending on the direction), indicating short autocorrelation distances and high variability of the Øysand sand deposit. These values can be rounded off to 2 to 3 m in practice. Very high standard deviation of the horizontal scale of fluctuation (40 to 50%) were obtained, thus reinforcing the indication of variability of the Øysand sand. Such variability is typical for deltaic deposits, especially in the foreset unit where the layers are dipping at steep angles (Figure1). -Based on the results of the analyses and the calculation time it took to obtain the results, the autocorrelation function fitting (ACF) and the Simplified Vanmarcke methods seem to be more reliable than the maximum likelihood estimate (MLE). However, the MLE approach has a stronger and more fundamental mathematical basis. The challenge lies in how to solve large matrix properly and efficiently, not with the method itself. -Considering only 10 CPTU values in the spatial variability analysis in the horizontal direction (without the three outer CPTU's (CPT 11, CPT12 and CPT 13 in Figure 5 right) did not change the results obtained. -Neglecting doing the depth-wise "phasing" leads to a larger uncertainty in the calculated scale of fluctuation. To do reliable analyses of spatial variability, one needs to have access to large datasets comprising closely spaced data [44]. Abundance of data is rarely the case in practice. The paucity of data represents a major limitation for the use of the spatial variability statistical tools in practice.

Conclusion
The scale of fluctuation SOF is a convenient measure for describing the spatial variability of a soil property in a random field. It is a measure of the distance within which points are significantly correlated. Points separated by a larger distance than SOF will show little correlation, and practically no correlation will be observed when points are separated by a significantly larger distance than SOF. This paper presented a study of the spatial variability of the cone resistance in the medium dense Øysand sand in Norway. The corrected cone resistance q t derived from the cone resistance and pore pressure measurements was used to establish the autocorrelation structure. The q t data were first analyzed visually. A depth interval of 7 to 12 m was selected for further analysis. To study spatial variability numerically the scale of fluctuation was identified. Several autocorrelations functions were considered, and the single exponential function was found to be the one offering the best autocorrelation. The spatial variability in the vertical and horizontal directions was analyzed statistically, using three different approaches, i.e. auto-correlation fitting, maximum likelihood estimation and simplified Vanmarcke method. The results indicate short autocorrelation distances of 3.5 m or less in the horizontal direction, suggesting a very variable sand at the Øysand site. In the vertical direction the scale of fluctuation was less than one meter, as obtained for other soils.
The CPTU data at the Øysand presented some very interesting spatial variability aspects. The study of vertical and horizontal spatial variability of the corrected cone resistance q t led to the following observations.
-It was possible to establish a pattern of the corrected cone resistance with depth where the peaks and troughs of each CPTU could be "brought into phase" with a small adjustment of the depths (as illustrated by the surface in Figure 10). That aspect should probably be also considered at other sites with complex sand profiles where the data show a wide range of cone resistances (as in Figure 7). Neglecting to do such phasing will result in even greater spatial variability (and larger uncertainty). Such phasing further highlights the importance of understanding the geological history at a given site when assessing the spatial variability in soil properties. The results indicate very short autocorrelation distances, and high variability in the Øysand deposit. Very high standard deviations of the horizontal scale of fluctuation (COV of 40 to 50%) were obtained, reinforcing the variability of the measured values at the site. It is expected that such variability is typical for sands found in the foresets of fjord deltas. -The analysis results also illustrate the importance of understanding the geological setting, processes and history of a sand deposit to help interpret the result of the spatial variability analyses. -Based on the results of the analyses and the calculation time required to obtain the results, the ACF and the Simplified Vanmarcke methods seem to be more reliable at this time than the maximum likelihood estimate (MLE). However, improvement to this approach, mainly solving large matrices in a reasonable time, should gradually lead to good estimates. To study spatial variability, one needs large datasets, rather evenly and closely spaced. This is rarely the case in practice. This is a major limitation for the use of the elegant statistical approaches. It was not possible from this study to establish general trends that are applicable to other sand sites, except perhaps for the vertical scale of fluctuation one meter or less, which is a similar result to that obtained for other soils earlier.