Prediction of soil formation as a function of age using the percolation theory approach

: Recent modeling and comparison with field results showed that soil formation by chemical weathering, either from bedrock or unconsolidated material, is limited largely by solute transport. Chemical weathering rates are proportional to solute velocities. Nonreactive solute transport described by non-Gaussian transport theory appears compatible with soil formation rates. This change in understanding opens new possibilities for predicting soil production and depth across orders of magnitude of time scales. Percolation theory for modeling the evolution of soil depth and production was applied to new and published data for alpine and Mediterranean soils. The first goal was to check whether the empirical data conform to the theory. Secondly we analyzed discrepancies between theory and observation to find out if the theory is incomplete, if modifications of existing experimental procedures are needed and what parameters might be estimated improperly. Not all input parameters required for current theoretical formulations (particle size, erosion, and infiltration rates) are collected routinely in the field; thus, theory must address how to find these quantities from existing climate and soil data repositories, which implicitly introduces some uncertainties. Existing results for soil texture, typically reported at relevant field sites, had to be transformed to results for a median particle size, d50, a specific theoretical input parameter. The modeling tracked reasonably well the evolution of the alpine and Mediterranean soils. For the Alpine sites we found, however, that we consistently overestimated soil depths by 45%. Particularly during early soil formation, chemical weathering is more severely limited by reaction kinetics than by solute transport. The kinetic limitation of mineral weathering can affect the system until 1 kyr to a maximum of 10 kyr of soil evolution. Thereafter, solute transport seems dominant. The trend and scatter of soil depth evolution is well captured, particularly for Mediterranean soils. We assume that some neglected processes, such as Abstract 17 Recent modeling and comparison with field results showed that soil formation by chemical weathering, 18 either from bedrock or unconsolidated material, is limited largely by solute transport. Chemical 19 weathering rates are proportional to solute velocities. Nonreactive solute transport described by non- 20 Gaussian transport theory appears compatible with soil formation rates. This change in understanding 21 opens new possibilities for predicting soil production and depth across orders of magnitude of time scales. 22 Percolation theory for modeling the evolution of soil depth and production was applied to new and 23 published data for alpine and Mediterranean soils. The first goal was to check whether the empirical data 24 conform to the theory. Secondly we analyzed discrepancies between theory and observation to find out if 25 the theory is incomplete, if modifications of existing experimental procedures are needed and what 26 parameters might be estimated improperly. Not all input parameters required for current theoretical 27 formulations (particle size, erosion and infiltration rates) are collected routinely in the field; thus, theory 28 must address how to find these quantities from existing climate and soil data repositories, which 29 implicitly introduces some uncertainties. Existing results for soil texture, typically reported at relevant 30 field sites, had to be transformed to results for a median particle size, d 50 , a specific theoretical input generate soil depth and confirms decreasing production rates with age. A steady state for soils is not 40 reached before about 100 kyr.


Setting 47
The importance of quantifying chemical weathering and soil formation rates has as its basis their 48 relevance across a wide range of fields of study, from agricultural engineering (Montgomery, 2007ab;49 Lal, 2010) through climate change (Berner, 1992;Raymo, 1994 presence of humans at this boundary, and the many uses we have for, as well as burdens we place on, soil 56 (Hillel, 2005). 57 Concepts in soil formation have been developed since the 19 th Century (Darwin, 1881;Dokuchaev, 1948). 58 Soil is formed from bedrock by a complex interaction of biochemical and physical processes (Jenny,59 1941; Huggett, 1998). The role of physical fractures in soil formation, by allowing water penetration of 60 bedrock, was recognized as long ago as Gilbert (1877). Soil also forms in unconsolidated sediments 61 developed through processes such as alluvial deposition (White et al., 1996) tree-throw (Bormann et al., 62 1995), landslides (Smale et al., 1997;Trustrum and deRose, 1998), mining (Frouz et al., 2005), and 63 exposure of glacial sediments by glacial retreat (Mavris et al., 2010;Egli et al., 2014). In the latter, the 64 role of any gradual physical weathering processes is correspondingly reduced, since water penetration is 65 guaranteed. 66 The perception of the importance of chemical weathering to soil formation and geology in general has 67 been steadily climbing for decades (Anderson and Anderson, 2010), partly on the basis of the 68 experimentally determined proportionality between soil production rates, erosion and chemical partly through the increased attention to its role in past climate change (Berner, 1992;Raymo, 1994, 71 Algeo and Scheckler, 1998). In the latter, it is only the weathering of silicates that is relevant, as, e.g., 72 weathering of carbonates has no effect on the atmospheric carbon content (e.g., Berner, 1992). Here the 73 simplification offered by the representation of silicate weathering in the Urey reaction (Urey, 1952) 74 CaSiO 3 + CO 2 ↔ CaCO 3 + SiO 2 demonstrates the universality of its impact on climate; each mole of 75 silicate mineral weathered draws down one mole of carbon dioxide from the atmosphere, with only the 76 CO 2 present as a gas. The often-cited Urey reaction describes the weathering of wollastonite. A general 77 reaction for the more common plagioclase is (White at al., 1999): 78 Ca x Na (1-x) Al (1+x) Si (3-x) O 8 + (1+x)CO 2 + (2 + 2x)H 2 O ↔ xCa 2+ + (1-x)Na + + (3-x)SiO 2 + (1+x)HCO 3 -+ 79 (1+x)Al(OH) 3 . 80 The Urey representation by itself or the dissolution of plagioclase do not explicitly give form to the 81 known role of water fluxes in the geochemistry, which is typically, according to both experiment (Maher,82 2010) and theory (Hunt et al., 2015), the limiting factor. Identification of this particular limitation then 83 places the science of hydrology right at the center of study. 84 While one can focus on any one of a number of the individual components of the entire process of 85 weathering, such as: 1) the necessity for vegetation to deliver CO 2 fluxes above equilibrium values as 86 derived from atmospheric CO 2 levels, 2) the importance of fractures in facilitating water penetration into 87 bedrock, 3) the interactions of the various components of biota in the porous medium that contribute the 88 respiration necessary to produce the CO 2 , other research is currently providing grounds for focusing one's 89 attention on chemical weathering processes. For example, it has been shown (Eppes and Keanini (2007) 90 that crack formation in bedrock under subcritical conditions is chiefly a result of chemical weathering, 91 which changes mineral volumes, thus producing the stress that continues the fracturing. The particular 92 phenomenon of fracture formation, traditionally viewed as predominantly physical, consequently turns 93 out to be often, at its core, a chemical process. At the same time, however, theoretical and experimental 94 studies of chemical weathering indicate that the rate at which this process occurs is very often limited by 95 solute transport (Maher, 2010;Hunt and Ghanbarian, 2016;Braun, 2016;Yu and Hunt, 2017c), that is to 96 say, the rate at which reaction products may be removed from the weathering front. A quasi-universal 97 description of soil formation and production could answer a question first posed by Dokuchaev at the turn 98 of the 20 th Century (the translation into English, however, did not arrive until 1948; Dokuchaev, 1948), 99 regarding how to frame soil formation in terms of a set of consistent (soil forming) factors: Climate, 100 organisms, relief, parent material, and time (Yu et al., 2017). 101 In summary, while it is important to keep in mind that weathering phenomena may not be easily separable 102 into, e.g., chemical and physical processes, for interpretation it is still important to be able to identify the 103 rate-limiting input into the soil formation process. The evidence summarized briefly above places an 104 importance on the development of a model of soil formation, based on chemical weathering, that can 105 deliver specific predictions for soil depth in a general framework, with results that can be applied across 106 boundaries defining geological, climatological, and biogeochemical variability. 107 108

