Large-Scale Modeling of Solum Thickness Based on the Local Topographic Attributes of Ground and Bedrock Surfaces

This study was aimed to address the importance of neighborhood scale and using bedrock topography in the soil-landscape modeling in a low-relief large region. For this study, local topographic attributes (slopes and curvatures) of the ground surface (DTM) and bedrock surface (DBM) were derived at five different neighborhood sizes (3×3, 9×9, 15×15, 21×21, and 27×27). Afterward, the topographic attributes were used for multivariate adaptive regression splines (MARS) modeling of solum thickness. The results demonstrate that there are statistical differences among DTM and DBM morphometric variables and their correlation to solum thickness. The MARS analyses revealed that the neighborhood scale could remarkably affect the soil – landscape modeling. We developed a powerful MARS model for predicting soil thickness relying on the multi-scale geomorphometric analysis (R 2 = 83%; RMSE= 12.70 cm). The MARS fitted model based on DBM topographic attributes calculated at a neighborhood scale of 9×9 has the highest accuracy in the prediction of solum thickness compared to other DBM models (R 2 = 61%; RMSE = 19cm). This study suggests that bedrock topography can be potentially utilized in soil-related research, but this idea still needs further investigations.


INTRODUCTION
Soil-landscape modeling is an essential component of sustainable land management (SLM) (Khanifar et al., 2020). The topography is the primary pedogenic factor that affects the pattern of soil properties variability across the landscape by controlling hydro-geomorphic processes (Khanifar and Khademalrasoul, 2021). This relationship allows us to utilizing the topographic attributes to model soil properties. Geomorphometry is the knowledge that deals with quantitative analysis of the topographic surface, and in general mode focuses on the derivation of topographic attributes from the Digital Elevation Models (DEMs) (Pike et al., 2009). Various topographic attributes of DEM (or morphometric variables) were employed to model soil thickness (e.g. Thompson et al., 2006;Mehnatkesh et al., 2013), which is a fundamental variable in the soil quality assessment.
One idea is to utilize the topographic attributes of the bedrock surface in the soil properties prediction. The bedrock is defined as the integrated solid rock that lies under un-consolidated surface materials, like soil and gravel. The bedrock appears in some zones at the ground's surface, but in regions can be located at a depth of one thousand meters below the earth's surface (Shangguan et al., 2017). Predicting depth to bedrock is a critical research field in geophysical studies. Depth to bedrock data has special applications in various fields of geoscience. Depth to bedrock influences energy and hydrology cycles and has the potential to be applied as an input variable for mapping natural hazards like earthquake and landslide (Yan et al., 2020).
An important setting in soil-landscape modeling is the operational scale, which is less considered. Operational scale refers to the size of the space in which processes associated with a particular phenomenon operate (Bian, 1997); it is defined in the form of a locally moving neighborhood window in geomorphometry. Most geomorphometric analyzes are based on the neighborhood size of 3×3, but this scale cannot be ideal for all soil-landscape modeling because of the diversity of pedogenic processes (Wang et al., 2010;Khanifar and Khademalrasoul, 2021).
Compliance the operational scale of pedogenic processes with the neighborhood size that is applied for mathematical calculations of topographic attributes can be an influential factor in the model performance. This study was conducted in a multi-scale framework to demonstrate the importance of operational scale and evaluate the idea of utilizing bedrock topography in the modeling of solum thickness using MARS.

