Accurate simulation of transient landscape evolution by eliminating numerical diffusion: the TTLEM 1.0 model

. Landscape evolution models (LEM) allow studying how the earth surface responses to changing climatic and tectonic forcings. While much effort has been devoted to the development of LEMs that simulate a wide range of processes, 10 the numerical accuracy of these models has received less attention. Most LEMs use first order accurate numerical methods that suffer from substantial numerical diffusion. Numerical diffusion particularly affects the solution of the advection equation and thus the simulation of retreating landforms such as cliffs and river knickpoints with potential consequences for the integrated response of the simulated landscape. Here we test a higher order flux limiting finite volume method that is total variation diminishing (TVD-FVM) to solve the partial differential equations of river incision and tectonic displacement. 15 We show that the choice of the TVD-FVM to simulate river incision significantly influences the evolution of simulated landscapes and the spatial and temporal variability of catchment wide erosion rates. Furthermore, a 2D TVD-FVM accurately simulates the evolution of landscapes affected by lateral tectonic displacement, a process whose simulation was hitherto largely limited to LEMs with flexible spatial discretization. We implement the scheme in TTLEM, a spatially explicit, raster based LEM for the study of fluvially eroding landscapes in TopoToolbox 2.


Introduction
Landscape evolution models (LEMs) simulate how the earth surface evolves in response to different driving forces, including tectonics, climatic variability and human activity.LEMs are integrative because they amalgamate empirical data and conceptual models into a set of mathematical equations that can be used to reconstruct or predict terrestrial landscape evolution and corresponding sediment fluxes (Glotzbach, 2015;Howard, 1994).Studies that address how climate variability and land use changes will affect landscapes in the long term increasingly rely on LEMs (Gasparini and Whipple, 2014).
Landscape evolution is not always smooth and gradual.Instead, sudden tectonic displacements along tectonic faults can create distinct landforms with sharp geometries (Whittaker et al., 2007).These topographic discontinuities do not necessarily smooth out over time but may persist over long timescales in transient landscapes (Mudd, 2016;Vanacker et al., 2015).For example, faults may spawn knickpoints along river profiles.These knickpoints will propagate upstream as rapids or water falls (Hoke et al., 2007), thereby maintaining their geometry through time (Campforts and Govers, 2015).After an uplift pulse, the river will only regain a steady state when knickpoints finally arrive in the uppermost river reaches.Transiency is not limited to individual rivers but also affects entire orogens such as the Southern Alps of New Zealand where the landscape may never reach a condition of steady state due to the permanent asymmetry in vertical uplift, climatically driven denudation and horizontal tectonic advection (Herman and Braun, 2006).
Transient "shocks" and topographic discontinuities are inherently difficult to model accurately.Most of the widely applied LEMs use first-order accurate explicit or implicit fi-B.Campforts et al.: The TTLEM 1.0 model nite difference methods to solve the partial differential equations (PDEs) that are used to simulate river incision (Valters, 2016).These schemes suffer from numerical diffusion (Campforts and Govers, 2015;Royden and Perron, 2013).Numerical diffusion will inevitably lead to the gradual disappearance of knickpoints and will result in ever-smoother shapes.It has already been shown that numerical smearing decreases the accuracy of modeled longitudinal river profiles (Campforts and Govers, 2015).Here, we hypothesize that it is also relevant for the simulation of hillslope processes: hillslopes respond to river incision and inaccuracies in river incision modeling will thus propagate to the hillslope domain.
Whether and to what extent this occurs is still unexplored.
Tectonic displacement is similar to river knickpoint propagation; in both cases, sharp landscape forms are laterally moving.Numerical diffusion may therefore significantly alter landscape features when tectonic shortening or extension is simulated using first-order accurate methods.In principle, flexible gridding overcomes this problem through dynamically adapting the density of nodes on the modeling domain to the local rate of topographic change.However, models using flexible gridding have other constraints.They are more difficult to implement and impose the structure of the numerical grid on the natural drainage network since rivers must follow the grid structure.Furthermore, the output of flexible grid models is not directly compatible with most software that is available for topographic analysis.
Here we present TTLEM (TopoToolbox Landscape Evolution Model), a spatially explicit raster-based LEM, which is based on the object-oriented function library TopoToolbox 2 (Schwanghart and Scherler, 2014).Contrary to previously published LEMs, we solve the stream power river incision model using a flux-limiting finite volume method (FVM), which is total variation diminishing (TVD), in order to avoid numerical diffusion.Our numerical scheme expands on previous work (Campforts and Govers, 2015) by extending the mathematical formulation of the TVD method from one-dimensional to entire river networks.Moreover, we develop a two-dimensional TVD-FVM scheme to simulate horizontal tectonic displacement on regular grids, which enables simulation of three-dimensional variations in tectonic deformation.The objective of this paper is to evaluate TTLEM and assess the performance of the numerical methods for a variety of real and simulated topographic and tectonic situations.

Tectonic deformation
In its simplest form, tectonic processes are represented by their kinematics and the assumed vertical surface deformation field U (x, y, t) [L T −1 ].However, many tectonic configurations imply that displacements have both a vertical (uplift or subsidence) and a lateral (extension or shortening) component (Willett, 1999;Willett et al., 2001).The change in elevation of the earth surface over time due to lateral tectonic displacement, excluding vertical rock uplift (∂z/∂t) td , is then where u and v [L T −1 ] are the tectonic displacement velocities in the cardinal directions (horizontal u and vertical v).