Relationship of Chemical Weathering to Solute Transport 109
In order to understand properly the limits on chemical weathering placed by solute transport, it turns out 110 that it is necessary to use a modern solute transport theory for heterogeneous porous media, which 111 generates non-Gaussian predictions (Hunt and Skinner, 2008;Hunt et al., 2011;Ghanbarian et al., 2012). 112 Here, non-Gaussian means several things: 1) the solute velocity (or flux) is only proportional to the water 113 flow rate, not equal to it, 2) the solute velocity decays according to a power of the time, 3) the dispersivity 114 is linear in the time (Scher, et al., 1991). In the present case, condition 2) is most important, since the 115 predicted power-law decay of the solute velocity in time (Hunt and Ghanbarian, 2016) matches the 116 experimentally observed power in both soil formation (Egli et al., 2014;Friend, 1992) and for chemical 117 weathering (White and Brantley, 2003). Condition 1) is also relevant, since soil formation and chemical 118 weathering rates are also typically proportional to water fluxes, and thus to precipitation. This simple 119 statement is in accord with a great deal of data on soil production and/or chemical weathering (Maher,  solute transport in heterogeneous porous media is so often non-Gaussian, rather than Gaussian, (Cushman 122 and O'Malley, 2015) that the non-Gaussian condition is becoming recognized as the norm, rather than the 123 anomaly. It is also important that the non-Gaussian transport results used are derived from network 124 models (Hunt and Manzoni, 2016), which places a premium on identifying network properties not always 125 investigated in the field, such as a fundamental length scale, and which generate results fundamentally 126 different from continuum models of solute transport. 127 In treatments of solute transport limitations on chemical weathering by non-Gaussian solute transport, it 128 has been shown that the solute flux is an excellent proxy for reaction rates  as, aside 129 from a factor equal to molar density, the weathering rate (expressed in moles m -2 s -1 ) is equivalent to a 130 solute velocity. The more rapid the Urey reaction or the plagioclase dissolution reaction occurs, the more 131 likely it is that solute transport is the limiting factor over the entire period of soil formation. The 132 fundamental length scale brings in the most important single parameter from the parent material. The flow 133 rate brings in climate variables. 134 The theory of solute transport limited weathering ("percolation theory"; Hunt  formulas, as well as the numerical routine that have been derived in initial publications, with their closer 141 correspondence to chemistry, rather than geology, do not contain variability over time in such important 142 input parameters as flow rates, or particle size. It is known that in the field these quantities vary over time, 143 both systematically and randomly, and over all time scales, and some of these effects will ultimately need 144 to be considered. 145 This theoretical description has been thoroughly reviewed in some recent books Hunt 146 and Maonzoni, 2016), as well as in Reviews of Geophysics (Hunt and Sahimi, 2017). This theoretical 147 treatment so far does not incorporate directly changes in the medium associated with changing density, 148 changing particle size (network scale), or changing depth due either to deposition or erosion, but such 149 changes can, in principle, be introduced into the equation describing changes in soil depth as calculated 150 from the erosion and soil production rates, a topic deferred to later. 151 152

Relationship Between Chemical Weathering and Soil Formation 153
A range of studies over the past decade has shown empirically (Burke et al., 2006;Dixon et al., 2009;154 Egli et al., 2014), at least, that the processes of soil formation (and thus soil upbuilding and therefore 155 thickness) and chemical weathering are highly correlated. Since the "w" in Bw horizon indicates the 156 development of color or structure, or both (Soil Survey Division Staff, 2014), and thus weathering, using 157 the bottom of the Bw for calculating the soil depth would mean that such a correlation is partly a result of 158 definition. But some of these studies (Burke et al., 2006;Dixon et al., 2009) have also indicated that the 159 correlation between the two processes is rather loose. Sources of discrepancy can be that total soil 160 thickness may include an organic layer (O horizon), for example, which is deposited on top of the 161 weathered soil, or that aeolian dust deposition or alluvial deposits, that are disconnected from the parent 162 material, may be significant, neither of which is directly connected to soil production from weathering. At 163 many of the sites discussed in the present study the correlation has already been shown to be robust (Hunt 164 and Ghanbarian, 2016; Yu and Hunt, 2017bc). Furthermore, it has also been shown (Hunt and 165 Ghanbarian, 2016) that the time dependence of the soil formation rates of these particular studies 166 corresponds closely to the time dependence of chemical weathering from completely different sources 167 (White and Brantley, 2003), and over time scales from years to millions of years. Determination that 168 chemical weathering and soil formation rates are proportional across a range of field experiments lends 169 support to the assumption that such a proportionality holds more universally, but does not constitute a 170 proof that extension to all sites under all conditions will be valid. 171 172

Present Study 173
Our present study applies systematically this recent model of soil formation derived from the chemical 174 weathering depth as limited mainly by solute transport (percolation theory) to two suites of soils; i.e., 175 alpine and Mediterranean soils. What distinguishes this particular work is in its more systematic approach 176 to the analysis, an improvement which is appropriate for the greater richness of the data analyzed. Such 177 an approach also allows an increase in the corresponding richness of interpretation. improperly. The field data were drawn from alpine sites as well as from sites with Mediterranean climate. 184 Data collected include what is required to calculate (or estimate) such relevant parameters as the 185 infiltration rate, the erosion rate, and a characteristic particle size, d 50 . In order to incorporate any 186 temporal variability in length and time scales into predictions, it will be necessary to make some 187 straightforward extensions of the theoretical model and, possibly, to collect additional data as well. The Peclet number, P e , a ratio of diffusion to advection times (Saffman, 1959;Pfannkuch, 1963), was 199 developed to distinguish regimes where either advective (P e > 1) or diffusive (P e <1) solute transport 200 might dominate. If the ADE is valid, dominance of advection, P e > 1, implies a solute velocity equal to 201 the flow velocity; in the latter case, which reproduces a diminishing solute velocity with time, no 202 proportionality to the flow rate can be developed, since P e < 1 and the diffusion term dominates (Hunt and 203 Ewing, 2016). Although the ADE can be applied to predicting solute transport at the scale of a single pore 204 (e.g., Neuman and di Federico, 2003), scaling up the application of the ADE to larger length scales 205 produces serious problems. The only way the ADE can relate the observed reduction in reaction rates to a 206 diminishing solute transport capability in time (rather diffusion like) is to abandon the observed 207 proportionality to the flow rates. In other words, one cannot have it both ways with the continuum 208 approach. 209 The solute velocity is obtained from a known scaling relationship between transit time and system length 210 (Lee et al., 1999), plus the identification of the fundamental length and time scales (Hunt, 2017). The 211 Damköhler number, Da I , which is the ratio (Salehikhoo et al., 2013) of a solute advection time to a 212 reaction time under well-mixed conditions (Yu and Hunt, 2017a), measures the importance of solute 213 transport relative to reaction kinetics on chemical weathering rates. When Da I is larger than 1, the 214 transport limitation for the chemical weathering rate is expected to be valid; for smaller values of Da I , 215 reaction kinetics will dominate. With increasing time of soil evolution soil depth usually increases. 216 Because the Damköhler number is a ratio of transport time to reaction time and because the solute 217 transport time increases more rapidly than linearly in system length, Da I increases during soil evolution 218 from its initial value (Yu and Hunt, 2017a). At shorter time scales, when Da I is small and the kinetics of 219 the particular weathering reaction dominant at each specific site dominates, more structure is expected to 220 be visible in the time-dependence of the weathering rates. 221 It is necessary at the outset to be clear that the solute transport limitations discussed here arise from 222 advective solute transport, not from diffusion, as has been argued by Bandupadhyay et al., (2017). 223 Arguments based on the Peclet number as calculated from characteristic instantaneous flow velocities at 224 the scale of a single pore were used to justify this a priori (Hunt and Manzoni, 2016). However, this may 225 be criticized on the grounds that one should use a yearly average velocity, since diffusion may be relevant 226 throughout the year, even in the vadose zone. The structure of advective solute transport paths in 227 heterogeneous porous media is fractal, though not, in general, related to any structure of the medium, and 228 it is this fractality, which leads to the power-law decay in solute transport fluxes with time (Hunt and 229 Ewing, 2016). In applying percolation theory, it is assumed that the flow paths of least resistance can be 230 calculated using the critical percolation probability. As long as the optimal flow paths can be described 231 using such critical paths, their fractal dimensionality is as given in percolation theory, explaining why the 232 particular characteristics of the medium have reduced importance to flow path characteristics. 233 Percolation theory generates a suite of properties relevant to flow, diffusion, and dispersion, as well as to 234 structure, although it is sometimes necessary to choose which percolation results are appropriate, or 235 whether an alternate formulation, such as the effective medium approximation (EMA), may be more 236 suited to generating an accurate prediction in any specific case (Hunt and Sahimi, 2017). Nevertheless, 237 the general theory of percolation is best discussed in terms of its topological implications first. Also, 238 although its application need not be restricted to a regular grid, it is much easier to discuss under such 239 restrictions, so only regular grids, also known as lattices, are considered here in the theory review. 240 The best review for understanding the basics of percolation theory is by Stauffer and Aharony (1994), 241 from which much of the following discussion ultimately comes. Consider a square lattice, on which sites 242 may be occupied by either plastic or steel spheres of the same size. The choice, metal or plastic, at any 243 given site is randomly generated, though the conclusions do not change fundamentally if spatial 244 correlations are added. Nearest neighbor spheres touch at one point. If sufficiently many of the spheres 245 are metallic, a continuously connected infinitely long path through metallic spheres is produced. This 246 transition occurs, for any given medium, at a specific value of the fraction of spheres that is metallic, values are smaller than site percolation values on the same lattice. In order for a bond percolation problem 252 to be relevant to real systems, the (hydraulic or electric) conductance values connecting sites (i.e., pores) 253 must vary widely (Hunt and Sahimi, 2017). In a procedure called critical path analysis (CPA) the 254 percolation probability can then be used in a way that generates the connected path of lowest total 255 resistance (to flow), i.e., the preferential flow paths (Pollak, 1972;Hunt, 2001). This is particularly 256 relevant here, as Sahimi (1994) has emphasized how the topological structure of the critical paths will 257 dominate solute transport whenever they describe such flow concentration as are observed. In natural 3D 258 systems p c values tend to be around 15% (Scher and Zallen, 1970) or, as a consequence of correlations, 259 even less Garboczi et al., 1995). The small values of p c in bond percolation in natural media imply that 260 only a relatively small portion of the medium generates nearly all the fluid flow. It thus also means that 261 percolation theory may generate connected flow paths that look like what hydrogeologists refer to as 262 preferential flow. It is quite generally believed that such preferential flow paths are important to solute 263 transport (National Research Council, 1996). Here, and in previous work, we have made this assumption 264 regarding the relevance of such paths, as medium flow paths were assumed compatible with CPA. Then 265 the percolation descriptions of, e.g., the tortuosity or fractal dimensionality, of such paths should 266 correspond to observation (Hunt and Sahimi, 2017). 267 Once one has an infinitely large interconnected cluster of sites (or bonds), it is possible to define the 268 percolation backbone as that part of the infinite connected cluster which remains when all sites (bonds) 269 that can be connected to it through only one point have been removed (Stanley, 1977). All sites that 270 connect through two points are retained. Sites that connect only through one point are called dead-ends, 271 since they do not support flow. The remaining structure has mass fractal dimensionality D b , a quantity 272 which is relevant for solute transport (Lee et al., 1999). 273 Importantly, for a wide range of conditions, the value of D b remains the same, but it differs fundamentally 274 depending on the dimensionality of the system studied (Sheppard et al 1999). Flow in fractures as well as 275 in strongly layered (anisotropic) soils may be fundamentally 2D, but, more generally, 3D conditions are 276 observed. The value of D b does depend sensitively on the saturation condition, i.e., whether the medium is 277 fully saturated, or whether imbibition or drainage is occurring (Sheppard et al., 1999). The latter two 278 processes are known as invasion percolation (Wilkinson and Willemsen, 1983). Most other medium 279 variations, such as the particular lattice type chosen, or the particle size values, do not change the value of 280 D b . However, under certain long-range correlations in the positions of highly conducting bonds, such 281 paths can be straighter than in random percolation (Sahimi and Mukhopadhyay, 1996), for which values 282 of D b are closer to 1 (a value of D b = 1 corresponds rather closely to Gaussian transport). However, in the 283 present work, complications from such correlations are not addressed. Moreover, evidence for the 284 relevance of such correlations, though demonstrated for such frequently measured quantities as the 285 electrical conductivity (Hunt and Sahimi, 2017), has not yet been found for problems of chemical 286 weathering or soil formation. 287 288

Scaling, Chemical Weathering and Soil Production 289
In the following, the scaling relationships relating time, distance, and velocity of solute transport are 290 reproduced from known results for systems near the percolation threshold (Sheppard et al., 1999). What 291 can make their applicability universal, however, is the tendency for water flow in disordered media to 292 follow paths of least resistance, as defined in CPA. Then the solute transport is controlled by paths near 293 the percolation threshold, as quantified by using percolation theory. 294 In Figure 1, the concept of the percolation theory is schematically drawn and compared with a soil mass 295 balance. Soil depth depends on mass input and output. Soil production includes all mass and volume 296 changes due to the transformation of the parent material into soil (by chemical and physical weathering 297 processes, mineral transformation), the lowering of the bedrock/parent material -soil boundary 298 (Heimsath et al., 1997), but also atmospheric deposition and net organic matter input: 299 (1) 300 where TP Soil = the transformation of the parent material or rock into soil; according to Eq. (1), A = 301 atmospheric deposition and O = net organic matter input, G = organic matter decay (Zollinger et al.,302 2017)). Erosion (E) and chemical leaching (W) contribute to mass losses. Besides P Soil and the parameters 303 E and W, soil depth is strongly related to the hydrologic water balance, granulometry of the medium and 304 time (or velocity). 305 The specific role of percolation theory in describing solute transport is now discussed. When solute enters 306 a medium at a point (Lee et al., 1999), say a site on a grid that corresponds to a particle with a particular 307 mineralogy, the time, t, it takes for the solute to travel a distance x, is proportional to x Db . This 308 proportionality can be represented as an equation, if appropriate values of constants representing a 309 fundamental length scale, x 0 , and corresponding time scale, t 0 , can be identified, 310 The ratio of x 0 /t 0 ≡ v 0 , having units of velocity, is assumed to represent the pore-scale fluid flow rate, 312 which is the only externally imposed velocity in the physical system. We refer to v 0 as a pore-scale flow 313 rate. In the context of the following discussion, it will become clear that v 0 must be a value that 314 characterizes the mean of the local vertical flow rate. Determination of one additional parameter, either x 0 315 or t 0 , completes the parameter determination in Eq.
(2). 316 If the network representation is applied to a porous medium, essential suggestive choice for x 0 is a pore 317 separation (Hunt and Manzoni, 2016), which should be more or less equal to a particle diameter, because 318 this length scale defines the separations of the local connections between flow pathways. 319 In a highly disordered network, where particle sizes can vary widely, we have proposed (Yu and Hunt,320 2017bc) that the best choice for x 0 is d 50 , the median particle size. Although temporal mean values of v 0 321 are on the order of 1m yr -1 , the instantaneous flow rates are much higher, either during storms or seasonal 322 snowmelt. The larger instantaneous flow rates, which are limited primarily by the hydraulic conductivity, 323 can often be on the order of 0.1m day -1 , a value typical for saturated conditions (Blöschl and Sivapalan, 324 1996), and used to provide the rationale for neglecting diffusion. 325 While the instantaneous flow rates must be large enough to neglect diffusion, the time-averaged flow 326 rates must be small enough to limit reaction. Otherwise, reaction kinetics can be the limiting factor. In 327 order to diagnose the relative importance of these two factors, one can apply the Damköhler number, 328 defined as (Yu and Hunt, 2017a): 329 to flow rates over four orders of magnitude (Maher, 2010), as are soil production rates (Hunt, 2017), and 336 2) the known dominance of advective solute transport over diffusive transport (Hunt and Manzoni, 2016) 337 at typical instantaneous subsurface flow rates (Bloeschl and Sivapalan, 1996). The tendency for Da I to 338 increase with column length is evident in Eq.
(3). The pore volume is proportional to length, cancelling 339 the linear factor in L in the numerator, but the more rapid than linear increase in transport time with length 340 leaves a second factor in length to the 0.87 power, which is not cancelled in the denominator. As long as 341 Da I < 1, reaction rates are time-independent, if the only limiting factor considered explicitly is solute 342 transport. When Da I > 1, reaction rates decline over time. In the previously analyzed case of MgCO 3 343 dissolution (Yu and Hunt, 2017a), for all field conditions the value of Da I was never smaller than tens of 344 thousands. However, reaction rates of silicate minerals, even under well-mixed conditions, can be orders 345 of magnitude smaller (Stumm and Morgan, 1996) than for MgCO 3 , as treated by Salehikhoo et al. (2013), 346 and further analyzed by Yu and Hunt (2017a). 347 Inverting the relationship Eq. (2) for t(x) leads to 348 The time derivative of x(t) yields the solute velocity, argued above to be a proxy for a chemical 350 Eq. (5) can be rewritten in a form that depends only on the distance, x. Since the introduction of erosion 353 can, in principle, make it impossible to define a unique time for a given transport distance, it is more 354 useful to write Eq. (5) in a form that eliminates time from the equation: 355 Here, the soil production function, R, was equated with v, and v 0 was substituted for x 0 /t 0 as well. In the 357 absence of erosion, it is a reasonable hypothesis that the total solute transport distance is equal to the soil 358 depth. Then, the temporal derivative of the soil depth, v, is the soil production function (in units of depth 359 divided by time), as also given by Eq. (5), making the soil production function proportional to the 360 chemical weathering rate. The proportionality constant is equal to the ratio of the (bulk) density of the 361 chemically weathered material to its molecular mass. 362 Eq. (4) for the soil depth can hold only as long as erosion can be neglected, which is very rarely the case. 363 Since the soil production rate (Eq. 5) declines with age (or depth, Eq. 6), the period of time when erosion 364 can be neglected is always limited. 365 366

