Hybrid data-model-based mapping of soil thickness in a mountainous watershed

. Soil thickness plays a central role in the interactions between vegetation, soils, and topography where it controls the retention and release of water, carbon, nitrogen, and metals. However, mapping soil thickness, here defined as the mobile regolith layer, at high spatial resolution remains challenging. Here, we develop a hybrid model that combines a process-based 10 model and empirical relationships to estimate the spatial heterogeneity of soil thickness with fine spatial resolution (0.5 m). We apply this model to two examples of hillslopes (south-facing and north-facing, respectively) in the East River Watershed in Colorado that validates the effectiveness of the model. Two independent measurement methods—auger and cone penetrometer—are used to sample soil thickness at 78 locations to calibrate the local value of unconstrained parameters within the hybrid model. Sensitivity analysis using the hybrid model reveals that the diffusion coefficient used in hillslope diffusion 15 modelling has the largest sensitivity among all input parameters. In addition, our results from both sampling and modeling show that, in general, the north-facing hillslope has a deeper soil layer than the south-facing hillslope. By comparing the soil thickness estimated between a machine learning approach and this hybrid model, the hybrid model provides higher accuracy and requires less sampling data. Modeling results further reveal that the south-facing hillslope has a slightly faster surface soil erosion rate and soil production rate than the north-facing hillslope, which suggests that the relatively less dense vegetation 20 cover and drier surface soils on the south-facing slopes may influence soil characteristics. With only seven parameters for calibration, this hybrid model can provide a realistic soil thickness map at other study sites by with a relatively of sampling provides divergent convergent (depositional) low - land a Light Detection (LiDAR) based digital elevation model spatial-resolution mechanistic the soil of soil convergent soil The polynomial fitting and smoothing over time show the same best fitting results among the two hillslopes, which is m and around 2 kyr, respectively. Overall, smoothing the DEM over time provides the smallest RMSD with a relatively high efficiency compared to the other two approaches. Smoothing over time provides relatively constant and stable results for both 275 hillslopes. Therefore, in this study, we use the DEM smoothed over time at year 2 kyr to calculate curvatures. However, we still use the original lidar DEM for the mechanistic model of soil thickness.

