Interactive comment on “ Effect of changing vegetation on denudation ( part 2 ) : Landscape response to transient climate and vegetation cover

I am happy that the authors should proceed on the basis of the comments and reviews received so far. The reviews are admittedly mixed, but on the whole the reviewers feel that there is merit in this work. When addressing the the points that the reviewers identify, the authors should consider how they present and justify the model parameter selection, and whether any model calibration is possible or appropriate. There are also some suggestions for how the model results could be broadly compared with examples from the literature. As identified by a couple of reviewers, a conceptual model might help to summarise the model findings. Finally, I would also encourage the authors to


Model Description and governing equations
For this study, we use the open-source model framework Landlab (Hobley et al., 2017), which provides easily accessible methods for building a landscape evolution model.Landlab provides the computational environment to b uild an experimental set-up to test hypotheses and conduct sensitivity analyses of topography to different surface process parameterizations.
For our study, we chose a model-domain representative of each of the four Chilean areas using an area of 100km 2 which is implemented as a rectangular grid, divided into 0.01km 2 spaced grid cells.For simplification in the presentation of results, we present our results for the driest, northern most area (Parque Nacional Pan de Azucar) and for the Parque Nacional La Campana which is situated at 32°S Latitude (Fig. 1) and shows the highest values in analyzed basin metrics (Fig. 2), although the general behavior and results presented here are representative of the other two areas not shown.The topographic evolution of the landscape is a result of tectonic uplift and surface processe s, incorporating detachment limited fluvial erosion and diffusive transport of sediment across hillslopes (Fig. 3).These processes are linked to, and vary in their effectiveness due to surface vegetation density.Details of the implementation of these proc esses into Landlab are explained in the following subsections.This model setup is simplified in regards to hydrological parameters like, for example, soil moisture and groundwater and unsaturated zone flow.Also, the erosion and transport of material due to mass-wasting processes such as rockfalls and landslides are not considered.We argue that those processes do not play a major role in the basins we used for modelcalibration and that the processes acting continuously along hillslopes and channels have the largest impact on shaping our reference landscapes.Additional caveats and limitations of the modeling approach used are discussed in Section 5.4.
Main model parameters used in the model (and described below) are provided in Table 1.

Boundary and Initial Conditions, and Model Free Parameters
In an effort to keep simulations for the different EarthShape areas comparable, we minimized the differences in parameters between simulations.The exceptions to this include the surface vegetation cover and mean annual precipitation, which were varied between simulations.One of the main controls on topography is the rock uplift rate.We kept the rock uplift rate temporally and spatially uniform across the domain and at 0.2 mm/yr (Table 1).Studies of the exh umation and rock uplift history of the Coastal Cordillera, Chile, are sparse at the latitudes investigated here, but existing and in progress studies further to the north are broadly consistent with the rock uplift rate used here (Juez-Larré et al., 2010;Avdievitch et al., 2017).Furthermore, the thermochronometer cooling ages in the northern Coastal Cordillera suggest constant Cenozoic exhumation over >50 Ma at this rate.Thus, despite being located on an active plate boundary, existing observations suggest relatively slow, and temporally constant rock uplift of this region.
The initial topography used in our simulations was a random white-noise topography with <1 m relief.To avoid unwanted transients related to the formation of this initial topography we conduct equilibrium simulations for each set of the different climate and vegetation scenarios (see below), and run the model for 15 Myr until a topographic steady -state is reached.The equilibrium topography after 15 Myr was used as the input topography for subsequent experiments that impose transient forcings in climate, vegetation, or both (Fig. 4).The model simulation time shown in subsequent plots is the time since completion of this initial 15 Myr steady-state topography development.In the results section, we present these results starting with differences in the initial steady -state topographies (prior to imposing transient forcings) and Earth Surf.Dynam.Discuss., https://doi.org/10.5194/esurf-2018-13Manuscript under review for journal Earth Surf.Dynam.Discussion started: 21 February 2018 c Author(s) 2018.CC BY 4.0 License.
then add different levels of complexity by imposing either: (1) a single transient step -change for the vegetation cover (Fig. 4b); (2) a step change in the mean annual precipitation (Fig. 4d); (3) 100 kyr oscillations in the vegetation cover (Fig. 4a); (4) 100 kyr oscillations in the mean annual precipitation (Fig. 4c); or (5) 100 kyr oscillation in both the vegetation cover and mean annual precipitation (both Fig. 4a and 4c).This approach was used to produce a stepwise increase in model complexity for evaluating the individual, and then combined, effects of fluvial and hillslope processes to different forcings.
The magnitude of induced rainfall transient forcings where based upon the present -day conditions along the Coastal Cordillera study areas (Fig. 1b, c).The step change and oscillations in vegetation cover and mean annual precipitation imposed on the experiments were designed to investigate vegetation and precipitation change effects on topography over the last ~0.9Ma, the period during which a 100 kyr orbital forcing is dominant in Earth's climate.Given this timescale of interest, we impose a 10% magnitude change in the step-increase or decrease, or the amplitude change in oscillations for the vegetation cover.This magnitude of vegetation cover change is supported by dynamic vegetation modeling of vegetation changes over glacial-interglacial cycles in Chile (see companion paper by Werner et al. 2018-this journal) and to some degree elsewhere in the world (Allen et al., 2010;Prentice et al., 2011;Huntley et al., 2013).We assume that the present-day conditions of combined vegetation cover and mean annual precipitation along the north-south gradient of the coastal cordillera are directly linked (Fig. 1b, c), and therefore follow an empirical approach based on the present day mean annual precipitation vs. vegetation cover relationship in Chile (Fig. 5).We do this by associating each 10% change in vegetation cover (dV) with a corresponding change of mean annual precipitation (dP, Fig. 5) present in the study areas considered.
The boundary conditions used in the model were the same for all simulations explained above (Fig. 3).One boundary was held at a fixed elevation and open to flow outside the domain.The other three were allowed to increase in elevation and had a zero-flux condition.This design for boundary conditions is similar to previous landscape evolut ion modelin g studies (Istanbulluoglu and Bras, 2005) and provides a means for analyzing the effects of different vegetation cover and precipitation forcings on the individual catchment and subcatchment scale.

Vegetation Cover Dependent Geomorphic Transport Laws
The governing equation used for simulating topographic change in our experiments follows the continuity of mass.
Changes in elevation at different points of the model domain over time dz(x,y,t) depend on where z is elevation, x, y are lateral distance, t is time, U is the rock uplift rate , kd∇z the linear diffusive flux of sediment along hillslopes, and Dc the detachment capacity of channels (Tucker et al., 2001a).Our implementation of vegetation cover effects for the last two parameters in the Landlab model are described in more detail below.