Erosion and Soil Depth Evolution 367
When erosion, E, must be considered, one can construct an equation for the soil depth based on the 368 concept of mass balance, 369 Here, R is given by Eq. 6 and E is, for constant erosion rates, a parameter. In general, however, E is a 371 function of time. The term E in fact is equaled to denudation that includes besides the output of solid 372 material also chemical leaching of silicate particles. Compared to erosion, this type of leaching and, 373 therefore, loss is in most cases fully subordinate (Dixon and von Blanckenburg, 2012). In three 374 dimensions with moisture conditions corresponding either to wetting or full saturation, D b has the value 375 1.87 in a wide range of conditions (Sheppard et al., 1999), and has, up to now, been treated as universal. 376 Note that 1/1.87 = 0.53, which is very close to the exponent for diffusion, which explains in a single 377 relationship both the proportionality of weathering rates to flow rates, and the resemblance of their time-378 dependence to diffusion. 379 Substituting in R from Eq. (6) yields, 380 with I(t)/φ as the net infiltration rate that varies over time and φ = pore volume. 382 Owing to the fractional exponent of x introduced by R, Eq. (7) does not have an analytical solution, but it 383 may be readily solved numerically. One simply solves Eq. (4) for some initial (sufficiently small) time 384 and then calculates R from Eq. (6). Using the calculated value of R and the field value for E as well as an 385 arbitrary time step one can calculate the change in soil depth and add it to the initial value. Then one 386 calculates R from Eq. (6) using the new soil depth. This procedure is then simply followed until the total 387 time elapsed is equal to the age of the soil, or until the depth no longer changes in time, at which point a 388 steady-state soil depth has been generated. As long as no significant changes in parameters occur, steady-389 state conditions will then prevail. In steady state, dx/dt = 0, and the soil depth x may be obtained by 390 solution of E = R as, 391 At the opposite (short time) end of the time spectrum, an important complication can arise when chemical 393 weathering is not solute transport-limited. In particular, there is the theoretical possibility that for some 394 ranges of experimental, or field conditions, Da I < 1. In this case, a constant rate of weathering would 395 ensue, reaction kinetics provide a limitation which is unchanging in time. Since the limits imposed by 396 kinetics are stronger in this case, at least at short time scales, than those due to solute transport, data  Table S1), and Mediterranean (94 sites; Table S2). Most of 415 these data have already been published, but some of the Wind River Range, Wyoming, USA, are new. %!'" Other data have been accessed from the literature. All soils have developed from unconsolidated material. %!(" The sites are distributed on five continents, as shown in Fig. 2. %!)" To solve Eq. (8) the parameters"!"#$-"%"#$&'!&'#&"./0"( #)* "need to be known. Some of these parameters are %!*" more or less easily accessible (e.g. soil depth) and others need to be estimated (e.g. E(t)). Accompanying %#+" data relevant to assess actual evapotranspiration, AET, and soil particle sizes, which are necessary to %#!" make concrete predictions, have also been collected for each site within these regions. All data, together %##" with their sources, are given in the supplementary Tables S1-S6. %#$" %#%"

