Note

This page was generated from a jupyter notebook.

The Linear Diffusion Overland Flow Router

Overview

This notebook demonstrates the LinearDiffusionOverlandFlowRouter Landlab component. The component implements a two-dimensional model of overland flow, based on a linearization of the diffusion-wave approximation of the shallow-water equations.

Theory

Flow direction, depth, and velocity

The diffusion-wave equations are a simplified form of the 2D shallow-water equations in which energy slope is assumed to equal water-surface slope. Conservation of water mass is expressed in terms of the time derivative of the local water depth, \(H\), and the spatial derivative (divergence) of the unit discharge vector \(\mathbf{q} = UH\) (where \(U\) is the 2D depth-averaged velocity vector):

\[\frac{\partial H}{\partial t} = R - \nabla\cdot \mathbf{q}\]

where \(R\) is the local runoff rate [L/T] and \(\mathbf{q}\) has dimensions of volume flow per time per width [L\(^2\)/T]. The flow velocity is calculated using a linearized form of the Manning friction law:

\[\mathbf{U} = \frac{H^{4/3}}{n^2 u_c} \nabla w\]

\[w = \eta + H\]

Here \(\eta(x,y,t)\) is ground-surface elevation, \(w(x,y,t)\) is water-surface elevation, \(n\) is the Manning friction coefficient, and \(u_c\) is a characteristic scale velocity (see, e.g., Mariotti, 2018). Thus, there are two parameters governing flow speed: \(n\) and \(u_c\). The may, however, be treated as a single lumped parameter \(n^2 u_c\).

Rainfall and infiltration

Runoff rate is calculated as the difference between the rates of precipitation, \(P\), and infiltration, \(I\). The user specifies a precipitation rate (which is a public variable that can be modified after instantiation), and a maximum infiltration rate, \(I_c\). The actual infiltration rate depends on the available surface water, and is calculated in a way that allows it to approach zero as the surface-water depth approaches zero:

\[I = I_c \left( 1 - e^{-H/H_i} \right)\]

where \(H_i\) is a characteristic water depth, defined such that the actual infiltration rate is about 95% of \(I_c\) when \(H = 3 H_i\).

Numerical Methods

Finite-volume representation

The component uses an explicit, forward-Euler finite-volume method. The solution for water depth at a new time step \(k+1\) is calculated from:

\[H_i^{k+1} = H_i^k + \Delta t \left[ \frac{dH_i}{dt} \right]_i^k\]

The time derivative at step \(k\) is calculated as:

\[\left[ \frac{dH_i}{dt} \right]_i^k = R - I - \frac{1}{\Lambda_i} \sum_{j=1}^{N_i} \lambda_{ij} q_{ij}\]

where \(R\) is rainfall rate, \(I\) is infiltration rate (calculated as above), \(\Lambda_i\) is the horizontal surface area of the cell enclosing node \(i\), \(\lambda_{ij}\) is the length of the cell face between node \(i\) and its \(j\)-th neighbor, and \(q_{ij}\) is the specific water discharge along the link that connects node \(i\) and its \(j\)-th neighbor.

For a raster grid, this treatment is equivalent to a centered-in-space finite-difference arrangement. For more on finite-difference solutions to diffusion problems, see for example Slingerland and Kump (2011) and Press et al. (1986).

Time-step limiter

Because of the linearization described above, the flow model is effectively a diffusion problem with a space- and time-varying diffusivity, because the effective diffusivity \(D\) depends on water depth:

\[D = \frac{H^{7/3}}{n^2 u_c}\]

One numerical challenge is that, according to the Courant-Levy-Friedrichs (CFL) criterion, the maximum stable time step will depend on water depth, which varies in space and time. To prevent instability, the solution algorithm calculates at every iteration a maximum value for \(D\) using the current maximum water depth (or a very small minimum value, whichever is larger, to prevent blowup in the case of zero water depth). The maximum step size is then calculated as:

\[\Delta t_\text{max} = \alpha \frac{L_\text{min}^2}{2 D}\]

