The Landlab v1.0 OverlandFlow component: a Python tool for computing shallow-water flow across watersheds

<jats:p>Abstract. Representation of flowing water in landscape evolution models (LEMs) is often simplified compared to hydrodynamic models, as LEMs make assumptions reducing physical complexity in favor of computational efficiency. The Landlab modeling framework can be used to bridge the divide between complex runoff models and more traditional LEMs, creating a new type of framework not commonly used in the geomorphology or hydrology communities. Landlab is a Python-language library that includes tools and process components that can be used to create models of Earth-surface dynamics over a range of temporal and spatial scales. The Landlab OverlandFlow component is based on a simplified inertial approximation of the shallow water equations, following the solution of de Almeida et al.(2012). This explicit two-dimensional hydrodynamic algorithm simulates a flood wave across a model domain, where water discharge and flow depth are calculated at all locations within a structured (raster) grid. Here, we illustrate how the OverlandFlow component contained within Landlab can be applied as a simplified event-based runoff model and how to couple the runoff model with an incision model operating on decadal timescales. Examples of flow routing on both real and synthetic landscapes are shown. Hydrographs from a single storm at multiple locations in the Spring Creek watershed, Colorado, USA, are illustrated, along with a map of shear stress applied on the land surface by flowing water. The OverlandFlow component can also be coupled with the Landlab DetachmentLtdErosion component to illustrate how the non-steady flow routing regime impacts incision across a watershed. The hydrograph and incision results are compared to simulations driven by steady-state runoff. Results from the coupled runoff and incision model indicate that runoff dynamics can impact landscape relief and channel concavity, suggesting that, on landscape evolution timescales, the OverlandFlow model may lead to differences in simulated topography in comparison with traditional methods. The exploratory test cases described within demonstrate how the OverlandFlow component can be used in both hydrologic and geomorphic applications.                     </jats:p>

section_content: Numerical models of overland flow have a variety of applications. Examples include mapping urban flooding events (e. g., Dutta et al., 2000; Morrit and Bates, 2002; Maksimović et al., 2009; Kulkarni et al., 2014; Cea and Bladé, 2015), understanding the interactions between surface and subsurface water by way of soil infiltration (e. g., Esteves et al., 2000; Panday and Huyakorn, 2004; Kollet and Maxwell, 2006; Maxwell and Kollet, 2008; Shrestha et al., 2015), and exploring hydrogeomorphologic processes in natural landscapes (e. g., De Roo et al., 1996; Beeson et al., 2001; Francipane et al., 2012; Kim et al., 2013; Wang et al., 2014; Published by Copernicus Publications on behalf of the European Geosciences Union. Rengers et al., 2016 ). Yet to be deeply explored is how the details of hydrologic processes, specifically runoff generation, impact landscape evolution over centennial scales and longer. Pioneering work by Tucker and Bras (1998) and Sólyom and Tucker (2004) explored this problem, but many questions remain, including how hydrograph shape impacts erosion rates and topographic patterns. Models of landscape evolution share the same fundamental structure: all use numerical methods to model flow or transport of water and sediment across a representative mesh that is tessellated into discrete elements (e. g., Willgoose et al., 1991; Tucker and Slingerland, 1994; Willgoose, 1994; Braun and Sambridge, 1997; Tucker et al., 2001; Coulthard, 2001; Coulthard et al., 2002; Willgoose, 2005; Tucker and Hancock, 2010; Chen et al., 2014; Hancock et al., 2015). However, the complexity of the runoff mechanism varies. The representation of surface water flow in landscape evolution models (LEMs) is often simplified, as solving the shallow water equations in 2-D can be computationally intensive. Most models assume unidirectional steady-state water discharge, where surface water flux is modeled at each location as a product of drainage area and rainfall rate, or

where Q ss is the steady-state water discharge (L 3 T -1 ), P is the spatially averaged effective precipitation or runoff rate (L T -1 ) and A is drainage area (L 2 ). Discharge increases moving downstream with drainage area but only lasts for the duration of a precipitation event. If the precipitation rate is constant, the discharge rate at a given point in the domain will be constant for the duration of the storm event, creating a rectangular hydrograph (Fig. 1 ). In more physically based models, the steady-state assumption is replaced with non-steady runoff processes that simulate flowing water across a watershed. Figure 1 compares the steady-state discharge assumption to a non-steady method at one location in the watershed. The effective rainfall rate P is the same rate and duration for both the steady (Q ss ) and non-steady (Q h ) discharge simulations. The non-steady hydrograph (Q h ) lasts longer than rectangular steady-state hydrograph (Q ss ), as water takes time to flow across the landscape, a process controlled by the physical nature of the system. The simplifying assumption of steady-state discharge is made for two reasons: there can be significant differences between hydrologic timescales for individual flood and storm events (minutes to days) and geomorphic timescales of rock uplift and landscape evolution (thousands to millions of years) that may be complex to resolve. Additionally, computational power is often a limiting factor, as some processes in LEMs do not lend themselves to parallelization, so making assumptions about how water fluxes are calculated (e. g., Eq. 1) can speed up model processing time. 

Whereas many LEMs generalize surface water flow using steady-state assumptions, most physical models of runoff production simulate changing surface water discharge Image illustrating the differences between steady-state and non-steady hydrology and incision at a single point within a watershed. In this schematic, the effective precipitation rate (P ) is the same for both steady and non-steady cases. During the precipitation event, steady discharge (Q ss ) and incision rate (I ss ) are constant, driven by that effective precipitation rate and drainage area (A), erodibility (K), water surface slope (S), and stream power exponents (m sp, n sp ). In the non-steady case, a wave front begins to propagate and incise, producing time-varying discharge (Q h ), calculated using physical parameters such as water depth (h), water surface slope (S), and Manning's roughness coefficient (n). Nonsteady incision rate (I h ) is calculated using the time-varying discharge, erodibility, and water surface slope. At the end of the precipitation event, Q ss and I ss also end, while non-steady values Q h and I h continue until all water has completely exited the system at the outlet. 

through time, capturing the spatial and temporal variability of flowing water across a modeled landscape (e. g., Ogden et al., 2002; Downer and Ogden, 2004; Ivanov et al., 2004; Hunter et al., 2007; Moradkhani and Sorooshian, 2009; Devi et al., 2015). Surface water runoff is one of many physical processes and parameters explored in these models. Some of these runoff models have been paired with erosional models at the watershed scale (e. g., Aksoy and Kavvas, 2005; Francipane et al., 2012; Coulthard et al., 2013; Kim et al., 2013). However, there are a limited number of studies that integrate a physically based distributed runoff method into a LEM framework; the steady-state discharge assumption (Eq. 1) is often used instead. 