River incision
Detachment-limited fluvial erosion (∂z/∂t) fluv is calculated with the stream power law (SPL) (Howard and Kerby, 1983): The equation is solved on a dendritic stream network domain , where x refers to the distance from the outlet.A [L 2 ] is catchment area and proxy for the local discharge, and K [L 1-2m T −1 ] is an erodibility parameter that depends on local climate, hydraulic roughness, lithology and sediment load.m and n are the area and slope exponents: their values reflect hydrological conditions, channel width and the dominant erosion mechanism.K, m and n are interdependent and it is usually impractical to constrain any of their values alone (Croissant and Braun, 2014;Lague, 2014).Thus, many studies provide estimates for the m/n ratio.For m/n ratios between 0.35 and 0.8, K values span several orders of magnitude between 10 −10 and 10 −3 m (1-2m) yr −1 (Kirby and Whipple, 2001;Seidl and Dietrich, 1992;Stock and Montgomery, 1999).

Hillslope processes
River incision drives the development of erosional landscapes by setting the base level for hillslope processes.Steepening of hillslope toes leads to increased sediment fluxes from hillslopes to the river system.Hillslope denudation (∂z/∂t) hill is equal to the divergence of the flux of soilregolith material (q s , [L 3 L −1 T −1 ]): Different geomorphological laws describe hillslope response to lowering base levels.The model of linear diffusion assumes that the soil-regolith flux is proportional to hillslope gradient ∇z (Culling, 1963): where D is the diffusivity [L 2 T −1 ] that parameterizes hillslope erodibility and determines rate of soil-regolith creep.Main controls on variations of D include substrate, lithology, soil depth, climate and biological activity.Values of D range between 10 −3 and 10 −1 m 2 yr −1 for slopes under natural land use (Campforts et al., 2016a;DiBiase and Whipple, 2011;Jungers et al., 2009;Roering et al., 1999;West et al., 2013).Linear hillslope diffusion produces convex upward slopes.Field evidence, however, suggests that the linear diffusion model is only rarely appropriate (Dietrich et al., 2013).Instead, hillslopes often tend to have convex to planar profiles because rapid, ballistic particle transport and shallow landsliding dominate when slopes approach or exceed a critical angle (DiBiase et al., 2010;Larsen and Montgomery, 2012).To account for this rapid increase of flux rates with increasing slopes, Andrews and Bucknam (1987) and Roering et al. (1999) proposed a nonlinear formulation of diffusive hillslope transport, assuming that flux rates increase to infinity if slope values approach a critical slope S c : 2 . (5)

Final model
In summary, TTLEM solves the following PDE, whereby an explicit distinction is made between the fluvial and hillslope domain.The fluvial domain is determined by cells having a contributing drainage area exceeding a critical drainage area (A c ): The detachment-limited incision model assumes that rivers incise directly into bedrock and instantaneously excavate all material entering rivers from adjacent hillslopes.Material fluxes on slopes mobilize either soil or regolith that have different bulk density than the bedrock.This is accounted for by multiplying the rock uplift rate with the density ratio between ρ r and ρ s [M L −3 ] representing the bulk densities of the bedrock and the regolith material, respectively (Perron, 2011).

Implementation and numerical schemes of TTLEM
We solve Eq. ( 6) using a set of numerical schemes that we implement in the software TTLEM (see also Fig. A1 in the Appendix).TTLEM is written in the MATLAB programming language and in C-code where this significantly improves performance (e.g., for the nonlinear hillslope diffusion algorithm of Perron, 2011).Integrating TTLEM into TopoToolbox (Schwanghart and Kuhn, 2010;Schwanghart and Scherler, 2014) provides access to efficient algorithms of digital elevation model (DEM) analysis, as well as numerous routines for visualizing and analyzing modeling outputs.In the following sections, we discuss the numerical schemes of TTLEM to solve the PDEs described in the previous section.
The section numbers correspond to the processes indicated in the model flowchart in the Appendix (Fig. A1).

Drainage network development
TopoToolbox provides a function library for deriving the drainage network and terrain attributes (Schwanghart and Scherler, 2014).The calculation of flow-related terrain attributes, i.e., data derived from flow directions, relies on a set of highly efficient algorithms that exploit the directed and acyclic graph structure of the river flow network (Phillips et al., 2015).Nodes of the network are grid cells and edges represent the directed flow connections between the cells in downstream direction.Topological sorting of this network returns an ordered list of cells in which upstream cells appear before their downstream neighbors.Based on this list, we calculate terrain attributes such as upslope area with a linear scaling, thus enabling efficient calculation (O(n)) at each time step even for large grids (Braun and Willett, 2013).DEMs of real landscapes frequently contain data artifacts that generate topographic sinks.Topographic sinks can also occur during simulations when diffusion on hillslopes creates "colluvial wedges" that dam sections of the river network.By adopting algorithms of flow network derivation from Topo-Toolbox, TTLEM makes use of an efficient and accurate technique for drainage enforcement to derive non-divergent (D8) flow networks (Schwanghart et al., 2013;Soille et al., 2003).Based on the thus-derived flow network, TTLEM uses downstream minima imposition (Soille et al., 2003) that ensures that downstream pixels in the network have lower or equal elevations than their upstream neighbors.

