terrainbento 1.0: a Python package for multi-model analysis in long-term drainage basin evolution

Models of landscape evolution provide insight into the development of specific field areas, create testable predictions of landform development, demonstrate the consequences of current geomorphic process theory, and spark imagination through hypothetical scenarios. While the last four decades have brought the proliferation of many alternative formulations for the redistribution of mass by Earth surface processes, relatively few studies have systematically compared and tested 5 these alternative equations. We present a new Python modeling package, terrainbento 1.0, that enables multi-model comparison, sensitivity analysis, and calibration of Earth surface process models. terrainbento provides a set of 28 model programs that implement alternative transport laws related to four model elements: hillslope processes, surface-water hydrology, erosion by flowing water, and material properties. The 28 model programs stem from 13 binary choices related to one of these 10 four elements—for example, the use of linear or non-linear hillslope diffusion. terrainbento is an extensible framework: model base classes that treat the elements common to all models (such as input/output and boundary conditions) make it possible to create a new model without re-inventing these common methods. terrainbento is built on top of the Landlab framework, such that new Landlab components directly support the creation of new terrainbento models. terrainbento is fully doc15 umented, has 100% unit test coverage including numerical comparison with analytical solutions for process models, and continuous integration testing. We support future users and developers with introductory Jupyter notebooks and a template for creating new terrainbento model programs. In this paper, we describe the package structure, process model theory, and software implementation of terrainbento. Finally, we illustrate the utility of terrainbento with a benchmark example highlighting 20 the differences in steady state topography between five different process models. 1 Geosci. Model Dev. Discuss., https://doi.org/10.5194/gmd-2018-204 Manuscript under review for journal Geosci. Model Dev. Discussion started: 12 October 2018 c © Author(s) 2018. CC BY 4.0 License.


Introduction
Computational models of long-term drainage basin and landscape evolution have a wide spectrum of applications in geomorphology, ranging from addressing fundamental questions about how climatic and tectonic processes shape topography, to performing engineering assessments of landform stability and potential for hazardous-waste containment (see, e.g., reviews by Coulthard, 2001;Pazzaglia, 2003;Martin and Church, 2004;Willgoose, 2005;Codilean et al., 2006;Bishop, 2007;Tucker and Hancock, 2010;Willgoose and Hancock, 2011;Pelletier, 2013;Temme et al., 2013;Chen et al., 2014;Valters, 2016).Although the basic principles of drainage basin evolution are reasonably well understoodsuch as the fundamental concept that erosion is driven by gravitational and water-runoff processes, the latter of which depend strongly on surface gradient and water flow -uncertainty remains concerning the appropriate forms of the governing transport laws for any particular set of materials and environmental conditions (Dietrich et al., 2003).This situation creates a need for comparative testing in order to gauge the overall performance of various mathematical formulations, to identify knowledge gaps in areas where numerical models perform poorly, and to assess which transport laws are most appropriate for various environmental conditions, timescales, and spatial scales.
Relatively few studies have systematically compared and tested alternative transport laws, and those that do usually address only a single, quasi-one-dimensional landform element, such as the shape of an idealized hillslope (Doane et al., 2018), the longitudinal profile of a particular stream channel (Tomkin et al., 2003;van der Beek and Bishhop, 2003;Loget et al., 2006;Valla et al., 2010;Attal et al., 2011;Hobley et al., 2011;Gran et al., 2013), or a fault scarp (Andrews and Hanks, 1985;Andrews and Bucknam, 1987;Hanks, 2013;Pelletier et al., 2006).Two-dimensional numerical models have focused on testing alternative formulations of soil production and transport (Herman and Braun, 2006;Roering, 2008;Petit et al., 2009;Pelletier et al., 2011), and studies that test two-dimensional numerical models that couple hillslope and channel processes are much more limited (Hancock and Willgoose, 2001;Willgoose et al., 2003;Hancock et al., 2010;Gray et al., 2018).Models that combine hillslope and channel processes -often referred to as landscape evolution models -can simulate the formation of three-dimensional landforms that arise from the interaction of multiple processes, and in principle comparative testing ought to be straightforward (Hancock et al., 2010).Yet the algorithms behind these numerical models commonly differ from one another in multiple ways, which makes one-to-one comparison difficult.For example, if two numerical model codes differ simultaneously in their treatments of hydrology, sediment transport, and material properties, diagnosing any differences in their performance requires disentangling each of these effects.Often research questions focus on a combination of geological process and boundary conditions; classic examples include morphologic dating of fault scarps and glacial moraines.In this way, boundary conditions become a core numerical model component.To help address this challenge, it would be useful to have a software framework in which an investigator could alter one "process ingredient" or boundary condition at a time and thereby conduct meaningful parameter studies, sensitivity analyses, calibrations, multimodel analyses, and comparisons with data.
Terrainbento is a Python-language software product designed to help meet this need.Terrainbento version 1.0 provides three resources for exploring alternative process formulations for landscape evolution.First, terrainbento 1.0 includes a collection of 28 distinct model programs for the long-term (order 10 4 -10 6 years) evolution of drainage basin topography; most of these programs vary from a simple "base" numerical model in just one or two particular process descriptions.We consider a systematic subset of the 2048 possible numerical models derived from 11 binary options (Sect.3.9 describes the basis for the subset).
Second, terrainbento takes advantage of Python class inheritance such that all common features of terrainbento model programs (such as input/output, and the handling of boundary conditions) are provided in a generic "Erosion-Model" base class from which specific programs are derived.This ErosionModel class enables modelers to craft and apply their own implementations without needing to reinvent the overarching software framework or the necessary utility functions.Terrainbento 1.0 builds on the Landlab toolkit (Hobley et al., 2017), using Landlab components to represent individual hillslope, hydrologic, and channel process components and taking advantage of Landlab to handle common tasks such as input and output management.
Finally, boundary conditions can have a profound impact on model behavior.Terrainbento has a set of extensible features called "boundary handlers" that can be used to implement many common and complex boundary conditions.
Earth's landscapes are incredibly diverse, and the scientific questions that they pose are equally extensive.No one model program, or even a general framework like terrainbento, can hope to encompass all of this diversity.Terrainbento 1.0 was originally created to address landscape evolution in a humidtemperate, soil-mantled, postglacial environment with moderate relief (order 10 2 m, on a horizontal scale of order 10 4 m) and relatively rapid erosion rates (10 −4 to 10 −2 m yr −1 ) over a timescale of order 10 4 years.The choices of algorithms and process laws among the constituent numerical models reflect this motivation.Nonetheless, through the model template, terrainbento provides a sufficiently generic platform that it can be readily adapted to address a range of other scales and environments.This paper presents and describes terrainbento version 1.0, including its basic structure, mathematical underpinnings, software implementation, and the 28 constituent model programs.
2 General structure of a terrainbento model program A terrainbento model program begins with a gridded representation of topography.Landlab's RasterModelGrid, HexModelGrid, RadialModelGrid, and irregular Delaunay-Voronoi grids are all supported in terrainbento version 1.0.The elevation, and regolith thickness if present, at each grid node evolves according to a specified set of erosion and/or sediment transport laws, which vary between model programs.Terrainbento is agnostic to time and space units but expects that a user will ensure consistency between the grid units, time step units, and model input parameters.
In this section, we start by outlining the governing equations in a generic form.We then examine the software framework that implements elements common to all terrainbento model programs.The subsequent section (Sect.3) then presents the collection of process laws and algorithms used to represent hillslope erosion, hydrology, water erosion, and material properties.Section 4 describes the handling of boundary conditions.The governing equations for all 28 model programs in terrainbento 1.0 are listed in Appendix B.

2.1 A note on terminology
The word "model" can have multiple meanings in scientific computing and indeed in science generally.Here we will use the term mathematical model to mean a set of governing equations, which in this case describe landscape evolution under a given set of assumed process dynamics, material properties, and boundary conditions.Under this definition, two mathematical models may have governing equations that are structurally quite similar, but which are nonetheless considered to be distinct models either because certain constants take on different values or because a term is included in one version but not the other.For example, as described below, water erosion is commonly treated as proportional to either hydraulic power or hydraulic stress.We consider these to be distinct mathematical models, despite the fact that the difference lies only in the choice of two exponent values in the governing equation.
Each mathematical model contains terms that represent individual processes (or closely related collections of processes), such as erosion by surface-water flow.The mathematical representation for an individual process will be referred to as a process law or rate law.By this definition, a mathematical model in terrainbento consists of a set of process laws embedded within an overall mass-conservation equation.
The term numerical model is used here to refer to a numerical algorithm that solves a particular mathematical model by marching forward in time from a given initial condition under given boundary conditions.The term model program will refer to a set of source code files that performs the calculations needed to implement a numerical model.In some cases in terrainbento 1.0, a single model program can be configured to implement many numerical models, depending on its input parameters.One can consider alternative models that require only a different set of model program parameter values alternative parametric models, whilst models that require a different program are different structural models.The combination of a model program plus the inputs that control this type of choice will be referred to as a model configuration.