The assumption of steady-state discharge in LEMs is not always reasonable. Steady-state hydrologic conditions are rarely achieved in larger catchments with long flow paths or in landscapes dominated by short-duration precipitation events. Additionally, the traditional steady-state model (Eq. 1) does not capture differences in basin organization or orientation, whereas discharge is known to be sensitive to these characteristics (Snyder, 1938). For example, watersheds with identical drainage areas but different shapes or orientations may have dramatically different hydrograph shapes that are not captured by the traditional steady-state assumption. 

Adding hydrologic variability to LEMs has also been shown to impact watershed morphology and landscape evolution. Previous work coupling spatially variable rainfall models with steady-state discharge in erosion models has illustrated impacts on landform morphology, including relief and drainage network organization (e. g., Anders et al., 2008; Colberg and Anders, 2014; Huang and Niemann, 2014; Han et al., 2015). Similarly, introducing storm and discharge variability into LEMs has implications for incision rates, channel profile form, and steepness in modeled landscapes (e. g., Tucker and Bras, 2000; Lague et al., 2005; Molnar et al., 2006; DiBiase and Whipple, 2011). Coulthard et al. (2013) integrated a semi-implicit hydrodynamic model into the CAESAR LEM and noted reduced sediment yields on decadal timescales of landscape evolution when using non-steady runoff. In another approach, Sólyom and Tucker (2004) estimated non-steady peak discharge as a function of storm duration, rainfall rate, and the longest flow length in a network. Incision rates were estimated using those peak discharge values. Their findings demonstrated that landscapes evolved with non-steady water discharge were characterized by decreased valley densities, reduced channel concavities, and increased relief when compared to landscapes evolved using steady-state runoff. 

To investigate the role of non-steady flow routing on landform evolution, a hydrodynamic model has been incorporated into the Landlab modeling toolkit. In this paper, we describe the fundamentals of the Landlab modeling framework, the theoretical background of the Landlab Overland-Flow component, based on a two-dimensional flood inundation model (LISFLOOD-FP; Bates and De Roo, 2000; Bates et al., 2010; de Almeida et al., 2012; de Almeida and Bates, 2013) and how this model was adapted to work in coupled geomorphic-hydrologic applications. This description of the new OverlandFlow component includes information on how to set up a model domain using a digital elevation model, how to handle boundary conditions, how Landlab components store and share data in "fields", and the validation against a known analytical solution. The OverlandFlow component is then used to route non-steady flow on one real and two synthetic watersheds. Model output demonstrates that the Over-landFlow component is sensitive to both catchment characteristics and precipitation inputs. Output hydrographs can be flashier or broader depending on changes in these parameters and model domain. Finally, the variable discharge from the OverlandFlow component is coupled to a detachmentlimited erosion component (DetachmentLtdErosion) to explore the feedbacks between hydrograph shape and shortterm (10-year) erosion patterns throughout a landscape.
section_name: Introduction

section_content: Landlab is a Python-language, open-source modeling framework, developed as a highly flexible and interdisciplinary library of tools that can be used to address a range of hypotheses in Earth-surface dynamics (Adams et al., 2014; Tucker et al., 2016; Hobley et al., 2017). The utilities in Landlab allow users to build two-dimensional numerical models (Fig. 2 ). This includes a gridding engine that creates structured or unstructured grids, a set of pre-built components that implement code representing Earth-surface or near-surface processes, and structures that handle data creation, management, and sharing across different process components. A diverse group of processes, such as uniform precipitation, detachment-limited incision, linear diffusion, crustal flexure, soil moisture, vegetation dynamics, and overland flow, is available in the Landlab library as process components. The Landlab architecture allows for a "plug-and-play" style of model development where process components can be coupled together. Coupled components share a grid instance and can operate on the data attached to the grid. 

Landlab offers several different grid types. However, because the core algorithm in the OverlandFlow component can only be applied to structured grids, only the RasterModel-Grid class is described here. The RasterModelGrid class can build both square ( x = y) and rectangular ( x = y) grids. OverlandFlow methods only operate on square grid cells and require x = y. Each grid type in Landlab is composed of the same topological elements: nodes, which are points in (x, y) space; cells, a polygon with area x y surrounding all non-perimeter or interior nodes; and links, ordered line segments which connect neighboring pairs of nodes and store directionality (Fig. 3 ). In the RasterModel-Grid library, each node has four link neighbors, each oriented in a cardinal direction. Each node has two "inlinks" connecting a given node to its south and west neighbors, and two "outlinks" connecting to the node neighbors in the north and east. The terms inlinks and outlinks are for topological reference only, as the direction of fluxes in a typical Landlab component is based on link gradients. Model data are stored on these grid elements using Landlab data fields. The data fields are NumPy array structures that contain data associated with a given grid element. To store and access data on these fields, data are assigned using a string keyword and are accessed using Python's mutable dictionary data structure. Data are attached to the grid instance using these fields and can be accessed using the string name keyword and updated by multiple Landlab components. For example, a field of values representing water depth at a grid node can be accessed using the following (Peckham, 2014; Hobley et al., 2017). Model boundary conditions are set within a Landlab grid object. Boundary conditions are set on nodes and links (Fig. 4 ). Node boundary statuses can be set to either "boundary" or "core". If a node is set to boundary, it can be further defined as an open, fixed gradient, or closed (no flux) boundary. In all RasterModelGrid instances, default boundary conditions are set as follows: perimeter nodes are open boundary nodes, while interior nodes are set as core nodes. Boundary conditions can also be applied to interior nodes (e. g., NO-DATA values on non-perimeter nodes in a digital elevation model can be set as closed boundaries). In OverlandFlow applications, open boundary nodes act as flow outlets, allowing water fluxes to move out of the model domain. Input rainfall is added to all core nodes, where water depths are updated at each time step to drive fluxes on grid links. 

There are three link boundary statuses: active, inactive, and fixed. Link boundary status is tied to the neighboring nodes. Once boundary conditions are set on the nodes, link (2) between any pair of boundary nodes of any type (Fig. 4 ). Core nodes and active links make up the computational domain of a Landlab model. 