Tectonic displacement
We implement a two-dimensional version of a flux-limiting total volume method to reduce numerical diffusion when simulating tectonic displacements on a regular grid.Equation (1) can be written as a scalar conservation law: where f (z) u = uz and f (z) v = vz are the flux functions of the conserved variable z.We refer to the Supplement of Campforts and Govers (2015) for a derivation of the differential form of Eq. ( 7), which can be converted to a numerical semiconservative flux scheme: where z k i,j is the elevation of the cell at row i and column j at time k × t. f represents the numerical approximation of the physical fluxes from Eq. ( 7).The incoming and outgoing www.earth-surf-dynam.net/5/47/2017/Earth Surf.Dynam., 5, 47-66, 2017 fluxes are approximated with a flux-limiting upwind method, which is TVD.A TVD scheme prevents the total variation of the solution to increase in time and hence prevents spurious oscillations that are associated with higher-order numerical methods (Toro, 2009).The flux limiter entails the method having a hybrid order of accuracy, being second-order accurate in most cases but shifting to first-order accuracy near discontinuities.Hence, the TVD-FVM method achieves two desirable properties: a higher order of accuracy than firstorder schemes and high numerical stability (Harten, 1983).TTLEM uses a staggered Cartesian grid for numerical discretization.The DEM grid centers represent the center of the computational cells, whereas the velocity fields (u and v) are located at the cell faces.
The numerical TVD fluxes are calculated following Toro (2009).In the following, we illustrate how to derive the flux over one out of the four cell boundaries: where f HI and f LO represent the high-and low-order fluxes, respectively: The low-order fluxes are solved with a first-order explicit upwind Godunov (1959) scheme: The high-order fluxes are solved with a Lax-Wendroff scheme (Lax and Wendroff, 1960): From Eqs. ( 10), ( 11) and ( 12) it follows that φ i+ 1 2 ,j represents the flux limiter, which is solved with the van Leer (1997) scheme: where r is a smoothness index calculated as The overall performance of the TVD-FVM is evaluated by comparing it with the first-order accurate upwind Godunov scheme (Godunov, 1959), which is not flux limiting Eq. ( 11).
In the remaining part of the text, we refer to this scheme as the first-order Godunov method (GM).

Numerical solution
TTLEM features a one-dimensional version of the fluxlimiting TVD-FVM to solve for river incision (Eq.2), which is written as a scalar conservation law: where f (z) represents the flux function of the conserved variable z, representing the river elevation.The method resembles the one described in the previous section but differs in that fluxes are calculated in one direction on a directed acyclic graph (Phillips et al., 2015).We refer to the Supplement provided by Campforts and Govers (2015) for a full derivation of this scheme.
In addition, we implement a first-order implicit FDM for the solution of the SPL detailed in Braun and Willett (2013).The method provides stable solutions regardless of the time step length, a property desired when simulating landscape evolution over long timescales and large spatial domains.Explicit schemes of river incision (both FDM and TVD-FDM), in turn, require time steps that satisfy the Courant-Friedrich-Lewy condition (CFL): where u max is the maximum velocity dictated by few river cells with high-drainage areas.Compared to these velocities, hillslope processes modeled by the linear diffusion equation are usually slow.Applying longer time steps for hillslope processes is a computational advantage that an implicit scheme increases even more (Pelletier, 2008).TTLEM thus uses two time steps: an outer time step ( t outer ) during which hillslope processes and the planform river network are calculated, and an inner time step ( t inner ) nested within the outer time step that is used to solve for river incision.Thus, while t outer should satisfy the CFL criterion for the explicit linear or nonlinear diffusion equation, the t inner is flexible and adheres to the CFL criterion of the explicit river incision method (Fig. A1).The adoption of implicit methods allows the relaxation of both time step constraints.However, TTLEM allows limits to be set to t outer and t inner , and it enables us to investigate the impact of the length of the time step on model outcomes (see Sect. 4.1.3).

Analytical solution
Ideally, numerical methods are benchmarked against analytical solutions.Albeit analytical solutions are available for spe-cific initial and boundary conditions only, they are accurate and grid-resolution independent, contrary to numerical solutions where model parameter values might depend on the grid resolution (Pelletier, 2010).We implemented an analytical solution for the SPL as an independent benchmark to compare the performance of the different numerical schemes of river incision under conditions where an analytical solution is available.First, we created an artificial DEM with topography in steady state between uplift and erosion (see Table 1).From this DEM, we extracted the drainage network and corresponding river elevations by selecting all cells exceeding 10 6 m 2 .Very short river profiles (< 10 km) are excluded from the analysis.Subsequently, we simulate landscape evolution using the numerical models documented in the previous sections assuming spatially invariant uplift rates.After each simulation, we obtain river elevations from the resulting DEMs and compare them with river elevations that we derived analytically using the pre-uplift, steady-state river profiles as input.Analytical solutions for the stream power law are based on the slope patch method of Royden and Perron (2013) that non-dimensionalizes the stream power law using a dimensionless height (λ) and transformed horizontal distance metric χ: where z x represents the dimensionless elevation along the river profile, h 0 is a reference length scale (set to 1 m) and A 0 is a reference value for the drainage area (set to 1 × 10 6 m 2 ).
To integrate over abrupt changes in the drainage area along the rivers, Eq. ( 19) is solved using the rectangle rule (Mudd et al., 2014).Steady-state river profiles appear as straight lines in this nondimensional coordinate system.The analytical slope patch solution then calculates the evolution of a dimensionless river profile in response to uplift.The method is detailed in the Appendix of Royden and Perron (2013) and is based on tracing individual patches that are initiated at the outlet of the drainage network and propagate upstream with a velocity dictated by upstream area and the parameters of the SPL (Eq.2).We applied the slope patch solution to the steady-state preuplift river profiles using the simulated uplift rates as input.We also assessed the accuracy of the numerical methods with the root mean squared error (RMSE): where z i,analytical and z i,numerical refer to the analytically and numerically calculated elevation of a river cell, respectively, and n riv is the total number of river cells.

