Colluvial deposits as a possible weathering reservoir in uplifting mountains

Colluvial deposits as a possible weathering reservoir in uplifting mountains Sébastien Carretier1, Yves Godderis1, Javier Martinez2, Martin Reich2,3, and Pierre Martinod2,1 1GET, IRD, CNRS, CNES, OMP, UPS, Université de Toulouse, 14 avenue E. Belin, F-31400, Toulouse, France 2Department of Geology, FCFM, University of Chile, Santiago, Chile 3Andean Geothermal Center of Excellence (CEGA), FCFM, University of Chile, Santiago. Correspondence to: S. Carretier (sebastien.carretier@get.omp.eu)


Introduction
Mountain building steepens the relief and drives river incision and hillslope erosion. Fresh minerals are exposed by erosion and are transformed into grains of different sizes that are subjected to weathering by meteoric water. The weathering of silicates, in particular, consumes atmospheric CO 2 over timescales of many thousands to millions of years (Walker et al., 1981). This erosion-driven weath-20 ering process led to the debated "uplift" theory, in which mountains play a key role in regulating the global climate over geological time (Raymo et al., 1988). Soil column models have challenged this theory by predicting that above a certain erosion rate value, minerals do not stay in the regolith long 1 Earth Surf. Dynam. Discuss., https://doi.org/10.5194/esurf-2017-48 Manuscript under review for journal Earth Surf. Dynam. Discussion started: 28 August 2017 c Author(s) 2017. CC BY 4.0 License. enough to significantly weather, producing a hump in the weathering-erosion relationship (Dixon et al., 2009b;Ferrier and Kirchner, 2008;Hilley et al., 2010;Gabet and Mudd, 2009). Other models 25 argue that when the regolith vanishes at large erosion rates, weathering becomes significant in the fractured bedrock (Maher, 2011;Calmels et al., 2011;West, 2012), or that high reliefs consume more CO 2 than low reliefs during wetter periods (Maher and Chamberlain, 2014). Datasets from soil pits and riverine fluxes show a monotonic relationship between both the denudation rate and weathering rate in some cases (Millot et al., 2002;Dixon et al., 2009b;West et al., 2005), but also 30 evidence a possible maximum erosion rate above which the weathering rate decreases (Dixon and von Blanckenburg, 2012). Recent data from the Southern Alps in New Zealand have challenged the existence of this erosion rate limit by demonstrating that weathering was able to continue increasing at the highest erosion rates when rainfall is abundant (Larsen et al., 2014). In such regions, landslides constitute a significant weathering reservoir (Emberson et al., 2016a, b). Downstream from 35 the Andes and Himalaya, transported minerals may continue to weather significantly in the floodplain Bouchez et al., 2012;Moquet et al., 2016). As a result, the debate on the locus of weathering in mountains is still open and different weathering reservoirs from the hillslopes to plains may dominate at different stages of the mountain evolution. Until now, four main weathering reservoirs have been identified: soils (Dixon et al., 2009b), fractured bedrock (Calmels et al.,40 2011), basins (Bouchez et al., 2012), which also trap a considerable amount of organic carbon (Galy et al., 2015), and oceans (Oelkers et al., 2011). In this paper, we address the particular question of the relative contributions of in situ produced regolith and colluvial deposits in the weathering outflux of an uplifting mountain under a cooling climate. 45 None of the available models are able to discriminate between these weathering reservoirs. Moreover, few models (Vanwalleghem et al., 2013;Braun et al., 2016) account for the heterogeneity of erosion and weathering during relief adaptation to uplift which may control the overall evolution of the weathering rate of a rising mountain range (Anderson et al., 2012;Cohen et al., 2013;Carretier et al., 2014). None of these models can be used to trace the weathered material through its stochastic 50 displacement from the hillslopes to the basins.
We therefore developed a dynamical model (Cidre) that accounts for the heterogeneity of the erosion and weathering evolution under a scenario of climate change. This model uses a novel approach that couples the landscape evolution and moving clasts, which can be used to follow the weathered 55 material through different weathering reservoirs. By making a link between the weathering processes at the mineral, hillslope and river scales, we provide new insights into the particular effect of valley widening and the associated colluvial deposits on weathering rates in uplifted areas.
The elevation z (river bed or hillslope surface) changes on each cell (size dx) according to the balance between erosion [LT −1 ], deposition D [LT −1 ], sediment discharge per unit length from 70 lateral (bank) erosion q sl [L 2 T −1 ] and uplift U [LT −1 ] (e.g. Davy and Lague, 2009): and we define (Davy and Lague, 2009;Carretier et al., 2016) r = Kq m S n for river processes (2) h = κS for hillslope processes (3) 75 where K [T −0.5 ] and κ [LT −1 ] are lithology-dependent (different for bedrock or regolith/sediment) erosion parameters, S is the slope, q [L 3 T −1 ] is the water discharge per stream unit width, m and n are positive exponents, and D r = q sr ξq for river processes (4) D h = q sh dx 1−(S/Sc) 2 for hillslope processes (5) 80 where q sr and q sh are the incoming river and hillslope sediment fluxes (total q s = q sr + q sh ) per unit width [L 2 T −1 ], ξ is river transport length parameter [T L −1 ] and S c is a slope threshold. The deposition fluxes on a cell are a fraction of the incoming sediment. When the local q and S values are larger, less sediment eroded from upstream will deposit on the cell. The sediment leaving a cell is spread in the same way as water, i.e. proportionally to the downstream slopes. 85 Flowing water in each direction can erode lateral cells perpendicular to that direction. The lateral sediment flux per unit length q sl [L 2 T −1 ] eroded from a lateral cell is defined as a fraction of the river sediment flux q sr [L 2 T −1 ] in the considered direction (e.g. Murray and Paola, 1997;Nicholas and Quine, 2007): where α is a bank erodibility coefficient. α is specified for loose material (regolith or sediment) and is implicitly determined for bedrock layers proportionally to their "fluvial" erodibility such that α loose /α bedrock = K loose /K bedrock (K from Equation 2). If a regolith or sediment covers the bedrock of a lateral cell, α is weighted by its respective thickness above the target cell.