Vegetation Cover Influenced Diffusive Hillslope Transport
The change of topographic elevation in a landscape over time which is directly caused by hillslope-dependent diffusion can be characterized as: Landscape evolution models characterize the flux of sediment q s either as a linear or non-linear function of surface slope S (Culling, 1960;Fernandez and Dietrich, 1997).In order to keep the number of free parameters for the simulation to a minimum, we used the linear description of hillslope diffusion: Following the approach of (Istanbulluoglu and Bras, 2005), we assign the linear diffusion coefficient Kd as a function of surface vegetation density V, an exponential coefficient α, and a baseline diffusivity Kb, such that:

Vegetation Cover Influence on Overland Flow and Fluvial Erosion
The erosion and transport of material due to river incision is represented as a detachment -limited process where the amount of lowering of elevation (z) over time depends on the fluvial erodibility Kv, [L 1-2m t -1 ] (a function of soil type, climate, vegetation and the scaling of runoff to drainage area).The contributing drainage area to each node is represented by A [L 2 ], S is the corresponding channel slope, and m, and n are empirically derived constants (Howard and Kerby, 1983;Whipple and Tucker, 1999).The detachment-limited erosion law can be re-written in a form that directly links excess bed shear stress to detachment capacity of a river where: where ke represents the erodibility of the substrate, τb represents the bed shear stress acting on the bed surface, τc is the critical shear stress needed to erode the substrate, ρw is the density of water, g is gravitational acceleration, R is hydraulic radius and S is local channel slope.By combining the shear-stress formulation with Mannings equation and introducing two Mannings surface roughness values (nV for the influence of vegetation on surface roughness and n s for the influence of the soil surface to surface roughness, e.g., Istanbulluoglu and Bras, 2005) we can reformulate the hydraulic radius R (Wilgoose et al., 1991;Istanbulluoglu et al., 2004) and write the shear-stress formulation as a function of the Mannings coeffcients: By writing n v as a function of a basic Mannings coefficient for a reference vegetation cover n vr and V r , the vegetation cover at each individual node V and an empirical scaling parameter W, we arrive at: Combining equation 8 and 9 results in an equation for vegetation cover dependent shear stress at each node.The effects of surface vegetation cover on both diffusivity and fluvial erodibility are shown in Figure 6.

Model Evaluati on
Model performance was evaluated using the above equations and different initial vegetation covers and mean annual precipitation based on the steady-state predicted topography.Our focus in this study is on the general surface process response to different transient vegetation and climate conditions, rather than a calibrated modeling study of the Chilean study areas.Nevertheless, topographic metrics of relief, mean slope, and normalized steepness index (Ksn) were computed from the model results and compared to observed values from the 30 m SRTM DEM for each of the four areas (Fig. 2).

Results
Our presentation of results is s tructured around three groups of simulations.after the landscape has reached steady state (Section 4.3).For each g roup of transient simulations, we show the topographic evolution with help of standard topographic metrics and the corresponding erosion rates after the induced change.

Equilibrium Topographic Metrics
Topographic metrics from each of the four Chilean focus areas (Fig. 1a) were extracted for comparison to equilibriu m topographies predicted after 15 Myr of model simulation time.This comparison was done to document the model response to changing vegetation cover (with climate held constant) and changing vegetation cover and precipitation, and also to demonstrate the modeling approach employed throughout the rest of this study captures the general characteristics of different topographic metrics along the Chilean Coastal Cordillera.We refrain from conduc ting a more detailed modelobservation comparison for reasons previously mentioned.
Analysis of the digital elevation model for each of our four Chilean focus areas illustrates observed changes in catchment relief, slope, and channel steepness (Ksn) in relation to the surface vegetation (Fig. 7, red points) and latitude (Fig. 2).The general trend in the observed metrics shows a non-linear increase in each metric until a maximu m is reached for regions with 70% vegetation cover.Following this, all observed metrics show a decline towards the area with 100% vegetation cover.
The model predicted equilibrium topographies (Fig. 7a,b,c) from four different steady-state simulations with variable vegetation cover but a constant mean annual precipitation (900 mm/yr) show a nearly linear increase in all observed basin metrics with increasing vegetation cover and therefore do not reflect the overall trend observed in the DEM from the study areas.For example, basin relief and slope are both under predicted for simulations with V < 100% (Fig. 7a,b), and only the predicted maximu m relief for a fully-vegetated simulation resembles the DEM maximu m value.For the normalized channel steepness, only two observed mean values (for V = 10% and 70%) lie within the range of mean to maximu m predicted Ksn values (Fig. 7c).
The resulting equilibrium topographies from simulations with both variable mean annual precipitation an d vegetation cover (Fig. 7d,e,f) show an improved representation of the general trend of the DEM data.The vegetat ion cover and precipitation values used in these simulations come from the values observed in the Chilean areas (Fig. 1b, c; Fig. 5).In these simulations, the maximu m in the observed basin metrics is situated at values of V = 30% with a following slight decrease in the metric for V = 30% to V = 70%, followed by a steeper decrease in metrics from V = 70% to V = 100%.
Generally the model-based results tend to underestimate the basin relief and overestimate the bas in channel steepness (Fig. 7d,f).Variations in basin slope are captured for all but the non-vegetated state (Fig. 7e).
Although the above comparison between the models and observations demonstrates a range of misfits between the two, there are several key points worth noting.First, the model results shown are highly simplified in their setup (e.g.assuming similar rock uplift rate, identical lithology and constants), and assume the present day topography is in steady state for the comparison.Second, despite the previous simplifying assumptions, the degree of misfit between the observations and model are surprisingly small when both variable vegetation and variable precipitatio n, are considered (Fig. 7d,e,f).Finally (third), the general 'humped' shape curve observed in the Chilean areas is captured in the model predictions (Fig. 7d,e,f), with the notable exception that the maximu m in observed values occurs at a higher vegetation cover (V = 70%) than the model predictions (V = 30%).Explanations for the possible source of these differences are revisited in the discussion section.Without additional observations from the Chilean areas, reduction of the misfit between the observations and models is not tractable.

Transient Topography -Step Change
The evolution of our model topographies after a ind uced instantaneous disturbance (Fig. 4) of either only the surface vegetation cover (Fig. 8, green lines) or a step change in only the mean annual precipitation (Fig. 8, blue lines) is analyze d for changes in topographic metrics for either a pos itive disturbance (Fig. 8a,b,c) or a negative disturbance (Fig. 8d,e,f).
This scenario with changes in only vegetation or only precipitation was chosen to analyze sensitivity of drainage basins to changes in vegetation cover and precipitation and isolate the effects o f these specific transient forcings.Decoupling of vegetation and climate changes can occur via sudden disturbances, such as wildfires and lagged vegetation responses to climatic changes.Mean catchment erosion rates are also analyzed for their evolution a fter the disturbance (Fig. 9).For simplicity in presentation, results are shown for only two of the initial vegetation (V) and precipitation (P) values for vegetation covers of 10 and 70%, and precipitation rates that correspond to these vegetation covers (i.e.P(V=10%) or P(V=70%)) in the Chilean areas (Fig. 5).The results described below show a general positive correlation between all observed topographic metrics and surface vegetation cover and a negative correlation between observed topographic metrics and mean annual precipitation and therefore supports the data from our equilibrium topography simulations.The adjustment time until the system again reaches a new steady state varies between different simulations.