Hillslope processes
We implemented linear hillslope diffusion using the implicit Crank-Nicolson scheme (Pelletier, 2008).The scheme is unconditionally stable at large time steps.A numerical solution of the nonlinear hillslope equation, however, is more demanding.The maximum time step length of an explicit FDM sharply decreases as slopes approach the threshold gradient.To overcome this restriction, Perron (2011) developed Q-imp, an implicit solver that allows the increase of time step lengths by several orders of magnitude.Conversely, the per-operation computational cost of this algorithm is higher in comparison to the explicit solution, and the overall performance of this method is better than alternative solutions (Perron, 2011).Q-imp efficiently calculates hillslope diffusion even for high-resolution simulations.However, rapid incision during one time step may generate slopes along rivers that are greater than the threshold slope, a situation that Qimp cannot solve.An approach is thus needed that adjusts hillslopes to the threshold slope prior to calculating Q-imp.We assume that hillslopes instantaneously adjust to oversteepening by mobilizing the amount of material required to reduce the slope gradient to the threshold value S c (Burbank et al., 1996).We refrain from simulating individual landslides although we acknowledge that single high-magnitude low-frequency events may be relevant at the timescales of our simulations (Korup, 2006).Instead, our approach implicitly accounts for the combined effects of a large number and variety of landslides that effectively adjust slopes to a threshold slope.This threshold slope can be thought of as "an average effective angle of internal friction, which controls hillslope stability" (Burbank et al., 1996).We implement this hillslope adjustment using a modified version of the excess topography algorithm (Blöthe et al., 2015).In this algorithm, elevations z at time step t + 1 are calculated so that the absolute local gradient at each grid cell becomes less than or equal to S c .This is achieved by decreasing elevations at location i to the minimum elevation of all other locations j , to which we add an offset calculated as the product of the Euclidean distance i, j and S c : The equation above entails that z t+1 i at one location depends on all other grid cells and that the algorithm has a time complexity of O(N 2 ), which would render it unsuitable for frequent updating during LEM simulations.To avoid an excessively high computational load, we implement the algorithm using morphological erosion with a grayscale structuring element (see MATLAB function ordfilt2), which is a minimum sliding window with additive offsets calculated from the window size and S c .This significantly reduces run times since we calculate elevations at one location from the sliding window.However, this approach may retain gradients greater than S c at steep-and long-slope sections.We solve this by www.earth-surf-dynam.net/5/47/2017/Earth Surf.Dynam., 5, 47-66, 2017

Impact of numerical methods
We investigate how numerical schemes implemented in TTLEM affect simulated landscape evolution.As we focus on evaluating the schemes' performance, all simulations have synthetically generated landscapes as initial surfaces.Hence, our simulations are uncalibrated and results remain untested against an actual landscape: however, the chosen parameter values are within the range of previous studies (e.g., Gasparini and Whipple, 2014; Whipple and Tucker, 1999).We distinguish between the effects on simulated river incision on the one hand and on simulated tectonic displacement on the other.To investigate the accuracy and implications of river incision methods, we compare the explicit TVD-FVM with the first-order implicit FDM and further differentiate between the implicit FDM where no limitation is set on the time step and the implicit FDM where the CFL criterion limits the time step length.To investigate the accuracy and implications of river incision methods we compare an explicit first-order GM with the two-dimensional TVD-FVM.

One-dimensional river incision
The impact of numerical diffusion on propagating river profile knickpoints is most obvious in situations where an analytical solution is available.The first simulation illustrates such a situation, with an artificial river profile characterized by a major knickzone between 8 and 12 km from the river head (Fig. 1).We assume that the drainage area is increasing in proportion to the square of the distance and uplift equals zero.For this simplified configuration, an analytical solution for the SPL relies on the method of characteristics (Luke, 1972).Notwithstanding the relatively high spatial resolution of 100 m, the first-order implicit FDM suffers from considerable numerical diffusion when river incision is calculated over a time span of 1 My (Fig. 1).The TVD-FVM systematically achieves a much higher accuracy over a wide range of spatial resolutions and parameter values (Campforts and Govers, 2015).