where \(L_\text{min}\) is the length of the shortest link in the grid (which is just equal to node spacing, \(\Delta x\), for a uniform raster or hex grid). The stability factor \(\alpha\) is a user-controllable parameter that defaults to 0.2, and must be \(\le 1\).

If \(\Delta t_\text{max}\) is less than the user-specified “global” time step duration, the algorithm iterates repeatedly with time steps of size \(\Delta t_\text{max}\) (or the remaining time in the global step, whichever is smaller) until the global time step duration has been completed.

The component

Import the needed libraries, then inspect the component’s docstring:

import matplotlib.pyplot as plt
import numpy as np

from landlab import RasterModelGrid
from landlab.components.overland_flow import LinearDiffusionOverlandFlowRouter
from landlab.io import esri_ascii

Use the help function to get a description of the LinearDiffusionOverlandFlowRouter component. If you scroll down to the __init__ section, you will see a list of parameters.

help(LinearDiffusionOverlandFlowRouter)

Help on class LinearDiffusionOverlandFlowRouter in module landlab.components.overland_flow.linear_diffusion_overland_flow_router:

class LinearDiffusionOverlandFlowRouter(landlab.core.model_component.Component)
 |  LinearDiffusionOverlandFlowRouter(grid, roughness=0.01, rain_rate=1e-05, infilt_rate=0.0, infilt_depth_scale=0.001, velocity_scale=1.0, cfl_factor=0.2)
 |
 |  Calculate water flow over topography.
 |
 |  Landlab component that implements a two-dimensional, linearized
 |  diffusion-wave model. The diffusion-wave approximation is a simplification
 |  of the shallow-water equations that omits the momentum terms. The flow
 |  velocity is calculated using the local water-surface slope as an
 |  approximation of the energy slope. With this linearized form, flow velocity
 |  is calculated using a linearized Manning equation, with the water-surface
 |  slope being used as the slope factor. There are two governing equations, one
 |  that represents conservation of water mass:
 |
 |  ..math::
 |
 |      \frac{\partial H}{\partial t} = (P - I) - \nabla\cdot\mathbf{q}
 |
 |  where :math:`H(x,y,t)` is local water depth, :math:`t` is time, :math:`P`
 |  is precipitation rate, :math:`I` is infiltration rate, and :math:`\mathbf{q}`
 |  is specific water discharge, which equals depth times depth-averaged
 |  velocity. The other governing equation represents specific discharge as a
 |  function of gravity, pressure, and friction:
 |
 |  ..math::
 |
 |      \mathbf{q} = \frac{H^{4/3}}{n^2 U_c} \nabla w
 |
 |  where :math:`n` is the friction factor ("Manning's n"), :math:`U_c` is a
 |  characteristic flow velocity, and :math:`w` is water-surface height, which
 |  is the sum of topographic elevation plus water depth.
 |
 |  Infiltration rate should decline smoothly to zero as surface water depth
 |  approaches zero. To ensure this, infiltration rate is calculated as
 |
 |  ..math::
 |
 |      I = I_c \left( 1 - e^{-H / H_i} ) \right)
 |
 |  where :math:`H_i` is a characteristic water depth. The concept here is that
 |  when :math:`H \le H_i`, small spatial variations in water depth will leave
 |  parts of the ground un-ponded and therefore not subject to any infiltration.
 |  Mathematically, :math:`I \approx 0.95 I_c` when :math:`H = 3H_i`.
 |
 |  Examples
 |  --------
 |  >>> from landlab import RasterModelGrid
 |  >>> grid = RasterModelGrid((3, 3))
 |  >>> _ = grid.add_zeros("topographic__elevation", at="node")
 |  >>> olf = LinearDiffusionOverlandFlowRouter(grid, roughness=0.1)
 |  >>> round(olf.vel_coef)
 |  100
 |  >>> olf.rain_rate
 |  1e-05
 |  >>> olf.rain_rate = 1.0e-4
 |  >>> olf.run_one_step(dt=10.0)
 |  >>> grid.at_node["surface_water__depth"][4]
 |  0.001
 |
 |  References
 |  ----------
 |  **Required Software Citation(s) Specific to this Component**
 |
 |  None Listed
 |
 |  **Additional References**
 |
 |  None Listed
 |
 |  Method resolution order:
 |      LinearDiffusionOverlandFlowRouter
 |      landlab.core.model_component.Component
 |      builtins.object
 |
 |  Methods defined here:
 |
 |  __init__(self, grid, roughness=0.01, rain_rate=1e-05, infilt_rate=0.0, infilt_depth_scale=0.001, velocity_scale=1.0, cfl_factor=0.2)
 |      Initialize the LinearDiffusionOverlandFlowRouter.
 |
 |      Parameters
 |      ----------
 |      grid : ModelGrid
 |          Landlab ModelGrid object
 |      roughness : float, defaults to 0.01
 |          Manning roughness coefficient, s/m^1/3
 |      rain_rate : float, optional (defaults to 36 mm/hr)
 |          Rainfall rate, m/s
 |      infilt_depth_scale : float, defaults to 0.001
 |          Depth scale for water infiltration, m
 |      infilt_rate : float, optional (defaults to 0)
 |          Maximum rate of infiltration, m/s
 |      velocity_scale : float, defaults to 1
 |          Characteristic flow velocity, m/s
 |      cfl_factor : float, optional (defaults to 0.2)
 |          Time-step control factor: fraction of maximum estimated time-step
 |          that is actually used (must be <=1)
 |
 |  run_one_step(self, dt)
 |      Calculate water flow for a time period `dt`.
 |
 |      Default units for dt are *seconds*. We use a time-step subdivision
 |      algorithm that ensures step size is always below CFL limit.
 |
 |  update_for_one_iteration(self, iter_dt)
 |      Update state variables for one iteration of duration iter_dt.
 |
 |  ----------------------------------------------------------------------
 |  Readonly properties defined here:
 |
 |  vel_coef
 |      Velocity coefficient.
 |
 |      (1/(roughness^2 x velocity_scale)
 |
 |  ----------------------------------------------------------------------
 |  Data descriptors defined here:
 |
 |  rain_rate
 |      Rainfall rate
 |
 |  ----------------------------------------------------------------------
 |  Methods inherited from landlab.core.model_component.Component:
 |
 |  initialize_optional_output_fields(self)
 |      Create fields for a component based on its optional field outputs,
 |      if declared in _optional_var_names.
 |
 |      This method will create new fields (without overwrite) for any
 |      fields output by the component as optional. New fields are
 |      initialized to zero. New fields are created as arrays of floats,
 |      unless the component also contains the specifying property
 |      _var_type.
 |
 |  initialize_output_fields(self, values_per_element=None)
 |      Create fields for a component based on its input and output var
 |      names.
 |
 |      This method will create new fields (without overwrite) for any fields
 |      output by, but not supplied to, the component. New fields are
 |      initialized to zero. Ignores optional fields. New fields are created as
 |      arrays of floats, unless the component specifies the variable type.
 |
 |      Parameters
 |      ----------
 |      values_per_element: int (optional)
 |          On occasion, it is necessary to create a field that is of size
 |          (n_grid_elements, values_per_element) instead of the default size
 |          (n_grid_elements,). Use this keyword argument to acomplish this
 |          task.
 |
 |  ----------------------------------------------------------------------
 |  Class methods inherited from landlab.core.model_component.Component:
 |
 |  from_path(grid, path)
 |      Create a component from an input file.
 |
 |      Parameters
 |      ----------
 |      grid : ModelGrid
 |          A landlab grid.
 |      path : str or file_like
 |          Path to a parameter file, contents of a parameter file, or
 |          a file-like object.
 |
 |      Returns
 |      -------
 |      Component
 |          A newly-created component.
 |
 |  var_definition(name)
 |      Get a description of a particular field.
 |
 |      Parameters
 |      ----------
 |      name : str
 |          A field name.
 |
 |      Returns
 |      -------
 |      tuple of (*name*, *description*)
 |          A description of each field.
 |
 |  var_help(name)
 |      Print a help message for a particular field.
 |
 |      Parameters
 |      ----------
 |      name : str
 |          A field name.
 |
 |  var_loc(name)
 |      Location where a particular variable is defined.
 |
 |      Parameters
 |      ----------
 |      name : str
 |          A field name.
 |
 |      Returns
 |      -------
 |      str
 |          The location ('node', 'link', etc.) where a variable is defined.
 |
 |  var_type(name)
 |      Returns the dtype of a field (float, int, bool, str...).
 |
 |      Parameters
 |      ----------
 |      name : str
 |          A field name.
 |
 |      Returns
 |      -------
 |      dtype
 |          The dtype of the field.
 |
 |  var_units(name)
 |      Get the units of a particular field.
 |
 |      Parameters
 |      ----------
 |      name : str
 |          A field name.
 |
 |      Returns
 |      -------
 |      str
 |          Units for the given field.
 |
 |  ----------------------------------------------------------------------
 |  Static methods inherited from landlab.core.model_component.Component:
 |
 |  __new__(cls, *args, **kwds)
 |      Create and return a new object.  See help(type) for accurate signature.
 |
 |  ----------------------------------------------------------------------
 |  Readonly properties inherited from landlab.core.model_component.Component:
 |
 |  cite_as
 |      Citation information for component.
 |
 |      Return required software citation, if any. An empty string indicates
 |      that no citations other than the standard Landlab package citations are
 |      needed for the component.
 |
 |      Citations are provided in BibTeX format.
 |
 |      Returns
 |      -------
 |      cite_as
 |
 |  coords
 |      Return the coordinates of nodes on grid attached to the
 |      component.
 |
 |  definitions
 |      Get a description of each field.
 |
 |      Returns
 |      -------
 |      tuple of (*name*, *description*)
 |          A description of each field.
 |
 |  grid
 |      Return the grid attached to the component.
 |
 |  input_var_names
 |      Names of fields that are used by the component.
 |
 |      Returns
 |      -------
 |      tuple of str
 |          Tuple of field names.
 |
 |  name
 |      Name of the component.
 |
 |      Returns
 |      -------
 |      str
 |          Component name.
 |
 |  optional_var_names
 |      Names of fields that are optionally provided by the component, if
 |      any.
 |
 |      Returns
 |      -------
 |      tuple of str
 |          Tuple of field names.
 |
 |  output_var_names
 |      Names of fields that are provided by the component.
 |
 |      Returns
 |      -------
 |      tuple of str
 |          Tuple of field names.
 |
 |  shape
 |      Return the grid shape attached to the component, if defined.
 |
 |  unit_agnostic
 |      Whether the component is unit agnostic.
 |
 |      If True, then the component is unit agnostic. Under this condition a
 |      user must still provide consistent units across all input arguments,
 |      keyword arguments, and fields. However, when ``unit_agnostic`` is True
 |      the units specified can be interpreted as dimensions.
 |
 |      When False, then the component requires inputs in the specified units.
 |
 |      Returns
 |      -------
 |      bool
 |
 |  units
 |      Get the units for all field values.
 |
 |      Returns
 |      -------
 |      tuple or str
 |          Units for each field.
 |
 |  var_mapping
 |      Location where variables are defined.
 |
 |      Returns
 |      -------
 |      tuple of (name, location)
 |          Tuple of variable name and location ('node', 'link', etc.) pairs.
 |
 |  ----------------------------------------------------------------------
 |  Data descriptors inherited from landlab.core.model_component.Component:
 |
 |  __dict__
 |      dictionary for instance variables
 |
 |  __weakref__
 |      list of weak references to the object
 |
 |  current_time
 |      Current time.
 |
 |      Some components may keep track of the current time. In this case, the
 |      ``current_time`` attribute is incremented. Otherwise it is set to None.
 |
 |      Returns
 |      -------
 |      current_time

Example 1: downpour on a single cell

The first example tests that the component can reproduce the expected steady flow depth for rain on a single cell. The input water flow rate, in \(m^3 / s\), is:

\[Q_\text{in} = P \Delta x^2\]

The output flow rate is

\[Q_\text{out} = \frac{\Delta x}{n^2 u_c} H^{7/3} S_w\]

where \(S_w\) is the water surface slope. We can write the water-surface slope in terms of the water height of the (one) core node (which just equals \(H\), because the ground elevation is zero) and the water height at the adjacent open-boundary node, which is zero, so

\[S_w = \frac{H}{\Delta x}\]

We can therefore plug this into the equation for \(Q_\text{out}\) and solve for the expected equilibrium depth:

\[H = \left(\Delta x^2 P n^2 u_c \right)^{3/10}\]

Pick the initial and run conditions

# Process parameters
n = 0.01  # roughness coefficient, (s/m^(1/3))
vel_scale = 1.0  # velocity scale, m/s
R = 72.0 / (3600.0 * 1000.0)  # runoff rate, m/s

# Run-control parameters
run_time = 360.0  # duration of run, (s)
nrows = 3  # number of node rows
ncols = 3  # number of node columns
dx = 2.0  # node spacing, m
dt = 20.0  # time-step size, s
plot_every = 60.0  # plot interval, s

# Derived parameters
num_steps = int(run_time / dt)

Create grid and fields:

# create and set up grid
grid = RasterModelGrid((nrows, ncols), xy_spacing=dx)
grid.set_closed_boundaries_at_grid_edges(False, True, True, True)  # open only on east

# add required field
elev = grid.add_zeros("topographic__elevation", at="node")

# Instantiate the component
olflow = LinearDiffusionOverlandFlowRouter(
    grid, rain_rate=R, roughness=n, velocity_scale=vel_scale
)

# Helpful function to plot the profile


def plot_flow_profile(grid, olflow):
    """Plot the middle row of topography and water surface
    for the overland flow model olflow."""
    nc = grid.number_of_node_columns
    nr = grid.number_of_node_rows
    startnode = nc * (nr // 2) + 1
    midrow = np.arange(startnode, startnode + nc - 1, dtype=int)
    topo = grid.at_node["topographic__elevation"]
    plt.plot(
        grid.x_of_node[midrow],
        topo[midrow] + grid.at_node["surface_water__depth"][midrow],
        "b",
    )
    plt.plot(grid.x_of_node[midrow], topo[midrow], "k")
    plt.xlabel("Distance (m)")
    plt.ylabel("Ground and water surface height (m)")

Run the component forward in time, plotting the output in the form of a profile:

next_plot = plot_every
HH = []  # keep track of depth through time
for i in range(num_steps):
    olflow.run_one_step(dt)
    if (i + 1) * dt >= next_plot:
        plot_flow_profile(grid, olflow)
        next_plot += plot_every
    HH.append(grid.at_node["surface_water__depth"][4])

# Compare with analytical solution for depth
expected_depth = (dx * dx * R * n * n * vel_scale) ** 0.3
computed_depth = grid.at_node["surface_water__depth"][4]

print(f"Expected depth = {expected_depth} m")
print(f"Computed depth = {computed_depth} m")

Expected depth = 0.0037232911332721395 m
Computed depth = 0.0037082204564142748 m

plt.plot(np.linspace(0, run_time, len(HH)), HH)
plt.xlabel("Time (s)")
plt.ylabel("Water depth (m)")

Text(0, 0.5, 'Water depth (m)')

Example 2: overland flow on a DEM

For this example, we’ll import a small digital elevation model (DEM) for a site in New Mexico, USA, with 10 m cells.

# Process parameters
n = 0.1  # roughness coefficient, (s/m^(1/3))
uc = 1.0  # characteristic velocity scale (m/s)
R1 = 72.0 / (3600.0 * 1000.0)  # initial rainfall rate, m/s (converted from mm/hr)
R2 = 0.0 / (3600.0 * 1000.0)  # later rainfall rate, m/s (converted from mm/hr)
infilt_cap = 10.0 / (3600 * 1000.0)  # infiltration capacity, m/s (converted from mm/hr)

# Run-control parameters
heavy_rain_duration = 300.0  # duration of heavy rainfall, s
run_time = 1200.0  # duration of run, s
dt = 20.0  # time-step size, s
dem_filename = "../hugo_site.asc"

# Derived parameters
num_steps = int(run_time / dt)

# set up arrays to hold discharge, time, and rain rate
time_since_storm_start = np.linspace(0.0, run_time, num_steps + 1)
discharge = np.zeros(num_steps + 1)
rain_rate = np.zeros(num_steps + 1)
rain_rate[:] = R1
rain_rate[time_since_storm_start >= heavy_rain_duration] = R2

# Read the DEM file as a grid with a 'topographic__elevation' field
with open(dem_filename) as fp:
    grid = esri_ascii.load(fp, name="topographic__elevation", at="node")
elev = grid.at_node["topographic__elevation"]

# Configure the boundaries: valid right-edge nodes will be open;
# all NODATA (= -9999) nodes will be closed.
grid.status_at_node[grid.nodes_at_right_edge] = grid.BC_NODE_IS_FIXED_VALUE
grid.status_at_node[np.isclose(elev, -9999.0)] = grid.BC_NODE_IS_CLOSED

# display the topography
grid.imshow(elev, colorbar_label="Elevation (m)", cmap="pink")

It would be nice to track discharge at the watershed outlet, but how do we find the outlet location? We actually have several valid nodes along the right-hand edge. Then we’ll keep track of the field water__specific_discharge at the active links that connect to these boundary nodes. We can identify the nodes by the fact that they are (a) at the right-hand edge of the grid, and (b) have positive elevations (the ones with -9999 are outside of the watershed). We can identify the relevant active links as those connected to the outlet nodes that have active status (meaning they do not connect to any closed boundary nodes).

indices = np.where(elev[grid.nodes_at_right_edge] > 0.0)[0]
outlet_nodes = grid.nodes_at_right_edge[indices]
print(f"Outlet nodes: {outlet_nodes}")
print(f"Elevations of the outlet nodes: {elev[outlet_nodes]}")

links_at_outlets = grid.links_at_node[outlet_nodes]
links_to_track = links_at_outlets[
    grid.status_at_link[links_at_outlets] == grid.BC_LINK_IS_ACTIVE
].flatten()
print(f"Links at which to track discharge: {links_to_track}")

Outlet nodes: [1975 2051 2127 2203 2279 2355 2431 2507]
Elevations of the outlet nodes: [1661. 1660. 1661. 1662. 1665. 1669. 1672. 1676.]
Links at which to track discharge: [3849 4000 4151 4302 4453 4604 4755 4906]

# Instantiate the component
olflow = LinearDiffusionOverlandFlowRouter(
    grid, rain_rate=R, infilt_rate=infilt_cap, roughness=n, velocity_scale=vel_scale
)

def cfl():
    hmax = np.amax(grid.at_node["surface_water__depth"])
    D = hmax ** (7.0 / 3.0) / (n * n * uc)
    return 0.5 * dx * dx / D

q = grid.at_link["water__specific_discharge"]

for i in range(num_steps):
    olflow.rain_rate = rain_rate[i]
    olflow.run_one_step(dt)
    discharge[i + 1] = np.sum(q[links_to_track]) * dx

plt.plot(time_since_storm_start / 60.0, discharge)
plt.xlabel("Time (min)")
plt.ylabel("Discharge (cms)")
plt.grid(True)

grid.imshow(
    grid.at_node["surface_water__depth"], cmap="Blues", colorbar_label="Water depth (m)"
)

Voila! A fine hydrograph, and a water-depth map that shows deeper water in the channels (and highlights depressions and low-gradient spots in the topography).

References

Mariotti, G. (2018). Marsh channel morphological response to sea level rise and sediment supply. Estuarine, Coastal and Shelf Science, 209, 89-101.

Press, W. H., Vetterling, W. T., Teukolsky, S. A., & Flannery, B. P. (1986). Numerical recipes. Cambridge: Cambridge university press.

Slingerland, R., & Kump, L. (2011). Mathematical Modeling of Earth’s Dynamical Systems. Princeton University Press.

---

Generated by nbsphinx from a Jupyter notebook.