3 Component equations
section_name: Landlab modeling framework

section_content: Solving explicit two-dimensional hydraulic formulations can be computationally challenging. For example, the 1-D shallow water equation includes four terms:

where Q is water discharge (L 3 T -1 ); t is time (T); x is the location in space (L); A xs is cross-sectional area of the channel (L 2 ); g is gravitational acceleration (L T -2 ); h is water depth (L); z is the bed elevation (L); n is the Manning's friction coefficient (L -1/3 T) and R is the hydraulic radius (L). These terms represent, from left to right, local acceleration, advection, fluid pressure, and friction slope. To enhance stability, many solutions of the shallow water equations include numerical approximations that neglect terms from this solution. The simplest approximation, the kinematic wave model, neglects the local acceleration, advection, and pressure terms. A more complex approximation, the diffusive wave model, only neglects the local acceleration and advection terms (Kazezyılmaz-Alhan and Medina Jr., 2007). The Landlab OverlandFlow component adapts a twodimensional hydrodynamic algorithm to simulate flow at all points across the gridded domain. This algorithm, developed for the LISFLOOD-FP model, was incorporated into Landlab for modeling overland flow. Similar to the diffusive approximation, the LISFLOOD-FP algorithm assumes a negligible contribution from the advection term of the shallow water equations (Bates et al., 2010; de Almeida et al., 2012). Additionally, this solution assumes a rectangular channel structure and constant flow width, impacting the pressure and friction terms (A xs and R) in Eq. ( 2 ) (Bates et al., 2010). This formulation allows for a larger maximum time step than the more common diffusive approximation, enhancing the computational efficiency of the OverlandFlow component. The work of de Almeida et al. ( 2012 ) further stabilized this algorithm by introducing a diffusive term into LISFLOOD-FP, updating the Bates et al. (2010) algorithm to work on lower friction surfaces without sacrificing computational speed. 

To start the model, a stable time step is calculated. Stable time steps are set according to the Courant-Friedrichs-Lewy criteria which evaluate the ratio of time step size to grid resolution. If large time steps are used, areas of high slope are prone to wave oscillations, leading to a spatial "checkerboard" pattern of water depths. If time steps are very small, there are significant impacts on the computational performance of a model. To maximize the tradeoff between computational efficiency and stability of the de Almeida et al. (2012) solution, an adaptive time step (following Hunter et al., 2005) is used to keep the CFL condition valid:

where t max is the maximum time step that adheres to the CFL condition; α is a dimensionless stability coefficient less than 0. 7; x is the grid resolution (L); and √ gh max, the characteristic velocity of a shallow water wave, or the wave celerity (L T -1 ), is calculated using h max, the maximum depth of water in the modeling domain (L). When the OverlandFlow component is initialized, a thin film of water is set at all grid nodes to keep Eq. ( 3 ) valid. Flow stability and mass balance are controlled by the α value. On a case-by-case basis, α must be tuned to find the value that keeps the modeled flow stable while also reducing mass losses. Variables and parameters are defined in Tables 1 and 2. q x q x-1 q x+1

Figure 5. In the de Almeida et al. ( 2012 ) equation, flux information from neighboring links is used to calculate surface water discharge. In this sample one-dimensional grid, discharge is calculated in the horizontal (subscript x) direction on links. Here, discharge is calculated at location q x using the left neighbor (q x-1 ) and right neighbor (q x+1 ) flux values, following Eq. (4). 

To calculate water discharge at all grid locations, de Almeida et al. ( 2012 ) derived an algorithm, using the onedimensional Saint-Venant or shallow water equations, which simulates a flood wave propagating across the domain. This simplified algorithm calculates discharge at all points within the domain (for the full derivation, see de Almeida et al., 2012). The explicit solution follows the form

where q is water discharge per unit width (L 2 T -1 ), calculated on links, here given superscript t for the current time step and subscript x describing the location of links in space (Fig. 5 ). θ is a weighting factor between 0 and 1, given a default value of 0. 8, but it can be tuned by the user. Setting θ to 1 returns the semi-implicit solution of Bates et al. (2010), that is, removing the diffusive effects implemented by de Almeida et al. (2012). g is gravitational acceleration (L T -2 ); h f is the local maximum water surface elevation at a given time (L); t is the adaptive time step (T) (Eq. 3); S w is the dimensionless water surface slope; and n is Manning's friction coefficient (L -1/3 T) (Tables 1 and 2 ). Equation ( 4 ) is calculated as two one-dimensional solutions in a D4 (fourdirection) scheme: first calculated in the east-west direction (in the x direction) and then in the north-south direction (replacing x with y in Eq. 4). 

Water depth is calculated on nodes and updated at each time step as a function of the surrounding volumetric water fluxes (q • x) on both horizontal and vertical links:

where Q h(in) (L 3 T -1 ) are the summed water discharges moving into a given node and Q h(out) are summed water discharges moving out of a given node, following Fig. 3. Directionality of discharge is determined not by the orientation of inlinks or outlinks; instead, flow directions are determined by the water-surface gradient of each link. In this method,

water mass is conserved, as the flow moving out of a node is balanced by the flow moving into the nearest node neighbors. By default, this model assumes that all rainfall is spatially uniform and temporally constant, and all rainfall is converted to surface runoff. No infiltration or subsurface flow is considered within the model equations; however, the Overland-Flow component could be easily coupled with an infiltration component. Spatially or temporally variable rainfall could be generated by another process component or set manually by the user in a driver file. Effective rainfall depths are applied over the basin and added to the surface water depths at each time step. 

The de Almeida et al. ( 2012 ) equation is designed for urban flooding events and is most stable in low-to-zero slope environments. To adjust this component to work in steep mountain catchments, extra stability criteria were added to keep simulations numerically stable using the steep_slopes keyword flag. A similar criterion was implemented in the CAESAR-LISFLOOD model (Coulthard et al., 2013). This method reduces the calculated flow discharge as needed to keep flow regime critical to subcritical using the Froude number (Eq. 6), where subcritical flow is defined as F r < 1. 0. The Froude number is calculated as a function of wave veloc-Algorithm 1 Sample Landlab overland flow and erosion model 1: from landlab. components import OverlandFlow, DetachmentLtdErosion, SinkFiller #Import Landlab components and utilities 2: from landlab. io import read_esri_ascii 3: (grid, elevations) = read_esri_ascii(asc_file = "watershed_DEM. asc", name = "topographic__elevation") #Read in DEM and create grid 4: grid. set_watershed_boundary_condition(elevations, nodata_value = -9999. 0) #Set boundary conditions 5: effective_rain_rate_ms = 5. 0 × (2. 78 × 10 -7 ) #Convert rainfall from mm h -1 to m s -1 6: dle = DetachmentLtdErosion(grid) #Instantiate components and set parameters 7: of = OverlandFlow(grid, steep_slopes = TRUE, rainfall_intensity = effective_rain_rate_ms) 8: sf = SinkFiller(grid, routing = "D4") 9: sf. fill_pits() #Pre-process DEM and fill pits in D4 flow-routing scheme 10: elapsed_time = 0. 0 ity (u, calculated as q h f on all links) and wave celerity ( √ gh f ):

