**Research article**
11 Oct 2021

**Research article** | 11 Oct 2021

# A hybrid data–model approach to map soil thickness in mountain hillslopes

Qina Yan Haruko Wainwright Baptiste Dafflon Sebastian Uhlemann Carl I. Steefel Nicola Falco Jeffrey Kwang and Susan S. Hubbard

^{1},

^{1},

^{1},

^{1},

^{1},

^{1},

^{2},

^{1}

**Qina Yan et al.**Qina Yan Haruko Wainwright Baptiste Dafflon Sebastian Uhlemann Carl I. Steefel Nicola Falco Jeffrey Kwang and Susan S. Hubbard

^{1},

^{1},

^{1},

^{1},

^{1},

^{1},

^{2},

^{1}

^{1}Earth and Environmental Science Area, Lawrence Berkeley National Laboratory, Berkeley, CA, USA^{2}Department of Geosciences, University of Massachusetts Amherst, Amherst, MA, USA

^{1}Earth and Environmental Science Area, Lawrence Berkeley National Laboratory, Berkeley, CA, USA^{2}Department of Geosciences, University of Massachusetts Amherst, Amherst, MA, USA

**Correspondence**: Qina Yan (qinayan@lbl.gov)

**Correspondence**: Qina Yan (qinayan@lbl.gov)

Received: 26 Dec 2020 – Discussion started: 18 Jan 2021 – Revised: 30 Jul 2021 – Accepted: 31 Aug 2021 – Published: 11 Oct 2021

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 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 aspects of hillslopes (southwest- and northeast-facing, respectively) in the East River watershed in Colorado. 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 modeling has the largest sensitivity among all input parameters. In addition, our results from both sampling and modeling show that, in general, the northeast-facing hillslope has a deeper soil layer than the southwest-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 southwest-facing hillslope has a slightly faster surface soil erosion rate and soil production rate than the northeast-facing hillslope, which suggests that the relatively less dense vegetation cover and drier surface soils on the southwest-facing slopes influence soil properties. With seven parameters in total for calibration, this hybrid model can provide a realistic soil thickness map with a relatively small amount of sampling dataset comparing to machine-learning approach. Integrating process-based modeling and statistical analysis not only provides a thorough understanding of the fundamental mechanisms for soil thickness prediction but also integrates the strengths of both statistical approaches and process-based modeling approaches.

- Article
(8122 KB) -
Supplement
(5395 KB) - BibTeX
- EndNote

The soil layer is an element of the critical zone where water, carbon, nitrogen, and other elements exchange between air and plants and the subsurface. The thickness of a soil layer regulates the hydrologic response, including surface and base flow runoff, water partitioning, evapotranspiration, plant-available water, and water and nutrient residence time (Fan et al., 2019). It also determines hillslope stability (or landslide potential), channel initiation, drainage density, and other geomorphic processes (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). Therefore, an accurate soil thickness map can improve the estimation of water, carbon, nitrogen, and other elements dynamics for hydrologic and biogeochemical modeling (Carvalhais et al., 2014; Fan et al., 2019; Li et al., 2020; Patton et al., 2019; Pelletier et al., 2016). However, mapping soil thickness remains one of the key uncertainties in land surface process models because of the complexity of factors that affect soil thickness (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 (if there is no saprolite layer). Process-based geomorphologic models describe the soil thickness with mass conservation based on the balance between (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 bottom of the soil layer (Catani et al., 2010; Dietrich et al., 1995; Heimsath et al., 2001, 1997; Nicótina et al., 2011; Roering et al., 1999, 2001; Tesfa et al., 2009). These two processes are controlled by vegetation cover, topographic gradient, biogenic processes, and climate forcing. The surface soil transport is expressed as a diffusive-like process, either linear or nonlinear, with topographic gradient (Grant and Dietrich, 2017; Pelletier and Rasmussen, 2009; Temme and Vanwalleghem, 2016; Vanwalleghem et al., 2013). The soil production rate is usually expressed as relationships of exponential decay (Heimsath et al., 1997) or “bell-shaped” along soil depth (Heimsath et al., 2009; Pelletier and Rasmussen, 2009). However, previous studies have focused exclusively on erosional areas near ridges where the erosion rate keeps pace with soil production rate. In areas where the topography has convergent (mostly depositional) zones, which are commonly revealed in a Light Detection and Ranging (lidar) digital elevation model (DEM) for high spatial resolution modeling, these mechanistic models fail to capture the full soil thickness distribution. The missing part of soil thickness estimation in convergent areas underlines the need for a hybrid approach that couples mechanistic and empirical methods to map soil thickness.

Many studies have used curvature – defined as the second-order derivative of elevation – as a proxy for soil thickness (Patton et al., 2018; Tesfa et al., 2009) because 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). However, some studies show that 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., 2009). In addition, curvature is proven to show a weak correlation to soil thickness in catchments with high curvature variability (Patton et al., 2018). One of the reasons that curvature fails to be the dominant explanatory variable could be that curvature is sensitive to the DEM resolution even among adjacent hillslopes (Patton et al., 2018; Pelletier and Rasmussen, 2009). Moreover, because the relationship between curvature and other features are difficult to generalize or transfer from one study site to another site, it is not ideal to rely on curvature alone to estimate soil thickness. Despite this, the derivation of an empirical relationship may serve the needs for partial areas in a landscape such as zones with convergent topography (Patton et al., 2018).

Here we present a hybrid model that combines a process-based model with empirical relationships 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 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 such as curvature calculation with different DEM smoothing methods, sensitivity analysis of model parameters, and a machine-learning approach as a comparison with the hybrid model. This model was applied with a high-resolution DEM (i.e., 0.5 m) at two adjacent mountainous hillslopes in the East River watershed in Colorado, USA. Data from field observations at this site allow for calibration of the model. We first investigate the impacts of DEM resolution on the topographic 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.

## 2.1 Hybrid modeling approach

We introduce a hybrid approach that couples two methods – a mass conservation law (Dietrich et al., 1995; Roering et al., 1999, 2001; Yan et al., 2019) and 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 process, while the empirical relationship is applied to a convergent topography where deposition is the dominant process. In this hybrid model, seven parameters (Table 1) need to be calibrated for a specific hillslope area. In addition, 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. A diagram that highlights the workflow is shown in Fig. S1.