Positive Step Change in Vegetation Cover or Precipitation Topographic Analysis
A positive step change in vegetation cover (V) from V = 10% to V = 20% (solid green line Fig. 8a,b,c) leads to an increase of mean basin relief from 270m to 520m, mean basin slope from 11.2° to 15.9°, and mean basin channel steepness from 108m -0.9 to 222m -0.9 , which corresponds to a factor 1.9, 1.42, and 2.1 change, respectively.The adjustment time until a new steady state in each metric is reached is 3.1Ma.The corresponding positive change in mean annual prec ipitation (solid blue lines, Fig. 8a,b,c) leads to a decrease of mean basin relief to 176m, mean basin slope to 8.6° and mean basin channel steepness to 67m -0.9 .This corresponds to a decrease by factors of 1.5, 1.2 and 1.6, respectively.The adjustment time to new steady state conditions in this case are shorter and 1.1Ma (Fig. 8a,b,c).A second feature of these results is the brief increase and then decrease in basin average slope angles following the step change (Fig. 8b).
For simulations with V = 70% initial surface vegetation cover, a positive increase to V = 80% leads to an increase of mean basin relief from 418m to 474m, mean basin slope from 15.5° to 16.8° and mean basin channel steepness from 172m -0.9 to 199m -0.9 .This causes an increase in each metric by factors of 1.1, 1.1 and 1.2, respectively.The adjustment time to steady-state conditions is 1.9Ma (dotted green lines, Fig 8a,b,c).The corresponding positive change in mean annual precipitation leads to a decrease of relief to 268m, decrease in s lope to 11.9° and decrease of channel steepness to 105m -0.9 .This resembles a decrease by factors 1.5, 1.3, 1.6, respectively.Adjustment time in this case is 1.7Ma (dotted blue lines, Fig. 8a,b,c).The basin slope data shows similar behavior as the Vini = 10% simulations with an initial decrease and then increase for a vegetation cover step change and an initial increase and then decrease for a step change in mean annual precipitation.Comparison of the change in the topographic metrics for the low (V=10%) and high (V=70%) initial vegetation covers, the magnitude of change in each metric is larger when the step change occurs on a low, rather than higher, initial vegetation cover topography.

Erosion Rate Changes
Erosion rates show instantaneous reactions to positive disturbances in vegetation or precipitation.Generally, the model results show a negative relationship between increases in vegetation cover and erosion rates and a positive relationship between increases in precipitation and erosion rates.Altho ugh the reaction between the disturbances and changes in erosion rates are instantaneous, the specific maximu m or minimum is reached after some lag time and the magnitude and duration of non-equilibrium erosion rates varies between different simulation setups.
For initial vegetation cover of V = 10%, a change in vegetation cover (dV) of +10% leads to a decrease in erosion rates from 0.2 to 0.03mm/yr (factor of 5.7 decrease, Fig. 9a green line).The minimum erosion rate is reached 43.5kyrs after the step change occurs.Following this minimum in erosion rates, the rates increase until the steady -state erosion rate is reached after the adjustment time.An increase in mean annual precipitation corresponding to a vegetation cover of 10% (i.e.P(V=10%) to P(V=20%); Fig. 5) leads to an increase in erosion rates to a maximu m of 0.44mm/yr after 74.8kyrs (factor of 2.2 increase, Fig. 9a, blue line).For initial vegetation of V = 70% a vegetation increase of dV = +10% results in minimum erosion rates of 0.14mm/yr after 117.7kyrs (factor of 1.4 decrease, Fig. 9b, green line).A corresponding increase in precipitation for these same vegetation conditions leads to maximu m erosion rates of 0.44mm/yr after 107.5kyrs which is an increase by a factor of 2.2 (Fig. 9b, blue line).The previous results for a positive step change in vegetation or precipitation demonstrate that the magnitude of change in erosion rates is large for changes in precipitation rate than for vegetation cover changes, and in low initial vegetation cover sett ings (V=10%) that magnitude of change in erosion rates for changing vegetation is larger (compare green lines Fig. 9a with 9b).

Negative Step Change in Vegetation Cover or Precipitation Topographic Analysis
For negative step-changes in vegetation (green curves, Fig. 8d,e,f), the results show a sharp decrease in topographic metrics associated with shorter adjustment times compared to the positive step chan ge experiments (compare Fig. 8d,e,f with a,b,c).For step changes in precipitation (blue curves, Fig. 8d,e,f), the increase of topographic metrics happens slower and therefore with longer adjustment times.A negative step change in vegetation cover from V = 10% by dV = -10% leads to a decrease of mean basin relief from 269m to 35m, mean basin slope from 11.2° to 2.3° and mean basin channel steepness from 108m -0.9 to 11m -0.9 which resembles decreases by factors of 7.8, 3.8 and 9.6.leads to an increase in mean basin relief to 512m, mean basin slope to 15.8° and mean basin channel steepness to 223m - 0.9 .These increases reflect changes by factors of 1.9, 1.4 and 2.1 with an adjustme nt time of 4.9Ma (dotted green lines, Fig. 8d,e,f).Mean basin slope results (Fig. 8e) for a step change in vegetation illustrate a pulse-like feature of initially increasing slope values, followed by a decrease to lower slope values.In contrast, a negativ e step change in precipitation induce an initial decrease in slope, followed by a gradual increase in slope to a value higher than was initially observed before the change.
Simulations with initial vegetation cover V = 70% and dV = -10% show a decrease in mean basin relief from 418m to 356m, mean basin slope from 15.5° to 13.6° and mean basin channel steepness from 172m -0.9 to 144m -0.9 which resembles changes by factors of 1.2, 1.1 and 1.2 and an adjustment time of 2.1Ma (dott ed green lines, Fig. 8d,e,f).Corresponding negative changes in precipitation lead to increase of basin relief of 465m, basin slope to 16.4° and channel steepness to 195m -0.9 which resembles changes by factors of 1.1 for all three values.Adjustment time in this case is 2.2Ma (dotted blue lines, Fig. 8d,e,f).Behavior of mean basin slope after the step-change follows the V = 10% simulations but shows lower amplitudes of basin slope for both step-changes in vegetation and precipitation.