95
Following different authors, we assume that the regolith production rate follows a humped law, so that there is an optimum thickness at which the regolith production rate is at its maximum (Strudley et al., 2006) 100 where d 1 and d 2 [L] are the attenuation depths, k 1 is a non-dimensional coefficient, and w o [L T −1 ] is the regolith production rate for the exposed bedrock.
We also let w o depend on temperature T and precipitation rate P , following White and Blum (1995) and Dixon et al. (2009a) among others: where k w is a factor with the dimension of a weathering rate [L T −1 ], P [L T −1 ] is the amount of water entering the regolith (equal here to the rainfall rate), P o is the water flow reference value (1 m a −1 in this study), E a is the activation energy corresponding to the mineral that controls the weathering front advance, T o is a reference temperature (298 K), T is the local temperature expressed in 110 Kelvin and R is the gas constant.
Equations 7 and 8 are similar to the regolith production model tested by Norton et al. (2014), and were also used in Carretier et al. (2014). Nevertheless, except in rare studies (Wilkinson et al., 2005), no data fully support (or exclude) the humped function in Equation 7. The existence of an op-115 timum regolith thickness has been conceptually justified as resulting from water pumping by plants, or an optimum residence time of water within a porous soil to dissolve minerals (Gilbert, 1877 depend on groundwater discharge (e.g. Maher, 2010;Maher and Chamberlain, 2014;Hilley et al., 120 2010; Rempe and Dietrich, 2014). Lebedeva et al. (2010) and Braun et al. (2016) showed that in the absence of uplift and erosion, this type of model predicts that the regolith thickens as √ t. This means that the regolith production rate w varies as 1/B. Compared with the exponential trend of Equation  (Braun et al., 2016). Alternatively, if the weathering front advance is controlled by the rate of mineral fracturing, then the regolith production rate is predicted to be constant (Fletcher et al., 2006).
As we are analysing the effect of erosion on weathering, the 1/B depth attenuation for thick regoliths is not considered here. Nevertheless, we show one experiment that uses this 1/B attenuation to illustrate that the particular model for local regolith development is not crucial for our conclusions.

130
The dependence of the local regolith production on the temperature and precipitation rate is suggested by silicate weathering outfluxes on the soil and catchment scale (White and Blum, 1995;Gaillardet et al., 1999;Dupré et al., 2003;Oliva et al., 2003;Riebe et al., 2004;West, 2012), and experimental data (Brantley et al., 2008). Dixon et al. (2009a) showed that Equation 8 fits saprolite 135 data in the western Sierra Nevada Mountains in California. Nevertheless, Maher (2010Maher ( , 2011 and Maher and Chamberlain (2014) argued that in many mountainous situations, the weathering rate should essentially and linearly depend on the water flow in the soil, with a minor effect of temperature. In Equation 8, the linear dependency between the regolith production rate and runoff and the weaker dependence on temperature are consistent with that view, although our model clearly misses 140 the control of water flux partitioning between the surface and ground on the regolith development rate and pattern (Maher and Chamberlain, 2014;Rempe and Dietrich, 2014;Braun et al., 2016;Schoonejans et al., 2016). The drawback of Equations 7, 8 is that they are parametrical and not truly physically based. Nevertheless, given the lack of consensus, we assume the form of these laws and test the effect of varying their parameters. In particular, we test the difference between the humped 145 form and the exponential form (k 1 = 0) of this law (Heimsath et al., 1997).

Clast weathering
A clast has a specified radius r, with no particular limitation, between the size of a small mineral and a large cobble. Its probability to be detached, deposited or to pass through a cell depends on its size and on the associated fluxes calculated by Cidre on each cell (see Carretier et al., 2016). With 150 this algorithm, the spreading of the different clasts depends on the relative magnitude of diffusive and advective transport (Equations 2 to 5), while the mean population transport rate is determined by the transport discharge q s calculated in Cidre (Carretier et al., 2016 where V m is the molar volume [L 3 N −1 ], λ is a non-dimensional roughness coefficient defined by White and Brantley (2003) (see also White et al., 2008) and k m is a dissolution parameter depending on each mineral [N L −2 T −1 ] (e.g. Brantley et al., 2008). The other parameters are defined in Equation 8. The product λ4πr 2 m is the reactive surface. The second part of this equation comes from 165 experimental laws of mineral dissolution (e.g. Brantley et al., 2008). The first part, with the runoff dependence ( P Po ), accounts for the linear increase in the dissolution rate with water discharge (Maher, 2010). As for the regolith production law, we acknowledge that it is a simplification to assume a linear correlation between groundwater discharge in the soil and runoff.

170
If a clast is made up of different minerals, their proportions are specified at the beginning (χ mo ) and then evolve during their lifetime (χ m ). Modelling a complex mineralogical texture with a subdivision into different grains would be intractable in practice. In a simplified model, the main issue is to define a reactive surface for each mineral type. Given that minerals are spread into different grains within the clast, this surface is larger than the surface of a simple sphere made up of a particular 175 mineral (White and Brantley, 2003).
In order to take this reactive surface into account in the most simple manner, each mineral type is converted to an "equivalent" sphere with radius r m including the mineral and "virtual" vacuum.
In addition, this sphere surface is multiplied by the roughness coefficient λ to define the reactive 180 surface ( Figure 1). The sphere geometry is chosen for its simplicity. The "virtual" vacuum implies that the reactive surface of each mineral is larger than the surface of a "solid" sphere made up of this mineral only (even if λ = 1). λ is an adjusted factor that may account for the complex geometry of the crystals within the clast. This formulation also respects the fact that when smaller proportions of a mineral occur within the clast, its specific surface is larger (reactive surface over mineral mass).

185
At the beginning of the process, each "equivalent" sphere has the same radius as the clast ( Figure   1).
The total dissolution rate for the clast w c [L 3 T −1 ] is then The "solid" (real) volume lost by each mineral δv m is subtracted from its previous volume to calculate the new mineral volume v m . The new mineral radius r m is then calculated considering an "equivalent" sphere incorporating the solid and virtual vacuum of volume vm χmo : The sum of the new "solid" mineral volumes v m is the new "solid" clast volume. The new clast 200 radius is the radius of the largest "equivalent" sphere of its constitutive minerals ( Figure 1). Doing this, we assume that the largest mineral forms a mass that includes the other minerals.
This formulation was also designed to respect a basic mass balance: a clast of a given size constituted initially of one mineral (for example 100% of albite) evolves exactly in the same way as if 205 it was constituted of different proportions of the same mineral (for example 40% of albite + 50% of albite + 10% of albite). This equivalence is achieved by using the initial mineral fraction χ mo in Equation 13 to define each "equivalent sphere".
If a clast includes minerals of contrasting weathering rates, for example albite and quartz, the 210 rapid dissolution of albite ends up with a porous clast made up of a vacuum within the quartz. The true clast volume is therefore larger than the "solid" sphere corresponding to the mass of the quartz only, which is reproduced well by Equation 13. The porosity increases in these clasts, consistent with reality. Obviously, this approach supposes that the initial clast does not lose its cohesion and is not divided into different mineral grains, which can occur in nature.

215
This approach is probably less realistic when the clast size exceeds several centimetres. In that case, the advance of the annular weathering front may control the clast volume that is effectively being weathered (Lebedeva et al., 2010). This type of front could be simply introduced into the  Lebedeva et al. (2010). In the present form, the predicted 220 dissolved volume of such large clasts is therefore probably a maximum volume.
The weathering rate of the clasts does not depend on the depth in the regolith in the present version. The removal of water by plants, changes in the regolith porosity, pCO2, groundwater flow velocity, pH, sensitivity of the surface temperature variations, clay precipitation, etc., can modify the 225 weathering rate according to the depth. In the future, this could be accounted for by, for example, modulating the weathering rate by the same humped law (Equation 7) used for regolith production (Vanwalleghem et al., 2013).
Nevertheless, to validate our approach in a simple case, we simulated weathering on marine ter-230 races in Santa Cruz, California and found that the results agree with empirical observations (Supplementary Material).

Integration at the cell and grid scales
The philosophy developed here is to use the available limited information provided by the clasts present on some of the cells to estimate the chemical weathering outflux of the landscape. Weather- The dissolved chemical flux per layer w l [L 3 T −1 ] is the clast dissolution rate per clast volume wc vc , multiplied by the layer volume v l . This value is also weighted by the ratio vc vo between the current and initial clast volumes in order to take the fact that the volume of weathering material decreases 250 within the layer into account: The corresponding dissolved flux is also calculated for each mineral constituting the clasts, as 255 well as for some of their elements according to their stoichiometric proportions in the minerals. If a clast is completely dissolved, it remains in the regolith or in the deposit where it is trapped so that there is an integration layer around it which produces zero chemical flux. This allows to account for the depleted layers of the soil. A totally dissolved clast is killed as soon as it is detached from the regolith. The dissolved chemical flux of the entire cell w cell is obtained by summing w l ( Figure 2).

260
Finally, the total weathering outflux W [L 3 T −1 ] integrated over the model grid is weighted by the relative proportion of regolith or sediment that contain clasts:

Clasts revival
A clast is killed when detached after complete dissolution, or if it simply goes out of the model 265 grid. In both cases, the model offers the possibility to recycle the clast. It is put back into the same cell in which it was initially seeded, with the same characteristics, randomly, between depths 0 and B (regolith thickness) within the parent rock below the regolith (except if B = 0 in which case the revival depth is set at 10 m to avoid the handling of numerous clasts where weathering is null). The maximum depth B at which the clasts are repositioned is optimal in order to favour an equidistance 270 between the clasts within the regolith. Recycling a dead clast to its initial location also permits to densify the number of clasts where the exhumation is faster. By this approach, a limited list of clasts is handled, while optimising their distribution at depth to obtain the best estimate of the chemical outflux. This is particularly efficient and useful when the landscape is uplifting and exhumes clasts.

275
Assuming m = 0.5 and n = 1 (Whipple and Tucker, 1999) and using the scaling factors H for mountain height, L for mountain width, P for effective precipitacion rate (runoff) and U uplift rate, we obtain the non-dimensional form ( * ) of the mass balance equation 1: These numbers affect the morphology of the resulting topography at steady state. Smaller N riv , N depo and larger N hill and S c values yield topographies that are increasingly dominated by (diffusive) smooth and rounded hillslopes.
Following the same approach, the non-dimensional form of the regolith thickness variation  plementary Material) yields the number N reg = wop U . This non-dimensional number determines whether or not regolith exists at dynamic equilibrium. N reg > 1 produces a regolith-covered mountain, whereas N reg < 1 leads to a bare-bedrock mountain (Carretier et al., 2014).
Finally, the dimensional analysis of the clast dissolution rate leads to a non-dimensional number 295 N clast = τr τm , where τ r is a clast's residence time in the regolith at steady state (τ r = B/U ) and τ m is the characteristic dissolution time of the main mineral (here albite) defined as the time necessary to decrease the initial clast volume by a factor e 1 (see the Supplementary Material for the full expression including the model parameters). This number includes the clast size and the kinetic parameters associated with the particular mineral m. N clast is actually the Damköholer number (e.g. White and 300 Brantley, 2003;Hilley et al., 2010). This number indicates the weathering grade of a clast leaving the hillslopes. When N clast is large, a clast leaving the regolith is very depleted, whereas it remains fresh if N clast is small. The first situation has been called "supply" or "transport" or "erosion" limited weathering (e.g. Dixon et al., 2009a). The second situation has been called a "kinetically" limited regime (e.g. Ferrier and Kirchner, 2008). Hilley et al. (2010) identified N reg (" * " for them) and 305 N clast ("D i " for them) as key parameters controlling the weathering flux at the scale of a soil column.
It is worth noting that experiments sharing different model parameters but the same non-dimensional numbers give similar results (Supplementary Figure S1). Consequently, the complexity of this model 310 is actually reduced to seven non-dimensional numbers reflecting a great diversity of natural climatic, weathering and erosion situations.

Model parameters that matter
The number of parameters that matter in this contribution can be reduced to three, namely the valley widening parameter α, the Damköholer number N clast and the uplift-to-weathering number N reg .

315
The other four non-dimensional numbers S c , N riv , N depo and N hill affect the final relief, drainage density and hillslope roundness, and the response time for denudation to reach the uplift rate value. The main result of this contribution does not depend significantly on these parameters. On the contrary, α, N clast and N reg determine the time spent by the clasts in the different weathering reservoirs (regolith, colluvium, valley). Thus, we primarily vary the parameters included in these numbers to 320 study the different behaviours of the model with regards to the long-term trend of the weathering outflux.

Reference experiment WARM
We design a reference experiment, WARM, corresponding to a warm, wet and constant climate. An initial horizontal rough surface (σ = 0.5 m) is uplifted. Sediment can leave the southern boundary 325 but not the northern one (equivalent to a divide), and the two other sides are linked by periodic boundary conditions. The resulting half-mountain is 100 km wide and 150 km long (dx = 500 m), similar to the length of the Himalayan or Andes catchments. Rivers do not erode laterally (α = 0), uplift rate U (1 mm a −1 ) and the precipitation rate P (1 m a −1 ) and temperature T at z = 0 m (25 o C) are kept constant. We fix the erodibility parameters (K and κ) so that the maximum elevation 330 at steady state reaches a reasonable height ∼ 7000 m consistent with the Andes and Himalaya, on a time scale of ∼ 15 Ma: K = 1.5 10 −4 a −0.5 (bedrock) and 2.5 10 −4 a −0.5 (regolith or sediment) are within the range of previous estimates (Giachetta et al., 2015), κ=10 −4 m a −1 and S c = 0.84 (= tan 40) for both bedrock and loose material, and ξ = 0.1 a m −1 (Table 1).

335
In the regolith production law (Equation 8), E a = 48 kJ mol −1 is intermediate between albite and biotite, the minerals which control the weathering front advance. The reference temperature T o = 298.15 K. The humped attenuation parameters are from Strudley et al. (2006), with d 1 = 0.5 m, d 2 = 0.1 m and k 1 = 0.8. We acknowledge that these parameters are empirical and are not necessarily representative of chemical weathering. We come back to this point later. With these values, the regolith 340 production rate w is optimal (w op ) for B = 0.17 m. The parameter k w = 0.003 m a −1 is chosen so that N reg = 1.7, a value >1 implies that the hillslopes are mantled by a 0.55 m thick regolith at dynamic equilibrium.

Regolith-covered mountains
We run the simulation WARM until erosion balances uplift, which occurs when the maximum el-355 evation is ∼7000 m after ∼15 Ma of simulation. Figure 3 shows that during the adaptation of the topography and erosion to the imposed uplift, the mean regolith production rate increases. As predicted by the value of N reg > 1 (= 1.7), the mountain is covered by a ∼ 0.5 m thick regolith at dynamic equilibrium. The mean regolith thickness reaches a maximum in the early times of the surface uplift when erosion is still low in average. Then the regolith thickness decreases as the drainage The weathering rate is also plotted against the total denudation rate and compared to data in Figure 4, as well as the parametrical model of West (2012). Other experiments are plotted, in which 370 we vary the uplift (U / 10 to x 5), temperature (T / 1.4), precipitation (P / 5), clast size (r = 0.25-5 mm) and mineral roughness (λ x 160) and mineralogy (granitoid or pure albite). For experiments with larger N clast values than the reference experiment (mainly supply limited), the denudation and weathering rates fit the linear relationship observed for regolith-covered landscapes characterized by supply-limited weathering (e.g. Dixon et al., 2009a). For experiments with a smaller N clast value, 375 the weathering becomes progressively kinetically limited as in the reference experiment, and thus saturates at high denudation rates (Figure 4). Overall, Figure 4 shows that the range of weathering and denudation rates predicted by these different simulations through time fits a large range of data for regolith-covered mountains.

380
In the following, we explore the response of weathering in the case of a cooling and drying climate. In a first set of experiments (COOLING and OROGRAPHIC), the mountains become entirely bedrock at the end of the cooling period. These experiments represent an end-member model as pure bare bedrock mountains probably do not exist (Heimsath et al., 2012). Nevertheless, this case is useful to quantify the effect of colluvium temporally stored on valley borders. In a second time, The modelled cooling operates through a decrease in temperature at the mountain foot in four arbitrary steps of -2 o C at 3, 6, 9 and 12 Ma of the model time. Furthermore, a temperature gradi-390 ent of -6 o C km −1 is prescribed. Moreover, in order to account for the drying potentially associated with global cooling, we add a rainfall decrease of -5% per degree of cooling (Labat et al., 2004;Maher and Chamberlain, 2014). Figure 5 shows the response to this climate change for the COOL-ING experiment, which uses the same parameters as the previous WARM reference experiment.
During the mountain rise, progressive cooling and drying decreases the regolith production rate and 395 the clast weathering rate. Moreover, the erosion rate increases. During the first several millions of years, the weathering flux increases because there is a sufficient landscape area covered by regolith characterised by supply-limited weathering. Then, the regolith cover decreases dramatically and the weathering flux falls to zero everywhere in the bedrock mountain ( Figure 5). 400 We now allow valleys to widen by lateral erosion for the same experiment. The factor α controlling the widening rate of the valleys is poorly constrained. Nicholas (2013) used a slightly different equation for lateral erosion q sl = (E S l ) q sr , where S l is the lateral slope. He calibrated E between 1 and 10 in large alluvial rivers. With lateral slopes S l on the order of 0.01, α ranges between 0.01 and 0.1 for sediment. We use a lower reference value α = 0.001 for regolith or sediment, probably 405 better adapted to large pixels. Allowing rivers to erode laterally, the weathering rate follows the same initial evolution but then, it does not fall to zero ( Figure 5). Valley widening steepens the foot of the hillslopes on the borders of the valleys, which generates mass wasting and the deposition of fresh material on the valley borders ( Figure 5c). This fresh minerals weather before being removed by rivers. Colluvium reside long enough in valleys (99% of the clast residence times are smaller than 410 1500 years in the COOLING experiment with lateral erosion - Figure 6) to generate a significant and nearly constant weathering rate. As there is no remaining regolith on the hillslopes even at low elevations (erosion exceeds the regolith production rate), colluvium are the only loose material producing a weathering flux. The prescribed drops in rainfall have a limited impact on the weathering flux of the colluvium ( Figure 5). Indeed, a decrease in water discharge increases the colluvium residence 415 time (Figure 6), which counterbalances their lower weathering rate. Consequently the weathering flux reaches a steady-state when an equilibrium is reached between the rate of colluvium removal and their weathering rate. Dividing the lateral erosion parameter α by two and five only decreases the weathering flux by a quarter and a half, respectively ( Figure 5a). As soon as α is large enough for the valleys to widen and to drive mass wasting, the volume of the colluvial deposits depends

Cooling and mountains partially covered by regolith
Previous cooling models led to bare bedrock mountains. Nevertheless, soils almost always cover the bedrock at low elevations. In order to generate a more realistic regolith distribution, we only 440 increase the bedrock weatherability by increasing the parameter k w in the regolith production law (Equation 8). In the resulting COOLING REG. experiments, a regolith persists at low elevations (< 1500 m) when the climate is cooler and drier. In this case, the persistent regolith is able to sustain the weathering rate of the whole mountain at a significant value (Figure 9). Adding valley widening (α = 0.001) does not significantly modify the weathering flux evolution during the last 10 445 Ma. Nevertheless, colluvium still play a role in that case. Because the hillslopes steepen near valley borders, the area covered by regolith is reduced by half compared to the case without lateral erosion.
Yet, the weathering flux is similar with and without lateral erosion, which shows that colluvium account for half the weathering flux. This fraction is also observed for the OROGRAPHIC REG. experiments that use a Gaussian rainfall-elevation relationship ( Figure 10).

Other regolith production laws and pixel size
We made assumptions about the regolith production law. Yet, the form of the production law controls the spatial distribution of the regolith thickness and production rate in a mountain (Carretier et al., 2014;Braun et al., 2016). We thus test the robustness of our main result by assuming an exponential regolith production rather than a humped law. We do this by setting k 1 = 0 in Equation 7. In order 455 to compare experiments, we also decrease k w (Equation 8) so that the maximum regolith production rate for bare bedrock equals the optimum regolith production rate calculated using the humped version of the law (same N reg at z = 0 m - Table 1). Figure 11 compares the OROGRAPHIC ex-14 Earth Surf. Dynam. Discuss., https://doi.org/10.5194/esurf-2017-48 Manuscript under review for journal Earth Surf. Dynam. Discussion started: 28 August 2017 c Author(s) 2017. CC BY 4.0 License. periments using the humped law with the same experiments using the exponential law. For regolith thickness larger than the optimum thickness, the regolith production is faster for the humped law.

460
Consequently, the total volume of regolith produced with this law is larger. The weathering flux is thus greater with the humped law while some regolith remains.
We now assume that the regolith thickness decreases as 1/B for thicknesses larger than the optimum thickness (0.17 m), instead of exponentially (Braun et al., 2016). This different law produces a 465 much thicker regolith of several tens of meters in the early stage (∼ 1 Ma) of the mountain erosion.
This rapid regolith thickening generates a weathering peak, but which is only twice that produced with the humped law ( Figure 11). Indeed, the thick regolith is rapidly depleted, so that only the weathering of its deeper layer feeds the weathering outflux after several hundreds of thousand years.
Then, the weathering flux follows the same evolution as with the humped law.

470
Finally, in order to test the influence of mountain size and model resolution, we reduce the pixel size to 20 m ( Figure 12). The widening parameter α is multiplied by 5 in order to have a valley width larger than one pixel. The same transfer from the weathering reservoir on the hillslopes to the valley borders is observed. This suggests that this main outcome does not depend on the system size.

Discussion
Cidre does not model the precipitation of secondary minerals, or variations in the pH, pCO2 and changes in the chemical equilibrium related to the water-rock interaction (Oelkers et al., 1994;Brantley et al., 2008;Lebedeva et al., 2010), so that the predicted fluxes may be overestimated. Furthermore, the groundwater circulation is also neglected, although it can contribute significantly 480 to the weathering outflux (Calmels et al., 2011;Maher, 2011).
Nevertheless, our modelling approach presents several advantages: the model is at the scale of the whole landscape but also at the pedon scale (a denser clast distribution can be set in specific areas of interest); any mineralogical assemblage and clast size distribution can be studied; there is no need 485 to calculate the mountain weathering outflux W at each time step. If only a long-term trend of W is studied, it can be calculated at a low frequency, which is computationally efficient and allows this 3D approach to be applied to long time periods. Most importantly, 1-weathered material can be followed from the source to the sink and 2-the weathering outflux results from the stochastic residence times of the clasts in the hillslopes and in rivers. These two last points constitute the main differ-490 ences of our approach compared to previous pedon and landscape weathering models (e.g. Ferrier and Kirchner, 2008;Vanwalleghem et al., 2013;Braun et al., 2016).
All previous models founded on clast residence time in the regolith predict that weathering should be zero when the regolith disappears (e.g. Ferrier and Kirchner, 2008). Yet, documented catchment 495 weathering rates are significantly larger than zero (Dixon and von Blanckenburg, 2012). A simple explanation may be that there is always a sufficient fraction of hillslopes covered by soils to produce a significant weathering flux, even in fast-eroding mountains (Larsen et al., 2014;Heimsath et al., 2012). Alternatively, deep weathering within fractured bedrock may account for this difference. Calmels et al. (2011) showed that this deep weathering reservoir accounts for more than 1/3 of 500 the silicate weathering flux of a catchment in Taiwan. The significant contribution of deep groundwater weathering echoes the model proposed by Maher (2010). In this model, this is the ratio between the fluid residence time and the characteristic mineral dissolution time (our τ m ) which controls the weathering flux, rather than the ratio N clast between the clast's residence time and τ m . West (2012) argued also for a monotonic increase of chemical weathering with erosion due to water circulation in 505 the fractured bedrock, so that the weathering would not critically depend on the regolith thickness.
If alternatively the weathering layer corresponds to the vadose zone, then this layer may thicken and sustain the weathering flux in rapidly eroding hillslopes (Rempe and Dietrich, 2014;Maher and Chamberlain, 2014;Braun et al., 2016).

510
Our model points to another possible reservoir: the colluvium temporarily stored along valley borders. When the regolith thins, colluvium become the main locus of weathering, which prevents the weathering outflux of the catchment from dropping to very low values. This finding is supported by the correlation between the weathering rate and the volume of the landslides present in catchments in southwestern New Zealand (Emberson et al., 2016a). In another catchment in New Zealand, Em-  (Figure 6). These durations are consistent with residence times in the Andes for example, as determined from U series (Dosseto et al., 2008).
Thus, grains stay in the catchments long enough to yield a significant weathering flux. Nevertheless, our model does not account for the full stochasticity of landslides documented by Emberson et al. (2016a). In our model, colluvium are produced relatively continuously, so that differences in clast residence times are mainly due to progressive colluvium removal and differences in the initial distance from the outlet. In real landscapes, even those without significant lateral river erosion, there will still be colluvium storage on the hillslopes because the sediment production is stochastic. Thus, a better description should include the stochastic production of landslides (Gabet, 2007).  (Maher, 2011;Emberson et al., 2016b).
Despite its limitations, our model predicts weathering-erosion rates within the range of existing data (Figures 3, 8) and may explain this range in terms of topographic evolution. In cooling experiments, our model predicts weathering rates that initially increase with time due to supply-limited 545 conditions and increasing erosion, but then decline because the regolith is progressively stripped off.
This hump evolution is consistent with previous 1D models ( Gabet and Mudd, 2009;Dixon and von Blanckenburg, 2012). This is a remarkable similarity given the heterogeneity of the regolith thickness and denudation rate in the simulated landscapes during the relief growth. In our modellings, the hump in the weathering evolution results from the progressive stripping of the regolith via an in-550 creasing erosion rate and cooling climate, whatever the form of the regolith production law. The peak occurs for some optimum compromise between the area covered by regolith and its thickness distribution (Carretier et al., 2014). When accounting for colluvium, the weathering peak is followed by a nearly constant weathering flux in agreement with models assuming a constant weathering layer (West, 2012) or based on the residence time of the groundwater (Maher and Chamberlain, 2014).

555
However in our model, this sustained weathering rate does not result from a constant weathering layer but rather from a change in the weathering reservoir.
The contribution of colluvium to the weathering flux should depend on the ratio between river width and valley width, which itself depends on the lithology, uplift rate or flood distribution (e.g. 560 Brocard and der Beek, 2006). Thus, the colluvium reservoir cannot be considered as a general model for all catchments. Nevertheless, colluvium contribute significantly in all our cooling simulations, regardless of the rainfall pattern or catchment size. This suggests that fresh sediment that is temporarily stored along valleys, irrespective of the cause of their width (uplift variations, glaciations, etc.), may contribute to the long-term weathering fluxes trend. In particular, the weathering of these sediment 565 may be only slightly dependent on climate variations. An increase in water discharge fosters mineral weathering but at the same time decreases mineral residence times, so that the net weathering varia-17 Earth Surf. Dynam. Discuss., https://doi.org/10.5194/esurf-2017-48 Manuscript under review for journal Earth Surf. Dynam. Discussion started: 28 August 2017 c Author(s) 2017. CC BY 4.0 License. tion is negligible (Figure 6). The same balance may operate in foreland basins and may take part in the weathering stability over the last 12 Ma observed in the offshore sediment record (Willenbring and von Blanckenburg, 2010;von Blanckenburg et al., 2015).

570
In cooling experiments (Figure 7), the short weathering time of the grains implies that they are only partially weathered when they leave the mountain, so that their weathering in adjacent basins could also contribute to the total weathering outflux. Few data confirm this behaviour Bouchez et al., 2012;Moquet et al., 2016). The contribution of the basin still needs to be 575 analysed through the stochastic dynamics of the grains in the alluvial plain, a study within the scope of our model.

Conclusion
We designed a new model at the landscape scale which takes the weathering of distinct clasts in the regolith, colluvium and rivers into account. This model accounts for part of the stochasticity 580 of sediment transport, which is reflected by the distribution of the clast residence times in an uplifting mountain. The weathering model has limitations (no groundwater model, no dependence on PH or pCO2 or water-rock chemical disequilibrium, no precipitation of secondary phases). Nevertheless the model predicts a range of weathering rates consistent with the available data for a wide range of climatic and tectonic contexts. During the rising of a mountain and as the climate cools, 585 the weathering flux increases and then decreases, which is consistent with previous models. In addition, the dynamic adjustment of the topography, the tracing of weathered material and the stochastic transport of grains point to a possible significant contribution from colluvial deposits during cold periods. This weathering reservoir may contribute to a high and constant weathering flux in rapidly eroding mountains under cold conditions, in addition to deep weathering in fractured bedrock and 590 other potential reservoirs. This new model opens perspectives to study the weathering contribution of foreland basins during mountain growth and decline and the response of these reservoirs to cyclic climatic variations.
Acknowledgment. This paper is a contribution to the LMI COPEDIM (funding from IRD). M.R.     (Figures 5a and c). The clast residence time is defined as the time during which a clast weathers (e.g. Dosseto et al., 2006;Mudd and Yoo, 2010), either in the regolith on the hillslopes, or in the colluvium. Note that the distribution is cut at 5 ka, but residence times of several 10 ka are present at 1.5 Ma, although they constitute much less than 1% of the clasts. After 6 Ma, the regolith has completely disappeared, so that the residence times represent the time spent after their detachment from the bedrock, i.e. primarily within the colluvium and possibly sediment temporarily deposited by the rivers. This distribution of the residence times gives an estimate of the colluvium residence times along the valleys. 99% of the residence times are smaller than 1.5 ka after 4 Ma. This distribution increases slightly after 6 Ma because the rainfall decreases. The increase in residence time compensates for the rainfall decrease to produce the nearly constant weathering outflux seen in  Figure 10. Same as Figure 9 (cooling ending with a regolith at elevations < 1500 m), but the precipitations vary with elevation according to a gaussian curve with an elevation peak at 1300 m (orographic model 1 illustrated in Figure 8). The erodibility parameters (Equation 2) are also doubled to compensate for the smaller precipitation rate at high elevations and thus, to obtain comparable maximum elevations.