Drainage network
We assess the numerical accuracy of the entire drainage network with spatially and temporally constant values for all model parameter values (Table 1), assuming a fixed drainage network (see Sect. 3.3.2).We first create a steady-state artificial landscape (Fig. 2) on a 50 km × 100 km grid with a spatial resolution of 100 m that we initialize with uniformly distributed random elevation values between 0 and 50 m (Movie S1 in the Supplement).Our simulation uses   Dirichlet boundary conditions and inserts a spatially and temporally uniform vertical uplift of 1 km My −1 over a period of 150 My. t outer is set to 5 × 10 4 years.Following steady state, we impose four consecutive sinusoidal uplift pulses of equal magnitude on this artificial landscape over 1 My.Each uplift pulse has a wavelength of 0.25 My and an amplitude of 3×10 −3 m yr −1 (Fig. 3).We repeat the simulations with three different numerical schemes (implicit FDM without time step limitation, implicit FDM with time step limitation (CFL condition applied) and TVD-FVM), each at 22 different spatial resolutions (6.25, 12.5, 25, 50, 100, 150, ..., 950 m).Hillslopes are simulated using linear hillslope diffusion in combination with threshold slopes, a configuration typically used to simulate landscape evolution at geological timescales (e.g., Goren et al., 2014).The threshold slope is set to 0.8 m m −1 and hillslope diffusivity is 0.01 m 2 yr −1 .We record the CPU time required to run a 1 My simulation to assess computational performance.In order to facilitate the high-resolution run (at 6.25 m where the spatial domain covers 7950×15 950 cells), all model runs were executed on one computational node of the Flemish Super Cluster (VSC) using a single core (Broadwell, E5-2680v4) and 128 Gb RAM.We evaluate the numerical performance of the schemes and the impact of spatial resolution against an analytical solution (slope patch method) for the entire drainage network represented by all cells exceeding 1 km 2 (Fig. 2).
Figure 4 compares results obtained from the numerical methods and the analytical solution.The initial river profiles slightly differ depending on spatial resolution due to interpolation of the steady-state artificial landscape with a spatial resolution of 100 m.The results show that TVD-FVM and implicit numerical solutions converge at increasing spatial resolutions.Where the time step of the implicit scheme is unbounded by the CFL criterion, however, the solution deviates from those adhering to the CFL criterion.This illustrates that there is trade-off between numerical accuracy and numerical stability for an implicit scheme at long time steps.In addition, an implicit scheme at high spatial resolution and large time steps fails to converge to an analytical solution because uplift is modeled as a discrete stepwise function rather than a continuous function (e.g., the sinusoidal uplift history used here) that inserts artificial shocks in the solution.
The TVD-FVM is consistently more accurate than the implicit methods at all spatial resolutions, although the implicit FDM (CFL < 1) approaches the high accuracy of the TVD-FVM at very high resolutions (6.25 m) (Fig. 5a).At lower spatial resolutions (> 10 m) the numerical accuracy of the TVD-FVM is significantly higher compared to the accuracy obtained with the implicit methods at the cost of a slightly increased additional computation time.To achieve the same numerical accuracy as the TVD-FVM at 500 m spatial resolution (RMSE = 18.17, model runtime = 2.89 s), the implicit method (CFL < 1) would need to be evaluated at 150 m, which would take 12 times longer (model runtime = 36 s) (Fig. 5b).

River incision and catchment-wide erosion rates
We hypothesize that the diffusive nature of commonly applied first-order FDMs is not restricted to the simulation of river longitudinal profiles but has systematic consequences for other measures derived from LEM simulations.Such measures include catchment-wide erosion rates that consti-  tute the basis for model-field data comparison and model parametrization (Gasparini and Whipple, 2014;Moon et al., 2015).In order to investigate the sensitivity of LEMderived catchment-wide erosion rates to different numerical schemes of the river incision model, we use the steady-state artificial landscape described in the previous experiments (Sect.4.1.2).The simulation runs over 5 My with four consecutive uplift pulses of equal amplitude and a wavelength of 1.25 My with Dirichlet boundary conditions and a planform fixed drainage network.We use two spatial resolutions (100 and 500 m) and three different numerical methods (implicit FDM without time step limitation, implicit FDM with time step limitation (CFL condition applied) and TVD-FVM) to simulate river incision.The maximum length of t inner is set to 3×10 3 yr for all schemes to ensure that the implicit method converges at higher resolutions too (see Sect. 4.1.2).Hillslope response is simulated using a linear diffusion scheme in combination with a threshold slope (S c , see Fig. A2).
We compare differences in simulated erosion rates by randomly selecting > 200 catchments with drainage areas ranging between 1 and 50 km 2 (Fig. 7).We calculate the erosion rates for each time step by subtracting the elevation grid in the previous time step from the updated, current elevation grid.The sum of elevation differences within each catchment refers to the catchment-wide erosion rate integrated over the time step length.For each catchment, we then derive the difference between erosion rates calculated by the different numerical schemes and summarize them using the RMSE statistics (O TVD−FDM ): where ε i,TVD and ε i,FDM refer to the catchment-wide erosion rates simulated with the TVD-FVM and FDM, respectively, to simulate river incision, and nb t is the total number of discrete time steps of the simulated erosion record.We rank the catchments in increasing order of O TVD−FDM for each simulation to investigate variations in catchmentwide erosion rates.Figure 6 shows the results for the catchments at the 10, 50 (median) and 90 % percentiles.Ranks are derived separately for the model runs at 100 and 500 m since different catchments are randomly generated for both simulation runs.The percentiles shown in Fig. 6 therefore represent different catchments.
www.earth-surf-dynam.net/5/47/2017/Earth Surf.Dynam., 5, 47-66, 2017 For most catchments, we detect differences in catchmentwide erosion rates between the three numerical methods at a spatial resolution of 100 m.Generally, the amplitude of the response to a tectonic uplift pulse increases when using TVD-FVM: the use of a first-order implicit FDM without time step restriction results in a much smoother response in comparison to the TVD-FVM.The variations in response amplitude are significant: the majority of the catchments record amplitude reductions by more than 50 % when modeled with the implicit FDM without time step restriction.Time step restriction (and thereby sacrificing the main advantage of the implicit FDM) significantly reduces numerical diffusion so that most catchments display an erosional response comparable to that simulated by the TVD-FVM.However, this is only true for simulations with a 100 m spatial resolution.The advantage of a time-step-restricted implicit FDM over a nonrestricted implicit FDM disappears almost completely for a coarser grid resolution of 500 m.
Figure 7 shows that erosion rates diverge between the different methods with increasing distance to the outlet of the main river, while they are similar for larger catchments.A smaller effect of the numerical scheme on large catchment areas may partly arise from stronger averaging of local variations in catchment erosion rates.In addition, catchments at a large distance from the outlet -and thus likely with smaller catchment areas -will experience upstream migrating knickpoints only after several model time steps.If catchments are far from the fault zone, knickpoints will then be significantly smoothed by a first-order accurate implicit FDM, which will ultimately affect the response of the catchment.Again, spatial resolution matters: a larger grid size not only results in larger differences on average but also in larger differences between small and large catchments (Fig. 7).
The differences in catchment response relate to the differences in simulated erosion rates within the catchments.Figure 8 illustrates the spatial difference in erosion rates calculated with the two numerical methods during the final step of the model run (after 5 My).This figure shows that spatial differences are significant and form a systematic banded pattern related to the upslope migration of the erosion waves of the individual uplift pulses.