Erosion Rates
The positive step-change results (Fig. 9a, b) indicated that erosion rates reach their minimum or maximu m after a lag time after the change, and show significant differences in the magnitude and duration of non -equilibrium conditions depending on if vegetation or precipitation were changing.Simulations with a decrease from V = 10% to V = 0% (Fig. 9c) show a sudden increase in erosion rates to a maximu m value of 3.5mm/yr which is an increase from steady state conditions by a factor of 17.7 which is reached after 19.5kyrs (green line, Fig. 9c).A step decrease in precipitation for this corresponding vegetation difference (i.e.P(V=10%) to P(V=0%) leads to a smaller, and protracted (longer adjustment time) decrease in erosion rates to 0.03mm/yr after 50.1kyrs.These conditions cause a factor of 5.6 decrease (blue line, Fig. 9c).
Simulations of V = 70% with a vegetation change of dV = -10% show an increase in erosion rates to 0.27mm/yr which is a factor of 1.4 increase after 126.3kyrs (Fig. 9d).For the corresponding decrease in precipitation the data show a decrease in erosion rates to 0.15mm/yr after 124.5kyrs.This resembles change by factor of 1.2 (blue line, Fig. 9d).

Transient Topography -Oscillating
In addition to simulations where a transient step-change in either surface vegetation density or mean annual precipitation was conducted, we set up two distinct sets of simulations with an oscillating transient signal with a period of 100kyrs.
This period resembles the eccentricity driven part of the Milankovitch cycle that is dominant in Earth's climate over the last 0.9Ma.

Oscillating Surface Vegetation Cover, Constant Precipitation Topographic Analysis
The topographic evolution in simulations with a constant precipitation (10 and 360mm/yr for V=10%, and V=70%, respectively) and oscillating vegetation cover show a different response than the previous step change experiments.The differences depend on the initial steady-state vegetation cover prior to the onset of 100kyr oscillations.All observed basin metrics (Fig. 10) show an initial oscillating decrease in values until a new dynamic steady-state is reached where the amplitude in oscillations is less than in the preceding initial adjustment period.Simulations with V = 10% (Fig. 10a) show a decline in the basin relief from 269m to 107m which resembles a decrease of the mean elevation of factor 2.5.The positive amplitude of oscillation is 9m, the negative amplitude is 8.3m but time intervals of negative amplitudes are longer shows a similar behavior with a decrease of the mean slope from 11.2° (prior to the onset of oscillation s) to 6.0° (Factor 1.6, Fig. 10b).However, before this new equilibrium is reached, the slopes show an increase in mean slope for the first two periods of vegetation oscillation which then declines towards the new long -term stable dynamic equilibrium which is reached after approximately 500kyrs.Local maxima of mean basin slope coincide with local minima in basin relief.
For the V = 70% simulations, the reaction is significantly smaller with no change in mean slope for the new dynamic equilibrium and amplitudes of both positive and negative of 0.16°.Mean b asin channel steepness (Fig. 10c) reflects the behavior of mean basin elevation.For V = 10% simulations the mean channel steepness decreases from values of 108m - 0.9 to 40m -0.9 (factor 2.7 change) with a positive amplitude of 3.7m -0.9 and a negative amplitude of 6.1m -0.9 .For V = 70% simulations the response is again only minor, compared to the lower initial vegetation cover simulations with a change of mean channel steepness from 186m -0.9 to 167m -0.9 and positive amplitudes of 1.1m -0.9 and negative amplitudes of 0.9m - 0.9 .Like the elevation data, the steepness data shows a distinct oscillating pattern with slow increases to local maxima and rapid decreases to local minima which coincide with maxima/ min ima of elevation data.Taken together, the previous observations demonstrate a larger change in topography for oscillations in poorly vegetated areas compared to those with higher vegetation cover.Furthermore, the magnitude of topographic change that oscillations in vegetation impose on topography are largest in the first ~500 kyr after the onset of an oscillation, and diminis h thereafter.

Erosion Rates
The erosion history for simulations with oscillating vegetation cover (Fig. 11) demonstrate large variations in the erosion rate that depend on the average vegetation cover of the oscillation.Furthermore, pronounced differences in the amplitude of erosion occur if the vegetation cover is above or below the mean of the oscillation (Fig. 4a).More specifically , simulations with V = 10% show a pattern of a small decrease in erosion rates (from 0.2 to 0.03mm/yr) when vegetation cover increases above the mean cover, in contrast to a large increase in erosion rates up to 3.3mm/yr when vegetation cover decreases below the mean of the oscillation (Fig. 2, Fig. 11).Maximum erosion rates decline over multiple periods of oscillation until they reach a dynamic steady-state with maximu m rates of 1.2mm/yr at 760kyrs after the onset.Time periods of higher erosion rates (>0.2mm/yr) have a mean duration of 28kyrs, whereas periods of lower erosion rates (<0.2mm/yr) have a mean duration 72kyrs.For simulations with high vegetation cover (V = 70%) the maximu m and minimum erosion rates are 0.28 and 0.15mm/yr, respectively.The magnitude of maximu m and minimum erosion rate are not significantly time-dependent and are reached at each local vegetation cover minimum.The mean duration of periods with higher erosion rates (> 0.2mm/yr) is 55kyrs whereas the duration for periods with lower rates (< 0.2mm/yr) is 45kyrs .
These results demonstrate that areas with low vegetation cover experience not only larger amplitudes of change in erosion rates, but also an asymmetric change whereby decreases in erosion rates are lower magnitude than the increases in erosion rates.

Oscillating precipitation, Constant Vegetation Topographic Analysis
The evolution of topographic parameters for simulations with oscillating mean annual precipitation and constant surface vegetation cover (V=10 or 70%, Fig. 12) show a less extreme and smaller temporal change in erosion rate to variations in precipitation compared to the previous discussed effects of oscillating vegetat ion cover (Fig. 11).In Fig. 12a, the mean basin relief results for V = 10% and oscillating precipitation show small variations (+4.9m to -3.8m) in relief around a mean of 269m, which is similar to the mean relief prior to the onset of oscillations at 5000kyrs.For simulations with V = 70% the change in relief is slightly more pronounced with adjustment to a new mean of 380m from 418m in steady -state conditions.This change in mean relief equates to factor of 1.1 change.The evolution of topographic slope (Fig. 12b) for V = 10% simulations shows an adjustment to a new dynamic equilibrium from 11.2° to 10.6° (Factor 1.05) with a negative amplitude of 1.2° and positive amplitude of 0.9°.For V = 70% the mean slope values do not significantly change fro m steady-state to transient conditions and the amplitudes of oscillation are 0.6° (positive) and 0.7° (negative).Mean channel steepness (Fig. 12c) for V=10% shows an adjustment from 108m -0.9 to 110m -0.9 (Factor 1.01).The amplitude of oscillation is 4 m -0.9 for both negative and positive amplitudes.For Vini = 70% simulations the topography adjusts from 171m -0.9 to 152m -0.9 (Factor 1.1) with amplitudes of 4.5m -0.9 for both positive and negative changes.Thus, although Figure 12 illustrates changes in topographic metrics that result from oscillations in precipitation occurring around vegetation covers of 10 and 70%, these changes are significantly smaller than those predicted for constant precipitation, but oscillating vegetation conditions (Fig. 10).

Erosion Rates
Predicted erosion rates from simulations with constant surface vegetation cover and oscillating mean annual precipitation indicate different amplitudes of change around the mean erosion rate depending on the vegetation cover.For simulations with V = 10% (Fig. 13, blue solid line) erosion rates oscillate symmetrically around the steady -state erosion rate of 0.2mm/yr.The maximu m and minimum erosion rates of 0.42 and 0.01mm/yr, respectively, result in no change in the mean rate.In contrast, predicted rates with a higher vegetation cover of V = 70% (Fig. 13, blue dotted line) demonstrate an asymmetric oscillation in rates around the mean, whereby the maximu m in rates (0.43mm/yr) has a larger difference above the mean rate, then do the minimums in the oscillation (0.12mm/yr).For this higher vegetation cover scenario, a gradual increase in the mean erosion rate from 0.2 to 0.25 mm/yr as time progresses is evident.Furthermore, the maximu m and minimum erosion rates decline over several oscillation periods to v alues of 0.38 and 0.15mm/yr, respectively.Taken together, these results indicate that oscillations in precipitation impact erosion with different magnitude depending on the amount of vegetation cover.Areas with low vegetation cover demonstrate the highes t and symmetric oscillation of erosion rates due to changes in precipitation whereas in areas with high vegetation cover the effect of negative changes in precipitation is dampened by vegetation.

Discussion
The previous results highlight predicted topographies with different sensitivities to changes in either surface vegetation cover or mean annual precipitation.The previous simulations, were conducted to isolate the magnitude of effect each parameter has on topography and erosion.In the following, we synthesize the previous results and then build upon them to discuss the, over longer time-scales, more common scenario of synchronous variation in both precipitation and vegetation cover.

Interpretation of Steady-State Simulations
Landscapes in a topographic steady-state show distinct features in topographic metrics that are widely used to estimate catchment-averaged erosion rates and therefore the leading processes of erosion within a landscape (DiBiase et al., 2010).
In most studies, focusing on comparing in-situ measured 10 Be erosion rates with topographic metrics, this is done in catchments with low variations in precipitation to focus on distinct topographic controls on soil erosion and transport processes.By conducting simulations with equal soil properties and assuming that the basic processes of sediment erosion and transport do not change between different climate-settings we can reproduce (Fig. 7) variations in topographic metrics over different climates seen in other studies (Langbein and Schumm, 1958;Walling and Webb, 1983) in steady-state landscapes with homogeneous erosion rates.Comparison of simulations with homogeneous precipitation and changing values of vegetation cover (Fig. 7a,b,c) to simulations with both changing precipitation and vegetation cover (Fig. 7d,e,f) we see that we are only able to reproduce a similar trend with a distinct peak in topographic metrics when both variable precipitation and vegetation cover are considered.From this, we conclude that modern model-based landscape evolution studies that aim to compare areas with different climates should incorporate vegetation dynamics in their simulations.
Misfits between the predicted and Chilean observed topographic metrics (Fig. 7d,e,f) present when the vegetation and precipitation both vary likely stem from the simplicity of the model setup used and the likelihood of differences of the rock uplift rate and lithology's present in these areas.

Interpretation of Step-Change Experiments
Our analysis shows that changes in vegetation-cover typically have a higher magnitude of impact on topographies for lower values of initial vegetation cover, compared to simulations with high initial vegetation cover (Fig. 8, 9).In those settings the influence of vegetation cover outweig hs the influence of precipitation in cases of negative and positive directions of the step change.The reasons for this is due to a higher impact of changes in vegetation on erosivity and diffusivity (parameter Kv, Kd; equation 4, 8) then changes in precipitation and therefore changes in runoff have on overall erosion values.
Furthermore, a negative step change in vegetation cover impacts the topographic metrics a factor of two more than do positive step change changes (Fig. 8d,e,f).This response is interpreted to be due to the non-linear reaction of diffusivity and fluvial erodibility to changes in vegetation cover (See Fig. 6).Negative changes in vegetation cover lead to a higher overall change in diffusivity and erodibility which leads to a higher sens itivity of equations 4 and 8 to negative stepchanges compared to positive step-changes.
Model results for the topographic metrics and erosion rates also indicate a difference in the adjustment times of the system until a new steady state is reached when either precipitation or vegetation cover changes (Figs. 8, 9).For simulations with topographic metrics or erosion rates is three times higher than the adjustment time for changes in precipitation.
Simulations with a negative step-changes in vegetation cover show an adjustment time which is lower by a factor of 18 compared to negative changes in precipitation.This difference in adjustment t ime again is a result of the non-linear behavior of erosion parameters Kd and Kv which influence how effective a signal of increasing or decreasing erosion can travel through a river basin (Perron et al., 2012).High values of Kd and Kv are associated with lower adjustment times and are a result of negative changes in vegetation cover.The influence of changing precipitation on adjustment time behaves in a more linear fashion and therefore mostly depends on the overall magnitude of change.Therefore, positive step-changes in vegetation cover decrease Kd and Kv which leads to higher adjustment times than the corresponding changes in precipitation.
An increase and then decrease, or decrease and then increase, in predicted slope and erosion rates is observed for both the positive and negative s tep changes experiments (Fig. 8b,e; and Fig. 9).This non-linear response in both positive and negative step changes in precipitation and vegetation cover is also manifested in the subsequent oscillation experiment s , but most clearly identifiable in the step change experiments.The explanation for this behavior is as follows.A positive step change in vegetation cover (Fig. 8b) leads to a decrease in fluvial capacity because increased vegetation cover increases the Manning's roughness (parameter nv, equation 8).The effect of changing the Manning's roughness varies with the location in the catchment and influences which processes (fluvial or hillslope) most strongly influence slopes and erosion rates.In the upper part of catchments where contributing areas (and discharge) are low, this increase in Manning's roughness causes many areas to be below threshold conditions such that fluvial erosion is less efficient, and hillslope diffusion increases in importance's and lowers s lopes.In the lower part of catchments, where contributing area and discharge are higher, changes in the Manning's roughness are not large enough to impact fluvial erosion because these areas remain at, or above, threshold conditions for erosion.With time, the lower regions of the catchments that are at or above threshold conditions propagate a wave of erosion up to the higher regions that are below threshold conditions.
The propagating wave of erosion eventually leads to increase in slope angles, essentia l due to the response time of the fluvial network to adjust to new Manning's roughness conditions.
In contrast, a positive step change in mean annual precipitation leads to an initial increase in fluvial shear stress which initially causes headward incision of river channels and leads to wave of erosion that propagates upstream and increas es channel slope values (Fig. 8b, see also e.g.Bonnet and Crave, 2003).The increase in channel slopes leads to an increase in the hillslope diffusive flux adjacent to the channels that then propagates upslope.Eventually, this increase in hillslope flux leads to a decrease in hillslope angles, and an overall reductions in mean catchment slopes after the systems reaches equilibrium.
Negative step-changes in vegetation cover or precipitation (Fig. 8e, green curves) shows the opposite behavior of the previous positive step change description.A negative step change in vegetation cover leads to an initial increase of fluvial erosion everywhere in the catchment because the Manning's roughness decreases everywhere.This catchment wide decrease in Manning's roughness leads to fluvial incision everywhere in the catchment and an increase in mean slope.
However, eventually hillslope processes catch up with increased slopes near the channels and with time an overall reduction of slope occurs.Negative changes in precipitation (Fig. 8e, blue curves) lead to an initial decrease in fluvial erosion which leads to an increase in the significance of hillslope processes such that slope ang les between channel and ridge decrease as hillslope processes fill in channels.With time, the fluvial network equilibrates to lower precipitation conditions by increasing slopes maintain equilibrium between erosion and rock uplift rates.Thus, the contrasting behavior of either initially increasing or decreasing slopes and erosion rates, followed by a change in the opposite direction of this initial change highlight a complicated vegetation -climate induced response to changes in either parameter.This non-linear behavior, and the timescales over which these changes occur, suggest that modernsystems that experienced past changes in climate and vegetation will likely be in a state of transience and the concept of a dynamic equilibrium in hillslope angles and erosion rates may be difficult to achieve in these natural systems.
Previous studies have inferred relationships between mean catchment erosion rates derived from cosmogenic radionuclides and topographic metrics (e.g., DiBiase et al., 2010;DiBiase and Whipple, 2011).However, the previous discussion of how topographic metrics change in response to variable precipitation and vegetation suggest that empirical relationships between erosion rates and topographic metrics contain a signal of climate and veget ation cover in the catchment.We illustrate the effect of step changes in climate and vegetation on the new steady -state of topographic metrics in Figure 14.In this example, the new steady state conditions in basin relief and mean slope after a modest (+/ -10%) change in vegetation or precipitation (triangles) differ from the initial steady -state condition (circles).These changes in topographic metrics when the new steady -state is achieved occur despite the rock uplift rate remainin g constant.Thus, differences in mean slope and relief can occur solely due to changes in climate or vegetation and are not necessarily linked to variations in erosion rate.The change in relief or slope is most pronounced for catchments with initially low (e.g.10%) vegetation cover.

Interpretation of Oscillation Experiments
The results from the 100kyr oscillating vegetation and precipitation conditions shows that oscillating vegetation cover without the corresponding oscillations in precipitation leads to adjustments of t opographic features, to a new dynamic equilibrium after approximately 1.5Ma (Figs. 10,11).The results indicate that the magnitude of adjustment depends on the initial vegetation cover, whereby simulations with 10% initial vegetation cover (solid lines, Fig. 10) show the largest changes from the initial (pre-oscillation) conditions to the new dynamic steady -state.Simulations with 70% initial vegetation cover (dashed lines, Fig. 10) show only minor adjustment to a new dynamic steady -state and lower amplitudes of oscillation.This is also represented by the mean basin erosion rates which show a significant peak for the first period of oscillating vegetation cover with erosion rates being 16 times higher than steady -state erosion rates for simulations with 10% initial vegetation cover whereas the peak erosion rate for 70% vegetation cover simulations is only higher by a factor of 1.4 in the first period of oscillation.The previously described response of topographic metrics and erosion rates to oscillating vegetation are due to processes described in the previous step change experiments.For example, the asymmetric oscillations in topographic metrics for V=10% (Fig. 10) are due to the superposition of positive, then negative changes described in section 5.2.Variations in the imposed Manning's roughness, and relative strengths of fluvial vs hillslope processes in different parts of the catchments at different times causes the topographic metrics and erosion rates to have a variable amplitude and shape of response from the symmetric oscillations imposed on the topography (Fig. 4a).
Simulations with oscillating precipitation and constant vegetation cover however show a less pronounced shift to new equilibrium conditions and in general lower amplitudes of oscillat ion in both topographic metrics and erosion rate (Figs. 12, 13).This difference in the response of the topographic metrics and erosion rates in Figures 12 and 13 diffusivity and fluvial erodibility to changes in vegetation cover compared to the linear response to changes in precipitation.

Coupled Oscillations in Both Vegetation and Precipitation
The previous sections present a sensitivity analysis of how step changes or oscillations in either vegetation cover or precipitation influence topography.Here we present a step -wise increase towards reality by investigating the topographic response to changes in both precipitation and climate at the same time.The amplitude of change prescribed for both precipitation and vegetation is based upon the present empirical relationship observed in the Chilean study areas for initial vegetation covers of 10 and 70%, and mean annual precipitations for 10 and 360mm/yr (Fig. 5).As with the previous experiments, oscillations in parameters were imposed upon steady -state topography that developed with the previous values, and a rock uplift rate of 0.2mm/yr.
Figure 15 shows the evolution of topographic metrics for simulations with combined oscillations in precipitation and vegetation.The variation in topographic metrics resembles those described for simulations with constant vegetation cover and oscillating climate by showing little to no significant adjustment towards new dynamic steady -state conditions.The amplitudes of oscillation are dampened from those of previous results because of the opposing effects of changes in precipitation and vegetation cover (e.g.compare blue and green curves in Figs. 8 and 9).
However, inspection of the predicted erosion rates (Fig. 16) for the combined oscillations indicates a significant (~0.1; -~0.15mm/yr), and highly non-linear response.The response between the 70% and 10% vegetation cover scenarios are very different such that for heavily vegetated areas (P(V=70%)) erosion rates typically increase during an oscillation, whereas for the low vegetation cover conditions (P(V=10%)) erosion rates initially show a decrease, and then a n increase and decrease at a higher frequency.
To better understand this contrast in the response to combined precipitation and vegetation changes, the first cycle of the imposed oscillation is shown in Figure 17.After an oscillation starts, the 10% initial vegetation cover simulations show a decline in erosion rates with the minimum erosion rate correlated with highest values of both vegetation cover and mean annual precipitation (compare top and bottom panels).This part of the response is interpreted a s resulting from the hindering effect of increased vegetation on erosion rates outweighing the impact of higher values of precipitation on erosion rates (Fig. 17).After values of vegetation cover and precipitation start to decline, erosion rates show a ve ry rapid increase to values of ~0.3mm/yr.This increase in erosion rates is due to an increase in both Kv and Kd (Fig. 3b, equations 4, 5) which outcompetes the effect of precipitation decrease.
Following this, a sudden drop in erosion rates to 0mm/yr occurs and lasts for 3kyrs due to the onset of hyper arid conditions at minimum precipitation.After this low in erosion rates, they increase again to 0.3mm/yr as precipitation and vegetation cover increase while the effect of increased precipitation outweigh s the effect of the non-linear decrease in Kv and Kd (Fig. 3b, c; equations 4, 5).Finally, at the end of this complex cycle a decrease in erosion rates occurs (Fig. 17b ) while vegetation and precipitation are increasing (upper panel) because the effect o f vegetation increases Kv/Kd and outweights the effect of increasing precipitation.
Lastly, a clearly different behavior in erosion rates occurs for settings with higher vegetation cover (e.g.P(V=70%) ) compared to the previous lower vegetation cover scenarios.As the vegetation cover and precipitation increase (Fig. 17a) in the first half of the 100kyr cycle, the erosion rates increase to values of approximately 0.35mm/yr.This is due the increase in precipitation which outcompetes the decline in erosivity/diffusivity parameters Kd and Kv.Following this, when vegetation cover and precipitation decrease in the second half of the cycle, little to no change occurs in the erosion rates.This near static behavior in erosion rates while precipitation an d vegetation cover decrease is due to an equilibrium between the negative effect on erosion rates for decreasing precipitation and the positive effect on erosion rates for decreasing vegetation cover.
In summary, the non-linear shape of the vegetation dependent erosivity (Kv) and hillslope diffusivity (Kd) in combination with linear effects of mean annual precipitation on erosion rates, exert a primary control on the direction and magnitude of change in catchment average erosion rates.Despite a simple os cillating behavior in precipitation and vegetation cover, a complex and non-linear response in erosion rates occurs.The implications of this are large for observational studies of catchment average erosion rates and suggest that the direction and magnitud e of response observed in a setting is highly dependent on the mean vegetation and precipitation conditions of the catchment, as well as what time the observations are made within the cycle of the varying vegetation and precipitation.Furthermore, these re sults highlight the need for future modeling studies (and motivation for our ongoing work), to investigate the response of catchment topography and erosion rates to more realistic climate and vegetation change scenarios, as well as a broader range of initial vegetation covers and precipitation rates than those explored here such that the threshold in behavior between th e two curves shown in figure 17b can be understood.
This could be achieved by using simulation results from state-of-art dynamic vegetation models (e.g Smith et al., 2014; see also companion paper by Werner et al., 2018) as inputs into the landscape evolution model or by full coupling between both model types.Vegetation simulations could, e.g.benefit from simulated changes in soil depth, wh ich can crucially determine plant water stress, provided by the landscape evolution model.Coupling between Landlab and the dynamic regional to global vegetation model LPJ-GUESS (Smith et al., 2014;Werner et al., 2018) is envisioned.

Model Restrictions and Caveats
Like other previous work on this topic (Collins et al., 2004;Istanbulluoglu and Bras 2005) the model setup used in this study was intentionally simplified to document how different vegetation and climate related factors impact topography over long (geologically relevant) timescales.We acknowledge that future model-studies should address some of the restrictions imposed by our approach to evaluate their significance for the results presented here.Future work should consider a transport-limited fluvial model or a fully-coupled alluvial sedimentation and transport model.Addition of this could bring new understanding in to how vegetation not only influences detachment limited systems, but also influences sedimentation and entrainment of material.This added level of complexity could however limit (due to computational concerns) the temporal scales over the investigation can be conducted.Future studies could improve upon this work by considering a more in-depth parameterization of how vegetation related processes (e.g.root depth and density, plant functional type) influence topographic metrics and erosion rates.Due to the long -timescales considered here, mean annual precipitation rather than a stochastic distribution of precipitation was imple mented.Future work should evaluate how stochastic distributions in precipitation and extreme events in arid, poorly vegetation settings, impact these results.
Regarding the vegetation and water-budget, a more sophisticated model of evapotranspiration and infiltration as a function to surface plant cover and plant functional traits such as rooting depth would improve model predictions and is a priority goal of future research within this project.Improvements will come from planned coupling of surface p rocess model with the dynamic vegetation model LPJ-GUESS (Werner et al., 2018, this issue).
We also acknowledge here that the transient forcings we have chosen for driving our model are simplistic and could be improved by a higher-fidelity time-series of climate over the last millennia.We choose a 100kyr, eccentricity driven, periodicity because it is widely recognized that the eccentricity cycles are a main control in driving Earths glacial cycle over the last 0.9Ma.While this approach is reasonable for a sensitivity analysis such as we've conducted, it prohibits a detailed comparison to observations without additional refinement.Our results suggest that a shorter periodicity, which would resemble other periodicities in the Milankovitch cycle (e.g., 41ky rs, 23kyrs) or shorter time scale climate variations, such as Heinrich events (see Huntley et al., 2013) would lead to smaller magnitudes of adjustment to new dynamic equilibria, because of short timespans in high-/low-erosive climate conditions within one period.

Conclusions
The results from our model-based experiments in comparison to observations from topographic analysis from four different areas in the Coastal Cordillera in Chile show that the interactions of vegetation cover and mean annual precipitation on the evolution of landscapes is a complex system with competing effects.Main conclusion which emerge from this study are: (I) vegetation cover in general has a hindering effect on eroding surface processes but the magnitude on which changes in vegetation cover affect these processes is a function of the initial state of the system.Changes in systems with higher initial values of vegetation cover have a less pronounced effect than changes in systems with lower initial vegetation cover.
(II) In comparison to the Coastal Cordilleras of Chile, the relationship between precipitation and surface vegetation cover shows a distinct shape: For a 10% increase in surface vegetation cover, the corresponding increase in mean annual precipitation is smaller in areas of lower vegetation cover and increases for areas with higher vegetation cover.This has an effect on transient topographies by shifting the equilibrium of vegetation and precipitation effects on erosion rates.
(III) Following our step-change simulations, our model results show different behaviors for changes in vegetation -cover and mean annual precipitation.While increases in mean annual precipitation have an increasing effect on erosion rates and therefore a long-term negative effect on topographic metrics, an increase in vegetation cover hinders erosion, and leads to higher topographic metrics.The magnitude of these changes is highly dependent on the initial vegetation cover and precipitation before the step-change.
(IV) Simulations with either oscillating vegetation cover or oscillating precipitation show adjustments to new dynamic mean values around which the basin metrics oscillate.The magnitude of adjustment is highly sensitive to initial vegetation cover, where simulations with 10% initial cover show higher magnitudes than simulations with 70% cover, for oscillating vegetation.Oscillating precipitation leads to lower-/no adjustments but an oscillation of basin metrics around the initial mean values with generally lower amplitudes compared to simulations with oscillating vegetation cover.
(V) Simulations with coupled oscillations of both vegetation cover and precipitation show only small magnitudes of adjustments of topography metrics to new dynamic equilibriums similar to simulatio ns with a oscillation in only precipitation.However corresponding erosion rates show a complex pattern of rapid increases and decreases which results from a interplay of competing effects of hindering of erosion by vegetation and aiding of erosion by prec ipitation.
Taken together, the above findings from this study highlight a non -linear and highly variable behavior in how variations in vegetation cover impact erosion and topographic properties.The complexity in how vegetation cover and precipitation changes influence topography highlights the need for future work to consider both of these factors in tandem, rather singling out either parameter (vegetation cover or precipitation) to understand potential transients in topography.
Earth Surf.Dynam.Discuss., https://doi.org/10.5194/esurf-2018-13Manuscript under review for journal Earth Surf.Dynam.Discussion started: 21 February 2018 c Author(s) 2018.CC BY 4.0 License.compared to positive amplitudes.For simulations with V = 70% th e reaction and adjustment to the new dynamic steadystate is less pronounced with a decline in relief from 410m to 407m (Factor 1.01) with positive and negative amplitudes in dynamic steady-state of 1.6m.Analysis of mean basin slope for the model topograp hies with low (V=10%) vegetation positive step-changes (Fig. 8a,b,c) the adjustment time for changes in vegetation cover to reach a new equilibrium in Earth Surf.Dynam.Discuss., https://doi.org/10.5194/esurf-2018-13Manuscript under review for journal Earth Surf.Dynam.Discussion started: 21 February 2018 c Author(s) 2018.CC BY 4.0 License.
, compared to the oscillating vegetation cover experiments(Figs 10, 11) is due to a generally higher impact of changes in vegetation cover on parameters which guide erosion rates and therefore adjustment to topographic metrics, compared to the calibrated, corresponding changes in precipitation in our model domains.Especially for simulations with low initial vegetation cover the effect of changing vegetation shows larger magnitude effects because of the non -linear response of Earth Surf.Dynam.Discuss., https://doi.org/10.5194/esurf-2018-13Manuscript under review for journal Earth Surf.Dynam.Discussion started: 21 February 2018 c Author(s) 2018.CC BY 4.0 License.

Figure 1
Figure 1 Overview of the geographic location, precipitation, and vegetation cover of the Coastal Cordillera, Chile studies areas used for model setup up.A) Digital topography of the areas considered and corresponding to the EarthS hape (www.earthshape.net)focus areas where ongoing related research is located.B) Observed present day mean annual precipitation from the WorldClim and CHELS A datasets used as model input.B) Present day maximum surface vegetation cover from MODIS data.