If the steep_slopes flag is set when initializing Overland-Flow, restrictions are imposed to keep flow conditions critical to subcritical, a reasonable assumption for steep mountain catchments (Grant, 1997). Specifically, if the water velocity calculated by the component drives the Froude number greater than 1. 0, water velocity is reduced to a value that maintains a Froude number equal to 1. 0 for that given time step. This prevents water from draining too quickly, creating oscillating flow depths in steep reaches.
section_name: The deAlmeida OverlandFlow component

section_content: To illustrate the flexibility of the OverlandFlow component, we present an example in Sect. 7, in which water discharge calculated by the OverlandFlow component is used in the erosion component. Specifically, we explore a case where incision rate is solved explicitly and depends on local water discharge and water surface gradient (e. g., Howard, 1994; Whipple and Tucker, 1999, 2002; Pelletier, 2004). This equation follows the form

where I is the local incision rate (L T -1 ); K is a dimensional erodibility coefficient, where the units depend on the positive dimensionless stream power coefficient m sp, whereas the value of m sp is correlated with the other dimensionless stream power coefficient n sp. Q is total water discharge on a node at a given time step (L 3 T -1 ); S w max is the local maximum water surface slope, which is dimensionless, and β is the optional threshold, below which there is no incision (L T -1 ) (Tables 1 and 2 ). By default, m sp and n sp have set values of m sp = 0. 5 and n sp = 1. 0 that can be adjusted by the model user. This erosion formulation is implemented with the Landlab DetachmentLtdErosion component. This solution allows for only the local detachment of material and assumes that transport rate is much larger than sediment supply rate; therefore, no deposition is considered here. For simplicity, no threshold (β) is applied in the following applications.
section_name: DetachmentLtdErosion component

section_content: To use the coupled Landlab OverlandFlow and Detach-mentLtdErosion model, the user interacts with a driver file (Fig. 2 ). A simple Landlab driver file can run a model using fewer than 20 lines of code (Algorithm 1). There are four parts to running the coupled OverlandFlow-DetachmentLtdErosion model: (1) creating a domain using RasterModelGrid, either explicitly or using a digital elevation model (DEM) in the ArcGIS ASCII format;

(2) setting boundary conditions on the domain;

(3) initializing the components; and (4) coupling them using the Landlab field data structures.
section_name: OverlandFlow model implementation in Landlab

section_content: To set up a grid instance, the user can create a rectangular grid by passing the number of rows, number of columns, and grid resolution ( x) as keywords to the RasterModelGrid object. After Landlab and RasterModel-Grid are imported, this can be accomplished in one line of code: grid = RasterModelGrid((number_of_node_rows, number_of_node_columns), x). In this method, only an empty instance of the grid is created, so elevation data must be assigned to grid nodes by the user. 

An alternative method is to read in gridded terrain data from other file types. The original intent of Bates et al. (2010) was to develop a new flood inundation algorithm that could work easily with the growing availability of terrain data collected by satellite, airborne, or terrestrial sensors. Landlab's input and output utilities include functionality to read in data from an ASCII file in the Esri ArcGIS format (Algorithm 1, Line 3). In this method, elevation data are read in and automatically assigned to a Landlab data field called to-pographic_elevation, set using the name keyword.
section_name: Initializing a grid: user-defined or DEM

section_content: Node boundary conditions are set throughout the grid in a Landlab OverlandFlow model to delineate the modeling domain (Algorithm 1, Line 4). For flow to move out of a watershed or system, an open boundary must be set at the outlet(s). If the node location of the outlet is unknown, there is a utility within the grid (set_watershed_boundary_condition; Algorithm 1, Line 4) that will find a single outlet and set it as an open boundary, in addition to setting all NODATA nodes to closed boundaries across the DEM or model domain. For landscapes with multiple potential outlets, such as urban environments, which are not discussed here, the user would have to manually identify and set nodes to open boundary status. 

The de Almeida et al. ( 2012 ) equation uses neighboring link values when calculating water discharge (Fig. 5 ). By default, links on the edge of the watershed are set to inactive status and are assigned a value of 0, meaning no input from outside of the watershed for the simulation. If the user wants to simulate an input discharge on these links, an alternative method is the set_nodata_nodes_to_fixed_gradient method. If this method is called, the user can manually update discharge values on links with fixed link boundary status outside of the OverlandFlow class. Fixed links are accessed through their IDs using the RasterModelGrid class (grid. fixed_links). In this method, the user can set a discharge value per unit width (L 2 T -1 ) on all fixed links. This method is advised if the user has a known input discharge they want to force at the watershed or domain edge.
section_name: Boundary condition handling

section_content: Landlab components have a standard initialization signature and take the grid instance as the first keyword (Algorithm 1, Lines 6-8). Any default parameters are also in the component signature and can be updated when the component is called. These parameters can be adjusted according to the physical nature of the landscape being tested. For the Over-landFlow component, Eq. ( 4 ) parameters Manning's n and discharge weighting factor θ can be adjusted. To keep the time step equation (Eq. 3) valid, an initial thin film of water is set across the model domain using the keyword h_init (Ta-ble 2). A steady, uniform precipitation rate can also be passed as a system input using the rainfall_intensity parameter (Algorithm 1, Line 7). Additionally, a stability criterion flag for steep catchments can be set (steep_slopes = TRUE, as described in Sect. 3. 1. ). In the DetachmentLtdErosion component, stream power exponents m sp and n sp, threshold β, and erodibility parameter K are also set by passing arguments to the component on instantiation.
section_name: Initialize OverlandFlow and DetachmentLtdErosion

