The Use of Various Rainfall Simulators in the Determination of the Driving Forces of Changes in Sediment Concentration and Clay Enrichment

: Soil erosion is a complex, destructive process that endangers food security in many parts of the world; thus, its investigation is a key issue. While the measurement of interrill erosion is a necessity, the methods used to carry it out vary greatly, and the comparison of the results is often di ﬃ cult. The present study aimed to examine the results of two rainfall simulators, testing their sensitivity to di ﬀ erent environmental conditions. Plot-scale nozzle type rainfall simulation experiments were conducted on the same regosol under both ﬁeld and laboratory conditions to compare the dominant driving factors of runo ﬀ and soil loss. In the course of the experiments, high-intensity rainfall, various slope gradients, and di ﬀ erent soil surface states (moisture content, roughness, and crust state) were chosen as the response parameters, and their driving factors were sought. In terms of the overall erosion process, the runo ﬀ , and soil loss properties, we found an agreement between the simulators. However, in the ﬁeld (a 6 m 2 plot), the sediment concentration was related to the soil conditions and therefore its hydrological properties, whereas in the laboratory (a 0.5 m 2 plot), slope steepness and rainfall intensity were the main driving factors. This, in turn, indicates that the design of a rainfall simulator may a ﬀ ect the results of the research it is intended for, even if the di ﬀ erences occasioned by various designs may be of a low order.


Study Area and Soil Properties
The field experiments were conducted near Gerézdpuszta, situated in the Koppány Valley adjacent to the floodplain of the Koppány River ( Figure 1). This area is located in the east of the Somogy Region, southwest Hungary. Owing to the loess-like deposits, the area is prone to severe soil erosion; however, a further impact on the soil is the reason why conventional tillage has become less popular in recent decades [26]. The hillslopes that are less eroded due to the lack of tillage are characterized by Cambisols, while mostly Regosols and Colluviums can be found in the cultivated areas [27]. The climate is moderately warm and wet, with an average annual temperature of 10.0-10.2 °C and average annual precipitation of 605-700 mm. Most of the hillsides are characterized by agricultural cultivation, almost half of which is situated on slopes steeper than 12% [28]. The consequent tillage depths have resulted in hard plough-pan with extremely low porosity and infiltration.
The field experiments and the sampling for the laboratory investigations were conducted on recently tilled (seedbed preparation) Regosol, which had a silty clay loam texture, as well as low CaCO3 and organic carbon contents (Table 1). Owing to the loess-like deposits, the area is prone to severe soil erosion; however, a further impact on the soil is the reason why conventional tillage has become less popular in recent decades [26]. The hillslopes that are less eroded due to the lack of tillage are characterized by Cambisols, while mostly Regosols and Colluviums can be found in the cultivated areas [27]. The climate is moderately warm and wet, with an average annual temperature of 10.0-10.2 • C and average annual precipitation of 605-700 mm. Most of the hillsides are characterized by agricultural cultivation, almost half of which is situated on slopes steeper than 12% [28]. The consequent tillage depths have resulted in hard plough-pan with extremely low porosity and infiltration.
The field experiments and the sampling for the laboratory investigations were conducted on recently tilled (seedbed preparation) Regosol, which had a silty clay loam texture, as well as low CaCO 3 and organic carbon contents (Table 1). Field rainfall simulations were carried out on plots that were 3 m in length and 2 m in width, in triplicate. Altogether, 6 plots were investigated: 3 on the gentler parts of the slope (with an incline of 7-8%), and 3 on steeper parts (17-18%) (see Appendix A). The plots on the same parts of the slope were placed 2-3 m from each other ( Figure 1). The soil samples for laboratory simulations (from the topmost 0-25 cm tilled layer) were collected from the site at which the in situ simulations were conducted. In the laboratory, the collected samples were loaded into the flume in the state in which they had arrived. Initially, the flume was covered by woven geotextiles to ensure free drainage through the flume bottom to the collector taps. The soil was compacted solely by the energy of two successive (simulated) rainfall events (30 mm) to achieve the 20 cm soil depth required for the investigations. To exclude the effects of splash erosion, we covered the soil surface by a polyethylene mesh during these compaction pretreatments.