Figure 2
Figure 2 Normalized basin metrics for study-areas derived from 90 m S RTM digital topography from the study areas shown in Figure 1a Dots represent cumulative mean values calculated for all locations using 5-8 representative catchments in each area.Dotted lines represent linear interpolation for normalized slope, relief and channel steepness (k sn) values.Note the gradual increase, then decrease in all metrics around study area at ~32°S .

Figure 3
Figure 3 Example model setup used in simulations in this study.Boundary conditions and parameterizations used in the models are labeled.Blue colors represent low elevations, brown colors represent highe r elevations.Additional details of parameters used are specified in Table1.

Figure 4
Figure 4 Transient forcings in vegetation (a,b) and precipitation (c,d) considered in model experiments.S imulations were run for 15 Myr prior to the runtime show in the figure.All transients imposed started a runtime of 5 Myr.A) Variations in vegetation cover imposed in the oscillating experiment conditions for initial ve getation cover of 10 and 70%.B) Positive and negative step change parameterizations for vegetation cover.C) Variations in mean oscillating annual precipitation.D) Positive and negative step changes in mean annual precipitation used in experiments.The i nitial precipitation amounts prior to the transient oscillations or step changes correspond to the observed precipitation corresponding to the vegetation cover in each of the observed study areas (Fig. 1).S ee also Figure 5.