section_content: To couple the OverlandFlow and DetachmentLtdErosion components, values for water discharge (Q h ), water surface slope (S w ), and topographic elevation (z) are shared as data fields through the RasterModelGrid instance (e. g., Algorithm 1, Lines 14-15). At each time step, the water discharge and surface water slope fields are updated by the Overland-Flow component (Eq. 4). These new values are used to calculate an incision rate in the DetachmentLtdErosion component (Eq. 7). At each grid location, topographic elevation (z) is reduced according to the incision rate. Changes in topographic slope caused by erosion throughout the landscape will drive changes in surface water slope (S w max ) and discharge (Q h ) in the next iteration of the OverlandFlow component.
section_name: Coupling using Landlab fields

section_content: To validate the OverlandFlow component, we compared model output against an analytical solution for wave propagation on a flat surface, following Hunter et al. (2005). This test case propagates a wave over a flat horizontal surface (a slope of 0), given a uniform friction coefficient (n) and constant, single-direction velocity (u). (For the full derivation, see Hunter et al., 2005; Bates et al., 2010; de Almeida et al., 2012). The analytical solution is

Solving for the leftmost boundary of the modeling domain

All analytical solution tests were modeled across a rectangular RasterModelGrid instance with dimensions of 800 m by 6000 m. The water depth boundary condition (Eq. 9) is applied to the left edge of the domain through time, whereas the top, right, and bottom edges of the grid are set to closed boundary status to keep flow moving uniformly to the east and contained within the computational domain. All input flow remains on the surface of the domain, as no infiltration is considered. Although not illustrated here, mass was conserved in all analytical test cases. Grid set up and test parameters are described in Table 3.
section_name: Analytical solution

section_content: Following Bates et al. (2010), the behavior of OverlandFlow was modeled across a range of grid resolutions. Velocity and surface roughness were held constant throughout all runs (n = 0. 03 s m -1/3, and u = 1. 0 m s -1 ) and θ was set to 1. 0 (Bates et al., 2010, Fig. 2 ). Wave fronts were plotted at model time t = 3600 s. Four grid resolutions were tested: x = 5, 10, 25, and 50 m. These tests envelop a range of resolutions, including the 10 and 30 m dataset resolutions of the United States Geological Survey National Elevation Dataset (USGS-NED, 2017), as well as 30 m datasets from the European Environmental Agency's digital elevation model over Europe (EU-DEM, 2013). Larger grid resolutions ( x > 50 m) are not shown here, as at those coarser grid resolutions the Over-landFlow component becomes sensitive to the initial thin film of water (h_init) that is used to keep the time step (Eq. 3) valid. h_init was set to 1 mm in all test cases described here. The smallest time step over the duration of the x = 50 m test case can be compared to the published value of Bates et al. (2010). Time steps will decrease with increasing water depth, per Eq. ( 3 ). The minimum time step from the Over-landFlow component tests was 7. 25 s, identical to the value provided by Bates et al. (2010). 

In all grid resolution tests, the OverlandFlow predicted wave fronts closely approximate the analytical solution (Fig. 6a ). At the front of the wave, the predicted water elevations from OverlandFlow better approximate the analytical solution as grid resolution increases (Fig. 6b ), as noted by Bates et al. (2010) for the semi-implicit (θ = 1. 0) solution in LISFLOOD-FP. Figure 6 demonstrates that, with only a minor sensitivity at the leading edge of the wave front, the Landlab OverlandFlow model can effectively operate on a wide range of grid resolutions.
section_name: Sensitivity to grid resolution