Basic ingredients and governing equation
Topography in a terrainbento mathematical model is represented as a two-dimensional field of elevation values, η(x, y, t).The general governing equation describes the rate of change of η as the sum of two terms: one representing erosion (or deposition) of mass by water-driven processes and one representing gravitational ("hillslope") processes: where E W is the rate of erosion (or deposition, if negative) by water-driven processes such as channelized flow, and E H is the rate for gravitationally driven processes such as soil creep and shallow landsliding (the subscript H stands for "hillslope," recognizing that gravitational processes will tend to be most important on hillslopes).Water erosion is assumed to depend on local slope gradient, S, water discharge, Q (which in many terrainbento model programs will be treated using drainage area as a surrogate, as discussed below and in Sect.3.6), and material properties.Erosion or accumulation by gravitational processes is assumed to be a function of gradient, material properties, and (in some model programs) soil thickness.While many terrainbento model programs only treat the evolution of the topographic surface η, two types of model program treat a more complex set of layers.The first type adds an explicit mobile-regolith layer on top of bedrock.The second type treats two different bedrock lithologies in a userdetermined spatial configuration.

Soil-tracking model programs
Several of terrainbento's model programs explicitly track a layer of regolith, defined here as unconsolidated and potentially mobile sediment, such as soil or alluvium (see Sect. 3.7.1).Here, for simplicity we will refer to this material as soil, keeping in mind that our operational definition is more general than the one commonly used by soil scientists (Soil Survey Staff, 1999).The land surface height is the sum of bedrock elevation, η b , and soil thickness, H : (2) Here too the term "bedrock" is used in its broadest possible sense and may include, for example, cohesive sedimentary material such as glacial till or saprolite.The time rate of change of soil thickness is the difference between the rate of soil production and erosion, where P is the rate of soil production from bedrock, and E WHS denotes the total rate of soil erosion (or accumulation, if negative) resulting from water-driven and gravity-driven transport processes.Similarly, the rate of change of bedrock surface height is the sum of the soil production rate (scaled by any density contrast between rock and soil) and the rate of bedrock incision by running water, E WR : In model programs that honor both a soil layer and two different lithologies, the surface elevation is in which case the height of the bedrock surface is Where the top layer exists, it lowers as a result of water erosion and (if soil is tracked) rock-to-soil conversion.This can be expressed mathematically as where δ L is a spatially varying function equal to 1 where L 1 > 0 and 0 elsewhere (here P is considered to be zero in non-soil-tracking model programs).The rate of change of elevation of the top of L 2 is given by which means that the lower layer L 2 is vulnerable to erosion and weathering wherever the top layer is missing (for example, having been eroded through).Note that the BasicHySa model program allows for the simultaneous water erosion of soil and rock, as discussed below.

Process formulations
Each of the 28 model programs in the terrainbento 1.0 collection has four elements, reflecting the treatment of hillslope processes, surface-water hydrology, erosion by running water, and material properties.The possible formulations for each of these elements are constructed around a set of binary choices.Each choice represents a decision about how a particular element might be formulated.For example, the downhill soil transport rate could be represented as either a linear or nonlinear function of local topographic gradient, while the lithology could either be treated as uniform or divided into two distinct types as discussed in Sect.2.2.2.The binarychoice design makes it possible to test the behavior of one alternative process element at a time.The binary options that form the basis for the terrainbento 1.0 constituent model programs are listed in Table 1.In Table 1, option B in each row usually represents a more sophisticated choice than option A: one that may bring more realism, but generally involves more parameters.
Each of terrainbento's model programs uses Landlab components to implement the numerical algorithms behind channel erosion, hillslope processes, and water-flow routing.
The components used are briefly identified by name in the following descriptions.The software architecture that supports this component-based approach is then discussed further in Sect. 5. Further information about Landlab and its component-modeling capability is provided by Hobley et al. (2017).

Model domain options
Terrainbento supports all types of Landlab model grids.Many options for the creation of synthetic topography and instantiation from a DEM are available.These options are described in the user manual.