Rainfall Simulators Used in the Study
Two rainfall simulators were used in the experiment, the Shower Power (SP) and the "ELTE" simulator. SP 02 is a frequently used [29,30], portable field unit that can be employed at the plot-scale with alternating nozzles (Figure 2a,b). The device has a plot of 3 × 2 m; however, in this instance, the irrigated plot size was 12 m 2 to exclude border effects [30]. The drop forming unit, at a height of 3 m, comprises an alternating axis equipped with two adjacent 80100 Veejet nozzles spaced at 2 m [31]. This distance ensures complete overlap between the two nozzles to gain a uniform spatial intensity and drop size distribution pattern. The intensity can be set to 30-100 mm h −1 , depending on the axis alternation frequency. The average kinetic energy of the simulator is 24 J m −2 mm −1 , ensuring a reasonable drop size distribution and kinetic energy simulation of an intense natural rain [31]. Field measurement results are presented in Table A1. The properties affecting the momentary intensity and drop spectrum during a simulation (pressure regulation, evaporation, wind, pump properties) may well vary in the field, causing anomalies in the preset intensity and kinetic energy [32]. Thus, it was the real, independently measured rainfall intensities that had been collected (an average for each single run) that were used for calculations instead of the preset ones. The ELTE simulator [33,34] is a point-scale laboratory device equipped with a fixed 1/2 HH 40 WSQ full-cone nozzle (Figure 2c,d). It was set at a height of 3 m using a constant pressure of 20 kPa, which created a standard intensity of 80 mm h −1 , with an average kinetic energy of 17 J m −2 mm −1 , well representing the effects of an intense natural rainfall [35]. The ELTE simulator has a sample/monolith carrier flume measuring (lwd) 100 × 50 × 20 cm. The steepness of the monolith can be regulated over a range of 0-100%. On the bottom of the flume, under the geotextile cover, a collector system is fitted to drain the leaching water through four taps. In the present investigation, each tap was closed, modelling the quasi impermeable plough-pan layer detected in the field and thus increasing/maximizing the comparability of the research design of the field and laboratory measurements. This research design does not take the effect of splash out (runoff and sediment loss due to raindrop reflex from the surface, which leave the flume across the borders) into account, causing lower runoff and sediment yield values.

Experimental Design
A total of 18 simulation runs were conducted: 12 in the field and 6 in the laboratory ( Table 2). As the soil was the same, the impact of rainfall intensity, slope gradient, and soil surface parameters (roughness, moisture content, and crust state) on soil loss could be examined. The notion of the comparability of the erosion observed in the 2 experiments (comparing an in situ soil to a transported disturbed soil sample) was based on conclusions reached by Thomaz and Pereira [36], who found that the hydro-erosional responses of soils in a conventional tillage system or unstructured soils are similar. The ELTE simulator [33,34] is a point-scale laboratory device equipped with a fixed 1/2 HH 40 WSQ full-cone nozzle (Figure 2c,d). It was set at a height of 3 m using a constant pressure of 20 kPa, which created a standard intensity of 80 mm h −1 , with an average kinetic energy of 17 J m −2 mm −1 , well representing the effects of an intense natural rainfall [35]. The ELTE simulator has a sample/monolith carrier flume measuring (lwd) 100 × 50 × 20 cm. The steepness of the monolith can be regulated over a range of 0-100%. On the bottom of the flume, under the geotextile cover, a collector system is fitted to drain the leaching water through four taps. In the present investigation, each tap was closed, modelling the quasi impermeable plough-pan layer detected in the field and thus increasing/maximizing the comparability of the research design of the field and laboratory measurements. This research design does not take the effect of splash out (runoff and sediment loss due to raindrop reflex from the surface, which leave the flume across the borders) into account, causing lower runoff and sediment yield values.