^{a} We defined *h* as the distance along the normal direction to the land surface, which gives ${e}^{-h\mathrm{cos}\mathit{\theta}/{h}_{o}}$, where *θ* is the slope of the land surface in degrees (Pelletier and Rasmussen, 2009). In this case, *h*_{o} is adjusted to include cos *θ* when referring to other literature. ^{b} Heimsath et al. (1997). ^{c} Heimsath et al. (2000).
^{d} Heimsath et al. (2005). ^{e} Dietrich et al. (1995).
^{f} Patton et al. (2019). ^{g} Fernandes and Dietrich (1997). ^{h} Kilinc and Richardson (1973).

### 2.1.1 Mass conservation method

The mass conservation equation of soil thickness that combines soil surface transport and soil production processes can be expressed as follows:

where *h* is soil thickness [L], *t* is time [T], *P* is soil production
rate [L/T], *q*_{d} and *q*_{s} is the volume flux of sediment transport per unit width resulting from diffusion-driven and overland flow-driven erosion, respectively [L^{2}/T]. The diffusion-driven soil transport rate ∇⋅*q*_{d} 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., 1999, 2001):

where *K*_{d} is the diffusion coefficient [${L}^{\mathrm{2}}/T$] and *S*_{c}=1.25
(Roering et al., 1999) is a critical surface
slope. The equation implies that when the slope (∇*η*) is far less
than *S*_{c}, the relationship between diffusion flux and slope is almost
linear; when the slope approaches *S*_{c}, *q*_{d}
increases rapidly.

The overland flow-driven soil transport rate ∇⋅*q*_{s} is
expressed as the spatial divergence of stream power (Yan et al., 2019):

where *K*_{s} is the soil erodibility coefficient [${L}^{\mathrm{2}-\mathit{\alpha}}/T$], *S* is the slope along-flow direction [–], and *β* is an empirical constant for surface erosion, where *β*=1.68 (Yan et al., 2019). *H*_{w} is surface water depth [*L*], which is expressed in a 2-D diffusive form (Lal, 1998):

where *D*_{h} is a diffusion coefficient controlled by water depth, land surface gradient, and Manning's coefficient (Lal, 1998; Yan et al., 2019).

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 “belly-shape”) 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 soil production rate:

where *ρ*_{r} and *ρ*_{s} are bedrock and soil bulk density,
respectively, [$M/{L}^{\mathrm{3}}$]. The mean value of the soil bulk density at
sampling sites is 0.948 [g/cm^{3}], and the bedrock bulk density for
weathered shale from our deep samples is estimated to be 1.26 [g/cm^{3}].
*P*_{o} [$L/T$] and *h*_{o} [*L*] are empirical constants for soil production,
and *θ* is the slope of the land surface in degree [–] because the
direction of soil thickness is perpendicular to the land surface.