Tectonic displacement
We test the performance of the two-dimensional version of the flux-limiting TVD-FVM to simulate tectonic displacement.A synthetic DEM forms the initial surface for a simulation of a constant lateral tectonic displacement with neither fluvial incision nor hillslope diffusion.Theoretically, this should result in a laterally displaced landscape that, apart from this displacement, remains unchanged in comparison to the initial state.We compare the flux-limiting TVD-FVM with a first-order accurate upwind GM simulating a tectonic displacement in two directions (u = v = 10 mm yr −1 ) over a time span of 1 My. Figure 9 illustrates that the explicit GM strongly smooths the resulting DEM whereas the twodimensional TVD-FVM scheme produces a DEM that is very similar to the initial DEM, with reduced amounts of numerical diffusion.
In order to quantify the amount of numerical diffusion (D N [L 2 yr −1 ]) introduced by the GM and the TVD-FVM method, we test a range of different model configurations and calculate the numerical diffusivity, D N , corresponding to the observed smoothing.D N is the diffusivity required to transform the initial DEM (DEM ini ) to the final DEMs produced at the end of the simulations (DEM fint ).The optimum amount of diffusion is determined by minimizing the misfit function H with a sequential quadratic programming method (Nocedal and Wright, 1999).H is given by where nb px is the number of pixels in the DEM.We find that numerical diffusivity of the GM exceeds commonly used values of hillslope diffusivities as soon as spatial resolution exceeds 90 m (Fig. 10a).The two-dimensional TVD-FVM decreases numerical diffusion by a factor of 5-60 compared to the GM (Fig. 10b).The accuracy increases for both schemes with increasing resolution and increasing CFL numbers.However, the gain in accuracy with increasing spatial resolution is higher for the TVD-FVM than for the GM.Our analysis shows that the explicit FDM performs best with a CFL criterion close to one where additional required iterations within a given time interval are at a minimum (Gulliver, 2007).