Experimental Design
A total of 18 simulation runs were conducted: 12 in the field and 6 in the laboratory ( Table 2). As the soil was the same, the impact of rainfall intensity, slope gradient, and soil surface parameters (roughness, moisture content, and crust state) on soil loss could be examined. The notion of the comparability of the erosion observed in the 2 experiments (comparing an in situ soil to a transported disturbed soil sample) was based on conclusions reached by Thomaz and Pereira [36], who found that the hydro-erosional responses of soils in a conventional tillage system or unstructured soils are similar. A gentle slope of 5% was set in the laboratory, whereas all values <10% were classified as gentle in the field. The steep slope in the laboratory measurements was 12%, while slopes >10% were considered as steep in the field. Due to the dynamic nature of variable changes in surface roughness and moisture content under the effects of precipitation, we conducted 8 dual experiments, resulting in a total of 16 simulations. The first simulation runs (FG1; FG3; FG5; FS1; FS3; FS5; LG1; LS1) modeled the seedbed (recently tilled, relatively dry) soil with a rough surface, and the second runs (FG2; FG4; FG6; FS2; FS4; FS6; LG2; LS2) were conducted over the sealed and smoothed soil with the evolved crust and higher soil moisture content.
In the remaining two laboratory simulations (LP1; LP2), we simulated extremes of moisture content (modeling "drought" and "inland inundation") on a nearly flat surface (2% slope). A wide moisture content interval was necessary because the soil moisture content is an important factor controlling runoff [37], and its combined effect with the slope gradient could present an increased impact [38].
Four degrees of intensity (55, 70, 80, and 100 mm h −1 ) were planned and set for the field simulations. To measure the real amount of rainfall reaching, we placed the soil surface collector dishes next to the plot borders. Owing to the changes in evaporation, wind, and fluctuation in water pressure at the nozzles, the intensities as actually measured differed from those that had been set. For calculations, we used the measured amount and intensities. The time was recorded during all simulated rainfalls. All experiments were divided into 4 time periods: (1) time to ponding-in this period the infiltration exceeded rainfall intensity. The length of this period is proportional to the initial infiltration capacity; (2) ponding period (the length of which refers to the water storage capacity of the surface, or, in other words, the surface roughness); (3) runoff; and (4) the post-rain runoff period (runoff after the rain stops)-the length of this period also reflects the volume of water stored on the surface and indicates the degree of spatial connectivity. As the length of the ponding period describes the surface storage capacity, it is therefore an indirect measure of microtopography. During sample collection, all runoff was measured and collected, and runoff volumes were related to time intervals. After precipitation stopped, the runoff was collected into a separate bucket, and its volume was measured. The times at which ponds appeared and the runoff started and ended, as well as the length of runoff time after precipitation ended, were also recorded [29,30].

Runoff and Soil Loss Characterization
The total amount of soil loss was weighed after drying for 48 h at 60 • C. The sediment concentration (SC) of each sub-sample was calculated as the dry sediment mass per runoff volume (g L −1 ). The PSD of each sub-sample was determined in triplicate with a laser diffractometer (Horiba LA-950). For the PSD measurements, we selected the Mie theory of light scattering because it is more precise in the determination of clay fractions than Fraunhofer diffraction [39]. The refractive indices (RI) based on the Horiba software calculation were set to 1.63 (real RI) and 0.3 (imaginary RI).
The constant runoff rate was calculated using linear regression on the measured cumulative runoff volumes with time (practically during the second part of period 3) on the basis of the method described by Jakab et al. [30], in which the steepness of the fitted linear function provides the constant runoff rate (Table A2). The runoff coefficient was determined as the percentage of the constant runoff rate related to rainfall intensity. The clay enrichment CE values of each simulation were calculated as the median clay content values of the sub-samples in relation to the median clay content of the in situ soil. The upper limit for the clay fraction was set at 8 µm [40].