Under the assumption of steady-state conditions, which have been observed in
several mountainous area studies (Dietrich et 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 Sect. 3. By adopting this assumption,
the soil thickness (*h*) can be solved from the mass conservation equation
(Eq. 1) as follows:

The soil thickness value *h* can be directly solved here because ∇⋅*q*_{d} and ∇⋅*q*_{s} are independent of soil thickness.
One issue that arises here is that there is no real number for soil
thickness if $\mathrm{\nabla}\cdot {q}_{\mathrm{d}}+\mathrm{\nabla}\cdot {q}_{\mathrm{s}}<\mathrm{0}$ where there is a
continuous depositing of surface soil. Therefore, we introduce Patton's
method (Patton et al., 2018) in depositional
areas to overcome this drawback in the following section.

### 2.1.2 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:

where $\stackrel{\mathrm{\u203e}}{h}$ is the spatial mean value of soil thickness and *a* is a
constant value which is determined by having a negative linear relationship
with the standard deviation of curvature. In our model, we take *a* as an
independent parameter instead of being calculated based on curvature, which
adds one more degree of freedom to the model.

### 2.1.3 Investigation of the lidar DEM smoothing range for curvature

In the empirical equation (Eq. 7), the curvature (∇⋅∇*η*) is the only spatial variable that determines the soil thickness.
Further, the topographic curvatures are very 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 to reproduce the DEM by (1) smoothing
the DEM over space, (2) polynomial fitting of the DEM, and (3) smoothing the
DEM over time. The smoothed DEM is for calculating curvature in the imperial
method only (Eq. 7), and the rest of all other calculations (Eqs. 1–6)
still use the original lidar DEM as the input.

Smoothing of the DEM over space is done by replacing the value of a 2-D grid
cell with the mean value of its surrounding neighbors. The size of its
neighbor cells follows 3Δ*x* (8 neighbors), 5Δ*x* (24
neighbors), 7Δ*x* (48 neighbors), … , (2*N*+1)Δ*x* (${\left(\mathrm{2}N+\mathrm{1}\right)}^{\mathrm{2}}-\mathrm{1}$ neighbors), respectively, where Δ*x* is the original resolution (i.e., 0.5 m) and *N* is an integer; following this, a moving window replaces the value of every single 2-D grid cell in the 0.5 m lidar. The polynomial fitting of the DEM is achieved by fitting a second-order polynomial to grid cells within a specified radius and repeating this at each grid cell within the study area. For example, the elevation is fitted by $\mathit{\eta}={a}_{\mathrm{1}}{x}^{\mathrm{2}}{y}^{\mathrm{2}}+{b}_{\mathrm{1}}{x}^{\mathrm{2}}y+{b}_{\mathrm{2}}x{y}^{\mathrm{2}}+{c}_{\mathrm{1}}{x}^{\mathrm{2}}+{c}_{\mathrm{2}}{y}^{\mathrm{2}}+\mathrm{d}xy+{e}_{\mathrm{1}}x+{e}_{\mathrm{2}}y+f$,
where *a*_{1}, *b*_{1}*b*_{2}*c*_{1}, *c*_{2}*d*, *e*_{1}, *e*_{2}, and *f* are
parameters to fit the polynomial curve. The curvature value is $\mathrm{\nabla}\cdot \mathrm{\nabla}\mathit{\eta}=\mathrm{2}{c}_{\mathrm{1}}+\mathrm{2}{c}_{\mathrm{2}}$. The DEM smoothing over time is
performed by discretizing the diffusion equation $\frac{\partial \mathit{\eta}}{\partial t}={K}_{\mathrm{d}}\mathrm{\nabla}\cdot \mathrm{\nabla}\mathit{\eta}$, which gives ${\mathit{\eta}}_{n}={\mathit{\eta}}_{\mathrm{1}}+\mathrm{\Delta}t\times {K}_{\mathrm{d}}{\mathrm{\Sigma}}_{i=\mathrm{1}}^{n-\mathrm{1}}\mathrm{\nabla}\cdot \mathrm{\nabla}{\mathit{\eta}}_{i}$, where subscript *i*, 1, and *n* denote the time step. A
longer time period (higher *n* value) results in a smoother topography, so
the DEM is smoothed into different levels with different *n* values. To the
authors' knowledge, the smoothing over time approach is new and original in this
study.

### 2.1.4 Combine the mass conservation method with the empirical method

For the mass conservation equation, the steady-state assumption is used for
the soil thickness estimation in the assumption that the soil production
balances the physical erosion, as used in other studies (Pelletier and
Rasmussen, 2009; Pelletier et al., 2016; Dietrich et al., 1995). Therefore,
the mass conservation method with the steady-state assumption can be used to
solve the soil thickness at erosional sites but has a limitation at
depositional sites (Eq. 6, Dietrich et al., 1995). Patton's method is better
adapted to depositional sites. However, it can provide negative values of
soil thickness at erosional sites where there are zones with high
negative-curvature values. In addition, in a low gradient and divergent area, if the soil transport rate is assumed as a linear relationship with curvature
(i.e., $\mathrm{\nabla}\cdot {q}_{\mathrm{d}}=-{K}_{\mathrm{d}}\mathrm{\nabla}\cdot \mathrm{\nabla}\mathit{\eta}$, and $\mathrm{\nabla}\cdot {q}_{\mathrm{s}}=\mathrm{0.0}$), then Eq. (6) can be further simplified in that the soil thickness has a natural logarithm relationship with curvature (i.e., $h=-m\times \mathrm{ln}\left(\mathrm{\nabla}\cdot \mathrm{\nabla}\mathit{\eta}\right)+C$, where *m* and *C* are constant parameters that can be calibrated from field sampling data. However, Patton's method (Patton et al., 2018) always assumes a linear relationship between soil thickness and curvature. This may be why his empirical relationship does not work very well in the erosional areas. However, this can be compensated for by using the mass conservation method.

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 negative but close to zero. The threshold value used to separate
these two methods is represented as *E*_{thre} [$L/T$] (*E*_{thre}≥0):

where the threshold, *E*_{thre}, is a condition of the soil erosion rate and
equal to or larger than zero value. If $\mathrm{\nabla}\cdot {q}_{\mathrm{d}}+\mathrm{\nabla}\cdot {q}_{\mathrm{s}}>{E}_{\mathrm{thre}}$ at a 2-D grid cell, then this cell must be an erosional
site; if $\mathrm{\nabla}\cdot {q}_{\mathrm{d}}+\mathrm{\nabla}\cdot {q}_{\mathrm{s}}\le {E}_{\mathrm{thre}}$, then this
cell can be either a depositional site (if $\mathrm{\nabla}\cdot {q}_{\mathrm{d}}+\mathrm{\nabla}\cdot {q}_{\mathrm{s}}\le \mathrm{0})$ or a slightly erosional site (if $\mathrm{0}<\mathrm{\nabla}\cdot {q}_{\mathrm{d}}+\mathrm{\nabla}\cdot {q}_{\mathrm{s}}\le {E}_{\mathrm{thre}}$). In most of areas, a divergent
topography corresponds to erosional areas and vice versa for depositional
areas. However, here we use the transport rate instead of the curvature as the
criteria to choose between the two methods because there are possibly sites
that are convergent but erosional where overland flow erosion is stronger
than the diffusive deposition. In other words, areas where $\mathrm{\nabla}\cdot {q}_{\mathrm{d}}+\mathrm{\nabla}\cdot {q}_{\mathrm{s}}\le \mathrm{0}$, it must be a convergent area and
undergoing deposition, but if it is a convergent area, it is unnecessary to
be $\mathrm{\nabla}\cdot {q}_{\mathrm{d}}+\mathrm{\nabla}\cdot {q}_{\mathrm{s}}\le \mathrm{0}$. In addition, we assign
*E*_{thre}≥0 instead of equal to zero, aiming to provide more
flexibility to switch between the two methods. Overall, *E*_{thre} is
supposed to be very close to zero. We conducted a grid search to calibrate
the parameters and discuss the corresponding posterior distribution of
parameters in Sect. 4.3.

## 2.2 Sensitivity analysis of the model parameters

We introduce the Morris one-step-at-a-time (OAT) method for a sensitivity
analysis of parameters used in the hybrid model. Given the uncertainty of
the input parameters (Table 1), we applied the Morris 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 *k* parameters {x} in that $\left\{{x}_{i}|i=\mathrm{1},\mathrm{2},\mathrm{\dots},k\right\}$, the output from the combined model is
*f*(**x**). The Morris method first
generates a randomly assigned set of parameters in a discrete parameter
space and then changes one parameter at a time.
The (*k*+1) simulations are required to complete one path, which has one fixed parameter but changes all other parameters. This path is repeated for randomly generated parameter sets,
with the total run being *r*(*k*+1) where the number of paths is *r*. By
normalizing parameters and uniformly spacing from 0 to 1 with (*p*−1)
intervals, the elementary effect is given as follows:

where *τ* is a scaling factor, {x} is a
randomly selected parameter set, and the fixed increment $\mathrm{\Delta}=p/\left(\mathrm{2}\left(p-\mathrm{1}\right)\right)$
(Campolongo et al., 2007).

Here, we use iTOUGH2 (Wainwright et al., 2014)
to generate sets of parameters. Then we sample *r*(=20) elementary
effects for each parameter. The standard deviation (*σ*) and absolute
of the mean value ($\left|\mathit{\mu}\right|$) of elementary effect of each
parameter can be used to investigate the importance and nonlinearity in the
combined model. To avoid the influence of non-monotonicity when some 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 hillslope diffusion coefficient,
soil erosion coefficient, and the saprolite-to-soil weathering capacity.

## 2.3 Random forest regression

In addition to the hybrid modeling approaches, we use a machine-learning approach random forest (Breiman, 2001) to predict the soil thickness for comparison purposes to the hybrid method. The features to train the model are topographic and land cover data obtained from remote sensing. 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.

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, 50, and 90 m) (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 (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.

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 test bed aiming to improve predictions of hydrology-driven biogeochemical activities. There is a focus on this watershed because it hosts a wide spectrum of vegetation cover and
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 22 000 km^{2} (Hubbard et al., 2018). This ∼126.5 km^{2} watershed has a continental climate with long, cold winters and short, cool summers, with monthly mean temperature ranging from −9.2 to 9.8 ^{∘}C. The annual precipitation, based on 3 years of monitoring (2015–2018), is between 659 and 750 mm, the majority of 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 thickness
(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. S2) spanning a total of 0.4 km^{2} and located with an elevation range from ∼2755 to 2922 m at ∼38.93^{∘} N latitude and ∼106.95^{∘} W longitude. The differences in slope and aspect drive
differences in hillslope energy balance, snowmelt 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). Based on
the field campaigns in 2019, soils in this area are composed of 15 % sand,
47 % silt, and 38 % clay. Several remote sensing products available at
this site were used in this study including 0.5 m lidar DEM and topographic
and vegetation features (Table S1).

The last glacial advance and retreat in the upper Colorado River basin is dated between 16.1 and 20.8 ka (Brugger, 2010). Glacial deposits are mapped at many locations throughout the watershed (Gaskill, 1991; Fig. S3), but they are rather isolated and have a limited spatial extent, including in the area analyzed in this study. A former study at the same site analyzed 40 hand-augured soil cores and showed progressive changes in color and texture among soil, weathered zone, and unweathered bedrock with depth (Wan et al., 2019). Further, among the total of five wells drilled at the site, none of them reported the presence of a glacial deposit (Tokunaga et al., 2019; Wan et al., 2019). It is likely that the glacial legacy scoured the valley to bedrock and that the glacial retreat reset the “clock” for soil formation mostly via in situ bedrock weathering and (to a lesser extent) via colluvial deposition.

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). Here, we chose two independent methods to measure the soil thickness of the two hillslopes. We used an auger to drill and sample 78 points within the two 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).