%#&"
3.2.1 Soil Depth %#'" Soil depths are determined mostly through a process of excavation and measurement."We used datasets %#(" where information was available about soil horizons designation and thickness. In addition, data was %#)" collected (where available) about soil density, coarse fragment content and grain size distribution. where FE stock (calculated as kg m -2 ) = the fine earth stock summed over all soil horizons, FE i = the %$*" proportion or concentration of fine earth, !z i (m) = the thickness of layer i and ! (t m -3 ) = soil density, M %%+" = soil mass including coarse fragments (calculated as kg m -2 ). %%!" " %%#"

Dating and Erosion
Sites were selected where numerical indication about the surface and its soils was available (e.g. 14 C, %%&" 10 Be, etc.) or where soils could be correlated to terraces, moraines or geologic events that were dated.

%%'"
Erosion rates are sensitive to local vegetation, relief and slope angle, aspect, climate, and topographic %%(" curvature, as well as substrate. Not all of these influences can be quantified. In situ measurements of %%)" erosion or denudation were available only for a few sites (e.g. the sites 'Sila', several sites in the %%*" European Alps and some of the Wind River Range). Otherwise, present-day erosion rates had to be %&+" time that is related to major climate changes. The transition from the Pleistocene to the Holocene was 459 accompanied by distinctly higher rates. Raab et al. (2018) showed that the erosion rates can be up to 10 460 times higher when climate very distinctly changes (transition Pleistocene to Holocene). These 461 fluctuations, which occurred during the evolution of many of the soils considered, are shown in Fig. 3 and  462 had to be considered for sites where erosion values were derived from the previously mentioned maps 463 (mostly alpine sites). We assumed that the maps provide an average and relatively reliable erosion rate 464 value for undisturbed sites for the entire Holocene. The erosion rate of soils having an older age had to be 465 corrected using the average trend given in Fig. 3. 466 The experimental procedures to determine the particle size fractions were mostly determined using the 469 techniques of wet sieving, for particles greater in diameter than 32µm, and at smaller sizes, by 470 hydrometry or by using X-ray techniques. Using these data it was possible to determine directly d 50 . Since 471 the fraction of coarse sand is less relevant for soil formation when compared to smaller fractions d 50 was 472 calculated for the fractions < 500µm. A mechanical disintegration of the rock material into small units Previously published data, for example, often give only the percentages of the three fundamental size 482 classes: sand, silt, and clay. Since it is necessary to know d 50 in order to make a concrete prediction of soil 483 depths, we developed a regression routine over well-characterized soils for this purpose, using the 484 percentages of sand, silt, and clay to calculate a mean diameter for the input, and the observed d 50 as an 485 output. Since we expected, at the very least, key differences in regression parameters between Alpine and 486 Mediterranean sites, these regressions were performed separately. The results are given in Figures 4a and  487 4b. Once these relationships are established, we can use the appropriate regression relationship to 488 generate automatically a reasonable median particle size for any soil in these two environments, as long 489 as the sand, silt, and clay fractions are available. Note that the particle size distribution is a function of 490 height in the soil column. This result implies that the soil texture is changing over time. Specific results 491 indicate that the studied soils either become finer over time, or do not change perceptibly (Fig. 5). The analysis, one should expect that those soils most severely impacted by diminishing particle sizes could be 496 a factor 10 0.53 = 3.4 thinner than what the predictions using a constant value of d 50 would generate. Even 497 use of a median grain size in such soils, rather than the final value, could still overpredict the depth by a 498 factor roughly as high as the square root of 3.4, or about 80%. 499

Infiltration 501
Infiltration is the parameter that is most difficult to obtain accurately, while it is also demonstrated below 502 to be the parameter, to which the predicted soil depth responds most sensitively. Strictly speaking, the 503 infiltration rate required here, the fraction of precipitation actually relevant for soil weathering, is what 504 penetrates to the bottom of the soil layer, and is the difference between precipitation and (actual) 505 evapotranspiration plus whatever surface water runs on to the site less the amount that runs off. 506 Infiltration is seldom measured. Some combination of stream flow data, base and storm, can be used to 507 estimate mean regional infiltration values, as long the fraction can be excluded (as "weathering 508 inefficient") that travels exclusively by overland flow. It is furthermore difficult to estimate the local 509 variability in this variable. Evapotranspiration data, (surface) runoff and infiltration rates were obtained During the cold (glacial) periods, climate in Central Europe and the northern USA was much colder (4 to 529 14 °C) and often drier (25 to 75 % of present-day's precipitation). In southern Europe and mid to southern 530 USA, the climate was distinctly colder and precipitation rates varied from slightly lower to higher, 531 depending on the area. 532 The present-day hydrologic mass balance was used for soils that started to form during the Holocene. For 533 older soils, climatic conditions also of the Pleistocene had to be considered. Because the percolation 534 theory approach uses the hydrologic mass balance, precipitation, overland flow and evapotranspiration 535 had to be estimated also for periods prior to the Holocene. Even though the climatic conditions were cold 536 to very cold at the alpine sites, it does not mean that no weathering has occurred during this period. very cold conditions, a high rate of soil formation is possible, provided that water is available. The 539 hydrologic mass balance is relevant for the percolation theory approach -temperature is therefore at first 540 sight less important. The present-day and averaged precipitation, evapotranspiration and infiltration rates 541 over the entire soil formation period are given in the supplementary Tables S1 and S2. 542

Damköhler Number Estimations 544
Our present procedure for determining Da I is less an independent calculation, and more of a scaling 545 estimate. We use the above referenced formula from Yu and Hunt (2017) as our starting point. In that 546 calculation, all the parameters were given in the publication describing the experiments (Salehikhoo et al., 547 2013). These experiments were performed on MgCO 3 , a mineral with dissolution rate that can easily be 548 orders of magnitude faster (Stumm and Morgan, 1996) than for the silicate minerals predominantly 549 responsible for the soil formation in the studies considered here. In fact, Fig. 13.9 from Stumm and 550 Morgan (1996) shows that, for dolomite, the reaction rate ratio can vary from around 3.5 orders of 551 magnitude at pH near 8 to more nearly 5 orders of magnitude at pH = 6. This is the range of pH values 552 expected for conditions with, or without, carbonate present as a solid phase. Such a contrast makes it 553 important to estimate a value of Da I relevant for the field conditions relevant to this study. 554 Under field conditions, parameters for surface area and equilibrium solution concentration are simply not 555 available over the entire range of sites investigated, as the surface area depends on the particle 556 morphology, and the equilibrium concentrations on their mineralogy. Consequently, we do not attempt to 557 represent any distinction in their values from those of the experiments of Salehikhoo et al. (2013). 558 Further, although surface area is a strong function of particle size; in the alpine regions studied here, at 559 least, their relatively coarse particle sizes make the soils a fairly reasonable physical analogue to the 560 media of Salehikhoo et al. (2013) where the typical particle size was 400µm, allowing us to estimate this 561 parameter by its experimental value as well. We thus concentrate on the contrasts between field and 562 experimental values for flow rates and well-mixed reaction rates. In the case of MgCO 3 , the reaction rate 563 used was R m = 10 − 9 mol m -2 s -1 . Four orders of magnitude slower for silicates leaves R m = 10 -13 mol m 564 -2 s -1 . The slower flow rates in the field, i.e., on the order of meters per year, instead of experimental 565 values measured in meters per day, will tend to mitigate the impact of the slower reaction rates. 566 If we keep unknown parameters including the equilibrium concentration and the Brunauer-Emmet-Teller 567 BET surface area at the same values as those in the MgCO 3 experiment of Salehikhoo et al. (2013), but 568 substitute in a well-mixed reaction rate, of 10 − 13 mol m -2 s -1 for R m , and a yearly average infiltration rate 569 of I = 0.7m yr -1 for the flow rate, it is found that Da I = 1 at about 11cm. For smaller length scales, 570 chemical weathering should be limited by reaction kinetics, instead of transport. 571 Note that the calculated value of Da I is sensitive to the particular choice of particle length scale; smaller 572 lengths increase the numerator but decrease the denominator, since the total surface area increases with 573 diminishing particle size. Thus, the impact of reaction kinetics in limiting chemical weathering diminishes 574 rapidly with diminishing particle sizes. Except for the Alpine sites, most sites have d 50 values well under 575 500µm, and such a complication from a kinetic-limited regime is much less likely. However, considering 576 only the range of particle sizes observed in the field, we estimate that there is probably close to an order 577 of magnitude uncertainty in our estimation of Da I in the alpine sites, and incorporating variability 578 expressed in the geographic regions investigated would probably tend to increase Da I from our calculated 579 value, even for alpine sites. A similar range of uncertainty may result from considerations of effects of 580 variations in pH on R m , while field pH conditions here may tend to reduce Da I from our calculated value. 581 Note that R m and d 50 are quantities for which some guidance for estimating the site-to-site uncertainty, at 582 least, exists. 583 584

Computations 585
Programming and calculations were done using the software R (R 3.3.2.). In order to facilitate 586 reproducibility, all R code and the required input data is made available under 587 https://github.com/curdon/soilDepth. 588 589 590 4. Results and Discussion 591

Alpine Sites 592
Soil depth and, therefore, also soil mass increase with time (Fig. 6). The rate of change, however, 593 decreases considerably with time. The factor time is a very strong predictor of soil mass but it does not 594 provide any further process specific details. Our model, which reduces chemical weathering of solute 595 transport in soils to a few physical parameters such as infiltration rate, particle size, erosion rate and 596 surface age, consistently overestimates soil depths at alpine sites by a factor approximately 1.45 (Fig. 7a). 597 The Mediterranean sites are shown in Fig. (7b), but are discussed below, except to note that their depths 598 are modeled quite precisely. Although soil mass and soil depth are two different parameters, their overall 599 trend over time is similar. 600 At short time intervals, the overestimation of alpine soil depths is particularly severe, as is also apparent. 601 Noting that the slope of the early (in time) data appears also to be steeper than predicted, we suggest that 602 one of the problems in applying Eq. 8 to find the soil depths of alpine sites is that at early times, the 603 chemical weathering is more severely limited by reaction kinetics than by solute transport. Indeed, our 604 most basic calculation (section 3.3) of the complications of alpine sites in calculating Da I , suggested that 605 for distances less than about 10 cm (corresponding to times of less than about 100 years), chemical 606 weathering is kinetic-limited. Comparing the predictions of a constant rate of soil development in the case 607 of kinetic-limited weathering leads to a linear dependence of soil depth on time at early times, consistent 608 with the trend of the data in Fig. (7a), for which a number of sites lie well below the predicted depths at 609 times less than 100 years. Further, the trend of these data is nearly linear, as evidenced by the 610 approximate thickness of 1 cm at ten years to about 10 cm at 100 years. The later slowing in the increase 611 of soil thickness with increasing age is compatible with a cross-over to transport-limited soil formation at 612 later times. Thus, the linear increase in depth with increasing time, while predicting a thinner soil than 613 what would be due to transport limitations at short times, would predict a deeper soil than is observed at 614 longer times. Soil depth is, however, reasonably well predicted for the Wind River Range sites. No data 615 for very young soils (< 100 yr) were available. The Wind River Range soils have, in contrast to the other 616 alpine sites, a small loess layer. The relatively small grain size of loess leads to slower infiltration rates, 617 longer contact time with minerals so that transport is limiting from the beginning of soil formation (what 618 can be better predicted with our model). 619 An alternate, but related, framework to address the overprediction of soil depths at short times is 620 discussed in the literature (e.g., Dixon and von Blanckenburg, 2012), and is based on the description of 621 the trend curves, the interrelation of erosion and mineral weathering and the designation of certain 'speed 622 limits' in soil production. In regions with solid bedrock as soil parent material, discussions of the possible 623 relevance of a "humped" soil production function can help provide an introduction to this topic, though 624 discussion of a humped soil production function also brings in additional complications. 625 Gilbert (1877) was the first to propose the existence of a "humped" soil production function, with rates of 626 bedrock to soil conversion diminishing both for thinner soils and thicker ones and, thus, a maximum at 627 intermediate depths. Such a soil production function can lead to two modes of stability depending on 628 initial conditions, one with, and one without soil mantle, making this also an environmentally relevant 629 model. There are two proposed causes, of which one is relevant here. The first owes to Gilbert (1877) and Wilkinson , 2007). This decrease is due to the gradually advancing consumption of fresh and reactive 642 minerals. The texture of the parent rock, and thus the available surface area, is an important factor in 643 controlling the reaction rates (Taboada and García, 1999). If however a moraine or other sediment layer is 644 present, not a humped but rather an exponential, power law or logistic function over time is observed 645 (Egli et al., 2014). Both types of function usually predict a maximum soil production rate under thin soils, 646 and an inverse relationship between soil thickness and soil production. Based on this, some speed limits 647 were derived. Transport limitation and kinetic limitation are often considered processes that determine the 648 speed of silicate weathering. A transport-limited weathering regime exists then when the supply of water, 649 inorganic and organic acids and ligands is large relative to the supply of silicate (West et al., 2005). The 650 supply of fresh minerals limits the weathering rates. Old and flat topographies typically are in this 651 category. A kinetic limitation is then given when the main control on chemical weathering, silicate 652 weathering ω depends on the kinetic rate of mineral dissolution W, the supply of material (e.g., by 653 erosion) e, and the time t available for reaction (West et al., 2005) 654 Mineral dissolution (W) depends on environmental conditions such as temperature and runoff (or 656 precipitation). The weathering rate of a mineral or rock decreases with the time that the mineral spends in 657 the weathering environment, as ! ⋉ ! with ω V = instantaneous volumetric weathering rate, λ = 658 (erosion) exponent (White and Brantley, 2003;West et al., 2005). Apparently, young surfaces or surfaces 659 exhibiting a high erosion rate (leading to a continuous rejuvenation of the surface and soils) are in this 660 category. 661 Our model overestimates soil production rates for young soils. The kinetic limitation of mineral 662 weathering, that is particularly effective at young sites, is most likely the cause for this. This seems to be 663 the case up to about 100 years but sometimes up to 1000 or more years of soil formation. Alpine soils 664 often have a coarse soil structure (Egli et al., 2001) that gives rise to a high hydraulic conductivity, 665 minimizing the opportunity for overland flow, and decreasing the reactive surface area. The model 666 predicts soil production rates (Fig. 8) for very young soils of about 10 -2 m yr -1 . In terms of soil mass 667 (using an average soil coarse fragment content of young soils of about 60% and soil density about 1.4 t m -668 3 ; Fig. 6), this would lead to a maximum soil production rate of ca. 5600 t km -2 yr -1 . Particles > 2 mm 669 (coarse fragments) are usually considered 'chemically inert' for plant growth, although they can very well We can also consider maximum soil production rates in the present and related models. An exponentially 678 declining soil production function, such as those considered by the Heimsath group, introduce an obvious 679 "speed limit" equal to the value of the pre-exponential factor, as found by substituting a soil depth of zero 680 into the exponential function. The inverse power-law formulation here generates a maximum conceivable 681 speed of soil formation equal to the infiltration rate, I. But that value is obtained for a soil thickness of 682 only a single particle. For typical infiltration rates (on the order of 1 m yr -1 ), I is about four hundred times 683 greater than the largest soil formation rates observed to date, around 2400 m/Myr (Larsen et al., 2014), 684 which are found in areas of precipitation with more nearly 10 m yr -1 in New Zealand. Clearly, if such a 685 thin soil could be defined, most precipitation rates around the world would effectively "cleanse" the 686 bedrock of the soil, rather than forming additional soil. For a soil with a thickness of one thousand 687 particle diameters, however, according to Eq. (6), with x/x 0 = 1000, the soil production rate has already 688 dropped by a factor 1000 -0.87 , or ca. 400 to about 2500 m Myr -1 . So, if layers thinner than some hundreds 689 of particles are not considered soils, the excessive "speed limits" implied by the model do not disqualify it 690 generally from being relevant, but may help to explain limitations on its relevance at very short time 691 scales, when weathering reactions are more likely to be kinetic-limited. Consider that a thickness of 10 692 cm for particles of 500 µm size represents a total of only 200 particles. Finally, since percolation theory is 693 a statistical theory, for soils much thinner than a few hundred particles, its quantitative predictions will 694 not be accurate anyway. This implies that for, say media less than 10 particles thick, percolation concepts 695 would not be relevant at all. In the range 10 to a few hundred particles, percolation predictions of soil 696 formation rates may be too high, especially when weathering may be kinetic-limited. 697 The minor over-prediction at larger time scales must also be addressed. Within our theoretical framework, 698 the three most likely candidates to explain this discrepancy are: 1) the systematic decrease in particle size 699 over time (Fig. 5), 2) the possibility that our estimations of infiltration are inaccurate, 3) the use of a too 700 small erosion rate. We consider these possible factors separately. 701 In constructing Fig. 7a, we used site-to-site variability in a median particle size, which incorporates 702 information on this distribution throughout the soil column, including the original material. This 703 particular value represents likewise an approximation to a dominant particle size over the temporal 704 evolution of the soil. For some sites, the implication of the vertical cross-section is that the median 705 particle size did not change over time. For others, it apparently decreased, in some cases by as much as a 706 factor 10. Trying to reconstruct a history of particle size evolution for each of the 147 alpine soils in the 707 database is unpractical, but we did construct an amalgamated history of a typical such site and an extreme 708 site. A variability for a typical site was determined by a power-law regression on the entire data set 709 represented as a function of time, and yielded d 50 = 165 µm t -0.123 , where t is measured in years. Use of 710 this function of t for d 50 for all sites yielded an R 2 value of 0.2 instead of 0.39, and did not visually 711 improve the agreement in the trend. A similar investigation of the extreme of using the smallest particle 712 size as a function of time was even less satisfying. Thus, it appears that temporal variability in median 713 particle sizes averaged across a wide range of sites is less important information than the site-to-site 714 variability in the vertically averaged profile at the observed time, at least for this study. If the Wind River 715 sites (many of them having very old surfaces ages) are excluded, then the predicted and modeled values 716 correlate better (R 2 = 0.48). One difficulty therefore might be due to less precise estimation of the 717 hydrologic mass balance over very long time periods. We, therefore, cannot exclude the possibility that 718 the discrepancy between theory and observation would be reduced if the temporal history of each soil 719 could be reconstructed by careful analysis of the vertical profile and environmental conditions. In this 720 context, it may be pointed out that Yu and Hunt (2017bc) used current values of the soil particle sizes, 721 corresponding to the smallest in the profile, and did not generate an overestimation. 722 Consider next the effects of uncertainty in AET. The global average actual evapotranspiration (AET) rate 723 is roughly 65% of precipitation (Lvovitch, 1973;Schlesinger and Jasechko, 2014), and can be 724 considerably higher. A 10% underestimation of AET, when its actual value is 80% can lead to a nearly 725 50% overestimation of the net infiltration rate I. 726 We expect that erosion rates are less likely to be the cause. While erosion rates higher than those used in 727 the model could also produce a reduction in predicted soil depth at longer times, we found no support for 728 choosing higher erosion rates. Further, the largest effect of a higher erosion rate would be to shorten the 729 time scale at which steady-state conditions (unchanging depth) are approached. Cursory examination of 730 the data suggests that the given suite of erosion rate values is in relative accord with this time scale. 731 The magnitude of the observed scatter at any given time is reasonably well accounted for by the 732 variability in predicted values that arises from the variations in the input parameters. However, 733 comparison of all individual predicted and observed depths yields a somewhat unimpressive R 2 value of 734 about 0.39. Thus, we conclude that, in Alpine soils, it is possible that for small scale predictions local soil 735 rearrangement from such processes such as bioturbation need to be incorporated into our model of soil 736 production. For somewhat larger scales, a statistical correspondence could suffice. 737 738