Statistical Analyses
To obtain an overall picture of the dataset, we produced box and whisker plots. Differences between the parameters of the various laboratory and field samples were calculated using independent samples Mann-Whitney U tests [41], while in other cases when more than 2 groups were compared, an independent samples Kruskal-Wallis test [42] was used. It is difficult to determine the combined impact of the most influential independent variables (e.g., SC) in a complex system. However, in many cases, multivariate statistical tools can assist in resolving such problems. For example, variance analysis has been employed in a study on the effect of rain intensity, slope length, and gradient on SC [25]. As in other studies [43][44][45], principal component analysis (PCA) [46] was used to determine how SC influences response parameters measured in the soil.
In the present study, to explore how soil parameters affect the SC, primary (slope gradient (%), rain intensity (mm h −1 )) and derived properties (constant runoff/infiltration (%), ponding period length (s), constant runoff rate (mm h −1 )) of the field and laboratory measurements were selected as input parameters for PCA. Since they are all related to hydrology and runoff formation and, furthermore, independent of soil loss, they are suitable for the measurement and prediction of the likely susceptibility of a given soil to erosion, and thus, also, the sole chosen explanatory parameter, SC. The median diameter of the soil that was lost (µm) was also selected as a response parameter, since it also might reflect SC.
Prior to PCA, we conducted Bartlett's test [47] to investigate the applicability of the input correlation matrices (field and laboratory) to PCA. Both the field (χ 2 (15) = 56.63, p = 9.47 × 10 −7 ) and the laboratory datasets (χ 2 (15) = 85.87, p = 5.8 × 10 −12 ) were found to be statistically suitable for the application of PCA. Thus, the variability of the variables input to PCA was compressed into principal components on the basis of their linear relationships, as presented by Hatvani [48,49]. It should be noted that the observations' principal components (PCs) are referred to as PC scores. The elements of the eigenvectors of the empirical correlation matrix will be referred to as loadings, with these measuring the relationship of the coordinates and the PCs with the Pearson correlation coefficient. Loadings outside the ±0.7 interval are considered as meaningful.
During PCA, we took the PCs into account on the basis of their scree plots [50] and their eigenvalues, which had to be above 1 [51]. Thus, taking into account the previous considerations, the input variables' time series were practically reduced to vectors with uncorrelated coordinates (PC scores) using the first two PCs.
All computations were performed using STATISTICA 10, IBM SPSS 26, R, and visualized in MS Excel 2016.

Variations in Ponding and Runoff Related Properties
Ponding and time to runoff periods did not differ significantly in the lab and the field (U(N Field = 12, N Lab = 6) = 49, p = 0.25) and (U(N Field = 12, N Lab = 6) = 56, p = 0.067) at α = 0.05. However, the duration of runoff and post-rain runoff in the field and the lab were significantly different (U(N Field = 12, N Lab = 6) = 69 and 11, respectively, p < 0.018) (Figure 3).
Water 2020, 12, x FOR PEER REVIEW 8 of 17
In the field, the time to runoff periods were shorter than 2 min. The time to ponding varied from 3 to 77 s, the ponding period varied from 9 to 73 s, and the time to runoff varied from 34 to 119 s. Shorter ponding times corresponded to more intense precipitation (Table A3). The time intervals of the before-runoff periods were shorter, and the post-rain runoff period was lengthened by the inhibited infiltration due to crusting and the improved connectivity on the surface. On the crusted surfaces, the duration of the ponding periods was only half that of the tilled surfaces, except for the following dual experiments: FG1-2, FS1-2, and FS3-4. In these three exceptions, the tilled and crusted surfaces entered the ponding stage for the same duration. The post-rain runoff on crusted surfaces was longer by approximately 1-9 min when compared with that on the tilled surface in the same dual experiment. The effect of slope steepness was only notable during the post-rain runoff period, which was longer on gentle slopes. The median of this closing period was 465 s on gentle slopes and 198 s on steep slopes.
In the laboratory, the runoff in three experiments (LG1-2 and LS1) began after 6-15 min; two of these plots were recently tilled, and one had a crusted surface on a slope of 5%, and which also experienced the highest rainfall intensity of 141 mm h −1 . The length of the ponding period decreased with an increasing slope gradient, and the time to runoff decreased with the appearance of the crust. No runoff was detected after rain in experiments where the surface was crusted (LG2, LS2; LF2).
The runoff coefficients did not differ between the slope gradients or locations: H(18) = 0.711, p = 0.701 (Figure 4). The constant runoff rate was significantly higher in the laboratory (U(NField = 12, NLab = 6) = 66, p = 0.003) than in the field. This difference, however, was due to the inclusion of the results of the extreme moisture and dry experiments (Figure 4). It mostly depended on the intensity and moisture content/surface conditions, especially in the field. Asterisks refer to * crusted, ** wet, and *** dry surface conditions. F-field, L-laboratory; P-plain (2% slope); G-gentle slope; S-steep slope; for details, see Table 2.
In the field, the time to runoff periods were shorter than 2 min. The time to ponding varied from 3 to 77 s, the ponding period varied from 9 to 73 s, and the time to runoff varied from 34 to 119 s. Shorter ponding times corresponded to more intense precipitation (Table A3). The time intervals of the before-runoff periods were shorter, and the post-rain runoff period was lengthened by the inhibited infiltration due to crusting and the improved connectivity on the surface. On the crusted surfaces, the duration of the ponding periods was only half that of the tilled surfaces, except for the following dual experiments: FG1-2, FS1-2, and FS3-4. In these three exceptions, the tilled and crusted surfaces entered the ponding stage for the same duration. The post-rain runoff on crusted surfaces was longer by approximately 1-9 min when compared with that on the tilled surface in the same dual experiment. The effect of slope steepness was only notable during the post-rain runoff period, which was longer on gentle slopes. The median of this closing period was 465 s on gentle slopes and 198 s on steep slopes.
In the laboratory, the runoff in three experiments (LG1-2 and LS1) began after 6-15 min; two of these plots were recently tilled, and one had a crusted surface on a slope of 5%, and which also experienced the highest rainfall intensity of 141 mm h −1 . The length of the ponding period decreased with an increasing slope gradient, and the time to runoff decreased with the appearance of the crust. No runoff was detected after rain in experiments where the surface was crusted (LG2, LS2; LF2).
The runoff coefficients did not differ between the slope gradients or locations: H(18) = 0.711, p = 0.701 (Figure 4). The constant runoff rate was significantly higher in the laboratory (U(N Field = 12, N Lab = 6) = 66, p = 0.003) than in the field. This difference, however, was due to the inclusion of the results of the extreme moisture and dry experiments (Figure 4). It mostly depended on the intensity and moisture content/surface conditions, especially in the field.