At this study site, we used and compared both auger and CPT measurements to estimate soil thickness. The CPT measurements provide a vertical profile of soil resistance for a soil column. We tested the accuracy of the CPT measurement and found that the CPT (i) shows the largest change in resistance when entering weathered bedrock and (ii) can be stopped very sharply only in the presence of a boulder, in which case the resistance is so strong that the measurement was deemed suspicious and repeated nearby. Because the CPT may not clearly identify the potential presence of moraine deposits, we also visually inspected the soil and saprolite materials extracted by the auger. From the auger, the transition zone from soil to the saprolite or bedrock is based on the material size and color of retrieved samples (Fig. S4). When the auger reaches the bedrock shale, it cannot penetrate easily. We believe our measurement is relatively accurate and efficient, which provides a consistent assessment of soil thickness over space in comparison to other existing methods. Figure 1b–e 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. S5). The CPT and auger data are mostly in agreement. 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.

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. Next, 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.

## 4.1 Model parameterization of curvature with different smoothing techniques of DEM

The topographic curvature is the key variable for estimating the soil thickness for the empirical approach. However, curvature is an inherently resolution-dependent topographic feature that is derived from a DEM. A 0.5 m DEM can provide “noises” for the results of curvature. The goal here is to determine the optimal DEM resolution for curvature to match with the sampling data, and the smoothing methods provided here are only for the calculation of curvature. To investigate the optimal resolution of the DEM for calculating curvature, we tested three approaches as explained in Sect. 2.1.3: smoothing the DEM over space, polynomial fitting of the DEM, and smoothing the DEM over time. For smoothing over space, the elevations of original 0.5 m DEM were replaced by utilizing the mean elevation of the surrounding adjacent cells at the range of 1.5, 2.5, … , and 13.5 m. For the polynomial fitting, the diameter is also chosen over the same range and intervals as the way of smoothing the DEM. For smoothing over time, the evolution period of the topography is taken as 0, 0.25, 0.5, … , 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 parameter increasing linearly from the minimum value to the maximum (the same datasets as in the Morris OAT method). We chose the set of parameters which gives the smallest RMSD.

Different smoothing methods and study sites correspond to their own optimal DEM resolution for curvature (Fig. 2). For example, in this study the north-facing hillslope shows the best fitting with 4.5 m DEM smoothing over space, but the south-facing 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 3 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 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.

## 4.2 Global sensitivity analysis