determines hillslope stability (or landslide potential), channel initiation, drainage density, and other geomorphic processes 30 (Dietrich et al., 1995). Moreover, soils hold the largest reservoir of organic carbon in the terrestrial ecosystem and function as a reservoir of other elements' accumulation, sequestration, and biogeochemical reactions (Grant and Dietrich, 2017;Tokunaga et al., 2019). An accurate soil thickness map can improve the estimation of soil water-carbon cycling for hydrologic modeling, as well as the parameterization of biogeochemical models. However, despite the importance of mapping the spatial distribution of soil thickness, it remains one of the key uncertainties in hydrologic and surface process models because complex factors 35 determine soil thickness at local scale (Jackson et al., 2017;West et al., 2012;Pelletier et al., 2018;Tesfa et al., 2009).
The soil layer, here defined as the mobile regolith layer, extends from the land surface to the top of the saprolite layer or bedrock. Soil thickness dynamics include: (1) soil transport (i.e., erosion and deposition) on the land surface and (2) soil production resulting from the bedrock-to-soil or saprolite-to-soil weathering at the base of the soil unit. The surface transport 40 is controlled by vegetation cover, topographic gradient, biogenic processes, and climate forcing; and is expressed as a diffusivelike process, either linear or nonlinear, with topographic gradient (Grant & Dietrich, 2017;Pelletier & Rasmussen, 2009;Temme & Vanwalleghem, 2016;Vanwalleghem, et al., 2013). The soil production rate is controlled by bedrock or saprolite physical, chemical, and biogeochemical processes, as well as soil thickness itself. Empirical relationships for soil production rates exist for both exponential decay (Heimsath et al., 1997) and "belly-shape" depth-dependent rate along soil depth 45 (Heimsath, et al., 2009;Pelletier & Rasmussen, 2009). Process-based geomorphologic models describe the soil thickness with mass conservation based on the balance between the two processes (Catani, et al., 2010;Dietrich et al., 1995;Heimsath et al., 2001;Heimsath et al., 1997;Nicótina, et al., 2011;Roering, et al., 1999Roering, et al., , 2001Tesfa et al., 2009). However, these studies have focused exclusively on erosional areas near ridges where the erosion rate keeps pace with soil production rate. In areas where the topography is a mixture of divergent (erosional) and convergent (depositional) zones (e.g., a low-land area)-which are 50 commonly revealed in a Light Detection and Ranging (LiDAR) based digital elevation model (DEM) for high spatialresolution modeling-these mechanistic models fail to capture the soil thickness distribution. The missing part of soil thickness estimation in convergent areas underlines the need for a hybrid approach to map soil thickness.
Many studies have used curvature as an empirical proxy for soil thickness (Patton, et al., 2018;Tesfa et al., 2009), because 55 hillslope curvature has an inverse (either linear or nonlinear) relationship with soil production rate based on mass conservation laws (Dietrich et al., 1995;Heimsath et al., 1997;Jackson et al., 2017). Soil production rates also decrease exponentially or with a 'hump' at depth (Heimsath et al., 2001(Heimsath et al., , 2009. However, studies focusing on surface curvature may not be sufficient for predicting soil thickness. For example, curvature is a secondary or less significant variable than other topographic or land cover features in predicting soil thickness (Roering, et al., 2010;Shangguan, et al., 2017;Taylor et al., 2013;Tesfa et al., 60 because curvature and other features vary under the different resolutions of DEM, it is difficult to generalize or transfer the results from one study site to another site. Still, the derivation of an empirical relationship may serve the needs for zones with 65 convergent topography where deposition is the dominant process (Patton et al., 2018).
Here we present a hybrid model that combines a process-based model with empirical relationships, field observations, and data-driven statistical analysis of model parameters and simulation results to explore the fundamental mechanisms of soil thickness and estimate the spatial variability. The advancement embodied in this hybrid model is that it generalizes the 70 calculations needed to predict soil thickness and is therefore applicable to various sites. In the methodology section, we introduce our hybrid modeling approach and relevant concepts, a sensitivity analysis of model parameters, and a machinelearning approach for soil thickness prediction to compare with the hybrid model. This model was applied with high resolution DEM (i.e., 0.5 m) at two adjacent mountainous hillslopes in the East River watershed in Colorado, U.S.A. Additional data available at this site allowed for validation of the model. We first investigate the impacts of DEM resolution on the topographic 75 curvature, an essential variable in this hybrid model, and then discuss the sensitivity of parameters to determine the importance of each parameter. The spatial maps of soil thickness, surface transport rate, and soil production rate are then presented. A Random Forest machine-learning model is used to correlate and predict soil thickness based on topographic and vegetation features and is compared with the hybrid model.

Methodology
We seek to estimate soil thickness, which is controlled by erosion and deposition on the surface and soil production at the base of soil unit. In this study, we couple two methods-one using a mass conservation law (Dietrich et al., 1995;Roering et al., 1999Roering et al., , 2001Yan et al., 2019) and another using an empirical relationship (Patton et al., 2018)-to overcome the limitations of each individual method. The mass-conservation method is suitable for a divergent topography where erosion is the dominant 85 process, while the empirical relationship is applied to convergent topography where deposition is the dominant process.
Further, we synthesize methods that are used to investigate the impacts of DEM resolution on topographic curvature, which is a key input variable for soil thickness in the empirical relationship. We also introduce the Morris method for a sensitivity analysis of parameters used in the hybrid model. A Random Forest model is applied to predict soil thickness for comparison purposes to the hybrid method. 90

Mass conservation method
The mass conservation equation of soil thickness that combines soil surface transport and soil production processes can be expressed as: where ℎ is soil thickness [ ], is time [ ], is soil production rate [ / ], $ and % is the volume flux of sediment transport 95 per unit width resulting from diffusion-driven and overland flow-driven erosion, respectively [ & / ]. The diffusion-driven soil transport rate • $ is the outcome of combining wind erosion, biogenic disturbance, soil creep, and rain drop splash. On steep slopes, the following nonlinear slope-dependent transport law is often used for topographic analysis and numerical experiments and has been successfully demonstrated by field studies and laboratory experiments (Andrews and Bucknam, 1987;Perron, 2011;Roering et al., 1999Roering et al., , 2001: 100 where $ is the diffusion coefficient [L 2 /T], and 0 = 1.25 (Roering et al., 1999) is a critical surface slope. The equation implies that when the slope ( ) is far less than 0 , the relationship between diffusion flux and slope is almost linear; when the slope approaches S 1 , $ increases rapidly.

105
The overland flow-driven soil transport rate • % is expressed as the spatial divergence of stream power (Yan et al., 2019): where % is the soil erodibility coefficient [ &+5 / ], is the steepest slope to the eight surrounding grid cells [-], and and are empirical constants for surface erosion. 4 is surface water depth [ ], which is expressed in a 2-D diffusive form (Lal, 110 1998): where " is a diffusion coefficient controlled by water depth, land-surface gradient, and Manning's coefficient (Lal, 1998;Yan et al., 2019).

115
The soil production rate function has been developed to calculate only the development of mobile soil (it excludes the immobile saprolite layer) (Heimsath et al., 1997). The calculation of the soil production rate makes use of two assumptions: (1) production rate exponentially decreases with soil thickness, and (2) production rate has a humped relationship (or a 'bellyshape') along soil depth (Dietrich et al., 1995;Heimsath et al., 2001;Pelletier and Rasmussen, 2009). Because the humped relationship is not yet well validated by field observations, in this study we rely only on the exponential model to represent 120 soil production rate: The landscape also evolves along with the soil thickness in that: 125 where is the land-surface elevation [ ], U is a constant uplift rate [L/T], and • $ and • % are the surface soil transport rates as described above.
Under the assumption of steady-state conditions, which have been observed in several mountainous area studies (Dietrich et 130 al., 1995;Vanwalleghem et al., 2013), we can solve soil thickness from the mass-conservation equation. The validity of this assumption for the site studied here will be presented in Section 3. By adopting this assumption, the soil thickness (ℎ) can be solved from the mass-conservation equation (Eqn. 1) as: One issue that arises here is that there is no finite value for soil thickness if • $ + • % < 0 where there is a continuous 135 depositing of surface soil. Therefore, we introduce Patton's method (Patton, 2019) to overcome this drawback in the following section.

An empirical approach for depositional areas
Based on field measurements among five mountainous watersheds, Patton et al. (2018) found that soil thickness has a linear relationship with curvature: 140 where ℎ C is the spatial mean value of soil thickness, and is determined by a negative linear relationship with the standard deviation of curvature. In our model, we take as an independent parameter instead of being calculated based on curvature, which adds one more degree of freedom to the model.

145
Patton's method is best adapted to convergent zones where deposition is the main process. Indeed, it can provide negative values of soil thickness at divergent zones with high negative-curvature values where erosion is the main process. However, the negative soil thickness values predicted with this method can be compensated for by using the mass-conservation method.
Therefore, we combine the two methods by applying the mass-conservation method to most erosional sites, and then applying the empirical method to all depositional sites and a potentially small portion of erosional sites where the curvature values are 150 negative but close to zero. The threshold value used to separate these two methods is represented as #"=> [L/T] ( #"=> ≤ 0):

Investigation of the DEM smoothing range for curvature
In the empirical equation, the curvature ( • ) is the only spatial variable that determines the soil thickness, and soil thickness is sensitive to the local topographic curvature (Pelletier and Rasmussen, 2009). Further, the topographic curvatures are 155 sensitive to the resolution of the DEM. If the grid size of the DEM is large, then the topography could be over-smoothed, thus underestimating the actual curvature. On the other hand, if the grid size of elevation is small, then there could be temporary pits or burrows in the topography, which can result in large local curvature values that do not represent the soil production processes. To identify a reasonable range of DEM resolution for calculating • , we explored three approaches: smoothing the DEM over space, polynomial fitting of the DEM, and smoothing the DEM over time. 160 Smoothing of the DEM over space is done by calculating the mean value of a moving window surrounding each grid cell. The smoothed DEM has the resolution reduced from the original to 3 , 5 , 7 , …, (2 + 1) , where is an integer.
The polynomial fitting of the DEM is achieved by fitting a 2 nd order polynomial to points within a specified radius and repeating this at each grid cell within the study area. For example, the elevation is fitted by = are parameters to fit the polynomial curve. Therefore, We generated additional curvature distributions by varying the patch radius used to perform the polynomial fit. The DEM smoothing over time is performed by discretizing the diffusion equation where subscript , 1, and denote the time step. A longer time period (higher value) results in a smoother topography, so the DEM is smoothed into different levels with different values. 170

Sensitivity analysis of the model parameters
Given the uncertainty of the input parameters (Table 1), we applied the Morris one-step-at-a-time (OAT) method to quantify parameter sensitivity (Campolongo, et al. 2007;Morris, 1991). The Morris method provides global sensitivity indices over the parameter space at a relatively limited computational cost (Wainwright, et al., 2014). With a set of parameters { } in that { A | = 1,2, … , }, the output from the combined model is ({ }). The Morris method first generates a randomly assigned set 175 of parameters in a discrete parameter space, and then changes one parameter at a time. (k + 1) simulations are required to complete one path, having a set of parameters, and changing all the parameters. This path is repeated for randomly generated parameter sets, with the total run being ( + 1), when the number of paths is . By normalizing parameters and uniformly spacing from 0 to 1 with ( − 1) intervals, the elementary effect is given as: where is a scaling factor, { } is a randomly selected parameter set, and the fixed increment Δ = /q2( − 1)r (Campolongo et al., 2007).
Here, we sample elementary effects of the standard deviation ( ) and absolute of the mean value (| |) to investigate the importance and nonlinearity of each parameter in the combined model. To avoid the influence of non-monotonicity when some 185 effects cancel each other out (Campolongo et al., 2007), we use the absolute mean value instead of the mean value in this study. Our work focuses on the sensitivity of how soil thickness is dictated by landscape characteristics such as the soil diffusion coefficient, soil erosion coefficient, and the saprolite-to-soil weathering capacity.

Random Forest regression
In addition to the hybrid modeling approaches, we use a machine-learning approach Random Forest (Breiman, 2001), to predict 190 the soil thickness based on topographic and land-cover features obtained from remote sensing data. The Random Forest analysis also helps to identify important predictors for soil thickness. Random Forest is capable of averaging the generated regression trees from bootstrapped subsampled data. This algorithm is nonparametric by assuming no specific data distribution ahead of time (Hastie, 2001). We use the 'sklearn' Python package for this study.

195
The Random Forest algorithm input dataset comprises topographic and land cover features, including aspect, gradient, curvature, topographical flow accumulation, normalized difference vegetation index (NDVI), canopy water content, topographical position index (TPI) etc. (Brodrick et al., 2020;Chadwick et al., 2020;Goulden et al., 2020) at 10 m resolution and smoothing with a moving window at different sizes (e.g., 30 , and 50 , and 90 ) (See Table S1 for the full list of features). Because the field observations have limited sampling points (78), we use the leave-one-out cross-validation method 200 (Efron, 1982) where the number of folds equals the number of instances in the dataset. The Random Forest algorithm is applied once for each instance, using all other instances as a training set and using the selected instance as a single-item test set.

Study site description and data sampling
Our study site is within the East River watershed, which represents a typical headwater mountainous watershed in the Upper Colorado River Basin in Colorado. The East River watershed is a testbed aiming to improve predictions of hydrology-driven 205 biogeochemical activities. There is a focus on this watershed because it hosts a wide spectrum of vegetation cover, features various hydrologic and geomorphologic processes. It is a headwater to the Colorado River that supplies 1 in 10 Americans with water for municipal usage and irrigation for over 2.2 million ha (Hubbard et al., 2018). This ~126.5 km & watershed has a continental climate with long, cold winters and short, cool summers, with monthly mean temperature ranging from -9.2 to 9.8 o C. The annual precipitation, based on three years of monitoring (2015)(2016)(2017)(2018), is between 659 and 750 mm, the majority of 210 which falls as snow, followed by mid-to late-summer monsoonal rainfall.
We consider opposite facing hillslopes to determine whether the variability in meteorological forcing drives differences in the soil characteristics (Pelletier et al., 2018). Our study sites (Fig. 1a) include northeast-and southwest-facing hillslopes (for simplicity, referred to as north-facing and south-facing hillslopes) connected to a floodplain (Fig. S1) spanning a total of 0.4 215 km & and located with an elevation range from ~2755 to 2922 m at ~38.93 O N latitude and ~106.95 O W longitude. The differences in slope and aspect drive differences in hillslope energy balance, snow-melt timing, and vegetation seasonal dynamics. The south-facing slope shows a thinner snowpack resulting from episodic melt events scattered throughout the winter and spring; whereas the north-facing slope develops a thicker seasonal snowpack that primarily melts during a large snowmelt event occurring over several weeks in spring (Hinckley et al., 2014). Several remote sensing products available at 220 this site were used in this study including, 0.5 m lidar DEM, topographic and vegetation features (Table S1).
Soil thickness was measured at 78 locations across the north-and south-facing slopes. Many studies dig soil pits or use augers to distinguish the contact layer between mobile and immobile regolith layers (Catani et al., 2010;Heimsath et al., 1997;Patton et al., 2018;Pelletier and Rasmussen, 2009), and one study used a sharpened copper-coated steel rod graduated at 5 cm intervals 225 vertically into the ground until refusal (Tesfa et al., 2009). Here, we chose two independent methods to measure the soil thickness of the two hillslopes. We used an auger to drill and sampled 78 points within the delineated hillslopes, and used a dynamic cone penetrometer technique (CPT) (Vanags, et al., 2004) to measure two transects along each side of the hillslopes, sampling 54 locations in total (Fig. 1a).

230
The bottom soil boundary was estimated in two ways. From the auger, the transition zone from soil to the saprolite or bedrock is based on the material size and color of retrieved samples. Based on a number of surface soil samples being analysed, the soil in this area is composed of 15% sand, 47% silt, and 38% clay based on field campaigns in 2019. The weathered shale (or saprolite) is between pebble-to cobble-size with light brown color. When the auger reaches the bedrock shale, it cannot penetrate. From the CPT-inferred vertical profile of the soil resistance, we used the sharp transition from low to high values 235 to separate the two layers. Figures 1b-1e show the relationship between soil thickness estimated from auger, CPT, and local elevation. There is a high variation in soil thickness from local to hillslope scales. To fully take advantage of all the sampling data, we used auger data to fit values for CPT (Fig. S3) since more locations are sampled from auger drilling than CPT. The CPT and auger data are mostly in agreement. However, for soil thicknesses less than ~0.5 m, the CPT data are slightly higher than the auger, and for soil thickness larger than ~0.5 m, the CPT data are slightly lower.

4 Results and Discussion
In this section, we first investigate the scaling issues from DEM resolution with curvature and then discuss the sensitivity of all the parameters used in this hybrid model. Further, we show both the soil thickness predictions based on the hybrid model with the optimal curvature and the soil thickness results from the Random Forest machine-learning approach. Finally, we discuss the relationship between soil thickness, surface transport rate, and weathering rate as determined by the hybrid model. By coupling soil thickness with landscape evolution, we found that the soil thickness reaches a dynamic steady-state after approximately 20 kyr at this study site (Fig. S2), which is consistent with other studies in mountainous areas (Dietrich et al., 1995;Vanwalleghem et al., 2013). This implies that the current soil thickness in the East River Watershed may already have reached a steady-state condition.

Model parameterization of curvature with different smoothing techniques of DEM 255
The topographic curvature is the key variable for estimating the soil thickness as introduced in the empirical approach.
However, curvature is an inherently resolution-dependent topographic feature that is derived from a DEM. To investigate the optimal resolution of the DEM for calculating curvature, we tested three approaches as explained above: smoothing the DEM over space, polynomial fitting of the DEM, and smoothing the DEM over time. For smoothing over space, the original 0.5 m DEM was resampled to 1.5 m, 2.5 m, …, and 13.5 m by utilizing the mean elevation of the surrounding adjacent cells. For the 260 polynomial fitting, the diameter is also chosen over the same range and intervals as the way of smoothing the DEM. For smoothing over time, we use a fixed diffusion coefficient, $ , and the evolution period of the topography is taken as 0 kyr, 0.25 kyr, 0.5 kyr, …, 3.25 kyr. We compare the root-mean-square deviation (RMSD) between the sampling results and simulation results of soil thickness as an indicator of the optimal resolution for calculating the curvature (Fig. 2). The sets of parameters used for each spatial or temporal resolution are determined by comparing with the field sampling data, with each 265 parameter increasing linearly from the minimum value to the maximum (the same data sets as in the Morris OAT method).
We chose the set of parameters which gives the smallest RMSD.
Different smoothing methods and study sites corresponds to their own optimal DEM resolution for curvature (Fig. 2). For example, in this study the north-facing hillslope shows best fitting with 4.5 m DEM smoothing over space, but the south-facing 270 hillslope corresponds to 8.5 m. Other studies that calculate curvature with various spatial smoothing constraints show the highest accuracy in model predictions with smoothing ranging between 3m and 10 m (Patton et al., 2018;Tesfa et al., 2009).
The polynomial fitting and smoothing over time show the same best fitting results among the two hillslopes, which is 8.5 m and around 2 kyr, respectively. Overall, smoothing the DEM over time provides the smallest RMSD with a relatively high efficiency compared to the other two approaches. Smoothing over time provides relatively constant and stable results for both 275 hillslopes. Therefore, in this study, we use the DEM smoothed over time at year 2 kyr to calculate curvatures. However, we still use the original lidar DEM for the mechanistic model of soil thickness.

Sensitivity analysis
We apply the Morris OAT method to investigate the sensitivity of the seven parameters (Table 1) in the hybrid model. We compared the "absolute of the mean elementary effect," | |, and the standard deviation (Fig. 3). Each dot is a value from one of the sampling sites. Each color represents one parameter, but at different sampling locations. In general, the parameters in the mass-conservation model have more significant impact on soil thickness than the parameters in the empirical model. 285 The diffusion coefficient, $ , is the most important factor and thus should be carefully calibrated. It represents the soil diffusive-like process such as soil creeping and biogenic activities, which implies that the diffusive process is the most important transport mechanism for the soil mantled hillslopes rather than the soil erosion from overland flow (Dietrich et al., 1995;Nicótina et al., 2011;Roering et al., 1999Roering et al., , 2001.

290
The potential soil production rate, represented as J ( J = : , : $ * 0;%< $ ; ), plays a small role compared to the normalized soil depth ℎ ; in the mechanistic model. The two parameters from the empirical method, and ℎ C , are used for soil depositional areas. The relationship of is nearly linear (since s is close to zero), but when the soil thickness reaches an upper limit (2 m in our model), this causes a nonlinear increase in soil thickness and hence a rapid increase of .

Hybrid data-model soil thicknesses predictions and their comparison to measurements
A spatial map of soil thickness that compares modelling and observation results is shown in Fig. 4. In general, the soil thickness of the north-facing hillslopes is greater and has higher variation than that of the south-facing hillslopes. Specifically, based on the sampling data, the soil thickness of the south-facing side ranges from about 0.2 m to 1 m, with a mean value of 0.55 m, 305 while on the north-facing side, the soil thickness ranges from 0.15 m to 1.5 m, with a mean value of 0.67 m. We use the sampling data from both auger and the CPT to calibrate the seven parameters (Table 1) for the south-facing and north-facing hillslopes separately (Fig. S4). The calibration follows the same sets of parameters generated in the Morris OAT method by increasing each parameter linearly from the minimum value to the maximum in a constant interval. The RMSEs of the southfacing and north-facing hillslopes are 0.20 m and 0.195 m, respectively. One of the reasons for the mismatch between modeling 310 and field observations (Fig. 4b) could be that parts of the hillslopes are on terraces, which may have fluvial deposits on ancient floodplains before they into terraces (Yan et al., 2018).
The difference in soil thickness between the two hillslopes is evident, controlled by insolation due to the topographic aspect.
The air temperatures and potential evapotranspiration rates produce significantly different microclimates that determine the 315 structure of different ecosystems and surface process regimes. Weathered shale under the soil layer appears to be mechanically weaker on the north-facing slope because of a higher thickness of the saprolite layers (Fig. 4b), which results in less resistance during excavations in the field sampling. This implies that the thicker soil depth results from a higher soil water content associated with a longer snowmelt period in the north-facing hillslope. The thicker soil depth in turn provides a higher water-storage capability, higher concentration of fine particles, and more biomass, which leads to a positive feedback (Pelletier et 320 al., 2018;Roering et al., 2010).

Predictive value of landscape characteristics for estimating soil thickness
We evaluate the correlation between soil thickness from sampling and other landscape characteristics obtained from remote sensing data (Fig. 5). Among about 18 topographic features with various spatial resolutions, we generated ~50 topographic matrices (Table S1). The top five topographic matrices, correlating with soil thickness higher than 25%, from high to low are 330 topographic position index (TPI), curvature, slope degree, topographic wetness index (TWI), and elevation. Other factors such as NDVI (normalized difference vegetation index), leaf mass area, and canopy liquid water do not show obvious correlation with soil thickness. This is consistent with other studies that commonly use the environmental variables such as TWI, elevation, and curvature as the most highly correlated variables in the geostatistical interpolation of soil thickness (Hengl, et al., 2004;Kuriakose, et al., 2009;Shangguan et al., 2017;Taylor et al., 2013). Among the five highest correlated features, TPI and 335 curvature have the highest correlation (-0.69 and -0.87, respectively), which implies that the local relief of ridges and valleys are scalable with the size of the corresponding local curvature. Since we have a limited number of soil samples (78), we use the five most highly correlated features as a collection of 345 topographic metrics to perform the Random Forest modeling. We use the Random Forest model with leave-one-out cross validation to predict the soil thickness and compare with the hybrid modeling results (Fig. 6). The result shows that the hybrid model (RMSE = 0.196 m) outperforms the Random Forest model (RMSE = 0.225 m) by ~13%. The comparable performance between Random Forest and the hybrid model suggests that (1) the correlations of soil thickness with topographic metrics are the major driving factors, and (2) the hybrid model is applicable to other sites given similar soil types and topography, without 350 requiring many datasets given the application of process-based laws. To improve the results from machine learning, one may need to collect more data from additional sampling points. However, the advantage of this hybrid model is that there are only seven parameters to calibrate, and fewer sampling points are required. The extension to other watersheds is also easier with the hybrid approach. Finally, we note that the hybrid modeling approach not only provides the soil thickness distribution, but also other outputs from the model, including the surface soil transport and soil production rates, that may be useful (discussed 355 below).

Relationship between soil thickness, soil surface transport rates, and soil production rates in two aspects
One of the advantages of this hybrid model compared with machine-learning models is that we can obtain the spatial map of surface soil transport rates (Fig. 7) and soil production rates (Fig. 8). The probability density functions (PDFs) show that in general, the south-facing hillslope has a thinner soil thickness and a faster erosion rate than the north-facing hillslope (Fig. 7b).
On the south-facing hillslope, where solar radiation is sufficient, soil moisture is a limiting factor that controls vegetation cover 365 of critical zone processes. In contrast, on the north-facing hillslope where solar radiation is limited, the energy becomes a limiting factor in the coupled hydro-bio-geomorphological processes (Pelletier et al., 2018). The higher potential evaporation of south-facing hillslopes results in lower mean soil moisture and thinner vegetation cover, thus triggering feedbacks that result in loss of soil, fewer steep slopes, but a higher surface soil erosion rate. Moreover, on the south-facing hillslope, the lower interception of rainfall due to thinner vegetation cover also increases sediment mobilization by raindrop impacts. 370

375
The south-facing hillslope corresponds to a higher actual soil production rate than the north-facing hillslope, which is consistent with transport rates in that the south-facing hillslope has a higher erosion rate. The actual soil production rate is controlled by the soil thickness and potential production rate, ; (Eqn. 5). A thinner soil layer (south-facing) results in a faster actual soil production rate. A high evapotranspiration and rapid vertical infiltration result in lower ; . The potential production rate on 380 the south-facing hillslope is slower than the north-facing hillslope (Table 1) due to the microclimate differences between the two hillslopes. On south-facing hillslopes, more of the incoming precipitation is lost to evapotranspiration, which results in less water available for runoff and infiltration . In addition, the south-facing slopes experience brief periods of rapid vertical transport following snowmelt events and are drier overall than north-facing slopes (Hinckley et al., 2014). 385 Pelletier et al. (2013) uses an energy-based variable, Effective Energy and Mass Transfer (EEMT), which is a function of precipitation, temperature, and vapor pressure deficit, to quantify the potential rate of bedrock breakdown into soils. Their study also suggests that under comparable climatic conditions north-facing hillslopes have a higher EEMT, leading to higher ; . However, the thicker soil thickness offsets the impact of ; in that thicker soil thickness results in slower actual soil https://doi.org/10.5194/esurf-2020-110 Preprint. Discussion started: 18 January 2021 c Author(s) 2021. CC BY 4.0 License. production rate. This is a possible reason that explains why the north-facing hillslope has slower actual soil production rate 390 even though it has a higher ; .

Conclusion
Soil thickness plays a central role in the feedbacks among surface-subsurface water flow, vegetation, soil production, drainage density, and topography, and these in turn control soil thickness. In this study, we developed a data-driven hybrid model approach to predict the spatial distribution of soil thickness. The hybrid model that we introduced in this study overcomes the 400 drawbacks of both mass-conservation laws and empirical relationships. To the authors' knowledge, this is the first study to use a hybrid approach to estimate soil thickness.
Our results show that this hybrid model provides slightly better accuracy than the Random Forest model (by ~13%) on soil thickness estimation. According to both the hybrid and Random Forest models, soil thickness is more strongly controlled by 405 topographic metrics than vegetational features. The sensitivity analysis of the input parameters (seven in total) show that the diffusion coefficient of hillslope erosion is the most sensitive parameter. We found that smoothing DEM over time has a high efficiency and provides the optimal curvatures that result in the least error.
This hybrid model is a flexible, generally applicable approach to predicting soil thicknesses. The hybrid model, with only 410 seven parameters for calibration, can provide a relatively realistic soil thickness map at other study sites making use of a relatively small number of samples. It can also provide additional output as compared to machine learning algorithms, including surface soil transport and soil production rates.
Based on field observation and the hybrid model simulation, the north-facing hillslope promotes deeper soil depth than the 415 south-facing hillslope as a result of the different insolation at different aspects. The model analysis suggests that the southfacing hillslope has a slightly faster actual surface transport rate and actual soil production rate than the north-facing hillslope.
The potential soil production rate is higher on the north-facing hillslope, caused by relatively denser vegetation cover, less solar radiation, and wetter surface soil material as fundamentally controlled by aspect. One should note that soil production rate here is different from bedrock weathering rate, which is controlled by water table dynamics . 420 The limitation of the hybrid modeling approach developed here is that it would fail in alluvial depositional sites where topography is controlled primarily by flooding events (Yan et al., 2018) and strong human intervention landscapes where the surface topography is intensively reshaped for farming and other purposes (Kuriakose et al., 2009;Yan et al., 2020). Integrating process-based modeling, inverse modeling, and statistical analysis provides a thorough understanding of the fundamental 425 mechanisms for soil thickness prediction. Although the example applications in this paper are at two hillslopes, there is no limitation on the size of the watershed that can be analyzed using this model framework.

Author contributions 435
QY performed the numerical and statistical calculations and produced the first draft of this paper. BD, SU, and QY conducted field sampling work. HW and BD supervised this study, contributed to the interpretation of this results, and reviewed several versions of the manuscript. HW also contributed significantly on the sensitivity analysis. JK contributed to the discussion and coding of curvature study and revised the manuscript. CIS supervised this study, particularly on the soil thickness evolution study, and revised the manuscript. SSH supervised this study and revised the manuscript. NF provided the remote sensing data 440 and revised the manuscript.