Discussion
Our analysis of numerical solvers focuses on three interrelated issues: numerical accuracy, spatial resolution and computational efficiency.Adopting highly simplifying assumptions allow us to benchmark the solvers against analytical solutions.Our focus is on testing an implicit, first-order accurate FDM against TVD-FVM.The implicit FDM has several desirable properties.It is unconditionally stable and tolerates time step lengths exceeding those prescribed by the CFL criterion.LEMs are often run over time spans of millions of years and the CFL criterion is dictated by a few grid cells with high upslope areas.Adopting an implicit scheme is therefore potentially interesting since it allows the decrease of the computation time while enabling simulations at high spatial resolutions.Our results, however, show that this major advantage vanishes if the aim of a LEM simulation is to capture transiency correctly.For CFL > 1 the implicit FDM introduces significant numerical smearing, and for CFL 1, the approach tends to insert an artificial shock wave of uplift because gradual uplift is approximated by a step function if time steps are (very) large.
For time step lengths approaching those prescribed by the CFL criterion, we show that computational gains by implicit FDM are marginal compared to TVD-FVM.The TVD-FVM code can be vectorized, i.e., it exploits single-instruction multiple-data parallelism to save CPU time.The implicit FDM requires a lower number of numerical operations but all stream network nodes need to be treated sequentially.Simulations at higher spatial resolutions increase the numerical accuracy and may balance the low accuracy of the implicit, first-order accurate FDM.Our results indicate that there is indeed a strong gain in numerical accuracy for all methods (Figs. 4 and 5) with increasing spatial resolution.However, to achieve the same numerical accuracy as the TVD-FVM, the implicit method with a CFL < 1 constraint requires the use of spatial resolution that is about 3 times higher, resulting in a computation time that is ∼ 12 times higher (Fig. 5).In summary, while a first-order implicit scheme is stable and accurate for long-term steady-state solutions (Braun and Wil- lett, 2013), it has severe shortcomings when simulating transient landscape evolution caused by knickpoint propagation in detachment-limited erosional basins.These shortcomings can, to a large extent, be avoided by using a TVD-FVM.
We also show that the impact of the numerical scheme used to simulate river incision is not limited to river profile development alone.Hillslopes adjust to local base level changes dictated by river incision.Hillslope denudation rates must therefore -at least partly -reflect the geometry and dynamics of a knickpoint and will respond differently to a diffuse signal that is the result of relatively slow, continuous uplift on the one hand and sharp discontinuity caused by a rapid base-level drop of major fault activity on the other hand.Our simulations show that, depending on the spatial and temporal resolution, catchment-wide erosion rates are more responsive to uplift when fluvial incision is calculated by TVD-FVM rather than by the first-order accurate implicit FDM.This is because first-order (explicit and implicit) FDMs fail to properly reproduce transient incision waves (Campforts and Govers, 2015) due to knickpoint smoothing.This also affects hillslope denudation since the drop in hillslope base level due to the passage of a knickpoint is smeared out in time when smoothing occurs.The response of catchmentwide erosion rates to uplift will therefore also be smoothed, resulting in significantly lower peak erosion rates.This effect will be most significant in upstream catchments that are The ratio between the amount of numerical diffusion for the first-order GM vs. the flux-limiting TVD-FVM.
far away from the base level since smoothing increases with time and knickpoint migration distance.
One might question the significance and necessity of numerical schemes that avoid diffusion of retreating knickpoints.Given the many assumptions and uncertainties that underlie many LEMs, numerical accuracy may seem like a problem of lesser importance.We argue that the simulations presented in this paper show that this is not the case and that it is indeed critical to simulate knickpoint retreat as accurately as possible.However, our analysis does not cover all situations wherein the accurate simulation of knickzones is important.Simulation of sharp knickpoints is also required in geomorphological and lithological settings where knickpoint retreat is caused by rock toppling, possibly triggered during extreme flood events (Baynes et al., 2015;Lamb et al., 2014;Mackey et al., 2014).Similarly, glacial incision often creates hanging valleys that are reshaped by migrating fluvial knickpoints after glacial retreat (Valla et al., 2010).In all of these cases simulation tools with a minimum of numerical diffusion are required to correctly quantify natural knickpoint diffusion and to study the underlying processes.
First-order numerical methods also inadequately simulate lateral tectonic displacement on a regular grid.The amount of numerical diffusion that is introduced by these methods will, in many cases, exceed natural diffusion rates, thus making accurate simulation of hillslope development impossible.A two-dimensional variant of the TVD-FVM reduces the amount of numerical diffusion to values well below natural diffusivity values, an effect that is especially apparent at high spatial resolutions.The two-dimensional TVD-FVM thus allows the accurate modeling of this process, which significantly impacts the evolution of topography and river networks (Willett, 1999), using a fixed grid.This was hitherto only possible with flexible spatial discretization schemes.
Although most LEMs use first-order accurate discretization schemes (Valters, 2016), the problem of numerical diffusion has been discussed in the broader geophysical community (Durran, 2010;Gerya, 2010).An alternative family of shock-capturing Eulerian methods are MPDATA advection schemes (Jaruga et al., 2015).These schemes are based on a two-step approach in which the solution is first approximated with a first-order upwind numerical scheme and then corrected by adding an anti-diffusion term (Pelletier, 2008).However, contrary to the TVD-FVM, the standard MPDATA scheme (Smolarkiewicz, 1983) is not monotonicity preserving (i.e., it is not TVD).Instead, MPDATA introduces dispersive oscillations in the solution if combined with a source term (such as uplift) in the equation (Durran, 2010).Adding limiters to the solution of the anti-diffusive step (Smolarkiewicz and Grabowski, 1990) renders the MP-DATA scheme oscillation free (Jaruga et al., 2015).However, by adding this additional correction, the method approaches the numerical nature of the TVD-FVM, which does not require further adjustments in any case.
Some of the weaknesses of the tested numerical solutions can be reduced by using LEMs that rely on irregular grid geometries.Irregular grids, for example, allow the simulation of tectonic shortening using a Lagrangian approach where grid nodes are advected with the tectonically imposed velocity field (e.g., Herman and Braun, 2006).In TTLEM the TVD-FVM solvers are implemented using a fixed grid, which has several advantages.First, input data such as topography, climate, lithology or tectonic displacement fields are typically available as raster datasets and thus require only minor modifications, whereas irregular grids require substantial preprocessing.Second, TTLEM output can instantly be analyzed and visualized using the TopoToolbox library (Schwanghart and Kuhn, 2010;Schwanghart and Scherler, 2014) or any other geographic information system.Thus, www.earth-surf-dynam.net/5/47/2017/Earth Surf.Dynam., 5, 47-66, 2017 while irregular grid geometries and flexible grids may have some advantages over rectangular grids, TTLEM's implementation of numerically accurate algorithms strongly reduces the shortcomings of rectangular grids while facilitating straightforward processing of model input and output.

Conclusion
Despite the growing interest in the development and use of LEMs, accuracy assessment of the numerical methods has received little attention.First-order accurate FDMs are the most commonly applied numerical methods.However, they introduce numerical diffusion and artificially smooth discontinuities that are inherent in transient landscapes.To overcome this problem, we developed the TVD-FVM.The TVD-FVM solves river incision more accurately than the first-order accurate FDMs with significant influences on the geometry of modeled river profiles and implications for catchment-wide erosion rates.Errors due to numerical diffusion depend on the spatial and temporal resolution as well as on the position of the catchment in the landscape.In addition, we introduce a two-dimensional version of the TVD-FVM that allows the simulation of lateral tectonic displacement with low numerical diffusion on a fixed computational domain.Our new numerical techniques are implemented in the open-access raster-based Landscape Evolution Model (TTLEM) contained within TopoToolbox.Together with numerical implementations of common hillslope process mod-els, TTLEM provides the community with a novel simulation tool for the accurate reconstruction, exploration and prediction of landscape evolution.In its current form, TTLEM is limited to uplifting, fluvially eroding landscapes.Further development will integrate other processes (e.g., glacial erosion) as well as the explicit routing of sediment through the landscape.
-surf-dynam.net/5/47/2017/calling the algorithm repeatedly until all slope values are less than or equal to S c .

