Quantitative Evaluation of the Spatial Variation of Surface Soil Properties in a Typical Alluvial Plain of the Lower Yellow River Using Classical Statistics, Geostatistics and Single Fractal and Multifractal Methods

Featured Application: Spatial variability has always been a hot topic in soil science. Through an example of calculation, we compare and analyse the results of four methods: classical statistics, geostatistics, fractal and multifractal. Given this research scale, different methods have different conclusions. Considering that classical statistics, geostatistics and the fractal method have different disadvantages in application, we think that the multifractal method has broad prospects in Abstract: The spatial variability of soil properties has always been a signiﬁcant research ﬁeld in geoscience. The types of soil properties cover a wide range, but most studies have focused on the spatial variability of soil physicochemical properties over the past decades. Studies on soil hydraulic characteristics are limited, and most of them are limited to the farmland scale. However, the spatial variability of regional soil properties (soil texture and hydraulic properties) is valuable for the study of sedimentation processes and soil water transport. Therefore, here, the spatial variation of six soil properties (sand, silt, clay content, bulk density, saturated water content and saturated hydraulic conductivity) in the typical alluvial plain area of the lower Yellow River is quantitatively studied, by using classical statistics, geostatistics and single fractal and multifractal methods. This study mainly quantitatively analysed the spatial variability of di ﬀ erent soil properties and compared four research methods. Although the coe ﬃ cient of variation, nugget coe ﬃ cient, single fractal dimension and multifractal spectral width can reﬂect spatial variability, diverse conclusions are drawn (on variability) if di ﬀ erent methods are used, and the di ﬀ erent soil properties show large disparities. These four methods show a di ﬀ erent variation order of soil properties, but there are some common conclusions based on analysis and judgment. In general, the silt content in the study area is stable, mainly originating from loess transported by Yellow River erosion, which is also reﬂected in the Kriging interpolation maps under the geostatistical models. The variation in bulk density and saturated water content is weak, and the spatial variability of sand and clay content is moderate. In addition, the saturated hydraulic conductivity ﬂuctuates violently. This may be related to the di ﬀ erences in local topography, human activity and the content of sand and clay, each of which signiﬁcantly a ﬀ ects the saturated hydraulic conductivity. Classical statistics has a limitation because it fails to corelate with spatial location. Due to the small sample capacity and calculation error of lag distance, the accuracy of geostatistics and single fractal dimensions needs to be improved. Multifractal spectral analysis does not need to consider the normality of data and can quantitatively represent local characteristics; therefore, its results have high reliability.