section_content: To test the Landlab OverlandFlow component with different roughness and resolution characteristics, a RasterModel-Grid instance with dimensions of 800 m by 6000 m was initialized with a resolution of x = 25 m. In order to evaluate the sensitivity to surface roughness (Manning's n), two analytical solution test cases were run on the domain. The first is a low friction test (n = 0. 01 s m -1/3, u = 0. 4 ms -1 ; Fig. 7a, c ) following the solution of Bates et al. (2010) and de Almeida et al. (2012, Fig. 2 ). In the second test, the friction value was increased by an order of magnitude, while velocity was unchanged (n = 0. 1 s m -1/3, u = 0. 4 m s -1 ; Fig. 7b, d ). The two Manning's n values in this test were selected to demonstrate model behavior across a range of conditions: n = 0. 01 s m -1/3 represents urban environments or man-made channel systems; n = 0. 1 s m -1/3 can be used in landscapes or channels characterized by dense brush and tree growth (Chow, 1959). To mirror previous tests using the LISFLOOD-FP model, Fig. 7 shows the water depth of wave fronts at three model times: t = 2700, 5400, and 9000 s. Each dashed line represents a changing θ value in Eq. ( 4 ), with θ = 1. 0 representing the semi-implicit solution of Bates et al. (2010). 

The smallest time step over the duration of the low friction model run (n = 0. 01 s m -1/3 ) can be compared to the published value of de Almeida et al. (2012). The minimum time step from the OverlandFlow component tests, sampled at t = 9000 s, was 8. 6 s, identical to the value provided by de Almeida et al. (2012). 

In all velocity-roughness conditions, the wave fronts predicted by the Landlab OverlandFlow component correlate well with the analytical solution defined using Eq. ( 9 ). In the low friction case (n = 0. 01; Fig. 7a, c ), the wave speed produced using Landlab OverlandFlow is slower than the predicted wave front speed. Increasing surface roughness (n = 0. 1; Fig. 7b, d ) leads to the predicted wave front overestimating the analytical solution. Overall, the close approximation of the modeled solutions to known analytical solutions, across a wide range of roughness values, demonstrates the sensitivity of the Landlab OverlandFlow component to different roughness coefficients and the flexibility of the component to work across a wide range of landscape conditions.
section_name: Sensitivity to surface roughness

section_content: The Landlab OverlandFlow component can be used in hydrology applications, routing precipitation across a real landscape DEM and estimating runoff for every point within a discrete RasterModelGrid instance. Discharge values can be calculated at every point in the watershed and updated at each time step. Updated water depths, driven by changing discharge, can be used to calculate shear stress following the depth-slope product:

Equation ( 10 ) calculates the bed shear stress τ (M L -1 T -2 ) as a function of fluid density ρ (M L -3 ), g gravity, h water depth, and S w surface water slope. Shear stress exerted on the bed can be used to estimate sediment transport driven by flowing water throughout the domain. 

Here, we illustrate a single storm routed across a DEM. In addition to water discharge, water depth and bed shear stress are calculated by the model at all grid locations. This implementation of the OverlandFlow component illustrates how hydrologists can use Landlab as a simplified distributed runoff model to estimate the flow of water and sediment resulting from a single storm on a real landscape.
section_name: Application: modeling OverlandFlow in a real landscape

section_content: To route runoff across a real landscape, a DEM can be read into Landlab and converted easily into a RasterModelGrid instance. The Spring Creek watershed is used in this example, as a preprocessed DEM for the watershed has been used before in Landlab applications (e. g., Adams et al., 2016; Hobley et al., 2017, Fig. 15 ). Spring Creek is a steep 27 km 2 watershed, located within Pike National Forest in central Colorado, USA (Fig. 8a ). This lidar-derived DEM has square cells with a resolution of x = 30 m (DEM data : Tucker, 2010). Using the set_watershed_boundary_condition utility, all NODATA nodes in the DEM are set to closed boundary status (Algorithm 1, Line 4). This method identifies the lowest elevation point along the edge of the watershed, the outlet, and sets it to an open boundary. 

The DEM was preprocessed using the Landlab SinkFiller component to ensure all surface water flow can be removed from the domain. This component fills pits in the DEM in a D4 routing scheme, where all nodes have at least one downstream neighbor in one of the four cardinal directions (Algorithm 1, Lines 8-9). If this step were to be skipped, flow may pond in "lakes" or "pits" in the domain, where flow cannot travel out of a given node location until the water surface elevation of the lake exceeds the bed elevation of one of the four neighboring nodes. 

To initiate flow across the domain, a single storm was routed across the watershed. A theoretical "base storm" (Table 4) was used as an example, with a constant, effective rainfall rate of 5 mm h -1 and a duration of 2 h. The storm event was spatially uniform across the domain and was estimated using NOAA precipitation data from a nearby site in Colorado (NOAA, 2014). For this storm, hydrographs were recorded at three points within the model domain. No infiltration or subsurface flow was considered in this test case. Water depths at every location in the watershed were used to calculate the shear stress, which can be used to make interpretations about the transport of sediment across the watershed as a result of the storm.
section_name: Methods: domain and parameterization

section_content: In order to illustrate the downstream movement of the flood wave, hydrographs were plotted at three locations within the channel. The three hydrographs correspond to the three starred locations on the watershed DEM in Fig. 8a : at the outlet (black line; Fig. 8b ), the approximate midpoint of the main channel (violet line; Fig. 8b ), and an upstream location in the main channel (lavender line; Fig. 8b ). In these hydrographs, peak discharge and time to peak increase as the sampling site nears the outlet (moving from lighter to darker color), demonstrating that the model behaves as expected. 

Water depths are variable at each point throughout the model run, changing as a function of discharge inputs, outputs, and effective rainfall rate at each time step (Eq. 5). Water depth values can be mapped across the domain at discrete time steps. In this example, water depth was plotted at the peak of the outlet hydrograph (Fig. 8c ). The scale in Fig. 8c emphasizes flow patterns in the channels, but water depth and discharge are calculated across the entire watershed, including on the hillslopes. These water depths can be used to calculate shear stress (following Eq. 10). Stress values were tracked at all points throughout the model run, and the local maximum value for each node was plotted in Fig. 8d. Shear stress (τ ) values can be used to interpret the size of particles that can be entrained and transported by surface flow. Greater τ values correspond to areas with greater water depths (e. g., channels), where more sediment transport would be expected in high flow conditions. 

In this example, we illustrate hydrographs across a real landscape and the resulting shear stress values. These results can be used to explore the processes controlling overland flow in a gauged landscape. Shear stress values can be used to estimate sediment transport rates and make interpretations about spatial patterns of erosion and deposition, as well as total sediment yields for particular storm events. These data can be used to explore landscape sensitivity to different rainfall events and runoff conditions.
section_name: Results and implications

section_content: The implementation of the OverlandFlow component in Landlab allows us to investigate the impact of storm characteristics on the resulting hydrograph and how these hydrographs drive erosion processes throughout the basin. Here, we demonstrate the abilities of this new component, how the component resolves the details of the storm hydrograph, and how these hydrographs compare to the traditional steadystate method used in LEMs. Additionally, in coupling this new component with the Landlab DetachmentLtdErosion component, these model results illustrate the erosion magnitudes and patterns in response to a hydrograph and allow us to make inferences about how this type of hydrodynamic model could impact long-term geomorphic evolution of similar watersheds.
section_name: Application: coupling with an erosion component in Landlab

section_content: To test the new Landlab OverlandFlow component, two synthetic watersheds were generated using the Landlab FlowRouter and StreamPowerEroder components (not described here; see Hobley et al., 2017). These basins were evolved to topographic or geomorphic steady state, where uniform rock uplift is matched by erosion at all grid locations, and topography is effectively unchanging through time. Two watershed shapes were modeled: a "square" watershed (Fig. 9a ) and a "long" watershed (Fig. 9b ) to evaluate how hydrograph shapes change with increasing maximum flow length, where the "long" basin has longer flow paths to the outlet when compared to the "square". Each watershed has a drainage area of approximately 36 km 2 at the outlet. The square basin has dimensions of 200 rows by 200 columns; the long basin has dimensions of 400 rows by 100 columns. Cells are square and have a resolution of x = 30 m. Each basin has an open boundary at the watershed outlet, located at the center node of the southernmost grid edge. The remaining southern nodes, along with the west, east, and north grid edges, were set to closed boundary status. 

To initiate flow and incision, three precipitation events were modeled across both watersheds. These storms were represented as spatially uniform across the model domain, and intensities were constant for the given storm duration. No infiltration or subsurface flow was modeled in these test cases. The base storm, following the example in the real land- scape, has a rainfall intensity of 5 mm h -1 falling over 2 h. To test the impacts of changing intensity and duration on model output, duration was extended compared to the base case (the "longer duration" storm; Table 4 ) and intensity was increased relative to the base storm (the "higher intensity" storm; Table 4 ). The storm with the longer duration maintained the 5 h -1 rainfall intensity, but duration was doubled to 4 h. In the higher intensity storm, rainfall rate was doubled to 10 mm h -1, while the base duration of 2 h was kept. 

Discharge was calculated at all grid locations during each model run. To capture the entire overland flow event, all simulations were run for 24 modeled hours, although flow had nearly stopped after 12 h of modeled time. A single base storm on the square watershed run for 24 modeled hours took approximately 80 s on a 2014 iMac with 4 GHz Intel Core i7 processors. 

The OverlandFlow results from the two test basins were coupled with the DetachmentLtdErosion component in Landlab to test the impact of non-steady hydrology on erosional patterns. At each time step, the DetachmentLtdErosion component calculated total incision depth at all points in the grid using Eq. ( 7 ). The initial condition for both test basins was topographic steady state, and so the predicted geomorphic "steady-state" incision rate was equal to the rock uplift rate applied in the model. Total incised depth for the hydrologic steady-state runs can be inferred from this steady-state incision rate. To test the erosional impact of non-steady hydrology, decadal simulations were run on each basin for the three precipitation events (Table 4 ). The known steady-state incision rate and depth can be compared to the predicted De-tachmentLtdErosion depth produced when coupled with the OverlandFlow component. In each basin, an annual precipitation rate of 0. 5 m yr -1 was set, and each simulation was run for 10 model years. Decadal-scale runs were selected, as they can be run quickly on a personal machine (on the order of hours), and the results can be used to make inferences about how erosion patterns would scale in long-term landscape evolution runs. Because of differences in intensity and duration, the base storm was run 500 times, assuming 50 storms per modeled year, while the longer duration and high intensity storms were run 250 times, assuming 25 storms per modeled year, to achieve 5 m total rainfall depth over 10 years. Cumulative incision depth at the end of each modeled run was saved at all points within the gridded terrain.
section_name: Methods: domain and parameterization

section_content: The hydrographs measured at the outlet of both the square and long basins are compared with the steady-state hydrographs (Fig. 10 ). The gray box represents the steady-state case, which produces the same discharge in both watersheds, as they have the same drainage area. In the non-steady method, hydrograph shapes are distinct between the different basins (Fig. 10a ). In the results from the base-case storm (Table 4), the hydrographs persist after precipitation and steadystate discharge end. In the case of the square basin, peak discharge exceeds that predicted by the steady-state case (∼ 50 m 3 s -1 ), a signal not seen in the long basin results. In the long basin, a singular peak discharge is not clear, and discharge values represented by the hydrograph are less than the predicted steady state at all time steps. Because flow in the long basin has to travel a greater distance from the upstream portion of the watershed, there is an elongated hydrograph with no clear peak discharge. 

As expected, the OverlandFlow component is also sensitive to changes in rainfall characteristics in both test basins. In the square basin, extending the duration of the storm (green line; Fig. 10b ) results in a higher overall peak discharge when compared to the base storm (light blue line; Fig. 10b ), as well as a longer overall hydrograph. The second peak in the longer duration hydrograph is due to the drainage organization in the square basin (Fig. 9a ), when flow from other tributaries reaches the outlet after the initial flood peak (see video in the Supplement). Increasing the rainfall intensity in the square basin (dark blue line; Fig. 10c ) increases peak discharge when compared to the base storm case. . OverlandFlow output for all storms described in Table 4. Hydrographs are taken from the active link upstream of the outlet node. Steady-state discharge is shown for each event, with the gray box representing the base storm in all cases. Panel (a) shows the base storm for both the square basin and the long basin; panel (b) compares outlet hydrographs from the base and longer duration storms in the square basin; panel (c) compares outlet hydrographs from the base and higher intensity storms in the square basin; panel (d) compares outlet hydrographs from the base and longer duration storms in the long basin; panel (e) compares outlet hydrographs from the base and higher intensity storms in the long basin. 

In the square basin, each storm has a clear hydrograph signature. These patterns are distinct from the long basin results. In the long basin, all three storm hydrographs have lower peak discharges than similar storms in the square basin (Fig. 10a ). The higher intensity storm run (mauve line; Fig. 10e ) has higher discharge values than both the base case and longer duration runs (Fig. 10d ), similar to what was seen in the square basin. However, the hydrograph shapes and discharge values are largely similar in all long basin cases, with longer, lower hydrographs that reflect the different travel time of water in the basin when compared to the square basin. 

To understand how non-steady hydrologic methods drive erosion in comparison to more traditional LEM methods, total incised depths for the three storm cases can be compared to predicted geomorphic steady-state incised depths after 10 modeled years. This application tests how the different hydrologic methods (steady vs. non-steady) impact morphology in LEM applications, following the work of Sólyom and Tucker (2004). The non-steady incision depth results demonstrate distinct patterns when compared to geomorphic steady state. Figure 11 shows that the coupled steady-state hydrology and stream power solutions predict higher incision rates than the non-steady method at all drainage areas. These patterns are clear in both the long watershed with a broad hydrograph and the square basin with a more peaked hydrograph. The depth of total incision in both basins is on the same order of magnitude, and the pattern of increasing incision depth downstream is also similar in both basins (Fig. 11a ). While the steady-state topography maintains the same land surface elevation, changing the hydrologic regime to nonsteady would lead to more relief in modeled landscapes, as the downstream will initially erode more rapidly than the upstream channels. In other words, the upstream locations will need to steepen more than the downstream locations in order to reach geomorphic steady-state incision rates throughout the landscape. Because the upstream locations must steepen more than the downstream locations in order to reach that geomorphic steady state, this will also lead to increased channel concavity on landscape evolution timescales. 

The pattern of increasing downstream incision is seen in all storm cases (Fig. 11b, c ). In both basins, total incised depth is least in the higher intensity storm, increases in the longer duration storm, and is greatest in the base case. The higher intensity storm exhibits a greater peak discharge in both basins, but there are fewer overall higher intensity and longer duration storms when compared to the base storm case to maintain the 5 m total rainfall depth over 10 years. Additionally, when calculating total incision using the stream power model, increases in discharge are less significant than the water surface slope due to the exponents m and n. While not explored here, changing the stream power exponents m and n will likely impact the steady and non-steady fluvial erosion results in this model, as would adding a threshold β to Eq. ( 7 ). 

Overall, these results suggest that, when compared to the OverlandFlow component, hydrologic steady-state predictions can over-or underestimate the peak of a hydrograph depending on basin orientation or shape (Fig. 10a ). As expected, the hydrodynamic algorithm from de Almeida et al. ( 2012 ) is sensitive to rainfall inputs, both with changes in duration and intensity (Fig. 10b-e ). This component can be applied across a range of timescales, used for predictions of overland flow for a single storm or multiple storms, and used efficiently with other process components in Landlab, as demonstrated by coupling to the DetachmentLtdErosion component. shows total incision results from the square basin for all three precipitation events after 10 years; and panel (c) shows total incision results from the long basin for all three precipitation events after 10 years. 

The patterns of erosion support earlier findings by Sólyom and Tucker (2004), which suggested that landscapes dominated by non-steady runoff patterns can be characterized by greater overall relief. Their results were generated using an incision rate controlled by the peak discharge. In contrast, the runs using the Landlab model were over shorter timescales, but these results were integrated over the entirety of the hydrograph, not just the peak discharge. These results suggest that, on longer timescales, watershed morphology would vary depending on the method used to calculate overland flow. Additionally, as the watershed morphology evolves in response to these spatial variations in incision rate, the hydrograph shape may change, impacting overall incision patterns and rates. The difference in patterns between steady and non-steady hydrology implies that flow patterns across a landscape during a runoff event, driven by non-steady hydrology, can have morphological significance over landscape evolution timescales.
section_name: Results and implications

section_content: The Landlab OverlandFlow model is flexible enough to be used in a number of scientific applications not discussed here. While the model does simulate surface flow over the entire domain, internally it makes no distinction between hillslope or channel processes, which can be problematic, as hillslopes make up the majority of a watershed area and supply sediment to the channels. If coupled with a hillslope sheet-wash component, OverlandFlow could be used to examine how non-steady channel processes interact with hillslope processes to sculpt watersheds across a range of spatial and temporal scales. Furthermore, these hillslope processes can be coupled with a fluvial transport-limited component and applied at event scales to explore sediment delivery from hillslopes to channels and how quickly sediment moves through a watershed. At landscape evolution timescales, evolved topographies resulting from more physically based hydrology and sediment transport components can be compared to traditional models to evaluate how physical parameters within the fluvial and hillslope models impact landscape relief and organization. 

Other opportunities include evaluating the impact of spatially variable parameters on model behavior. Spatial variability in rainfall could be explored with the development of new components that model orography or variability in storm cell size. Following the work of Huang and Niemann (2014), the OverlandFlow model can be used to explore patterns in runoff and erosion in response to changes in storm size, area, and location within a watershed. Spatially variable roughness could also be incorporated into the OverlandFlow component. A water-depth-dependent Manning's n method, similar to that of Rengers et al. (2016), could be implemented, where roughness at each grid node is calculated based on local wa-ter depths. Spatially variable roughness can also be used as input and set by the user based on field observations. Another potential application is coupling the Overland-Flow component to Landlab's ecohydrology components (Nudurupati et al., 2015). In this type of application, Over-landFlow could be used to calculate water depths across a surface. Surface water depths can be used to drive infiltration in the SoilInfiltrationGreenAmpt component. The Soil-Moisture component computes the water balance and rootzone soil moisture values. Soil moisture can drive changes in the Vegetation component, which simulates above-ground live and dead biomass. This coupled model would provide a more complete process ecohydrology model to be used in applications to understand how different flood events impact the succession of vegetation. 

Finally, the applications explored in this paper are on shorter timescales, ranging from event-to decadal-scale runs. An interesting future direction is exploring the OverlandFlow component in true landscape evolution runs (millennia or longer). Preliminary work modeling 10 3 to 10 4 years demonstrates that patterns seen in the decadal applications are clear; however, the full implications of hydrograph-driven erosion on longer timescales need to be further explored.
section_name: Future applications

section_content: This paper illustrates the theory behind the OverlandFlow component and how to use it as part of Landlab. Being part of the Landlab modeling framework comes with many advantages. The OverlandFlow component can make use of DEM input and output utilities and be coupled with other process components. Results from the real landscape application demonstrate that the OverlandFlow component can be used to route flow from observed rainfall events across a watershed DEM. This method can be used to estimate the grain sizes moved by real storm events and, in the future, could be coupled with other components and calibrated to understand the erosional response to flooding events. 

The OverlandFlow component can also be coupled to the DetachmentLtdErosion component to explore impacts of a hydrograph on erosion on decadal scales. In the synthetic landscapes explored here, the hydrograph results from the OverlandFlow component demonstrate a sensitivity to basin shape, precipitation duration, and intensity. The incision results predicted by using steady-state and non-steady water discharge are distinct in both the patterns and magnitudes of eroded depth and incision rates. Landscape evolution driven by non-steady runoff showed increasing incision rates moving downstream in the modeled watersheds. These results suggest that non-steady runoff could have important implications for predicting watershed relief and hypsometry in landscapes with different rainfall regimes and that choice of runoff method can have implications for both short-and long-term modeling results.
section_name: Conclusions

section_content: Acknowledgements. This research was supported by the National Science Foundation grants ACI-1147519 and ACI-1450338 (PI: Nicole M. Gasparini ), ACI-1148305 and ACI-1450412 ( PI: Erkan Istanbulluoglu ), ACI-1147454 (PI: Gregory E. Tucker ) and ACI-1450409 (PI: Gregory E. Tucker; co-PI: Daniel E. J. Hobley ), as well as the Tulane University Department of Earth and Environmental Sciences Vokes Fellowship ( Jordan M. Adams ). The authors are grateful to the topical editor Jeffrey Neal and reviewers Astrid Kerkweg, Katerina Michaelides, and Dapeng Yu, whose comments greatly improved the manuscript. 

Edited by: J. Neal Reviewed by: K. Michaelides and D. Yu

section_content: Code and data availability. The Landlab OverlandFlow and De-tachmentLtdErosion components are part of Landlab version 1. 0. 0. Source code for the Landlab project is housed on GitHub: http: //github. com/landlab/landlab. Documentation, installation instructions, and software dependencies for the entire Landlab project can be found at http://landlab. github. io/. A detailed user manual and driver scripts for the applications illustrated in this paper can be found at https://github. com/landlab/pub_adams_etal_gmd (Adams, 2016, GitHub Repository). The Landlab project is tested on recentgeneration Mac, Linux, and Windows platforms using Python versions 2. 7, 3. 4, and 3. 5. The Landlab modeling framework is distributed under a MIT open-source license. 

The Supplement related to this article is available online at doi:10. 5194/gmd-10-1645-2017-supplement.

section_content: Competing interests. The authors declare that they have no conflict of interest.