Mediterranean Sites 739
For the Mediterranean sites, predicted soil depths match quite well observed soil depths (Fig. 7b). There 740 is no obvious trend towards greater or smaller underestimation at longer time intervals. The scatter in 741 experimental data is also well captured by the scatter in predictions, as seen in Fig. (7b). However, direct 742 representation of the individual predictions as a function of observed depths reveals a relatively large 743 scatter, and an R 2 value of only about 0.5. All sites have a predominantly silicate geology. If the sites 744 having some carbonate additions are excluded, then this scatter is slightly reduced and R 2 increases to 745 0.59. Therefore, the diversity of the geology also seems to affect to a minor extent the performance of the 746 model. Furthermore, it is to be presumed that: 1) there is considerable uncertainty in the determination of 747 the parameters, which could be due to incompletely guessed surface flow routing, or 2) other influences 748 on soil depth that are important have been omitted. Since the mean depth as a function of age is relatively 749 well predicted, in the case of 2), we suggest that the best candidates for neglected processes should be 750 those that contribute to local reorganization of the soil, such as bioturbation, tree throw (Šamonil et al., 751 2017), eaolian deposits (after the start of soil formation), petrography, other hillslope processes or 752 unconsidered human impact. The Mediterranean soils have a typical soil density of about 1.6 t m -3 and a 753 soil coarse fragment content of about 20%. Production rates for soils having an age between 1 and 10 kyr 754 are in the range of about 70 to 380 t km -2 yr -1 . Very young soils (age of 10 to 100 yr) have rates of up to 755 1000 t km -2 yr -1 . Although the climate is warmer (but not moister), weathering and soil production are not 756 significantly faster than in alpine areas. 757 It is, furthermore, interesting to note from examination of Figures 7a and 7b that the model better reflects 758 the variability of the data of both systems, Mediterranean and alpine soils. Furthermore, soil erosion 759 (output) and soil production (input) appear to be in a similar range after about 100 kyr as seen in Fig. 8a  760 and 8b. We therefore hypothesize that soils tend to a quasi steady-state after about 100 kyr. We however 761 have to underline that soils are dynamic and that a change in the settings may cause a distinct shift away 762 from steady-state conditions. As already Phillips (2010) mentioned, steady state soil thickness is only 763 there applicable, where sufficient time has elapsed for regolith accumulation, and where effects of 764 processes other than weathering and surface removals on thickness are negligible. 765 766