We apply the Morris OAT method to investigate the global sensitivity of the
seven parameters (Table 1) in the hybrid model. For each parameter, we
calculate the absolute of the mean elementary effect, $\left|\mathit{\mu}\right|$, in that the higher number represents higher importance, and the standard
deviation of the elementary effect, *σ*, which represent the
nonlinearity effect or interactions with other parameters (Fig. 3). Each dot
represents an evaluation of one parameter at one sampling site. In general,
the parameters in the mass conservation model have higher $\left|\mathit{\mu}\right|$
values, meaning that they have a more significant impact on soil thickness
than the parameters in the empirical model. The diffusion coefficient,
*K*_{d}, is the most important factor (high $\left|\mathit{\mu}\right|$ value and
high nonlinearity (*σ*)) and thus should be carefully calibrated. It represents soil diffusive-like processes such as soil creeping and biogenic activities. The normalized soil depth (*h*_{o}) is also has a higher $\left|\mathit{\mu}\right|$ value, which suggests that it is a very important factor but more linear than *K*_{d} due to its relatively small *σ*. These factors imply that the diffusive process is the most important transport mechanism for hillslope soil erosion rather than the soil erosion from overland flow on the surface of a soil layer (Dietrich et al., 1995; Nicótina et al., 2011; Roering et al., 1999, 2001). They also imply that the normalized soil depth is the most important parameter for estimating the soil production rate at
the bottom of the soil layer. The two parameters from the empirical method, *a* and $\stackrel{\mathrm{\u203e}}{h}$, are used for soil depositional areas. The sensitivity of *a* and $\stackrel{\mathrm{\u203e}}{h}$ are nearly linear (since *σ* 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 in *σ*.

## 4.3 Hybrid data–model soil thickness predictions and their comparison to measurements

At erosional sites, the erosion from the land surface can be balanced out by the soil formation from the bottom; therefore, the soil thickness may reach a steady-state condition. By coupling soil thickness with landscape evolution, we found that the soil thickness reaches a dynamic steady state after approximately 25 kyr at this study site (Fig. S6), 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 have already reached a steady-state condition since the last glacial legacy. Here, we only focus on the steady-state condition at erosional sites because they are where we apply the mass conservation equation. For depositional sites, the soil gradually accumulates on the land surface; meanwhile, the soil weathers slowly at the bottom. Therefore, the soil thickness is supposed to continuously increase (Dietrich et al., 1995). Due to the complexity of soil depositional environments, such as expansion or compression of soils, we consider using an empirical relationship as appropriate for the soil thickness at depositional areas.

A spatial map of soil thickness that compares modeling predictions 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, soil thickness of the south-facing side ranges from about 0.2 to 1 m, with a mean value of 0.55 m, while on the north-facing side, the soil thickness ranges from 0.15 to 1.5 m, with a mean value of 0.67 m. We use the sampling data from both the auger and the CPT to calibrate the seven parameters (Table 1) for the south-facing and north-facing hillslopes separately. One explanation for the mismatch between modeling and field observations (Fig. 4b) could be that parts of the hillslopes are on terraces. These areas may have fluvial deposits on ancient floodplains before they were turned into terraces (Yan et al., 2018).

We use the sampling data from auger and the CPT to calibrate the model
parameters (Table 1) for the south-facing and north-facing hillslopes
separately. The calibration is performed using a grid search approach where
the model loops through the entire parameter space with six evenly
distributed values for each parameter. The range of each parameter is based
on the literature and shown in Table 1. The overland flow coefficient
*K*_{s} is not included in this process because it is not sensitive to soil thickness (Fig. 3). The ranges of the remaining six parameters are based on the literature and shown in Table 1. We created an ensemble of the soil thickness based on 6^{6}(=46646) combinations of parameter sets. We then calculate the root-minimum-square error (RMSE) between the simulated and measured soil thickness across the site. Each set of parameters corresponds to an RMSE, and the distribution of the RMSE as a function of the parameters is shown in Fig. S7. The global minimum can be defined as the set of parameters that provide the lowest RMSE. In this specific study site, we find a unique set of parameters that provide the global minimum with the grid search study. Moreover, when the threshold of RMSE decreases, the number of samples with RMSE below this threshold decrease and converges to the global minimum. This result indicates the absence of other potential strong minima, and we can therefore identify the set of the parameters that provides the global minimum. In addition, we have applied a Bayesian approach to estimate the posterior distribution of parameters, as well as the maximum a posterior estimates (see the Supplement for more information). We find that the maximum a posteriori estimates (Fig. S8) also correspond to the global minimum of the RMSE. Moreover, among the posterior distribution of each parameter, *P*_{o} and *E*_{thre} are closer to uniform in the south-facing hillslope than the north-facing hillslope. The reason for this may be that the north-facing hillslope has more sampling points, which provides a better estimation than the south-facing hillslope.

The difference in soil thickness between the two hillslopes is evident. This is controlled by insolation due to the topographic aspect. The air temperatures and potential evapotranspiration rates produce significantly different microclimates that determine the 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 the thicker saprolite layers (Fig. 4b), which results in less resistance during excavations in the field sampling. This implies that 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, a higher concentration of fine particles, and more biomass, which leads to a positive feedback (Pelletier et al., 2018; Roering et al., 2010).

## 4.4 Predictive value of landscape features for estimating soil thickness

We evaluate the correlation between soil thickness from sampling and other landscape features obtained from remote sensing data (Fig. 5). Among about 18 topographic and land cover 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 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 curvature have the highest correlation (−0.87), which implies that the local relief of ridges and valleys is 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 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 that (2) the hybrid model is applicable to other sites given similar soil types and topography, without 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 advantages of this hybrid model are that there are only seven parameters to calibrate and that fewer sampling points are required. The extension to other watersheds is also easier with the hybrid approach. This hybrid method also provides higher accuracy than Patton's method in this study site, particularly at very thin or thicker soil layers (Fig. S9). 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, as discussed below.

## 4.5 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 maps 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 in the coupled hydro–bio–geomorphological processes. In contrast, on the north-facing hillslope where solar radiation is limited, the energy becomes a limiting factor of critical zone 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 a higher surface soil erosion rate. Moreover, on the south-facing hillslope, the lower interception of rainfall due to thinner vegetation cover and more sand material from the soil texture (Fig. S10) also increases sediment mobilization by raindrop impacts.

The contribution of overland flow-driven soil transport (erosion and
deposition) to the soil thickness formation is minor compared to the
diffusion-driven soil transport. The water depth reaches steady-state after
about 6 d with a constant rainfall intensity (i.e., 363 mm/yr) in the
study area (Fig. S11a). We apply this spatial map of the water depth to
drive the soil transport by overland flow (Eq. 3). The soil transport rate
from overland flow mostly happens in the water pathways, which have no
ponding after snow melting and storm events. The ratio of the soil transport
rate from overland flow to the total soil transport rate is mostly less than
1 % in the two hillslopes (Fig. S11b). This minor impact from overland
flow also explains why the parameter *K*_{d} is not sensitive in this
hybrid model.

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. One
should note that the soil production rate here is different from bedrock
weathering rate, which is controlled by water table dynamics
(Wan et al., 2019). The actual soil production rate
is controlled by the soil thickness and potential production rate, *P*_{o}
(Eq. 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 *P*_{o}. *P*_{o} on 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
(Tran, et al., 2019). 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).

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 *P*_{o}. However, the
thicker soil thickness offsets the impact of *P*_{o} in that thicker soil
thickness results in slower actual soil production rate. This is a possible
reason why the north-facing hillslope has slower actual soil production rate
even though it has a higher *P*_{o}.

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 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 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 the lidar DEM over time has a higher efficiency than smoothing it over space to obtain the optimal topographic curvature values, which provides the least error between the modeling results and sampling soil thickness.

This hybrid model is a flexible, generally applicable approach to predicting soil thicknesses. The hybrid model, with only 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 south-facing hillslope as a result of the different insolation at different aspects. The model analysis suggests that the south-facing 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 dense vegetation cover, less solar radiation, and wetter surface soil material as fundamentally controlled by aspect.

The limitation of the hybrid modeling approach developed here is that it would fail in alluvial depositional sites (i.e., floodplains) 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 mechanisms for soil thickness prediction in hillslopes. Although the example applications in this paper are at two hillslopes, this hybrid model framework should have little limitation when analyzing soil-mantled mountainous hillslopes after calibration with sampling dataset.

The dataset and code used in the hybrid model are openly available at http://doi.org/10.5281/zenodo.4445383 (Yan, 2021).

Remote sensing products can be accessed at https://doi.org/10.15485/1618131 (Brodrick et al., 2020) and https://doi.org/10.15485/1617203 (Goulden et al., 2020).

The supplement related to this article is available online at: https://doi.org/10.5194/esurf-9-1347-2021-supplement.

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 the 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 the soil thickness evolution study, and revised the manuscript. SSH supervised this study and revised the manuscript. NF provided the remote sensing data and revised the manuscript.

The authors declare that they have no conflict of interest.

Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

We thank Stijn Wielandt for his help during the field campaigns with soil sample collection. We thank Ilhan Ozgen for helping with the overland flow simulation.

The Watershed Function Scientific Focus Area is funded by the U.S. Department of Energy, Office of Science, Office of Biological and Environmental Research under award no. DE-AC02- 05CH11231.

This paper was edited by Cornelia Rumpel and reviewed by Jon Pelletier and Nicholas Patton.

Andrews, D. J. and Bucknam, R. C.: Fitting degradation of shoreline scarps by a nonlinear diffusion model, J. Geophys. Res.-Sol. Ea., 92, 12857–12867, https://doi.org/10.1029/JB092iB12p12857, 1987.

Breiman, L.: Random Forests, Mach. Learn., 45, 5–32, https://doi.org/10.1023/A:1010933404324, 2001.

Brodrick, P., Goulden, T., and Chadwick, K. D.: Custom NEON AOP reflectance mosaics and maps of shade masks, canopy water content, Watershed Function SFA [data set], https://doi.org/10.15485/1618131, 2020.

Brugger, K. A.: Climate in the Southern sawatch range and Elk Mountains, Colorado, U.S.A., during the last glacial maximum: Inferences using a simple degree-day model, Arctic, Antarct. Alp. Res., 42, 164–178, https://doi.org/10.1657/1938-4246-42.2.164, 2010.

Campolongo, F., Cariboni, J., and Saltelli, A.: An effective screening design for sensitivity analysis of large models, Environ. Model. Softw., 22, 1509–1518, https://doi.org/10.1016/j.envsoft.2006.10.004, 2007.

Carvalhais, N., Forkel, M., Khomik, M., Bellarby, J., Jung, M., Migliavacca, M., Mu, M., Saatchi, S., Santoro, M., Thurner, M., Weber, U., Ahrens, B., Beer, C., Cescatti, A., Randerson, J. T., and Reichstein, M.: Global covariation of carbon turnover times with climate in terrestrial ecosystems, Nature, 514, 213–217, https://doi.org/10.1038/nature13731, 2014.

Catani, F., Segoni, S., and Falorni, G.: An empirical geomorphology-based approach to the spatial prediction of soil thickness at catchment scale, Water Resour. Res., 46, W05508, https://doi.org/10.1029/2008WR007450, 2010.

Chadwick, K. D., Brodrick, P. G., Grant, K., Goulden, T., Henderson, A., Falco, N., Wainwright, H., Williams, K. H., Bill, M., Breckheimer, I., Brodie, E. L., Steltzer, H., Williams, C. F. R., Blonder, B., Chen, J., Dafflon, B., Damerow, J., Hancher, M., Khurram, A., Lamb, J., Lawrence, C. R., McCormick, M., Musinsky, J., Pierce, S., Polussa, A., Hastings Porro, M., Scott, A., Singh, H. W., Sorensen, P. O., Varadharajan, C., Whitney, B., and Maher, K.: Integrating airborne remote sensing and field campaigns for ecology and Earth system science, Methods Ecol. Evol., 11, 1492–1508, https://doi.org/10.1111/2041-210X.13463, 2020.

Dietrich, W. E., Reiss, R., Hsu, M., and Montgomery, D. R.: A process-based model for colluvial soil depth and shallow landsliding using digital elevation data, Hydrol. Process., 9, 383–400, 1995.

Efron, B.: The Jackknife, the Bootstrap and Other Resampling Plans, edited by CBMS-NSF regional conference series in applied mathematics 1982, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1982.

Fan, Y., Clark, M., Lawrence, D. M., Swenson, S., Band, L. E., Brantley, S. L., Brooks, P. D., Dietrich, W. E., Flores, A., Grant, G., Kirchner, J. W., Mackay, D. S., McDonnell, J. J., Milly, P. C. D., Sullivan, P. L., Tague, C., Ajami, H., Chaney, N., Hartmann, A., Hazenberg, P., McNamara, J., Pelletier, J., Perket, J., Rouholahnejad-Freund, E., Wagener, T., Zeng, X., Beighley, E., Buzan, J., Huang, M., Livneh, B., Mohanty, B. P., Nijssen, B., Safeeq, M., Shen, C., van Verseveld, W., Volk, J., and Yamazaki, D.: Hillslope Hydrology in Global Change Research and Earth System Modeling, Water Resour. Res., 55, 1737–1772, https://doi.org/10.1029/2018WR023903, 2019.

Fernandes, N. F. and Dietrich, W. E.: Hillslope evolution by diffusive processes: The timescale for equilibrium adjustments, Water Resour. Res., 33, 1307–1318, https://doi.org/10.1029/97WR00534, 1997.

Gaskill, D. L.: Geologic map of the Gothic quadrangle, Gunnison County, Colorado, https://doi.org/10.3133/gq1689, 1991.

Goulden, T., Hass, B., Brodie, E., Chadwick, K. D., Falco, N., Maher, K., Wainwright, H., and Williams, K.: NEON AOP Survey of Upper East River CO Watersheds: LAZ Files, LiDAR Surface Elevation, Terrain Elevation, and Canopy Height Rasters, Watershed Function SFA [data set], https://doi.org/10.15485/1617203, 2020.

Grant, G. E. and Dietrich, W. E.: The crontier beneath our feet, J. Chem. Inf. Model., 53, 1689–1699, https://doi.org/10.1017/CBO9781107415324.004, 2017.

Hastie, R.: Problems for Judgment and Decision Making, Annu. Rev. Psychol., 52, 653–683, https://doi.org/10.1146/annurev.psych.52.1.653, 2001.

Heimsath, A. M., Dietrich, W. E., Nishiizumi, K., Finkel, R. C., Mass, A., and National, L. L.: The soil production function and landscape equilibrium, Nature, 388, 358–361, 1997.

Heimsath, A. M., Chappell, J., Dietrich, W. E., Nishiizumi, K., and Finkel, R. C.: Soil production on a retreating escarpment in southeastern Australia, Geology, 28, 787–790, https://doi.org/10.1130/0091-7613(2000)28<787:SPOARE>2.0.CO;2, 2000.

Heimsath, A. M., Dietrich, W. E., Nishiizumi, K., and Finkel, R. C.: Stochastic processes of soil production and transport: erosion rates, topographic variation and cosmogenic nuclides in the Oregon Coast Range, Earth Surf. Proc. Land., 26, 531–552, https://doi.org/10.1002/esp.209, 2001.

Heimsath, A. M., Furbish, D. J., and Dietrich, W. E.: The illusion of diffusion: Field evidence for depth-dependent sediment transport, Geology, 33, 949–952, https://doi.org/10.1130/G21868.1, 2005.

Heimsath, A. M., Fink, D., and Hancock, G. R.: The `humped' soil production function: eroding Arnhem Land, Australia, Earth Surf. Proc. Land., 34, 1674–1684, https://doi.org/10.1002/esp.1859, 2009.

Hengl, T., Heuvelink, G. B. M., and Stein, A.: A generic framework for spatial prediction of soil variables based on regression-kriging, Geoderma, 120, 75–93, https://doi.org/10.1016/j.geoderma.2003.08.018, 2004.

Hinckley, E. L. S., Ebel, B. A., Barnes, R. T., Anderson, R. S., Williams, M. W., and Anderson, S. P.: Aspect control of water movement on hillslopes near the rain-snow transition of the Colorado Front Range, Hydrol. Process., 28, 74–85, https://doi.org/10.1002/hyp.9549, 2014.

Hubbard, S. S., Williams, K. H., Agarwal, D., Banfield, J., Beller, H., Bouskill, N., Brodie, E., Carroll, R., Dafflon, B., Dwivedi, D., Falco, N., Faybishenko, B., Maxwell, R., Nico, P., Steefel, C., Steltzer, H., Tokunaga, T., Tran, P. A., Wainwright, H., and Varadharajan, C.: The East River, Colorado, Watershed: A Mountainous Community Testbed for Improving Predictive Understanding of Multiscale Hydrological-Biogeochemical Dynamics, Vadose Zone J., 17, 180061, https://doi.org/10.2136/vzj2018.03.0061, 2018.

Jackson, R. B., Lajtha, K., Crow, S. E., Hugelius, G., Kramer, M. G., and Piñeiro, G.: The Ecology of Soil Carbon: Pools, Vulnerabilities, and Biotic and Abiotic Controls, Annu. Rev. Ecol. Evol. Syst., 48, 419–445, https://doi.org/10.1146/annurev-ecolsys-112414-054234, 2017.

Joshua West, A.: Thickness of the chemical weathering zone and implications for erosional and climatic drivers of weathering and for carbon-cycle feedbacks, Geology, 40, 811–814, https://doi.org/10.1130/G33041.1, 2012.

Kuriakose, S. L., Devkota, S., Rossiter, D. G., and Jetten, V. G.: Prediction of soil depth using environmental variables in an anthropogenic landscape, a case study in the Western Ghats of Kerala, India, CATENA, 79, 27–38, https://doi.org/10.1016/j.catena.2009.05.005, 2009.

Kilinc, M. Y. and Richardson, E. V.: Mechanics of soil erosion from overland flow generated by simulated rainfall, Hydrol. Pap., 63, available at: http://hdl.handle.net/10217/61574 (last access: October 2021), 1973.

Lal, A. M. W.: Performance comparison of overland flow algorithms, J. Hydraul. Eng., 124, 342–349, https://doi.org/10.1061/(ASCE)0733-9429(1998)124:4(342), 1998.

Li, M., Foster, E. J., Le, P. V. V, Yan, Q., Stumpf, A., Hou, T., Papanicolaou, A. N., Wacha, K. M., Wilson, C. G., Wang, J., Kumar, P., and Filley, T.: A new dynamic wetness index (DWI) predicts soil moisture persistence and correlates with key indicators of surface soil geochemistry, Geoderma, 368, 114239, https://doi.org/10.1016/j.geoderma.2020.114239, 2020.

Morris, M. D.: Factorial sampling plans for preliminary computational experiments, Technometrics, 33, 161–174, https://doi.org/10.1080/00401706.1991.10484804, 1991.

Nicótina, L., Tarboton, D. G., Tesfa, T. K., and Rinaldo, A.: Hydrologic controls on equilibrium soil depths, Water Resour. Res., 47, W04517, https://doi.org/10.1029/2010WR009538, 2011.

Patton, N. R., Lohse, K. A., Godsey, S. E., Crosby, B. T., and Seyfried, M. S.: Predicting soil thickness on soil mantled hillslopes, Nat. Commun., 9, 3329, https://doi.org/10.1038/s41467-018-05743-y, 2018.

Patton, N. R., Lohse, K. A., Seyfried, M. S., Godsey, S. E., and Parsons, S. B.: Topographic controls of soil organic carbon on soil-mantled landscapes, Sci. Rep., 9, 6390, https://doi.org/10.1038/s41598-019-42556-5, 2019.

Pelletier, J. D. and Rasmussen, C.: Geomorphically based predictive mapping of soil thickness in upland watersheds, Water Resour. Res., 45, W09417, https://doi.org/10.1029/2008WR007319, 2009.

Pelletier, J. D., Barron-Gafford, G. A., Breshears, D. D., Brooks, P. D., Chorover, J., Durcik, M., Harman, C. J., Huxman, T. E., Lohse, K. A., Lybrand, R., Meixner, T., McIntosh, J. C., Papuga, S. A., Rasmussen, C., Schaap, M., Swetnam, T. L., and Troch, P. A.: Coevolution of nonlinear trends in vegetation, soils, and topography with elevation and slope aspect: A case study in the sky islands of southern Arizona, J. Geophys. Res.-Earth, 118, 741–758, https://doi.org/10.1002/jgrf.20046, 2013.

Pelletier, J. D., Broxton, P. D., Hazenberg, P., Zeng, X., Troch, P. A., Niu, G.-Y., Williams, Z., Brunke, M. A., and Gochis, D.: A gridded global data set of soil, intact regolith, and sedimentary deposit thicknesses for regional and global land surface modeling, J. Adv. Model. Earth Syst., 8, 41–65, https://doi.org/10.1002/2015MS000526, 2016.

Pelletier, J. D., Barron-Gafford, G. A., Gutiérrez-Jurado, H., Hinckley, E. L. S., Istanbulluoglu, E., McGuire, L. A., Niu, G. Y., Poulos, M. J., Rasmussen, C., Richardson, P., Swetnam, T. L., and Tucker, G. E.: Which way do you lean? Using slope aspect variations to understand Critical Zone processes and feedbacks, Earth Surf. Proc. Land., 43, 1133–1154, https://doi.org/10.1002/esp.4306, 2018.

Perron, J. T.: Numerical methods for nonlinear hillslope transport laws, J. Geophys. Res.-Earth, 116, F02021, https://doi.org/10.1029/2010JF001801, 2011.

Roering, J. J., Kirchner, J. W., and Dietrich, W. E.: Evidence for nonlinear, diffusive sediment transport on hillslopes and implications for landscape morphology, Water Resour. Res., 35, 853–870, https://doi.org/10.1029/1998WR900090, 1999.

Roering, J. J., Kirchner, J. W., and Dietrich, W. E.: Hillslope evolution by nonlinear, slope-dependent transport: Steady state morphology and equilibrium adjustment timescales, J. Geophys. Res.-Sol. Ea., 106, 16499–16513, https://doi.org/10.1029/2001JB000323, 2001.

Roering, J. J., Marshall, J., Booth, A. M., Mort, M., and Jin, Q.: Evidence for biotic controls on topography and soil production, Earth Planet. Sci. Lett., 298, 183–190, https://doi.org/10.1016/j.epsl.2010.07.040, 2010.

Shangguan, W., Hengl, T., Mendes de Jesus, J., Yuan, H., and Dai, Y.: Mapping the global depth to bedrock for land surface modeling, J. Adv. Model. Earth Syst., 9, 65–88, https://doi.org/10.1002/2016MS000686, 2017.

Taylor, J. A., Jacob, F., Galleguillos, M., Prévot, L., Guix, N., and Lagacherie, P.: The utility of remotely-sensed vegetative and terrain covariates at different spatial resolutions in modelling soil and watertable depth (for digital soil mapping), Geoderma, 193–194, 83–93, https://doi.org/10.1016/j.geoderma.2012.09.009, 2013.

Temme, A. J. A. M. and Vanwalleghem, T.: LORICA – A new model for linking landscape and soil profile evolution: Development and sensitivity analysis, Comput. Geosci., 90, 131–143, https://doi.org/10.1016/j.cageo.2015.08.004, 2016.

Tesfa, T. K., Tarboton, D. G., Chandler, D. G., and McNamara, J. P.: Modeling soil depth from topographic and land cover attributes, Water Resour. Res., 45, W10438, https://doi.org/10.1029/2008WR007474, 2009.

Tokunaga, T. K., Wan, J., Williams, K. H., Brown, W., Henderson, A., Kim, Y., Tran, A. P., Conrad, M. E., Bill, M., Carroll, R. W. H., Dong, W., Xu, Z., Lavy, A., Gilbert, B., Romero, S., Christensen, J. N., Faybishenko, B., Arora, B., Siirila-Woodburn, E. R., Versteeg, R., Raberg, J. H., Peterson, J. E., and Hubbard, S. S.: Depth- and Time-Resolved Distributions of Snowmelt-Driven Hillslope Subsurface Flow and Transport and Their Contributions to Surface Waters, Water Resour. Res., 55, 9474–9499, https://doi.org/10.1029/2019WR025093, 2019.

Tran, A. P., Rungee, J., Faybishenko, B., Dafflon, B., and Hubbard, S. S.: Assessment of spatiotemporal variability of evapotranspiration and its governing factors in a mountainous watershed, Water, 11, 243, https://doi.org/10.3390/w11020243, 2019.

Vanags, C., Minasny, B., and Mcbratney, A. B.: The dynamic penetrometer for assessment of soil mechanical resistance, available at: https://www.regional.org.au/au/asssi/ (last access: February 2020), 1–8, 2004.

Vanwalleghem, T., Stockmann, U., Minasny, B., and Mcbratney, A. B.: A quantitative model for integrating landscape evolution and soil formation, J. Geophys. Res.-Earth, 118, 331–347, https://doi.org/10.1029/2011JF002296, 2013.

Wainwright, H. M., Finsterle, S., Jung, Y., Zhou, Q., and Birkholzer, J. T.: Making sense of global sensitivity analyses, Comput. Geosci., 65, 84–94, https://doi.org/10.1016/j.cageo.2013.06.006, 2014.

Wan, J., Tokunaga, T. K., Williams, K. H., Dong, W., Brown, W., Henderson, A. N., Newman, A. W., and Hubbard, S. S.: Predicting sedimentary bedrock subsurface weathering fronts and weathering rates, Sci. Rep., 9, 17198, https://doi.org/10.1038/s41598-019-53205-2, 2019.

Yan, Q.: Qinayan/Soil-thickness: Soil thickness estimation (v1.0.0), Zenodo [data set], https://doi.org/10.5281/zenodo.4445383, 2021.

Yan, Q., Iwasaki, T., Stumpf, A., Belmont, P., Parker, G., and Kumar, P.: Hydrogeomorphological differentiation between floodplains and terraces, Earth Surf. Proc. Land., 43, 218–228, https://doi.org/10.1002/esp.4234, 2018.

Yan, Q., Le, P. V. V, Woo, D. K., Hou, T., Filley, T., and Kumar, P.: Three-Dimensional Modeling of the Coevolution of Landscape and Soil Organic Carbon, Water Resour. Res., 55, 1218–1241, https://doi.org/10.1029/2018WR023634, 2019.

Yan, Q., Kumar, P., Wang, Y., Zhao, Y., Lin, H., and Ran, Q.: Sustainability of soil organic carbon in consolidated gully land in China's Loess Plateau, Sci. Rep., 10, 16927, https://doi.org/10.1038/s41598-020-73910-7, 2020.