Introduction
The characteristics of soil particle distribution and hydraulic properties (hereafter "soil properties") are the basis for the study of soil moisture and solute transport [1]. However, due to soil heterogeneity, soil properties vary spatially. The determination of distribution characteristics of the soil properties at a certain scale is significant for reasonable sampling and parameter prediction.
Formerly, classical statistics and geostatistics were the main methods used to study spatial variation. Geostatistics establishes the relationship between lag distance and variables through the semivariance function, which is the main means of quantitative analysis of spatial variability and the basis of Kriging interpolation. However, geostatistics is frequently limited by the sample capacity, distribution and anisotropy, which always leads to errors in practical application. With the advancement of fractal theory, the validity and potential of applying fractals to analyse soil spatial variability have been verified by early research over the past decades; they are thus an important tool with which to quantify the spatial variability and scale transformation of soil properties [2]. The combination of fractals and geostatistics can realise the quantitative characterisation of soil property spatial variation [3]; that is, the fractal dimension increases with the distribution uniformity. Researchers all over the world have calculated the fractal dimension values of various soil properties and environmental variables of different soil types. Among them, the representative studies include soil organic matter [4,5], total nitrogen [6], pH value [7], iron manganese nodules [8], and so forth. Due to different soil types, research scales and sample capacity, the fractal dimension values are also varied, but they are all in the range of 1.7~1.96. A single fractal dimension is simple and convenient, but it can only represent the holistic spatial variability of soil properties [9]. With the development of multifractal theory, the multifractal spectrum is used to quantitatively evaluate the spatial variability of soil properties. Multifractal parameters can not only reflect the spatial variation of the soil properties, but also reveal the small-scale or partial characteristics embedded in the holistic space [10]. Zeleke and Si [11], Caniego et al. [2], Eghball et al. [12] and Liu et al. [13] have analysed the spatial variability of soil properties, such as electrical conductivity, organic matter content, soil pH value, soil nitrate, soil water and salt, according to multifractal theory and obtained good results.
In recent years, the aforementioned methods have been used to study the spatial variation of soil properties such as soil particle distribution [14][15][16], bulk density [17,18], saturated hydraulic conductivity [19][20][21][22] and unsaturated hydraulic conductivity [23,24], allowing researchers to make a series of achievements. However, multifractal research is mostly limited to the farmland scale [13,[25][26][27][28]. The study of multifractals when the scale increases gradually (with a sampling interval of approximately several hundred metres) are rarely reported in the literature. In addition, there is a less joint application of multiple methods and a lack of comparative study [29,30].
Yanlou township in Lankao county is in a former old channel of the Yellow River and the Yellow River flood area. Owing to the frequent sedimentary processes and human activities, the soil structure is complex and strongly changes spatially, which is typical of the alluvial plain area in the lower Yellow River area. Therefore, this study focuses on Yanlou township and comprehensively analyses the spatial variation characteristics of six soil properties in the research area by classical statistics, geostatistics and single fractal and multifractal methods. Meanwhile, the causes of the spatial distribution of soil properties were analysed according to the spatial variability. It is found that different research methods show different results and precision. We think that the multifractal method has advantages in the quantitative evaluation of spatial variability. This research can provide a reference for the study of spatial variation, whose scale extends from the field to the regional level.

Geographical Location
This study was conducted in Yanlou township, Lankao county, Henan province (114 • 57 23" E~115 • 00 02" E, 34 • 54 27" N~34 • 52 31" N), which is adjacent to the Yellow River in the northeast and is part of the lower Yellow River alluvial plain. The Quaternary Holocene strata in the study area are well-developed. Due to the repeated diversion and flooding of the Yellow River, deposits repeatedly appear, forming a large-scale alluvial fan comprising mainly loess mixed with silty fine sand and silty clay.
The study area was roughly square, with a side length of 4 km and a total area of 16 km 2 . The terrain was relatively flat, with a maximum gradient of no more than 1.25% and a maximum height difference of 6.9 m. Most of the research areas were farmlands and villages, and the planting mode was winter wheat rotation with summer maize. The study area was divided into 64 grids with a grid cell size of 500 × 500 m.

Sampling and Indoor Experiment
The centre point of each grid was the coordinate of the sampling point. After the grid of the study area is surveyed and determined, a real-time kinematic (RTK) instrument is used to determine the coordinates and elevation of each grid centre. The RTK precision of this instrument is ±(8 + 1 × 10 −6 × D) mm (plane precision) and ±(15 + 1 × 10 −6 × D) mm (elevation precision). To ensure that the sampling point was in the grid and avoided houses and streets, the actual sampling point deviated from the centre point ( Figure 1). Soil samples were taken from 20-25 cm below the surface to avoid wrapping plant roots.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 3 of 15 method has advantages in the quantitative evaluation of spatial variability. This research can provide a reference for the study of spatial variation, whose scale extends from the field to the regional level.

Geographical Location
This study was conducted in Yanlou township, Lankao county, Henan province (114°57′23″ E~115°00′02″ E, 34°54′27″ N~34°52′31″ N), which is adjacent to the Yellow River in the northeast and is part of the lower Yellow River alluvial plain. The Quaternary Holocene strata in the study area are well-developed. Due to the repeated diversion and flooding of the Yellow River, deposits repeatedly appear, forming a large-scale alluvial fan comprising mainly loess mixed with silty fine sand and silty clay.
The study area was roughly square, with a side length of 4 km and a total area of 16 km 2 . The terrain was relatively flat, with a maximum gradient of no more than 1.25% and a maximum height difference of 6.9 m. Most of the research areas were farmlands and villages, and the planting mode was winter wheat rotation with summer maize. The study area was divided into 64 grids with a grid cell size of 500 × 500 m.

Sampling and Indoor Experiment
The centre point of each grid was the coordinate of the sampling point. After the grid of the study area is surveyed and determined, a real-time kinematic (RTK) instrument is used to determine the coordinates and elevation of each grid centre. The RTK precision of this instrument is ±(8 + 1 × 10 −6 × D) mm (plane precision) and ±(15 + 1 × 10 −6 × D) mm (elevation precision). To ensure that the sampling point was in the grid and avoided houses and streets, the actual sampling point deviated from the centre point ( Figure 1). Soil samples were taken from 20-25 cm below the surface to avoid wrapping plant roots.  Several surface soil samples were taken by ring knives, and the samples were transported to the laboratory for air drying. The division of soil particle size, the symbols and the test methods of soil properties are shown in Table 1. The test methods in Table 1 are strictly in accordance with the Standards for Geotechnical Testing Method (CN), GB/T50123-2019 [31]. Based on the test results, the soil texture distribution of 64 samples is shown in Figure 1. The classification of soil texture refers to the Code for Investigation of Geotechnical Engineering (CN), GB50021-2001 [32]. The coordinates of sampling points and test results can be seen in Table S1.

. Mesh Generation
To study spatial variation by multifractal theory, it is necessary to divide the grid into different scales [13]. The study area was divided into three scales, corresponding to grids of 64, 16 and 4 ( Figure 2). When the grid grew larger, the arithmetic mean value of the soil property index of the points in the cell grid was taken as the soil property index of the large grid measuring points. Several surface soil samples were taken by ring knives, and the samples were transported to the laboratory for air drying. The division of soil particle size, the symbols and the test methods of soil properties are shown in Table 1. The test methods in Table 1 are strictly in accordance with the Standards for Geotechnical Testing Method (CN), GB/T50123-2019 [31]. Based on the test results, the soil texture distribution of 64 samples is shown in Figure 1. The classification of soil texture refers to the Code for Investigation of Geotechnical Engineering (CN), GB50021-2001 [32]. The coordinates of sampling points and test results can be seen in Table S1.

. Mesh Generation
To study spatial variation by multifractal theory, it is necessary to divide the grid into different scales [13]. The study area was divided into three scales, corresponding to grids of 64, 16 and 4 ( Figure  2). When the grid grew larger, the arithmetic mean value of the soil property index of the points in the cell grid was taken as the soil property index of the large grid measuring points.

Geostatistics and Fractal Method
Geostatistics is an effective method to study the spatial distribution structure characteristics of regional variables. Its basic tool is a semivariate function that can be estimated by the following formula [33]: where γ(h) is the variation function, Z(x) is the value of the regionalised variable at sampling point x, N(h) is the number of pairs with interval h and h is the interval, which is called the lag distance.
Variograms can reflect and describe many properties of regionalised variables, and it is an important tool to analyse their spatial variation. The fractal dimension, D, is a measure of complexity. Due to the complexity of soil structure and the local micro differences of internal soil factors, the soil property parameters are irregular and random, like random fractional Brownian motion [34]. The semivariate function in one dimension is:

Geostatistics and Fractal Method
Geostatistics is an effective method to study the spatial distribution structure characteristics of regional variables. Its basic tool is a semivariate function that can be estimated by the following formula [33]: where γ(h) is the variation function, Z(x) is the value of the regionalised variable at sampling point x, N(h) is the number of pairs with interval h and h is the interval, which is called the lag distance.
Variograms can reflect and describe many properties of regionalised variables, and it is an important tool to analyse their spatial variation. The fractal dimension, D, is a measure of complexity. Due to the complexity of soil structure and the local micro differences of internal soil factors, the soil property parameters are irregular and random, like random fractional Brownian motion [34]. The semivariate function in one dimension is: where Z(x) and Z(x + h) are the measured values at x and x + h, respectively; h is the interval (lag); and H is the power index. For Brownian motion, the power index H is 0.5, while for the variation in soil properties, the H value ranges from 0 to 1.0. With the increase in H, the variation in soil properties gradually weaken. The fractal dimensions of these Brownian motions are given by the formula D = 2 − H.
γ(h) and h are drawn in double logarithmic space according to Formula (2). The fractal dimension D is calculated from the slope m of the regression line.
The fractal dimension, D, represents the structure between samples. The smaller the value of D, the greater the difference in soil properties between samples; that is, the worse the degree of uniformity; on the contrary, the larger the value of D, the smaller the difference in soil properties between samples, and the better the degree of uniformity.

Multifractal Method
To further understand the role of local conditions in the formation of fractal bodies, researchers have proposed the multifractal theory that mainly discusses the probability distribution of a certain parameter. For the multifractal method, previous studies include Renyi [35] and Chhabra and Jensen [36]. The following is the calculation formula of multifractal parameters: where µ i (ε) is the mass probability, M i (ε) is the variable, ε is the scale, N(ε) is the number of variables under the scale, D(q) is the generalized fractal dimension and q is the order of statistical moment.
where µ i (q, ε) is the q-order mass probability, α(q) is the singularity index and f (α) is the singularity spectrum.
q-D(q) and α-f (α) are two sets of parameters for describing multifractals. In practical application, we often use the partition function method and Legendre transformation [37] to establish the relationship between two sets of parameters: where χ q (ε) is the partition function and τ(q) is the q-order mass index. The quality index corresponding to each q value can be obtained by calculating the slope of the fitting line between log(ε) and log[χ q (ε)] by the least-squares method. Additionally, D(q) and τ(q) have the following relationship: The relationship between α(q) and f [α(q)] and τ(q) is as follows: The value of q can range from negative to positive infinity. When q > 0 and D(q) decreases with q, the variable has multifractal characteristics [12]. When q >> 1, the large part of χ q (ε) is dominant, while when q << −1, the small part of χ q (ε) is dominant. The higher-order statistical moments can extract extreme values and enlarge their contribution. The relationship between α and f (α) can be used to analyse the local characteristics of multiple fractal bodies. The value of ∆α reflects the overall heterogeneity of variable distribution; the larger the value of ∆α, the greater the variability of the variable distribution. The asymmetry of the α-f (α) curve reflects the dominant power of the uneven variable distribution. If the curve is shaped like a "left hook", the spatial variation of the variable is dominated by smaller values; on the contrary, if the curve is shaped like a "right hook", it is dominated by the larger values.

Classical Statistical Features
Classical statistics can directly reflect the volatility of variables. In general, the coefficient of variation (Cv) (Cv = standard deviation/mean) is used to measure the degree of variation. When 0 ≤ Cv ≤ 30%, it can be called weak variation; 30% < Cv ≤ 100% can be called medium variation; and Cv > 100% can be called strong variation. It can be seen from Table 2 that the content of silt, bulk density and saturated water content have weak variation; clay content has medium variation; and sand content and saturated hydraulic conductivity have strong variation. However, classical statistics cannot establish the relationship between variables and spatial location; therefore, geostatistics and the fractal dimension are used for further analysis.

Geostatistical Analysis
According to the Kolmogorov-Simirnov normal distribution test, the indexes of sand content, clay content and saturated hydraulic conductivity in the study area do not conform to the normal distribution, but all the six indexes follow a lognormal distribution (Table S2). To facilitate comparison, all six indexes were transformed into logarithms, which were imported into the software GS + 9.0 for geostatistical analysis. GS9.0 + software allocates point pairs by establishing the lag class. When the active lag distance and lag distance interval were determined, the number of pairs was equal to the number of lag classes. Generally, the active lag distance was 80% of the maximum interval distance, 3120.86 m. Considering the distribution of sampling points, the minimum interval of 361.31 m was taken as the lag interval, and uniform distribution was selected for lag class. According to the parameter setting, there were eight pairs of points in this study. The results of isotropic analysis are shown in Table 3.
The nugget is caused by random factors such as sample experimental error [38], which reflects the variation of regional variables on a small sampling scale; the sill is the half square difference of the variation function when the variable is stationary, which is the overall characteristic of the variable. The larger the sill, the stronger the spatial heterogeneity; the smaller the nugget coefficient, the stronger the spatial autocorrelation of the system; otherwise, the change between samples is closely related to random factors. The range reflects the maximum autocorrelation distance of the variable. Table 3 reveals that there are four kinds of statistical models for the six soil properties: Gaussian, exponential, linear and spherical. The small nugget coefficient of the silt content indicates that the silt content in the study area has a strong spatial autocorrelation. Additionally, its range is the smallest, which shows that the scale range of its spatial continuity is the smallest [18]. The nugget coefficients of the other five soil properties were between 0.4 and 0.5, which can be considered to have moderate spatial autocorrelation. The order from the largest to the smallest is: K s > sand > γ > clay > θ s . The range values of the other five soil properties range from 1913.92 m to 5097.43 m. The sequence of the range values from the largest to the smallest is: clay > sand > θ s > K s > γ.
Based on the model fitted in Table 3, simple Kriging interpolation maps of soil properties (logarithm) were drawn (Figure 3), which can directly reflect the spatial distribution of soil properties. The spatial distribution of the six soil properties is varied, and there is no obvious regularity in the spatial distribution of bulk density, saturated water content and saturated hydraulic conductivity. However, two conclusions can be drawn from Figure 3: One is that the sand content gradually increases from north to south, while the clay content shows the opposite. Secondly, the colour of the silt content changes slightly, which is consistent with the extremely weak nugget coefficient of silt content, although in classical statistics, the coefficient of variation of silt content is not the lowest.

Single Fractals
The "Fractal" module (isotropy) in GS + 9.0, with the same geostatistical parameter settings as the previous section, was used to calculate the fractal dimension, D, and determination coefficient, R 2 , of six soil properties in the study area (Table 4).

Single Fractals
The "Fractal" module (isotropy) in GS + 9.0, with the same geostatistical parameter settings as the previous section, was used to calculate the fractal dimension, D, and determination coefficient, R 2 , of six soil properties in the study area (Table 4). According to Table 4, the fractal dimension order D is: K s > silt > θ s > γ > sand > clay. The fractal dimension, D, is different from the coefficient of variation, indicating that it is related to the spatial position. However, the determination coefficients of silt content, saturated water content and saturated hydraulic conductivity are less than 0.6, which led to the low accuracy.

Multifractal Analysis
To determine whether the six soil properties indexes in the study area have multifractal characteristics, it was necessary to construct a q-D(q) curve for analysis (Figure 4). After debugging, q was set to −1.4 < q < 1.4, and the step size was set to 0.2. The calculation results of the multifractal parameters can be seen in Table S3.  According to Table 4, the fractal dimension order D is: Ks > silt > θs > γ > sand > clay. The fractal dimension, D, is different from the coefficient of variation, indicating that it is related to the spatial position. However, the determination coefficients of silt content, saturated water content and saturated hydraulic conductivity are less than 0.6, which led to the low accuracy.

Multifractal Analysis
To determine whether the six soil properties indexes in the study area have multifractal characteristics, it was necessary to construct a q-D(q) curve for analysis (Figure 4). After debugging, q was set to −1.4 < q < 1.4, and the step size was set to 0.2. The calculation results of the multifractal parameters can be seen in Table S3. (e) (f) It can be seen from Figure 4 that the D(q) of the six soil properties in the study area gradually decreases with q, and it can be determined that the six soil properties have roughly multifractal characteristics. However, the variation in D(q), that is, the generalised fractal dimension of the six soil properties, is different. The order of ΔD(q) is sand > Ks > clay > silt > θs > γ, which is consistent with the variation coefficient. Further analysis of the multifractal characteristics requires the multifractal spectrum, represented by the α-f(α) curve ( Figure 5).   It can be seen from Figure 4 that the D(q) of the six soil properties in the study area gradually decreases with q, and it can be determined that the six soil properties have roughly multifractal characteristics. However, the variation in D(q), that is, the generalised fractal dimension of the six soil properties, is different. The order of ∆D(q) is sand > K s > clay > silt > θ s > γ, which is consistent with the variation coefficient. Further analysis of the multifractal characteristics requires the multifractal spectrum, represented by the α-f (α) curve ( Figure 5).
(e) (f) It can be seen from Figure 4 that the D(q) of the six soil properties in the study area gradually decreases with q, and it can be determined that the six soil properties have roughly multifractal characteristics. However, the variation in D(q), that is, the generalised fractal dimension of the six soil properties, is different. The order of ΔD(q) is sand > Ks > clay > silt > θs > γ, which is consistent with the variation coefficient. Further analysis of the multifractal characteristics requires the multifractal spectrum, represented by the α-f(α) curve ( Figure 5).   The shape of the multifractal spectrum curve can reflect the spatial distribution characteristics of the soil properties [25]. The spectral width, Δα, of the multifractal spectrum can reflect the heterogeneity of variables ( Table 5). The magnitude order of Δα of the six soil properties in the study area is Ks > clay > sand > silt > θs > γ. The research results of Δα, the variation coefficient and the single fractal dimension are not consistent and show that the principles of the multifractal method and geostatistics are different. The symmetry of the multifractal spectrum curve can reflect the distribution dominance of variables. According to Figure 5, the multifractal spectrum of sand, silt, clay content and saturated hydraulic conductivity is obviously "left hook". In the process of reducing the sampling scale of these four soil properties, the large data set tends to stabilise gradually, so the spatial variation is mainly dominated by the small value data. When |q| is small, the symmetry of the multifractal spectrum curve of bulk density and saturated water content is weak, so it is necessary to define Δf[Δf = f(αmin) − f(αmax)] to determine its symmetry. The results show that the multifractal spectrum is "right hook" and the variability is dominated by large-value data ( Table 5).   The shape of the multifractal spectrum curve can reflect the spatial distribution characteristics of the soil properties [25]. The spectral width, ∆α, of the multifractal spectrum can reflect the heterogeneity of variables ( Table 5). The magnitude order of ∆α of the six soil properties in the study area is K s > clay > sand > silt > θ s > γ. The research results of ∆α, the variation coefficient and the single fractal dimension are not consistent and show that the principles of the multifractal method and geostatistics are different. The symmetry of the multifractal spectrum curve can reflect the distribution dominance of variables. According to Figure 5, the multifractal spectrum of sand, silt, clay content and saturated hydraulic conductivity is obviously "left hook". In the process of reducing the sampling scale of these four soil properties, the large data set tends to stabilise gradually, so the spatial variation is mainly dominated by the small value data. When |q| is small, the symmetry of the multifractal spectrum curve of bulk density and saturated water content is weak, so it is necessary to define ∆f [∆f = f (α min ) − f (α max )] to determine its symmetry. The results show that the multifractal spectrum is "right hook" and the variability is dominated by large-value data ( Table 5). 1 α min -minimum value of α; 2 f (α min )-value of f (α), when α = α min . 3 α max -maximum value of α; 4 f (α max )-value of f (α), when α = α max . 5 ∆α = α max − α min ; 6 ∆f = f (α min ) − f (α max ).

Influencing Factors of Variability
Geostatistics shows that the distribution of silt content is stable, but the fluctuation of sand and clay content is obvious. From the perspective of structural factors, the study area is in the historical Yellow River floodplain. Its soil parent material mainly comes from fourth-century Holocene alluvial deposits, predominantly composed of silty sand. However, as the Yellow River flows through the study area, the local topographical differences result in different sedimentary separation. From the point of view of random factors, the content of clay in soil may be affected by human factors, such as farming and engineering construction, which cause variation in soil structural components.
The spatial variability of six surface soil properties in the study area is quite different, which indicates that soil properties show different development trends under the sampling interval. The soil particle distribution (sand, silt and clay content) and bulk density, as the basic physical properties of soil, largely determine the hydraulic properties of soil. The saturated hydraulic conductivity fluctuates as sharply as the content of sand and clay, which also reflects the strong effect of sand and clay on the saturated hydraulic conductivity. The variation in bulk density and saturated water content is small, and its multifractal characteristics are also weak. This is consistent with the results of Zeleke and Si [39], Guo et al. [20] and Guan et al. [25].

Comparison of Research Methods
Classical statistics can directly reflect the degree of variation of six soil properties, but it has significant limitations in analysing spatial variation. The spatial position and local characteristics of variable distribution are not considered in classical statistics, which leads to the inconsistency between the research results of the variation coefficient and geostatistics, and single and multifractal analyses. Geostatistics establishes the relationship between variable variance and lag distance, which is a powerful tool for analysing spatial autocorrelation and interpolation. The fractal dimension, D, is another important parameter in geostatistical analysis, which can quantitatively describe the exponential relationship between γ(h) and h. However, geostatistics has a large error due to the small sample capacity and lag distance calculation. For example, the geostatistical model of saturated hydraulic conductivity and the determination coefficient of the fractal value fitting line are small in this study, which led to the insufficient accuracy of the geostatistical analysis. Furthermore, the multifractal analysis considers spatial variation from another perspective, i.e., the geological process of multiple activities often produces self-similar fields, which refer to the geometric similarity maintained by the site when changing the measurement scale [40]. It establishes the relationship between the distribution of variables and the scale ε (sampling interval). It does not need the data to satisfy the normal distribution, which is a significant advantage of the multifractal analysis of spatial variation. Compared with the single fractal method, the multifractal spectrum parameters can quantitatively represent the degree of variation and help to screen the extreme value spatial distribution, which is an important indicator for local features in the region.

1.
The results of classical statistics show that the silt content, bulk density and saturated water content in the study area have weak variation; the clay content has medium variation; and the sand content and saturated hydraulic conductivity have strong variation. The variability of the six soil properties is quite different. The order of variability is sand > K s > clay > silt > θ s > γ, and the variation coefficient is consistent with the result of the generalised fractal dimension, ∆D(q) (q > 0).

2.
The geostatistical results (isotropy) after logarithmic transformation show that the nugget coefficient of the silt content in the study area is 0.002, so the silt content has a strong spatial autocorrelation. Its range is also the smallest, which indicates that the scale range of its spatial continuity is the smallest. The nugget coefficients of the other five kinds of soil properties are all between 0.4 and 0.5, which can be considered moderate spatial autocorrelation. The nugget coefficient can also reflect the spatial variability of soil properties to some degree, and the order from the largest to the smallest is K s > sand > γ > clay > θ s > silt.

3.
The single fractal (isotropy) shows that the relationship of the fractal dimension, D, is K s > silt >θ s > γ > sand > clay. According to the meaning of D, the larger the fractal dimension value, the smaller the variability, which shows that the variables after considering the spatial position factor (lag distance) have different statistical characteristics from the variables themselves. The fractal dimension, D, can be used as a simple and direct parameter to measure the relationship between the lag and variogram. However, due to the small sample capacity and lag distance calculation, the low accuracy leads to substantial limitations in practical application.

4.
The six soil properties in the study area have multifractal characteristics to different degrees. According to the multifractal spectrum curve, the relationship between the spectrum width, ∆α, and the soil properties shows that different soil properties have different variability. The relationship order of ∆α is K s > clay > sand > silt > θ s > γ. The symmetry of the multifractal spectrum curve can reflect the distribution dominance of variables. The content of sand, silt, clay and saturated hydraulic conductivity are obviously "left hook", and the spatial variation is mainly dominated by small-value data. ∆f values of bulk density and saturated water content are less than zero, which shows that the multifractal spectrum is "right hook", and its variability is mainly dominated by large-value data. The multifractal method does not need to consider the normality of the data. It has higher accuracy than geostatistics and single fractal analysis and has obvious advantages in the quantitative characterisation of local characteristics and extreme values.