Study Area and Data Collection
The study site with an area of 129560 km 2 is located in a region in eastern China (114° 05ʹ to 118° 16ʹ E and 32° 50ʹ to 35° 54ʹ N; Figure. 1). The land use of the major part of study site based on MODIS images (MCD12Q1; land cover product) is agricultural. The bedrock elevation map (DBM) of the study site represented in Figure (1) has a grid spacing of 90m. It is produced from the difference between the STRM DEM and depth to bedrock data that are having the same grid spacing. The data of Depths to bedrock were derived from www.globalchange.bnu.edu.cn, and more information about this dataset are discussed in (Yan et al., 2020). The lowest and highest bedrock elevation values are 307.80 and 1029.45m, which are nearby of these two locations the highest and lowest depth to bedrock, respectively ( Figure. (Patel et al., 2008). Solum thickness data have been derived from ISRIC's global database (Batjes et al., 2019).

Geomorphometry analysis
In this study, the local topographic attributes were only used in MARS modeling. The definition and formulas of local topographic attributes are given in Table (1). The topographic attributes of ground surface (DTM raster) and bedrock surface (DBM raster) were calculated based on Wood (1996) approach at five different neighborhood sizes (3×3, 9×9, 15×15, 21×21, 27×27). The Wood method is a generalized of the Evans algorithm (Evans, 1979) for larger operational scales.
Evans algorithm is approximates the partial derivatives by finite differences using only the 3×3 square-gridded moving neighborhood window. In this study, a Low Pass filter was used to eliminate local noises in the DEM and DBM. The effect of applying this filter on DEM is evident in the histograms presented in Figure (1). After geomorphometry analysis, the values of DTM and DBM morphometric variables were extracted for the sampling points' position.

Statistical analysis
Descriptive statistics of solum thickness and topographic attribute was performed in Statistica V.12 software. To assessing the statistical differences of topographic attributes across the groups of DTM and DBM, the Kruskal-Wallis test followed by a post hoc test (multiple comparisons of mean ranks) was used. The MARS algorithm was used to model the relationship between DTM and DBM topographic attributes and solum thickness at different neighborhood scales.

Multivariate Adaptive Regression Splines (MARS) is a non-linear and non-parametric regression
method that models the relationships between explanatory and response variables based on a set of spline functions called the basis function (BF) (Friedman, 1991). The forms of BFs are as follows: Where is called a knot and is one of the observed values of the explanatory variable ( ). The

MARS algorithm splits the range of predictor variables by knots into smaller regions and fits a
BF at each region. The general form of the MARS model is as follows: Where is the response variable predicted by the function ( ), is the number of terms, ° is a constant, is the mth BF's coefficient, and ℎ ( ) is individual BF or a product of two or more BFs. The MARS model is generated in two stages. In the first step, all possible BFs are consecutively added to the model, leading to a complex and over-fits model. In the second stage, the BFs with the least contribution are pruned. Finally, the optimal model is selected based on the generalized cross-validation (GCV), a measure of the goodness of fit. The prediction accuracy of the fitted models was evaluated using the coefficient of determination (R 2 ) and the mean square error (RMSE) (30% of dataset was selected for validation, N = 22). MARS analyses were performed using the Salford Predictive Modeler (SPM) software.

RESULTS AND DISCUSSION
The results of statistical analysis of the solum thickness are shown in in Table 2. The solum thickness of the studied soil profiles varies from 12 to 150 cm, with an average of 37.80. The Shapiro-Wilk test confirms that solum thickness does not have a normal distribution like topographic attributes. The solum thickness has high variability (CV = 80.80%) in the region.
According to the above BFs, when profile curvature is greater than 1.05337e-006, solum thickness will decrease by 2.73256e+007 times the amount in excess of the threshold value (BF1). Also, when profile curvature is less than the threshold, soil thickness will be 3.46253e+007 times the profile curvature in excess of 1.05337e-006 1/m (BF2). Some of BFs are produced by multiplying two BFs. BF12 indicates that the impact of maximal curvature has emerged through the interaction with the slope Aspect. The contribution of all BFs to the solum thickness values is shown in the 3-D plots presented in Figure 2b. Mehnatkesh et al. (2013) showed that the model based on slope gradient (local MV), catchment area (regional MV), wetness index, and sediment transport index (secondary MVs) could explain 76% of soil depth variability within a hilly area. In this study, the LMVs only have been utilized, but in many studies such as Gessler et al. (2000) and Thompson et al. (2006), it is observed that some MVs of other classes are also included in the soil depth models.
The DBM-9×9 model was able to explain 62% of the solum thickness variability in the region.
The statistical analysis of topographic attributes calculated at neighborhood scales of 9×9 and 27×27 are presented in Figure 3. The curvatures when are extracted from the topographic surface of the bedrock, a relatively broad range of positive and negative values can be observed, which representing high variability in bedrock surface bending between sampling positions. However, in the most DTM morphometric variables, the values range and the bending variability are remarkably lower. Up-scaling has also led to a dramatic decrease in the curvatures and slope gradient value range due to smoothing of the surface roughness, which results from the dilution of their potential information with other data on the large neighborhood size. The Kruskal-Wallis test followed by multiple comparisons post hoc test confirmed that the difference between some groups of slope gradient, maximum curvature, and minimum curvature was significant (at the 1 % level).

CONCLUSIONS
In this research, if the local topographic attributes were derived only at the basic neighborhood scale (3×3), which is used in most studies, just DTM and DBM models with lower performance would be generated. Increasing the neighborhood scale to an appropriate extent, in addition to reducing noise in the elevation data, can lead to improved modeling of soil properties. Since no procedure has been provided to find the optimal range of neighborhood size in each landscape, the use of multi-scale geomorphometric algorithms is appropriate. The modeling results confirm the idea of applying bedrock geometry in the soil-landscape modeling. It is recommended that a fine scale research be conducted to survey more the correlation between bedrock topography and soil properties in various landform units.