Variations in Sediment Concentration and Composition
SC values increased and differed significantly (H(18) = 11.075, p = 0.004) with the slope gradient, but were unaffected by the location or surface roughness ( ≥ 0.129; Figure 5). On the gentler (2-8%) slopes, SC values of 9-14 and 8-11 g L −1 were recorded, while on the steep (12-18%) slopes, SC values of 10-26 and 15-20 g L −1 were observed in the field and laboratory, respectively (Table A1). Figure 5. Comparison of the sediment concentration and clay enrichment in relation to the location, slope gradient, and surface properties. The most striking differences were caused by the slope gradient, and the SC increased while the CE decreased with an increasing slope gradient. The lowest SC values of 5.5 and 9.2 g L −1 were measured in the runoff from the experiments with a slope gradient of 2% slope and extremely wet or dry initial moisture contents, respectively. The SC was higher in the field when surfaces were crusted (Table 3). Table 3. Comparison of sediment concentrations (g L −1 ) among the different combinations of slope gradients and surface types. LF1 and LF2 were excluded from this comparison as their experimental setups were unique.

Slope Gentle Steep
Sediment concentration (g L −1 ) The difference in PSD between the in situ soil and soil loss was expressed by the CE ratio. As a general result, all soil loss had a CE > 1, except for the highest SC-related rain with the highest runoff ratio (CE = 0.94). Moreover, CEs were significantly higher in the laboratory (U(N Field = 12, N Lab = 6) = 61, p = 0.018) ( Figure 5, Table A1) compared to the field measurements. In contrast, CE did not vary significantly, neither with regards to the slope gradient (H(18) = 5.088, p = 0.079) nor the surface roughness (H(18) = 6.868, p = 0.076).

Regulating Properties of SC for Various Rainfall Simulations
By analyzing the linear relationships of the response variables, we found that the first two PCs (PC1 and PC2) were able to account for approx. 80% of the total original variance in both the laboratory and field experiments together (Table 4a). After correlating the PC scores with the explanatory variable, we found that the main influencing factors of SC differed between the field and laboratory (Table 4). The most influential parameters in the first PC for the field measurements were constant runoff/infiltration, constant runoff rate, particle size median of soil loss with a positive load, and ponding period with a negative sign indicating a response opposite to the others (Table 4a). Specifically, if the surface was smooth (short ponding period), the constant runoff rate increased, the flow could transport larger particles/aggregates, and the PSD median also increased. The first PC displayed a strong and significant (p < 0.01) linear relationship with the explanatory variable SC (Table 4b). The variance compressed into the second PC for the field measurements mostly accounted for the rainfall intensity, which had a smaller linear effect on SC (Table 4b). The loadings of slope gradient did not fall within the chosen ±0.7 interval in neither PC, indicating a secondary effect on SC in the field.
As the slope gradient increased in the laboratory experiments, the rainfall intensity, ponding period, and constant runoff rate decreased according to the first PC (Table 4a); however, the two most determining variables (constant runoff/infiltration and particle size median) in the second PC displayed a parallel change (Table 4a). Significant correlations were observed between the PCs and SC; it was exclusively related to the first PC (Table 4b). Thus, while the actual surface state and development acted together as influential processes on SC in the field, the first PC indicated that, in the laboratory, the constant runoff/infiltration and the particle size median of the soil loss were independent of SC.