Figure 5
Figure 5 Graphical representation of the observed precipitationvegetation relationship in the focus areas (Fig. 1) and how precipitation amounts were selected when perturbations in vegetation cover were imposed.Black dots represent vegetationprecipitation values used in the steady-state model conditions and prior to any transients.Red dots show how vegetation cover perturbations in +/-10% in the model simulations were used to select corresponding mean annual precipitation amounts.

Figure 6
Figure 6 Predicted values of hillslope diffusivity Kd (solid line) and fluvial erodibility Kv (dashed line) as a function of vegetation surface cover.Although absolute values can't be compared due to different units, the shape of the curves representing the different parameters show different sensitivities to changes in vegetation cover.Fluvial erodibility shows highest magnitude of change for vegetation cover values < 25% whereas hillslope diffusivity reacts in a more linearly with highest chan ge below < 65% vegetation cover.

Figure 7 S
Figure 7 S teady-state model predicted (shaded regions) and observed (red dots) topographic metrics from the study areas 913

Figure 9
Figure 9 Mean catchment-wide erosion rates after step-change disturbance in model boundary conditions.Blue lines represent erosion rates for models with changes in only precipitation, green lines represent erosion rates for models with changes in only vegetation cover.Panels a,b show evolution after positive step-change, panels c,d for models with negative step-change.Note that the direction of change (positive or negative) from the initial state is in opposite directions for precipitation and vegetation cover changes.This effect is manifested in the subsequent plots.

Figure 10
Figure 10 Evolution of topographic metrics for simulations with oscillating surface vegetation cover and constant precipitation corresponding to the initial vegetation cover prior to the transient i n vegetation cover.Panels a,b,c show mean basin relief, mean basin slope and mean basin channel stee pness (ksn), respectively.

Figure 11
Figure 11 Predicted mean catchment erosion rates for simulations with oscillating surface vegetation cover and constant precipitation.

Figure 12 Figure 13
Figure 12 Evolution of topographic metrics for simulations with oscillating mean annual precipitation and constant vegetation cover.The vegetation cover was held constant at the value corresponding to the precipitation rate prior to the onset of the transient at 5000 kyrs.Panels a,b,c show mean basin relief, mean basin slope and mean basin channel steepness (k sn), respectively.

Figure 14 SFigure 15
Figure 14 S hifts in mean basin slope/mean basin relief relationship for simulations with positive and negative step-changes in 954 This was done to evaluate if our implementation of the governing equations in Section 3.4 produced topographies within reason of present day topographies in the four Chilean areas.A more detailed model calibration is beyond the scope of this study, and not meaningful without additional observational constraints on key parameters such latitudinal variations in the rock uplift rate and erosivity.