Drainage area, flow direction, and flow accumulation
All terrainbento model programs calculate drainage area and surface-water discharge using the Landlab FlowDirectors and FlowAccumulator.Flow direction algorithms presently supported in Landlab include SteepestDescent/D4, D8, D ∞ , and Multiple Flow Direction (O'Callaghan and Mark, 1984;Freeman, 1991;Tarboton, 1997).Water routing across closed depressions is optionally handled using a lake-fill algorithm implemented by the Landlab DepressionFinderAndRouter component (the current version of which uses an implementation based on Tucker et al., 2001).Once flow directions and surface-water runoff are calculated, the contributing drainage area or surface-water discharge at a given grid cell i is calculated by adding up the area of all cells whose flow eventually passes through i, plus the area or discharge of i itself using the Landlab FlowAccumulator component.

Basic model program
The simplest of the component model programs in terrainbento is known as the Basic model program.It implements a discretized numerical solution to the following governing equation for land surface elevation η(x, y, t): where K is an erosion coefficient with dimensions of [L 1−3m T m−1 ], Q is the surface-water discharge, S is the gradient in the steepest downslope direction [-], m and n are the discharge and slope exponents, respectively, and Here the first term on the right-hand side represents channel and gully erosion, while the second term represents erosion or deposition by gravitational sediment movement.An example of a landscape simulated using the terrainbento Basic model program, using the common choice m = 1/2 and n = 1, is shown in Fig. 1a.Each of the five panels in Fig. 1 illustrates different terrainbento model configurations run with the same boundary conditions.The model runs that produced all of these example landscapes, including the input that specifies the  model run and slope-area diagrams like that shown in Fig. 2, can be found in the Jupyter notebook tutorials on GitHub (Barnhart et al., 2019b).
The second term on the right is the popular linear diffusion rule for hillslopes (Culling, 1963).The first term on the right represents channel incision and is based on the widely used stream-power formulation (Howard et al., 1994;Whipple and Tucker, 1999), in which the long-term average rate of channel downcutting is taken to be proportional to hydraulic power per unit bed area.A key assumption behind the Basic model program is that the erosion rate is limited by the capacity to detach and remove material, rather than by along-stream variations in the capacity to transport sediment.A common simplification, discussed further in Sect.3.6, is that drainage area, A, appears as a surrogate for effective water discharge.Often the choice m = 1/2 is made, reflecting the assumption that discharge per unit channel width scales as the square root of drainage area (Wohl and David, 2008).Similarly, n is commonly taken to be unity based on the derivation from stream power.The examples presented below use the exponent values m = 1/2 and n = 1 unless otherwise noted.
Rather than hard coding values for the water-erosion discharge exponent m and slope exponent n, terrainbento model programs permit varying these two exponents as parameters.Thus, within terrainbento 1.0 the same model program can be configured to represent either a stream-power or shear-stress representation of water erosion.Similarly, alternative forms of the fluvial erosion process law can be implemented by changing the value of m or n to reflect (for example) a different channel width-discharge scaling relation (Snyder et al., 2003b;Wohl and David, 2008) or different fluvial erosion mechanisms (Whipple et al., 2000a).For example, Fig. 1b shows an example of an alternative parametric model that uses the Basic model program with a value of m = 1/4.
Although Eq. ( 10) is rather simple, having just four parameters (K, D, m, and n), it represents a formulation that has been widely used in geomorphic modeling (e.g., Miller and Slingerland, 2006;Miller et al., 2007;Perron et al., 2009;Pelletier, 2010;Duvall and Tucker, 2015).The equations are commonly solved numerically on a regular or irregular grid.Discharge and drainage area are normally evaluated using a downslope routing algorithm in which the water output from one grid cell is passed to one or more downhill neighboring cells (see, for example, review in Tucker and Hancock, 2010).Despite its simplicity, this two-parameter numerical model has been shown to reproduce first-order properties of drainage basin topography, including dendritic drainage networks, concave-upward channel longitudinal profiles, and convex-upward hillslopes.
One arrives at the terrainbento Basic model program by choosing option A for each item in Table 1.In the following subsections, we review the various options that terrainbento offers for alternative treatment of hillslope pro-cesses, surface-water hydrology, channel incision, materials, and boundary conditions.

Hillslope processes
To simulate hillslope evolution processes in a soil-mantled landscape, we use components of varying complexity that treat soil transport as a diffusion-like process in which sediment flux is governed by topographic gradient.Terrainbento offers two alternative soil flux rules with which to simulate the downslope transport of soil and its dependence on topographic gradient: linear and nonlinear.
In addition, as discussed previously, terrainbento also allows for the option of explicitly tracking a dynamic soil layer.Our representation and treatment of soil evolution are simpler than in existing models focused on soil production and evolution (e.g., Cohen et al., 2010;Vanwalleghem et al., 2013).The dynamic soil option is provided to address the possibility that soil may become thin enough to limit flux, and this limitation may in turn influence the rate and style of landscape evolution.Inclusion of a dynamic soil layer requires an equation for soil production from the underlying lithology (P in Eq. 3) and, furthermore, that the flux law be modified to account for the local soil thickness such that soil flux goes smoothly to zero as thickness declines.

Continuity law for soil creep
The simplest forms of the so-called "geomorphic diffusion" equation (Dietrich et al., 2003) assume transport-limited conditions, in which the production rate of soil is always much greater than the transport rate; thus, the transport rate does not depend on soil availability or thickness.In this case, the hillslope term in the continuity Eq. ( 1) is where q h is the hillslope soil volume flux per unit width, φ is the porosity of the soil, and the ∇ operator represents differentiation in both horizontal directions (∇ = ∂/∂x + ∂/∂y).

Linear creep law
A variety of formulae exist for the soil flux, q h .The simplest and most common formula (Culling, 1963) treats the soil transport rate as a simple linear function of topographic gradient using a transport efficiency constant, D : where ∇η is the slope gradient.Using this flux rule with Eq. ( 11), the hillslope term in the continuity equation becomes where D, sometimes referred to as hillslope diffusivity, is equivalent to D /(1 − φ) and has dimensions of L 2 /T .This  simplest form of the evolution equation for soil creep on hillslopes results in convex-upward topography at steady state.

Nonlinear creep law
A more complex version of the creep law for soil-mantled slopes involves a nonlinear relationship between soil flux and topographic gradient.The nonlinear formulation captures accelerated creep and shallow landsliding as the gradient approaches an effective angle of repose for loose granular material.Several nonlinear creep-transport laws have been suggested in the literature.The most popular of these is the Andrews-Bucknam equation (Andrews and Bucknam, 1987), which performs reasonably well when compared with experimental and field data (Roering et al., 1999(Roering et al., , 2001;;Roering, 2008).One problem with the Andrews-Bucknam law, however, is that the flux diverges when the slope gradient, S, equals the threshold gradient S c and is undefined for S > S c .This property makes it challenging to incorporate into a land-  cazeau, 2005).Terrainbento uses a truncated Taylor series formulation for soil flux, which was derived by Ganti et al. (2012) for the Andrews-Bucknam law.The flux is given by where S = −∇η is topographic gradient (positive downhill), D is the transport efficiency factor, and S c is a critical gradient.The user specifies the number of terms N to be used in the approximation.The nonlinear flux rule results in convexup topography for shallow slopes and transitions to linear hillslopes for steeper slopes.An example terrainbento simulation using the nonlinear creep law is shown in Fig. 1c.

Linear depth-dependent creep law
For model programs that explicitly track a soil layer H (x, y, t), one needs to modify the creep law to incorporate a relationship between flux, q h , and local soil thickness.
Terrainbento uses an approach proposed by Johnstone and Hilley (2015), in which the flux decays exponentially as soil thickness approaches zero: where H 0 represents the soil thickness for which q h shrinks to (1 − 1/e) of its maximum value for a given slope gradient.
(Note that in the original formulation of Johnstone and Hilley (2015), D is treated as the product of H 0 and a transport coefficient with dimensions of length per time; here we lump them together as D.)

Nonlinear depth-dependent creep law
We can modify the nonlinear flux rule (Eq.14) to accommodate soil, again assuming an exponential velocity distribution in the subsurface (Johnstone and Hilley, 2015): This approach is somewhat similar to that used by Roering (2008) in a study that compared the predictions of a nonlinear, depth-dependent flux law with observed hillslope forms.

Soil production
Models that track a layer of soil must include an expression to specify the rate at which soil is produced from the under-lying parent material.The most commonly applied formula, and the one used by terrainbento's soil-tracking model programs, treats the rate of soil production from the underlying lithology as an inverse-exponential function of soil thickness (Ahnert, 1976;Heimsath et al., 1997;Small et al., 1999): where P 0 is the maximum production rate (with dimensions of length per time), and H s is a depth-decay constant on the order of decimeters.

Hydrology
Treatments of surface-water hydrology in landscape evolution models are commonly quite straightforward, reflecting the need for both simplicity and computational efficiency.
Erosion formulae normally require specification of water discharge or (less commonly) depth.The most common parameterization is to use contributing drainage area, A, as a surrogate for surface-flow discharge, Q.This effectively states that Q = rA, where r, a runoff rate per unit area, is equal to 1.This is the default option in terrainbento's model programs.Operationally, this means that the water-erosion law includes A as the effective discharge (see Sect. 3.6 below) and that the erosion law parameters embed information about climatic factors such as precipitation frequency and intensity, as well as material properties such as soil infiltration capacity (e.g., Tucker, 2004).

Variable source area hydrology
In vegetated, humid-temperate regions, storm runoff is commonly produced by the saturation-excess mechanism, in which rain falls on areas that have become saturated (Dunne and Black, 1970).Such areas tend to occur in locations with either gentle topography, large contributing area, or both.Because the source area for runoff generation is both limited in spatial extent and varies over time, the phenomenon has come to be known as variable source area hydrology (VSA).
Previous modeling studies have suggested that VSA can impact long-term landform evolution, as steeper upland areas tend to experience less intense and/or less frequent erosion and sediment transport by runoff (Ijjasz-Vasquez et al., 1993;Tucker and Bras, 1998).For this reason, terrainbento 1.0 includes a set of model programs that provide a relatively simple treatment of VSA.This treatment is based on the approach of O'Loughlin (1986) and Dietrich et al. (1993) and is similar to the TOPMODEL concept of Beven and Kirkby (1979).Each element on the landscape is considered to have an upper permeable soil layer of thickness H and saturated hydraulic conductivity K sat .The soil layer is assumed to overlie relatively impermeable material.From Darcy's law, the maximum shallow subsurface flow discharge when the soil is fully saturated is the product of conductivity, depth, and local hydraulic gradient, which is assumed to be equal to topographic gradient, S. The maximum subsurface discharge per unit contour width is therefore given by where T = K sat H is the soil transmissivity.Next, we consider a recharge rate, R, which represents the average rate of water input per unit area (dimensions of length per time).The total unit discharge is the product of recharge and drainage area per unit contour length, a: Using these two principles, the surface-water unit discharge, q, at any location is This threshold-based approach has been used, for example, in numerical models that explore how hillslope hydrology influences landform evolution (Ijjasz-Vasquez et al., 1993;Tucker and Bras, 1998).One drawback, however, is that the use of mathematical thresholds in numerical models can complicate the calibration process by creating "numerical daemons": sharp discontinuities in a numerical model's response surface (i.e., the N p -dimensional surface that describes a simulation output quantity as a function of its N p input parameters) (e.g., Kavetski and Kuczera, 2007;Hill et al., 2016).In this particular case, we can create a smoothed version of Eq. ( 20) without any loss of realism by positing that within any given patch of land there is actually a distribution of effective recharge rates.The use of a smooth version in this formulation has two benefits -it is physically more realistic and is less prone to numerical discontinuities.The simplest strictly positive probability distribution is an exponential function where p(R) is the probability density function of R, and R m is the mean recharge rate.The mean surface-water unit discharge can then be found by integrating as follows: where R c = T S/a is the minimum recharge needed to produce surface runoff.
It is useful to recast this in terms of an effective contributing area, A eff , defined as where x represents flow width (in a gridded digital elevation model, it would be natural to use cell width).Note that we have assumed that the contour width is approximated by the flow width.By this definition, the effective drainage area is always less than or equal to the actual drainage area, reflecting the fact that some of the water runs through the shallow subsurface rather than across the surface as overland (or channelized) flow.Where slope gradient is small or drainage area is large, the effective area approaches the actual area.If the surface is flat (S = 0), the exponential factor equals unity and A eff = A, reflecting the fact that no water can be conveyed by shallow subsurface flow.At locations where the gradient declines in the downslope direction, water previously carried by subsurface flow effectively "exfiltrates" from subsurface flow within the soil layer and is carried as surface-water discharge.Conversely, where S is large and/or A is small -as might be the case in steep headwater areas -the effective drainage area becomes much smaller than the actual area, indicating that most of the incoming water is traveling beneath the surface rather than contributing to overland flow.A final step is to note that one can collapse the various factors in Eq. ( 23) into a single parameter, α = T x/R m , the saturation area scale (dimensions of length squared).A high value of α represents soils that have a large capacity to carry subsurface flow relative to the recharge rate; a low value reflects a more limited subsurface flow capacity.
Seven of terrainbento's model programs implement variable source area hydrology by using A eff , as defined in Eq. ( 23), in place of drainage area, A (Table 2).One of these (BasicSaVs) also explicitly tracks a soil layer, and the timeand space-varying thickness of this soil layer, H (x, y, t), is used to calculate T = K sat H (x, y, t).An eighth model program (BasicStVs) also uses a stochastic treatment of precipitation; in this model program, the randomly generated precipitation rate p is used for R m in Eq. ( 23).
An example simulation with a terrainbento model program (BasicVs) that includes a variable source area component is shown in Fig. 1d.The only difference in formulation between this example and the Basic model program illustrated in Fig. 1a is that BasicVs calculates channel erosion using effective drainage area, A eff , as defined in Eq. ( 23), in place of total drainage area.The result is a drainage network bounded by steep, convex-upward ridges.These ridges are sufficiently steep that A eff A, so their erosion is dominated by soil creep.The bases of the hills represent locations where water emerges from the shallow subsurface to become surface flow that feeds the channel network.

Stochastic precipitation and runoff
Many landscape evolution models use an effective discharge approach, in which a single value of precipitation or runoff (either given explicitly or embedded in a lumped rate coefficient) is used as a surrogate for the full range of runoffproducing events (e.g., Willgoose et al., 1991;Kooi and Beaumont, 1994;Tucker and Slingerland, 1997) proach has the advantages of simplicity and computational efficiency, but also has limitations.For example, the appropriate effective discharge may vary in space and time (Huang and Niemann, 2006).One solution is to use a stochastic treatment of precipitation and/or discharge, in which events are drawn from a specified probability distribution (Tucker and Bras, 2000;Snyder et al., 2003a;Tucker, 2004;Lague et al., 2005).
In order to facilitate comparison between numerical models with deterministic and stochastic treatments of water discharge, terrainbento 1.0 includes a set of six model programs that each implement two stochastic precipitation algorithms available in the PrecipitationDistribution Landlab component.The aim of these algorithms is not to reproduce individual storm events, but rather to capture a spectrum of runoff and streamflow events of varying frequency and magnitude.The first of these two methods is a stochastic-in-time approach based on Tucker and Bras (2000).The second option uses deterministic time steps but stochastic precipitation intensity.
In the first option, a series of "storms" is generated based on a specified mean storm duration T r , mean inter-storm duration T b , and mean storm depth h r .The mean storm and inter-storm durations are generated from exponential distributions, after Eagleson (1978).For each individual storm, the mean storm depth is generated from a gamma distribution.The gamma shape parameter used to draw a random storm depth is equal to that specific storm's duration divided by the mean storm duration.The scale parameter is equal to the mean storm depth.The depth and duration of an individual storm are then used to calculate the rainfall intensity (Ivanov et al., 2007).
In the second option, in which the time step duration is fixed, the frequency of occurrence of rainfall is described using an intermittency factor, F , which is defined as the fraction of time rain occurs rain, and a mean event precipitation rate, p d .
Thus, the mean precipitation rate (averaged over wet and dry periods), p ma , is given as The probability distribution of precipitation rate, p, is simulated using a stretched exponential survival function, where c is a shape parameter and λ is a scale parameter.Use of the stretched exponential function is based on Rossi et al. (2016), who found that the function provides a good approximation for daily rainfall distributions in the continental US and Puerto Rico.While Rossi et al. (2016) found this distribution to be appropriate for mean daily rainfall, note that terrainbento is agnostic to the time units chosen by a user.Wilson and Toumi (2005) argued that theoretical considera-tions suggest c ≈ 2/3, while Rossi et al. (2016) found a mean value of c = 0.74 for weather stations in the continental US.
The shape parameter λ associated with a mean daily precipitation rate p d and shape factor c is given by where is the gamma distribution function.
To describe the frequency-magnitude spectrum probabilistically in terrainbento's stochastic model programs, time is discretized into a series of steps of duration δt/n ts , where δt is the "global" time step used for all other process components.During each step, an "event" with precipitation rate p is drawn at random from the cumulative distribution in Eq. ( 25).One of two approaches is then used to calculate the corresponding runoff rate, r.The first approach, which is the default used in five of the six stochastic model programs, assumes a mean soil infiltration capacity I m .The rate of runoff is calculated as This formulation is a smoothed version of the simple threshold approach r = max(p − I m , 0), which has been used in prior studies to represent infiltration-excess overland flow generation (e.g., Tucker and Bras, 2000).The smoothed version avoids the sharp discontinuity at p = I m and is arguably more realistic as it honors natural variability in soil infiltration capacity.The runoff rate approaches zero when p I m and approaches p when p I m .
The second approach uses the variable source area runoff generation formulation described in Sect.3.5.1 using p in place of recharge R m .This approach is used only in the model program BasicStVs (Table 2).

Water erosion
Several different expressions have been proposed as process formulations for long-term channel incision (and for erosion by surface water more generally).Terrainbento 1.0 was originally designed to address erosion of cohesive sediments (including glacial till) and clastic sedimentary rocks with a relatively high fracture density, both of which are prone to erosion by hydraulic detachment of sediment grains and fracture-bounded fragments ("plucking").This focus guided the choice of water-erosion laws in terrainbento 1.0.Each terrainbento model program uses one of two main types of erosion law: a simple area-slope detachment formula (sometimes referred to in the literature as the stream-power family of erosion laws (e.g.,Howard et al., 1994;Whipple and Tucker, 1999), and an erosion formula that accounts for sediment discharge, particle entrainment from the bed, and particle deposition onto the bed.Within these two broad categories, terrainbento model programs express several variations in form; for example, some include a threshold term, and in some of these the threshold increases with progressive incision depth.Each variation is presented and discussed in the sections below.Here, we start with a description of the simplest formulation, which serves as the default choice.
The area-slope (a.k.a., stream-power) family of erosion laws derives from the assumption that the erosion rate, E W , depends primarily on the hydraulic gradient, S, the water discharge, Q, and the channel width W : where k 1 is a coefficient that depends on material properties, channel geometry, and other factors, and ω c is a threshold below which no erosion occurs (in practice, the threshold is often assumed negligible, or its effects are taken to be subsumed in the exponents).The exponents µ and ν reflect the nature of the erosional processes; for example, Whipple et al. (2000a) argued that different values may be appropriate for abrasion-dominated and for plucking-dominated systems.
The discharge exponent µ also embeds information about channel geometry.Often, drainage area A is used as a surrogate for discharge.One limitation of Eq. ( 28) is that it does not allow for sediment deposition; for this reason, it is sometimes referred to as a detachment-limited law (a term first coined by Howard, 1994), reflecting the assumption that the rate of downcutting is limited by the rate at which material can be detached and removed.Despite the simplicity of Eq. ( 28), its various permutations have shown reasonable success when tested against field observations (Stock and Montgomery, 1999;Whipple et al., 2000b;Kirby and Whipple, 2001;Snyder et al., 2000;Lavé and Avouac, 2001;Tomkin et al., 2003;van der Beek and Bishhop, 2003;Duvall et al., 2004;Loget et al., 2006;Whittaker et al., 2007;Attal et al., 2008Attal et al., , 2011;;Yanites et al., 2010;Hobley et al., 2011;Gran et al., 2013).Landscape evolution models that use the generic stream-power approach are able to reproduce basic properties of erosional landscapes, such as dendritic channel networks with concave-upward longitudinal profiles (e.g., Howard, 1994;Whipple and Tucker, 1999;Tucker and Whipple, 2002).
One of the most commonly used versions of Eq. ( 28) is obtained by making the following assumptions: (1) channel width can be represented as a power function of drainage area, and (2) the erosion threshold is negligible.Under these conditions, the erosion law becomes where K is a coefficient that includes information about precipitation and hydrology as well as material properties and channel geometry.The change from ν to n reflects the possibility that slope-dependent information is included in K.This is the fluvial erosion law introduced earlier in the discussion of the Basic model program.If one also assumes that Q = rA, with r representing runoff rate, then Eq. ( 29) can be cast in the commonly used form where K = Kr m .With respect to the exponents m and n, if one assumes that (1) the rate of downcutting depends on stream power per unit surface area, (2) effective discharge is proportional to drainage area, and (3) channel width is proportional to the square root of discharge, then the exponent values become m = 1/2 and n = 1.
The simplicity of Eq. ( 29) -it has only one parameter if m and n are assumed to be accurate representations of process -together with its ability to reproduce common features of drainage basins and networks have led to its widespread use in landscape evolution studies, especially with the exponent values m = 1/2 and n = 1 (e.g., Duvall and Tucker, 2015).One might think of this particular configuration as the "model to beat": to justify a more complex formulation, one would ideally need to demonstrate that such a formulation performs distinctly better.Equation ( 29), which we will refer to as the simple streampower law, forms the default choice for water erosion in terrainbento's model programs.In the following subsections, we describe the variations and alternatives to simple unit stream power among the terrainbento 1.0 model programs.The complete governing equations for each of the terrainbento 1.0 model programs are given in Appendix A.

Erosion threshold
Bed-load sediment transport is well known to exhibit threshold-like behavior, in which the transport rate is negligible until a certain minimum hydraulic tractive stress is reached, at which point significant transport begins.Similar behavior applies to the erosion of highly cohesive sediment (e.g., Julien, 2010) and presumably also to bedrock (though the values of the operative thresholds for bedrock are not known).For this reason, models of landscape or longitudinal channel profile evolution often include a threshold term below which no erosion takes place.
Several terrainbento model programs include a threshold in the water-erosion law.In order to promote mathematically smooth behavior, acknowledge evidence for distributions of transport thresholds (Kirchner et al., 1990;Wilcock and McArdell, 1997;McEwan and Heald, 2001), and avoid numerical daemons associated with threshold-type equations (e.g., Kavetski and Kuczera, 2007), the basic thresholded erosion law in terrainbento uses an exponential smoothing function following Shobe et al. (2017).Terrainbento's thresholded erosion laws take the form Here ω represents the erosion rate that would occur in the absence of a threshold and is a function of slope gradient  and either drainage area or discharge.For example, for those model programs that add a threshold term to the area-slope erosion in Eq. ( 29), ω is defined as The factor ω c is a threshold with dimensions of length per time.The functional form of the smooth-threshold erosion function (Eq.31) is illustrated in Fig. 3.A constant threshold term is included in the water-erosion laws for five of terrainbento's constituent model programs (Table 2).Several others use a space-and time-varying threshold, as we describe next.

Incision depth-dependent erosion threshold
In a study of river incision into glacial deposits following ice recession in the US upper midwest, Gran et al. (2013) found evidence for an erosion threshold that increased with progressive incision depth.They attributed this to a downstream increase in median grain diameter resulting from enrichment of coarse gravel in bed material as the channel cuts through glacial deposits and the valley widens.In comparing alternative long-profile-evolution numerical models with the observed profile, they found that the best match was achieved when the erosion threshold was allowed to increase linearly as a function of cumulative incision depth.Inspired by the findings of Gran et al. (2013), terrainbento 1.0 includes the option to allow the erosion threshold ω ct to increase with erosion depth according to where D I is the cumulative incision depth at location (x, y) and time t, ω c is the threshold when no incision has taken place yet, and b (with dimensions of inverse time) sets the rate at which the threshold increases with progressive incision depth.As before, an exponential term is used to smooth the threshold such that the water-erosion rate approaches zero when ω ω c and asymptotes to ω − ω c when ω ω c (Fig. 3).The max function is included to prevent the threshold from decreasing in locations where hillslope processes produce net deposition (i.e., negative incision).

Shear-stress erosion law
Two important and commonly used measures of the erosional potential of streamflow are unit stream power and shear stress.The first represents the rate of energy dissipation per unit surface area, while the second represents the hydraulic traction force per unit area.Erosion rates in cohesive or rocky material tend to correlate strongly with both quantities (e.g.,Howard and Kerby, 1983;Whipple et al., 2000b), and both are widely used as the basis for long-term erosion laws.To support studies that compare and test these two approaches, terrainbento 1.0 allows one to configure the erosion law to represent bed shear stress rather than unit stream power.This is accomplished simply by changing the exponents on discharge (or drainage area) and channel gradient in Eq. ( 28).If one uses the Manning equation to describe channel roughness and assumes that channel width is proportional to the square root of discharge, the applicable exponent values are m = 3/5 and n = 7/10 (Howard and Kerby, 1983;Howard, 1994).Use of the Darcy-Weisbach roughness law leads to slightly different values of m = 1/3 and n = 2/3 (Tucker and Slingerland, 1997), which we use in the examples that accompany the terrainbento 1.0 documentation.In terrainbento 1.0, the choice of exponent values is set using an input file or keyword arguments, so separate code is not needed to implement the shear-stress option.Nonetheless, we consider the stream-power and shear-stress formulations to form distinct parametric models.

Sediment-tracking entrainment-deposition hybrid model program
The sediment-tracking model program computes changes in riverbed elevation resulting from competition between the entrainment of bed material into the water column and deposition from the water column onto the bed using a combined entrainment-deposition law discussed by Davy and Lague (2009).The governing equations, derived from mass balance, state that changes in channel bed elevation η over time are driven by bed material erosion E and bed material deposition D s : where E and D s are volume fluxes of bed material per unit bed area representing entrainment from the bed and deposition onto the bed, respectively, and φ is the porosity of bed material.Note that here our equation is different from Eq. ( 3) of Davy and Lague (2009) in that 1 − φ is not in the denominator of E. This discrepancy is due to a difference in whether bulk or sediment density is used to convert between mass and volume for E. Equation ( 34) is coupled with the conservation of sediment concentration c s in a water column of depth h: where x represents distance along the path of flow.The above states that sediment in the water column involves a balance between erosion, deposition, and the streamwise spatial gradient in fluvial sediment flux per unit width, q s .Again following Davy and Lague (2009), we assume that the time rate of change of sediment in the water column is negligible (as it is meant to represent an average over time) so that In other words, the sediment flux at a particular downstream point x is the integral of all the erosion minus deposition that has taken place upstream.
The erosion flux E may be written in a number of ways, but in general depends on water discharge Q (or drainage area as a proxy), bed slope S, and some parameter or set of parameters describing the erodibility of the channel bed.As with other terrainbento model programs, we treat m and n, the exponents of Q and S, as parameters.The entrainment term may also include a threshold, and that threshold may be constant or may vary with incision depth or with lithology.
Sediment deposition flux D s is a function of the concentration of sediment in the water column and the effective settling velocity V of the sediment particles.Adding that c s is the volumetric sediment flux divided by the volumetric water flux, the deposition flux may be written where Q is volumetric water discharge and Q s is volumetric sediment discharge (equal to q s times flow width).Importantly, V is the net settling velocity after accounting for upward-directed turbulence and sediment concentration gradients in the water column.Davy and Lague (2009) separate the latter effects into a dimensionless parameter d * such that D s = d * V Q s Q , but here for simplicity we combine both effects into an effective settling velocity V .
The entrainment-deposition model program provides greater flexibility than detachment-limited model programs in that it can freely transition between detachment-limited and transport-limited behavior, depending on the relative importance of the erosion and deposition fluxes.If the deposition flux is negligible relative to the erosion flux, behavior becomes detachment limited.In the opposite case, the numerical model becomes transport limited.The entrainmentdeposition model program is therefore uniquely able to treat landscapes that may exhibit both types of behavior at different points in space and time, at the cost of only a single extra parameter (V ) relative to basic stream-power-type model programs.For a full description of the entrainmentdeposition formulation and its implications, see Davy and Lague (2009).

Entrainment-deposition hybrid model program with fine sediment
In the entrainment-deposition approach proposed by Davy and Lague (2009), all material eroded from the channel bed is included in sediment flux and deposition calculations.While this fully mass-conservative approach is a useful general case, it neglects the fact that clay-and silt-sized sediment may have such a low settling velocity that they remain permanently suspended until and unless they enter a body of standing water.A simple modification to the entrainmentdeposition model program allows for the treatment of a scenario in which the finest fraction of eroded sediment is permanently removed from the landscape upon entrainment.In the general case, the change in Q s along the river is written where dx f is the width of flow.To account for permanently suspended fine sediment, represented as a fraction of total bed sediment F f , we simply exclude the fine sediment from the sediment flux and write such that the material incorporated into the sediment flux is reduced in proportion to the amount of fine sediment on the bed.This approach is simple and efficient, but would likely be limited in settings with very high proportions of fine sediment, as large concentrations of even very fine grains in the water column may inhibit further sediment entrainment (Davy and Lague, 2009).

Entrainment-deposition model program with bedrock and alluvium
One weakness of the entrainment-deposition model program described above is its limitation to a single type of bed material.For example, one can configure the parameters to represent erodible material, such as loose sediment, or resistant material, such as indurated bedrock, but not both at once.This limitation means that the basic form of the entrainment-deposition model program cannot honor the reality that many bedrock-incising rivers are blanketed by alluvium, nor can it be used to assess the relative contributions of sediment entrainment and bedrock erosion to channel morphology and sediment flux.One potential solution is to use the entrainment-deposition model program in conjunc-tion with a substrate layering system (i.e., a layer of sediment overlying bedrock), in which each layer is defined by its own erodibility factor and erosion threshold (e.g., Gasparini et al., 2004;Carretier et al., 2016).However, such an approach does not allow for the simultaneous erosion of sediment and bedrock, which can occur in real rivers when the alluvial cover is spatially discontinuous and/or intermittent in time.Some recent numerical modeling approaches allow for a smooth transition between alluviated and barebedrock beds and the simultaneous evolution of the sediment and bedrock surfaces (Lague, 2010;Zhang et al., 2015;Shobe et al., 2017).Lague (2010) Changes in sediment thickness are treated identically to the entrainment-deposition formulation (Eq.34), and changes in bedrock height are driven by bedrock erosion E r (there is no deposition of bedrock): Erosion and deposition of sediment are computed using the same approach as used in the more basic entrainmentdeposition formulation, with the addition of a factor that limits the rate of sediment entrainment, E s , as sediment availability declines: where K s is an entrainment coefficient for alluvium.Here H * is the bedrock surface roughness length scale.Large H * corresponds to a rough bedrock surface and vice versa.
The SPACE numerical model includes a similar formulation for the bedrock, whereby bedrock erosion becomes more efficient as sediment thickness declines: Here, the r subscripts denote bedrock parameters.Adding bedrock erosion to the entrainment-deposition formulation requires that eroded bedrock material be added to sediment flux calculations: where A( x) represents drainage area, which increases as a function of streamwise distance x.The factor F f indicates the proportion of the bedrock that is made up of fine sediment that goes into permanent suspension when entrained and is no longer included in model program calculations.Q s therefore only includes grains not considered "fine".
As demonstrated by Shobe et al. (2017), SPACE is capable of transitioning between detachment-limited and transport-limited behavior.In a further advance over basic entrainment-deposition numerical models, SPACE can simulate bare-bedrock channels, fully alluvial channels, and mixed bedrock-alluvial channels, allowing the transition between these states to be set by sediment flux and erosive power.SPACE enables the simulation of channels that may alternate between bedrock, bedrock-alluvial, and alluvial states in response to changing tectonic forcing, climate, or sediment supply conditions.For a full derivation and discussion of SPACE, as well as a development of steady-state analytical solutions, see Shobe et al. (2017).

How alternative hydrology formulations influence terrainbento's erosion laws
For those model programs that use variable source area hydrology, the drainage area factor in the water-erosion law is replaced by effective drainage area, A eff , as defined by Eq. ( 23).Models that use stochastic hydrology replace A with Q = rA using r as defined in Eq. ( 27).One model program, BasicStVs, combines stochastic runoff generation with variable source area hydrology.With this model program, as in the variable source model program more generally, the capacity to carry subsurface discharge is defined as where as before T is transmissivity, S is surface gradient, and x is flow width.Assuming that interception loss and leakage to deeper groundwater are negligible, the total discharge produced by a storm event with rainfall rate p is The surface discharge, Q, should then be the difference between these two quantities or zero if Q ss > Q tot .However, a simple "either-or" differencing formulation is somewhat unrealistic (given small-scale natural variability in T ) and if implemented numerically would risk creating numerical daemons in the response surface.To avoid these issues, the BasicStVs model program uses the exponentially smoothed formula The form of this equation is similar to that of the smooth-threshold erosion law illustrated in Fig. 3. Substituting the definitions of Q tot and Q ss above, The precipitation rate calculated for each stochastic event is used to calculate Q, which is then used as the discharge factor in the erosion law E W = KQ m S n .

Soil and alluvium
One of the binary options listed in Table 1 is the ability to explicitly track a dynamic soil layer.Models that use this option implement the depth-dependent form of the applicable soil-creep law (i.e., either the linear or nonlinear form).
When the dynamic soil option is used in combination with a sediment-tracking entrainment-deposition erosion law (model program BasicHySa), the SPACE numerical model described above is used in place of the simpler (singlematerial-type) entrainment-deposition law.In all other cases, the use of a dynamic soil layer does not directly influence the water-erosion law.
When dynamic soil is combined with variable source area hydrology (model program BasicSaVs), the actual soil thickness at each point H (x, y, t) is used to calculate transmissivity.

Multiple lithologies
With two-lithology model programs, the material-dependent parameters in the water-erosion equation, including the coefficient (K) and, if applicable, the threshold (ω c ), vary in space and time as a function of the local surface elevation, η, in relation to the elevation of the contact between lithologies 1 and 2, η C (x, y).If η > η C , lithology 1 is exposed at the surface; otherwise, the surface unit is lithology 2.
To acknowledge the fact that lithological contacts are not razor thin and to preserve smoothness in the numerical solution, we allow there to be a finite "contact zone" within which the two lithologies are both considered to influence the material erodibility.One might imagine this zone as representing a gradational transition from one unit to another, or alternatively an uneven contact surface.We define a weight factor w that defines the relative influence of each of the two lithologies: Here, w represents the influence of lithology 1, and 1 − w describes the influence of lithology 2. At each location, the channel erosion rate coefficient is calculated by applying this weight factor.For example, in the model program BasicRt, which uses the simple unit stream-power formula, the rate coefficient K is calculated as where K 1 and K 2 are the rate coefficients associated with each lithology, and W c is the contact-zone width.

Variable climate
As a simpler representation of variable climate than available in the PrecipChanger described in Sect.4.2, one model program (BasicCc) provides the ability to change parameter K linearly through time.The representation of change is as follows.At the beginning of a simulation run, K is assumed to be larger or smaller than its final value (K 0 ) by a factor f ; if f > 1, K starts out larger than K 0 (representing a more erosive climate) and declines through time and conversely if f < 1.K stops changing after a time period T s , whereupon it assumes its final value K 0 .Mathematically, this linear variation in K is where κ = (1 − f )K 0 /T s is the rate of change.

Pairwise process combinations
As noted earlier, the various process options described above can be arranged into a set of 11 binary choices (Table 1).Terrainbento 1.0 is designed to support experimentation and hypothesis testing among these (and other) alternative formulations.The number of possible unique combinations among this set of 11 options is unwieldy (2 11 , though some are not physically sensible).In creating the individual terrainbento model configurations, we used an approach that focuses on single and pairwise variations on the Basic (simplest) model program, which is the first entry in Table 2.The next 11 entries are model programs or model configurations that differ from Basic in just one element.The remaining entries represent pairwise combinations.Not all possible pairwise combinations are included.Instead, the pairwise process combinations selected represent those for which we thought there might be nonlinear interactions between the two process elements -in other words, those combinations for which we expected the whole to be greater (or less) than the sum of the parts.An example of such a nonlinear interaction that has been explored in the literature is temporal variability in water discharge in a river system where the erosion process is strongly thresholded (Tucker and Bras, 2000;Snyder et al., 2003a;Lague et al., 2005;DiBiase and Whipple, 2011) The particular list of model program choices in Table 2 is not meant to be exhaustive.The terrainbento software was designed to be easily extensible as needed for any given application so that, for example, if a researcher wishes to explore combinations that are not included in the present collection of model programs or to add a new process formulation, he or she can do so with relative ease.In the next section, we describe how the software is designed to promote extensibility.

Boundary conditions
Just as process representation influences simulation results, so do boundary conditions.Representing boundary conditions in a component-like fashion permits systematic and reproducible changes in boundary conditions through either boundary condition component choice or parameter choice.To support alternative boundary conditions, terrainbento 1.0 includes five boundary condition handler classes.These boundary condition handlers are similar in construction to Landlab components: they are Python objects, and they must have an __init__ method that takes as a first argument a Landlab model grid and a run_one_step method that takes as its only argument the time step duration dt.Four of these classes are called base-level handlers, reflecting the fact that they modify the elevations on the boundaries of the simulated terrain.The final class is the PrecipChanger, a boundary condition handler that either modifies precipitation statistics or the value of K, K r , K s , K 1 , or K 2 , depending on the model program.

Base-level handlers
Each of the four base-level handlers modifies the elevations of specific grid nodes.Before describing these base-level handlers it is worth reviewing the boundary condition types available to Landlab model grid nodes (Hobley et al., 2017).A boundary node can be open (with either a fixed value or fixed gradient), looped, or closed.A boundary node need not live on the edge of a rectangular grid -for example, many nodes may be closed boundary nodes if the simulated domain is a single watershed.All nodes that are not boundary nodes are called "core nodes".
The four base-level handlers were designed to capture the most common cases for boundary conditions in Earth surface processes modeling.The SingleNodeBaselevelHandler controls the elevation of a single, open, fixed-value boundary node, meant to represent a watershed outlet.The outlet lowering rate is specified either as a constant or through a user-supplied text file that specifies the elevation change through time.The NotCoreNodeBaselevelHandler moves either the core nodes or the not-core nodes at a constant rate through time or based on a text file.The CaptureNodeBase-levelHandler was designed to simulate drainage basin capture by a basin external to the simulated domain.It changes the boundary condition status of a single node from closed to fixed-value open at a user-defined time and lowers its elevation.
The final base-level handler is the GenericFunctionBase-levelHandler.It is similar to the NotCoreBaselevelHandler in that it either moves the core nodes or the not-core nodes.However, instead of taking a constant rate or time-elevation pattern as input, it requires that a user define a function of two arguments that returns an at-node field of uplift rate.The two required arguments are the model grid and the simulation integration time.As the model grid contains attributes x_of_node and y_of_node, this boundary condition handler thus permits a user to define the relative uplift rate as any function of space and time.

PrecipChanger
The final boundary condition handler was designed to implement the impacts of changing climate on the precipitation distribution and, by extension, the erodibility of material by water.For model programs with stochastic precipitation and uniform time steps, this method modifies the intermittency factor and the mean rainfall rate, whereas for model programs with an effective discharge it modifies the erodibility by water.This boundary condition handler does not presently support stochastic precipitation with stochastic event durations.
For "St" model programs that explicitly represent the intermittency factor F and mean rainfall rate p d (Sect.3.5.2), the PrecipChanger must only modify those parameters.For model programs with effective discharge, deriving a relation between K (or K r , K s , K 1 , or K 2 ), p d , and F requires defining a mathematical model of the underlying hydrology and tracing how variations in precipitation parameters influence the long-term erosion rate.We start by noting that drainage area serves as a surrogate for discharge, Q.We can therefore write an instantaneous version of the erosion law in the Basic model programs as This formulation represents the erosion rate during a particular daily event with daily average discharge Q, as opposed to the long-term average rate of erosion, E. It uses an instantaneous erosion coefficient K i .We next assume that discharge is the product of runoff rate, r, and drainage area: Combining these we can write This equation establishes the dependence of short-term erosion rate on catchment-average runoff rate, r.
Next we need to relate runoff rate to precipitation rate.A common method is to acknowledge the existence of a soil infiltration capacity, I c , such that when p < I c , no runoff occurs, and when p > I c , An advantage of this simple approach is that I c can be measured directly or inferred from streamflow records.
To relate short-term ("instantaneous") erosion rate to the long-term average, one can first integrate the erosion rate over the full probability distribution of daily precipitation intensity.This operation yields the average erosion rate produced on wet days.To convert this into an average that includes dry days, we simply multiply the integral by the wetday fraction F .Thus, the long-term erosion rate by water can be expressed as where f (p) is the probability density function (PDF) of daily precipitation intensity.By equating the above definition of long-term erosion E with the simpler definition in Eq. ( 52), we can solve for the effective erosion coefficient, K: In this case, what is of interest is the change in K given some change in precipitation frequency distribution f (p).Suppose we have an original value of the effective erodibility coefficient, K 0 , and an original precipitation distribution, f 0 (p).Given a future change to a new precipitation distribution f (p), we wish to know what is the ratio of the new effective erodibility coefficient K to its original value.Using the definition of K above, the ratio of old to new coefficient is Thus, if we know the original and new precipitation distributions, we can determine the resulting change in K.
We use a Weibull distribution for the precipitation intensity PDF after Rossi et al. (2016), where λ is the distribution scale factor.Its relationship with p d is defined as The above definition can be substituted into the integrals in Eq. ( 58).We are not aware of a closed-form solution to the resulting integrals.Therefore, the erosion model programs used for projection apply a numerical integration to convert the input values of F , c, and p d (the last of which can change over time) into a corresponding new value of K.
5 Software implementation

Overview
In creating a software product that manifests not one but rather dozens of potential model configurations, efficiency and reuse are key design considerations.To meet this goal, terrainbento 1.0 uses an object-oriented approach to its highlevel design.Each terrainbento model program is implemented as a Python class.The class that implements any particular terrainbento model program inherits from a common base class called ErosionModel.Here we describe the main functions of the base class, the typical structure of the derived class, and the use of a driver program to configure and execute a terrainbento model program.

Terrainbento base classes
Terrainbento contains three base classes to minimize duplicate code and maximize the extensibility of the modeling framework.The first of these, ErosionModel, handles common methods such as instantiation, run, output creation, and model finalization.These include creating the model grid, reading initial topography from a file, creating synthetic topography, calculating elevation change, writing NetCDF and xarray datasets of simulation output, and interfacing with the boundary condition handlers.All model programs except the "St" and "Rt" series inherit directly from the ErosionModel base class.
The stochastic and two-lithology model programs each have a sufficient number of specialized methods to justify having their own base classes, which are the StochasticEro-sionModel and TwoLithologyErosionModel, respectively.Both of these inherit from ErosionModel.The StochasticEro-sionModel handles setting up the stochastic rain generator; calculating precipitation, runoff, and water erosion; and keeping records of storm sequences.The TwoLithologyEro-sionModel handles setting up the lithology contact elevation and updating any fields that depend on the depth to the contact.

Basic model interface
The Community Surface Dynamics Modeling System (CS-DMS) has promoted the use of an interface standard known as the basic model interface (BMI) for geoscientific numerical models (Peckham et al., 2013).Although terrainbento does not yet fully implement a BMI, its model-control functions follow the conventions used by the Landlab toolkit, which themselves have a close parallel with the main BMI model-control functions.The terrainbento initialize method is fully compatible with the BMI method of the same name, which takes as an argument a string containing the name of a parameter-input file (terrainbento's version can alternatively accept a Python dictionary containing parameter name-value pairs).The terrainbento run_one_step method serves the same function as BMI's update, but accepts step size as an argument.Terrainbento's run_for is similar to BMI's update_until (the former takes a duration, whereas the latter takes an absolute time).

Derived classes and use of Landlab components
Two features make the process of writing a new model program in terrainbento relatively fast and efficient: the ability to inherit functionality from the terrainbento base classes and the use of process components in the Landlab toolkit to handle individual process laws.Having already discussed the base class, it is useful to say a few words about Landlab.The Landlab toolkit is a Python-language software library designed to support the efficient creation, exploration, and modification of two-dimensional numerical models of Earth surface processes (Hobley et al., 2017).Landlab accomplishes this by using a CSDMS-inspired plug-and-play method, in which the functionality needed for a numerical implementation of a single process is encapsulated in a standard-format process component.Process components are implemented as Python classes.Landlab also uses an objectoriented approach to grid creation and management so that a simulation grid is encapsulated as a Python object.Components normally interact with a grid object and share fields (arrays) of grid-linked data by creating and attaching the necessary fields to a common grid.More information about Landlab can be found in Hobley et al. (2017).
Terrainbento uses Landlab components to implement its process laws.Each terrainbento model program is implemented as a class that derives from the ErosionModel base class.The model program's __init__ method handles parameter retrieval and instantiates the necessary Landlab components.The model program's run_one_step method then advances each component in turn, normally by calling the component-level run_one_step.In addition to the definition of the model class, each terrainbento model program includes a short main function that allows the model program to be run in a stand-alone fashion (as opposed to being instantiated and run from an outside script, which can also be done).This simple design allows the main model program files to be quite short, often with between 100 and 300 lines, of which only 20-50 are "true" lines of code and the remainder are comments or built-in documentation.

Model and class naming scheme
The naming scheme for the classes that implement the individual terrainbento model programs starts with the name  2).For example, the BasicTh model program uses a threshold formula for water erosion, but is otherwise identical to the Basic model program.Model BasicRtTh uses a threshold and also implements two separate lithologies (here, Rt stands for "rock and till," a name that reflects the original motivation for this particular capability).

Input/output formats and semantics
Terrainbento 1.0 provides multiple options for model instantiation and the specification of parameter values and run-control options.Parameters can be listed in an ASCIItext input file using YAML format ("YAML Ain't Markup Language"), as in the example in Fig. 4, and passed to a from_file constructor.Alternatively, models can be instantiated using a dictionary or with parameters passed as keyword arguments.
If a user wishes to read in a digital elevation model (DEM) to use as the initial topography, the name of the DEM file is given as a parameter in the input file or dictionary.All input/output methods supported by Landlab's create_grid function are supported by terrainbento.For two-lithology (Rt) model programs, the user must specify the elevations of the contact between the two units at each grid node.
Gridded output is written in NetCDF format.The base name for the output files must be specified as an input parameter.When a terrainbento model program runs, output is written at regular intervals, with the frequency set by the user via an input parameter.One file is created for every output in- terval; these files are numbered sequentially.A terrainbento output file contains all of the grid fields used in that particular simulation, which is to say all the grid fields created by that model program's Landlab components plus any created in the main model program.
In addition to output in the form of NetCDF files, terrainbento supports the supply of one or more functions or classes, termed output writers that are run at output intervals.If writing and then post-processing the NetCDF files is not sufficient for a user's application, the user can define an output writer to suit the application.For example, if users wanted to make a diagnostic plot to monitor a simulation run as it progresses, they could define an output writer that does this.The interface constraints on output writers are minimal: a function must take only one argument, expected to be a terrainbento model instance; a class must take one argument at instantiation, also expected to be a model instance, and must have a run_one_step method that takes no arguments.Examples of output writer usage are presented in the Jupyter notebook "Introduction to terrainbento output writers" and in the coupled model notebooks.
Unique names are assigned to each terrainbento input parameter and each data field.Terrainbento 1.0 parameter and field names are listed in Table 4, together with their equivalent mathematical symbols.Terrainbento 1.0 follows the naming conventions used by Landlab (see Hobley et al., 2017).These conventions are loosely based on the CS-DMS Standard Names (Peckham et al., 2013), whose syntax uses an "object plus value" pattern (for example, to-pographic__elevation).Both Landlab and terrainbento 1.0 names seek a balance between brevity, information content, and consistency with the CSDMS Standard Names.Many of the terrainbento and Landlab names are shorter than their full Standard Name equivalents (which can be quite lengthy), but are designed to be similar enough to allow for one-to-one automated mapping.Examples of input-parameter names are shown in the input file example in Fig. 4. Similar principles apply to the field names, which are encoded in the NetCDF output files.The analytical solution unit tests represent model program verification.Tests of the ErosionModel base class include verification that the same random seed reproduces the same initial condition topography, that ErosionModel can work with different instantiation methods and Landlab grid types, and that ErosionModel is compatible with boundary condition handlers and output writers.For base classes like the StochasticErosionModel, we test random seed reproducibility and to ensure that the sequence of rain events generated matches the desired distribution.For each model program, we test at least two process end-members: a case with only hillslope processes and a case with only water-erosion process.Here, for the sake of illustration, we provide the example of an analytical solution for the Basic model program.

Unit tests and model program verification
The Basic model program has the following governing equation:

Symbol Name Dimensions
a Becomes a field rather than single-value parameter in Dd models.b Units depend on value of m.Here we use m = 0.5.The common simplification of Q = A also changes the units of K.
Table 5. Terrainbento field names and unit dimensions.

Symbol Name Dimensions
Given boundary conditions of a constant relative uplift rate U of the core nodes and taking Q = A, at steady state this equation becomes In the water-erosion-only end-member, D = 0 and the equation can be rearranged for a relationship between slope and drainage area: In unit tests we assert that the Basic model program run to steady state conforms to this slope-area relationship.With regard to the hillslope-process-only case (K = 0), the governing equation under conditions of constant uplift can be rearranged for an expression of elevation as a function of position within the domain.For a simulation domain of size L in the x dimension and only one row of core nodes in the y dimension, this analytical solution is Unit tests for all other model program can be found in the source code under the folder tests.
6 Support for users and potential developers Terrainbento includes eight Jupyter notebooks designed to introduce new users to terrainbento and demonstrate five benchmark examples.These notebooks are available on GitHub (Barnhart et al., 2019b).Three introductory notebooks go over the philosophy of terrainbento model programs and an introduction to using them, an overview of how to use the base-level handler methods, and examples of creating and using output writers.
The five example landscapes shown in Fig. 1 are benchmark examples in which a terrainbento model program is created from a Python dictionary and run to steady state with output saved to a compiled NetCDF.In each of these benchmark example notebooks a 2-D image of elevation and slope-area plot shows example model results.
As terrainbento was designed to be generic, it includes a template to support interested developers in building their own model programs within the framework.This template provides an example file with the skeleton of a terrainbento model program and extensive comments on the type of documentation and public functions required of new terrainbento model programs.Throughout the documentation we have made notes encouraging users and developers to make a GitHub issue if they have questions, find errors, or feel that the functionality should be expanded to meet research needs.

Conclusions
Terrainbento 1.0 is a model analysis package and collection of alternative model programs for long-term landscape evolution built using the Landlab framework.Terrainbento was designed to enable hypothesis testing among alternative numerical models of Earth surface processes.Terrainbento 1.0 focuses on 11 binary options for process formulation and 28 model programs that systematically explore these options.As boundary conditions can substantially influence simulation results, terrainbento also includes five boundary condition handlers, which permit the consideration of boundary conditions in a parameterized way.Integration between terrainbento and Landlab permits process component development within Landlab and the use of new components in terrainbento.Thus, while the process combinations available within terrainbento 1.0 are not exhaustive, its extensible design facilitates inclusion of additional terrainbento model programs using new and existing Landlab components.
Recent work has yielded a plethora of numerical models of Earth surface processes, yet comparison among these models has long been difficult and inconsistent.Terrainbento enables efficient numerical model intercomparison with standardized parameters, input/output, and handling of boundary conditions.Consistent, reproducible comparison among landscape evolution models using the terrainbento modeling package will support model evaluation and advance quantitative understanding of Earth surface dynamics.
Code availability.The terrainbento source code can be found in a publicly available GitHub repository distributed under an MIT license (Barnhart et al., 2019a); the user manual is built into the source code Docstrings and compiled into a Read the Docs web page, and Jupyter notebooks that introduce terrainbento usage and show example simulations are located in a public GitHub repository (Barnhart et al., 2019b).Pre-built binaries are available through conda-forge (https://anaconda.org/conda-forge/terrainbento).
Continuous integration and dependencies.Terrainbento 1.0 is tested with the continuous integration tools TravisCI and AppVeyor to ensure that it can be installed, and all tests pass on three operating systems (Windows, Ubuntu Linux, and Mac OSX) and three Python versions (2.7, 3.6, and 3.7).Installing terrainbento from source requires Python and setup tools.Running terrainbento additionally requires NumPy (version 1.11 or higher), SciPy, xarray, dask, six, pyyaml, pytest, and Landlab (version 1.7 or higher).Testing terrainbento additionally requires pytest.surface-water unit discharge L 2 T −1 q tot total unit discharge L 2 T −1 q s fluvial sediment flux per unit width L 2 T −1 q h hillslope sediment flux per unit width L 2 T −1 q ss maximum subsurface unit discharge L 2

B5 BasicCh
BasicCh uses a nonlinear law for hillslope erosion and transport:

B7 BasicVs
The BasicVs model program implements variable source area runoff using the "effective area" approach described in Sect.3.5.1: A eff = Ae −αS/A , (B16) The function δ(H ) is used to indicate that water erosion will act on soil where it exists and on the underlying lithology where soil is absent.To achieve this, δ(H ) is defined to equal 1 when H > 0 (meaning soil is present) and 0 if H = 0 (meaning the underlying parent material is exposed).Parameters: K, m, n, D, P 0 , H s , and H 0 .

B9 BasicRt
BasicRt modifies Basic by allowing for two lithologies, as described in Sects.2.2.2 and 3.7.2: where W c is the contact-zone width.Parameters: K 1 , K 2 , m, n, D, and W c (plus specification of η C (x, y)).

B10 BasicCc
BasicCc uses the same governing equation as Basic, but allows the parameter K to vary through time according to a linear function: Parameters: K 0 , m, n, D, f (factor by which K is larger (f > 1) or smaller (f < 1) than K 0 at t = 0), and T s (time at which K becomes constant).

B11 BasicStTh
The land surface evolution equation is

B12 BasicThVs
The BasicThVs model program implements variable source area runoff using the "effective area" approach plus a threshold on the water-erosion law: where W c is the contact-zone width.Parameters: K 1 , K 2 , m, n, D, ω c1 , ω c2 , and W c (plus specification of η C (x, y)).

B23 BasicChRt
This model program uses nonlinear hillslope transport and two lithologies:

B25 BasicSaVs
This Parameters: K, m, n, K sat , R m , D, H 0 , P 0 , and H s .

B26 BasicRtVs
BasicRtVs is a two-lithology model program that uses variable source area hydrology:

B27 BasicRtSa
This model program combines a dynamic soil layer and two lithologies:

Figure 1 .
Figure 1.Three-dimensional view of simulated topography illustrating the use of five different terrainbento model programs: (a) Basic, (b) Basic with m = 0.25, (c) BasicCh, (d) BasicVs, and (e) BasicRt.Each landscape was initialized with the same random noise field, has the same boundary conditions of the center core nodes uplifted relative to a fixed boundary, and was run to steady state (based on topographic change between 1000-year intervals exceeding 1 mm at no grid cell).Each domain is 1 × 1.6 km, has 10 m grid spacing, and is displayed with a vertical exaggeration of 5 to 30×.

Figure 2 .
Figure 2. Slope-area relationship for the example simulation in Fig. 1a.

Figure 3 .
Figure 3. Illustration of the functional form of the smooth-threshold erosion law (Eq.31) compared with the more traditional hardthreshold formulation.

Figure 4 .
Figure 4. Example of a terrainbento input file.
Terrainbento follows modern software engineering best practices by incorporating documentation and testing into the package source code.The terrainbento documentation Docstrings include simple examples showing, for each model program, the minimum requirements to instantiate a model program instance.For each model program, unit tests verify that all value and compatibility checks raise the correct errors and that all existing analytical solutions are reached.These unit tests ensure that any refactoring of the code, additions or improvements in later terrainbento versions, or changes to Landlab components do not change the results produced by the model.Terrainbento 1.0 has 100 % coverage, which means that all lines of code in the base classes, derived models, and boundary condition handlers are tested by unit or Docstring tests.
is a critical slope gradient.Parameters: K, m, n, D, S c , and N.B6 BasicStBasicSt uses a stochastic representation of precipitation, in which the rainfall rate p is a random variable.The evolution equation is∂η ∂t = −K Qm S n + D∇ 2 η .(B11)The discharge, Q, associated with a particular value of p isQ = p − I m 1 − e −p/I m .(B12)The probability distribution of p is given by a stretched exponential survival function:Pr(P > p) = exp − cand scale parameter λ.The relationship between λ and the mean rainfall rate p d is p d = λ (1 + 1/c) .(B14) Parameters: K, m, n, D, I m , p d , and c.

Table 1 .
Model program binary options.

Table 2 .
Summary of 28 individual model programs in the terrainbento 1.0 collection.
scape evolution model program, whereby other processes may produce gradients equal to or greater than S c .Some authors have addressed this problem with a modified form that avoids divergence at gradient S = S c (e.g.,Carretier and Lu- www.geosci-model-dev.net/12/1267/2019/Geosci.Model Dev., 12, 1267-1297, 2019 1274 K. R. Barnhart et al.: Terrainbento 1.0 Shobe et al. (2017)ickness and allowed progressively more bedrock erosion as sediment thickness H declined relative to median grain size D 50 .He tested both exponential and linear numerical models for the relationship between bedrock exposure and the ratio H /D 50 .Zhang et al. (2015)compared sediment thickness to a statistical description of the macroscale bedrock roughness to determine the probability of bedrock being exposed.The probability of bedrock exposure increased with declining sediment thickness and increasing bedrock surface roughness.In terrainbento we use the Landlab component developed byShobe et al. (2017), which expresses the Stream Power with Alluvium Conservation and Entrainment (SPACE) numerical model.They used an exponential expression describing increases in bedrock exposure as sediment thickness declines relative to bedrock surface roughness.The SPACE numerical model tracks topographic elevation η as well as bedrock surface elevation η b and sediment thickness H such that

Table 3 .
Public base class methods.Clean up single time step NetCDF files * Empty function intended to be overridden by child class.

Table 4 .
Terrainbento parameter names and unit dimensions.
This is a sediment-tracking (hybrid) model program with two lithologies: B16 BasicDdVsModel BasicDdVs uses variable source area hydrology and an erosion threshold that increases with progressive erosion Geosci.Model Dev.,12, 2019www.geosci-model-dev.net/12/1267/2019/ model program combines variable source area hydrology with a dynamic soil layer.Unlike other model programs with variable source area hydrology, here the actual soil thickness H (x, y, t) is used to calculate transmissivity.