Discussion
In this study, the soil type, high rainfall intensities, and slope gradients were the matching parameters during the experiments, while initial soil moisture content, surface state (crusting and smoothening), and the device used varied.
No differences were found in time to runoff in either the laboratory or the field. This may suggest that the disturbed soil transported in the laboratory yielded comparable infiltration values with that measured in the field, even though the high variance makes the comparability difficult. In contrast, a larger plot size may trigger poorer surface flow connectivity, resulting in higher ponding capacity in the field, as reported by Langhans et al. [52]. The runoff coefficient values demonstrated the impact of the crust on runoff in both locations, as roughness affected the drainage network [53] and a sealed surface inhibited infiltration [54], although, without the measurement of surface roughness and micro-topography, this is just an indication of its veracity. Initial soil moisture content could also affect the initiation of runoff; it begins earlier on a wet surface, as reported by Jin et al. [7]. Gómez and Nearing [55] also observed delayed time to runoff owing to higher surface roughness, but they did not observe differences in the runoff coefficient. This might have been due to the changes in surface connectivity, which could drastically affect the runoff volume even on rough surfaces [56]. Above and beyond these processes, variations in research design (e.g., differences in kinetic energy, splash loss due to a lack of fences in the laboratory) would surely have contributed to the difference.
As the infiltration-runoff equilibrium of the complex system is affected by several properties, it is impossible to provide a clear explanation of the differences. The point is that the method used will impact the results.
One of the main goals of this study was to test the comparability of field and laboratory SC and CE results. Concerning the constant runoff and constant runoff rate measured in this study, we found that the values from the experiments were similar when conducted on similar slope gradients (Figures 3  and 4). However, the research design (border effect in splash erosion (that is, unmeasured soil loss due to the lack of fences) in the laboratory might decrease SC. At the same time, in the present study, the CE was higher in the laboratory ( Figure 5). Zemke [57] also demonstrated the enrichment of fine particles in soil loss on small plots under simulated rainfalls. Presumably, this was due to the smaller plot size, in which deposition of the coarse fraction within a plot is inhibited, as also suggested by Jakab et al. [58], or is simply the result of the higher kinetic energy of the applied rainfall. In contrast, the similar CE from tilled and crusted surfaces of the present study was in contrast to the results of high-intensity small plot experiments on silty loams by Ding and Huang [11].
Concerning the overall erosion process, in the field, both the SC and CE influenced the surface properties, as also reported by Koiter et al. [5]. In the laboratory, the velocity of the moving water and its influencing factors (slope and rainfall intensity) were the main variables driving SC, and the available soil particles and amount of water determined the CE. The above difference, however, is an artefact of the different research designs, including, among others, the plot size, the border effect, and the differences in the drop spectra and kinetic energy of the rain. To achieve more precise and detailed results, additional measurements are needed with a greater degree of standardization in the conduct of the experiments (e.g., using the same nozzles in the field and the laboratory). Moreover, given that the above findings are the results of experiments carried out on a single soil, more investigations on other soil types are required before the results can be generalized.

Conclusions
Results gained using the two simulators agree with each other and reflect in similar ways on the soil and surface properties. Therefore, both simulators are suitable for calculations concerning the susceptibility of soils to erosion, and their results are comparable. The slope gradient was found to be an effective regulator only in the laboratory. Rainfall intensity was also more effective in the laboratory than in the field simulations. This, then, suggests that soil-related properties had a prominent role in driving sediment concentration in the field, whereas in the laboratory, slope and rainfall intensity were found to be driving factors independent of soil-regulated sediment concentrations. However, these findings are based on a limited number of simulations and investigate only one type of soil, and some may conclude that different rainfall simulators behave differently, resulting in the overestimation of the role of various sub-processes. Thus, the research design of a rainfall simulator may affect the results even though the differences thus evoked seem to be small in magnitude, probably beneath the sampling error caused by natural spatial heterogeneity of the soil. Ideally, in the future, an optimal research design will be to use paired scenarios for various simulation purposes. To do so, many comparative measurements are needed.