Figure 1 .
Figure 1.Solution of the linear one-dimensional stream power law for a synthetic knickzone over a time span of 1 My.The analytical solution is obtained with the method of characteristics.The spatial resolution is 100 m.Table 1 lists other model parameter values.

Figure 2 .
Figure 2. A synthetic steady-state landscape produced as the testing environment to verify and compare the different numerical schemes implemented in TTLEM.Model runtime was 150 My, while uplift rate was assumed to be spatially uniform over the area (block uplift) and fixed to 1 km My −1 .Other model parameter values are listed in Table 1.Dynamic landscape evolution is presented in Movie S1.The gray lines indicate the drainage network for which the solution has been calculated analytically as a benchmark solution.The blue line indicates the river profile for which model results at different resolutions are plotted in Fig. 4.

Figure 3 .
Figure 3. Uplift imposed on the steady-state landscape shown in Fig. 2 to investigate the impact of different numerical schemes.

Figure 4 .
Figure 4. Comparison between different modeled resolutions for the river profile indicated in blue in Fig.2.The green line is the analytical "true" solution, obtained with the slope patch method ofRoyden and Perron (2013).The full red line represents the firstorder accurate implicit solution when the CFL < 1, and the dotted blue line represents the first-order accurate implicit solution when the time step is left free.The implicit solutions where CFL < 1 are simulated with a time step equal to the time step used for the TVD-FVM.

Figure 5 .
Figure 5. (a) Performance of the different numerical schemes where the RMSE is calculated between the analytical and numerical methods.(b) CPU time required to perform the model runs at the indicated resolutions.

Figure 6 .
Figure 6.Temporal variation in simulated catchment-wide erosion rates using different numerical methods to simulate river incision.The black lines represent simulations where a flux-limiting TVD-FVM is used, the blue lines represent the first-order accurate implicit FDM without constraints on the time steps, and the red lines represent the first-order accurate FDM with an inner time step calculated with the CFL criterion.(a) Simulations performed at a spatial resolution of 100 m.(b) Simulations performed at a spatial resolution of 500 m.Here, a median filter with a window of three time steps is applied to the simulated erosion rates to eliminate spikes that might occur at low resolutions.

Figure 7 .
Figure 7. Spatial variation of differences between simulated erosion rates calculated with a flux-limiting TVD-FVM for simulating river incision and a first-order accurate implicit FDM.Here, we compare methods that are both run with an inner time step constrained with the CFL criterion (see text).O TVD−FDM is thus calculated between the black and red lines from Fig. 6.The left column represents simulations run at a spatial resolution of 100 m, the right column at 500 m.(a, b) Location of the randomly selected catchments with an area > 1 and < 50 km 2 .Colors refer to the O TVD−FDM between the two simulations.(c, d) Differences between the schemes increase with increasing distance from the river outlets and are inversely correlated with the catchment area.

Figure 8 .
Figure 8. Spatial pattern of erosion rates during one model time step when simulating landscape evolution with the flux-limiting TVD-FVM vs. the first-order accurate implicit FDM.(a) Simulation at a resolution of 100 m where the time step of the implicit method is not constrained.(b) Simulation at a resolution of 100 m where the time step of the implicit method is constrained with the CFL criterion.(c) Simulation at a resolution of 500 m where the time step of the implicit method is not constrained.(d) Simulation at a resolution of 500 m where the time step of the implicit method is constrained with the CFL criterion.

Figure 9 .
Figure 9. Impact of numerical schemes when simulating horizontal shortening on a fixed grid.The simulations are performed at a spatial resolution of 50 m and a CFL of 0.5.(a) Extract from synthetically produced DEM from Fig. 2. (b) Horizontal shortening in two directions simulated with a two-dimensional explicit first-order GM.The green lines represent transects of the theoretically unchanged topography after lateral displacement.The red lines represent transects of the topography produced with the GM.(c) Horizontal shortening in two directions simulated with a two-dimensional explicit flux-limiting TVD-FVM (represented by black lines).

Figure 10 .
Figure 10.(a) Amount of numerical diffusion (D N ) introduced in the system when simulating lateral tectonic displacement in two directions as a function of raster resolution.The gray zone indicates the range of naturally observed diffusion rates.(b)The ratio between the amount of numerical diffusion for the first-order GM vs. the flux-limiting TVD-FVM.

Figure A1 .
Figure A1.Schematic representation of the TTLEM model flow.The numbered methods correspond with the paragraphs from Sect. 3 in the main text.

Figure A2 .
Figure A2.Hillslope response to river incision.(a) Standard SRTM DEM (30 m) included in TopoToolbox representing the Big Tujunga region.The dotted gray line indicates the location of the transect shown in panel (g).(b) Resulting topography after 500 000 years using four different descriptions for hillslope evolution.(c) Linear diffusion over all slope values (lin in panel g).(d) Threshold landscape where no slopes exceed the threshold slope (sc in panel g).(e) Linear diffusion combined with immediate adjustment to a threshold slope (lin and sc in panel g).(f) Nonlinear diffusion combined with immediate adjustment to a threshold slope (non-lin in panel g).(g) Elevation profiles of the different model runs compared with the initial profile.Model parameter values are listed inTable 1.

Table 1 .
Model parameters used for the TTLEM simulations.
Table 1 lists other model parameter values.