Sensitivity analysis 767
Owing to the simplifications made to run the model, it is necessary to know which effect a change in 768 these parameters has on the results. The influence of the parameters on soil depth and soil production (Eq. 769 8) was checked by a sensitivity analysis. A convenient way to express the sensitivity of a system to 770 changes in parameters is the normalized sensitivity coefficient, Q i,m , that is defined by 771 where C i is the dependent variable (in this case, it can be either soil production or soil depth) and P is the 773 value of parameter m (Furrer et al., 1990). This formula expresses a partial derivative based on the percent 774 change in the dependent variable caused by a 1% change in the value of the parameter m. Sensitivity 775 coefficients are negative for parameters in processes that lead to a reduction and positive for parameters 776 that cause an increase in variables of interest (Table 1). Processes that lead to a deeper soil exhibit 777 positive values: a greater particle size, a higher infiltration rate and, not surprisingly, a longer soil 778 evolution. Due to the fact that the rate of soil depth increase decreases with soil evolution, an increase in 779 age has a negative effect on the soil production rate, both for the alpine and Mediterranean sites. The 780 strongest influence of a change in the infiltration rate on soil depth and production rate is particularly 781 noteworthy. The estimation of the water fluxes and infiltration rates becomes particularly difficult for old 782 surfaces (having an age > 20 kyr). Although the paleoclimate studies give some hints about the 783 hydrological conditions over time, they are all not precise enough and only enable a rough, regional 784 estimate at its best (Minnich, 2007 in soil depth. However, the sensitivity even for the infiltration rate I is still < 1. 787

Most closely related model 788
We have discussed the exponential soil production model in several contexts above. We have also 789 mentioned some of the difficulties of models of chemical weathering based on the advection-dispersion 790 equation under conditions of full saturation. However, there is one related model due to Braun et al.

791
(2016), that is rather closely related to ours, and that deserves separate mention. These authors investigate 792 regolith production from chemical weathering in the context of unsaturated flow, which, in turn, is 793 addressed using Philip infiltration theory (Philip, 1967(Philip, , 1968 conditions show the same time-dependent scaling and proportionality to flow rates as observed in the 808 field. Thus, the percolation scaling theory yields the same time-dependence for unsteady flow under 809 unsaturated conditions as for solute transport under either saturated or wetting conditions, but the validity 810 of the Philip infiltration model explanation is restricted to unsaturated conditions and is incompatible with 811 identical results obtained under saturated conditions. 812

Conclusions 813
We applied a theoretical model of soil formation rates as limited by chemical weathering to prediction of 814 the evolution of soil depths at about 250 sites in two climate regimes, Mediterranean and alpine. The 815 chemical weathering was assumed to be limited by advective solute transport, a result which made 816 relevant the rate of the mean vertical flow through the soil column over time. These rates were estimated 817 using the difference between precipitation and evapotranspiration, but without accounting for the effects 818 of topography on the divergence of surface water flow. The model requires, as a second input parameters, 819 a characteristic particle size, which we assumed could be equated with d 50 , equivalent to a fundamental 820 soil network spatial scale. In most cases this value had to be obtained from a regression function relating 821 observed soil sand, silt, and clay fractions to observed d 50 values in the same climatic region. The other 822 two input parameters are the erosion rate and the time. The output depth, generated from the model was 823 then compared with the observed depth on a site-by-site basis. We found that the model matched soil 824 depths in the Mediterranean climate zone, but overpredicted depths in alpine regions by approximately 825 45%. The results suggest that it is possible to model soil formation rates with our analytical formula and 826 generate the major trends. Further, the range of predicted soil depths, at any given time is generally in 827 accord with our predicted variability. If this variation could be traced to an accurate representation of the 828 actual variability in the net infiltration rate I, erosion rate E, and the characteristic particle size d 50 at these 829 times, one would expect a very high R 2 value. Unfortunately, our two direct comparisons between 830 predicted and observed soil depth yielded R 2 values of 0.39 for the Alpine soils and 0.50 for the 831 Mediterranean soils. Such relatively low values of R 2 indicate the presence of a significant discrepancy. 832 While there are a number of possible implications, including that the theoretical model leaves out 833 important local redistributions of soils, such as due to bioturbation, we think that the first place to check 834 for a source for the discrepancy is in the estimation of the parameters. Systematic overestimation of soil 835 depths at later times will arise from exclusion of, or only partial accounting for, the fining of particle size 836 with time. But systematic overestimation of soil depth can also arise from random errors in actual 837 evapotranspiration (AET) estimation, when AET is significantly larger than 50% of the precipitation. For 838 example, if AET is 80% of the total precipitation, an underestimation of AET 10% can make an almost 839 50% overestimation of the value of I, while an overestimation of AET of the same magnitude leads to a 840 100% underestimation of I. Beyond systematic errors, neglect of surface water routing will fail to 841 generate any patchiness of soil depth due to differences in local infiltration rates. 842 In the future, we would like most importantly to improve the model performance with respect to its ability 843 to generate the actual variability in soil depths of a given age, not just the envelope of such soil depths. 844 We think that a logical means to address this issue is to solve for the variations of soil depth along 845 hillslopes and within individual basins. To accomplish this, we must first develop a protocol that yields a 846 more accurate value of I, which is possible only through integration with models that accurately describe 847 surface water routing. Although it would be useful as well to couple with hillslope subsurface flow 848 models, including dependences on surface topography, this particular improvement has a lower priority.

849
A goal is to be able to apply our soil production model at regional and larger scales using maps of 850 relevant data. Describing the subsurface hydrologic processes requires another level of input data, 851 including hydraulic conductivity variability, more difficult to access and predict. It will also be important 852 to incorporate into the numerical routine information regarding any known time dependence of 853 parameters, such as precipitation or erosion. With these advances, it is hoped that we can meet the 854 ultimate goal of the research, to transform mapped data regarding climate, erosion rates, and soil 855 characteristics and so forth into a reliable and compatible representation of predicted soil depths.