"""Landlab component for overland flow using a local implicit solution to the
kinematic-wave approximation.
Created on Fri May 27 14:26:13 2016
@author: gtucker
"""
import numpy as np
from scipy.optimize import newton
from landlab import Component
from landlab.components import FlowAccumulator
[docs]
def water_fn(x, a, b, c, d, e):
r"""Evaluates the solution to the water-depth equation.
Called by scipy.newton() to find solution for :math:`x`
using Newton's method.
Parameters
----------
x : float
Water depth at new time step.
a : float
"alpha" parameter (see below)
b : float
Weighting factor on new versus old time step. :math:`b=1` means purely
implicit solution with all weight on :math:`H` at new time
step. :math:`b=0` (not recommended) would mean purely explicit.
c : float
Water depth at old time step (time step :math:`t` instead
of :math:`t+1`)
d : float
Depth-discharge exponent; normally either 5/3 (Manning) or 3/2 (Chezy)
e : float
Water inflow volume per unit cell area in one time step.
This equation represents the implicit solution for water depth
:math:`H` at the next time step. In the code below, it is
formulated in a generic way. Written using more familiar
terminology, the equation is:
.. math::
H - H_0 + \alpha ( w H + (w-1) H_0)^d - \Delta t (R + Q_{in} / A)
.. math::
\alpha = \frac{\Delta t \sum S^{1/2}}{C_f A}
where :math:`H` is water depth at the given node at the new
time step, :math:`H_0` is water depth at the prior time step,
:math:`w` is a weighting factor, :math:`d` is the depth-discharge
exponent (2/3 or 1/2), :math:`\Delta t` is time-step duration,
:math:`R` is local runoff rate, :math:`Q_{in}` is inflow
discharge, :math:`A` is cell area, :math:`C_f` is a
dimensional roughness coefficient, and :math:`\sum S^{1/2}`
represents the sum of square-root-of-downhill-gradient over
all outgoing (downhill) links.
"""
return x - c + a * (b * x + (b - 1.0) * c) ** d - e
[docs]
class KinwaveImplicitOverlandFlow(Component):
r"""Calculate shallow water flow over topography.
Landlab component that implements a two-dimensional kinematic wave model.
This is a form of the 2D shallow-water equations in which energy slope is
assumed to equal bed slope. The solution method is locally implicit, and
works as follows. At each time step, we iterate from upstream to downstream
over the topography. Because we are working downstream, we can assume that
we know the total water inflow to a given cell. We solve the following mass
conservation equation at each cell:
.. math::
(H^{t+1} - H^t)/\Delta t = Q_{in}/A - Q_{out}/A + R
where :math:`H` is water depth, :math:`t` indicates time step
number, :math:`\Delta t` is time step duration, :math:`Q_{in}` is
total inflow discharge, :math:`Q_{out}` is total outflow
discharge, :math:`A` is cell area, and :math:`R` is local
runoff rate (precipitation minus infiltration; could be
negative if runon infiltration is occurring).
The specific outflow discharge leaving a cell along one of its faces is:
.. math::
q = (1/C_r) H^\alpha S^{1/2}
where :math:`C_r` is a roughness coefficient (such as
Manning's n), :math:`\alpha` is an exponent equal to :math:`5/3`
for the Manning equation and :math:`3/2` for the Chezy family,
and :math:`S` is the downhill-positive gradient of the link
that crosses this particular face. Outflow discharge is zero
for links that are flat or "uphill" from the given node.
Total discharge out of a cell is then the sum of (specific
discharge x face width) over all outflow faces
.. math::
Q_{out} = \sum_{i=1}^N (1/C_r) H^\alpha S_i^{1/2} W_i
where :math:`N` is the number of outflow faces (i.e., faces
where the ground slopes downhill away from the cell's node),
and :math:`W_i` is the width of face :math:`i`.
We use the depth at the cell's node, so this simplifies to:
.. math::
Q_{out} = (1/C_r) H'^\alpha \sum_{i=1}^N S_i^{1/2} W_i
We define :math:`H` in the above as a weighted sum of
the "old" (time step :math:`t`) and "new" (time step :math:`t+1`)
depth values:
.. math::
H' = w H^{t+1} + (1-w) H^t
If :math:`w=1`, the method is fully implicit. If :math:`w=0`,
it is a simple forward explicit method.
When we combine these equations, we have an equation that includes the
unknown :math:`H^{t+1}` and a bunch of terms that are known.
If :math:`w\ne 0`, it is a nonlinear equation in :math:`H^{t+1}`,
and must be solved iteratively. We do this using a root-finding
method in the scipy.optimize library.
Examples
--------
>>> from landlab import RasterModelGrid
>>> rg = RasterModelGrid((4, 5), xy_spacing=10.0)
>>> z = rg.add_zeros("topographic__elevation", at="node")
>>> kw = KinwaveImplicitOverlandFlow(rg)
>>> round(kw.runoff_rate * 1.0e7, 2)
2.78
>>> kw.vel_coef # default value
100.0
>>> rg.at_node["surface_water__depth"][6:9]
array([ 0., 0., 0.])
References
----------
**Required Software Citation(s) Specific to this Component**
None Listed
**Additional References**
None Listed
"""
_name = "KinwaveImplicitOverlandFlow"
_unit_agnostic = False
_info = {
"surface_water__depth": {
"dtype": float,
"intent": "out",
"optional": False,
"units": "m",
"mapping": "node",
"doc": "Depth of water on the surface",
},
"surface_water_inflow__discharge": {
"dtype": float,
"intent": "out",
"optional": False,
"units": "m3/s",
"mapping": "node",
"doc": "water volume inflow rate to the cell around each node",
},
"topographic__elevation": {
"dtype": float,
"intent": "in",
"optional": False,
"units": "m",
"mapping": "node",
"doc": "Land surface topographic elevation",
},
"topographic__gradient": {
"dtype": float,
"intent": "out",
"optional": False,
"units": "m/m",
"mapping": "link",
"doc": "Gradient of the ground surface",
},
}
[docs]
def __init__(
self,
grid,
runoff_rate=1.0,
roughness=0.01,
changing_topo=False,
depth_exp=1.5,
weight=1.0,
):
"""Initialize the KinwaveImplicitOverlandFlow.
Parameters
----------
grid : ModelGrid
Landlab ModelGrid object
runoff_rate : float, optional (defaults to 1 mm/hr)
Precipitation rate, mm/hr. The value provided is divided by
3600000.0.
roughness : float, defaults to 0.01
Manning roughness coefficient; units depend on depth_exp.
changing_topo : boolean, optional (defaults to False)
Flag indicating whether topography changes between time steps
depth_exp : float (defaults to 1.5)
Exponent on water depth in velocity equation (3/2 for Darcy/Chezy,
5/3 for Manning)
weight : float (defaults to 1.0)
Weighting on depth at new time step versus old time step (1 = all
implicit; 0 = explicit)
"""
super().__init__(grid)
# Store parameters and do unit conversion
self._runoff_rate = runoff_rate / 3600000.0 # convert to m/s
self._vel_coef = 1.0 / roughness # do division now to save time
self._changing_topo = changing_topo
self._depth_exp = depth_exp
self._weight = weight
# Get elevation field
self._elev = grid.at_node["topographic__elevation"]
# Create fields...
self.initialize_output_fields()
self._depth = grid.at_node["surface_water__depth"]
self._slope = grid.at_link["topographic__gradient"]
self._disch_in = grid.at_node["surface_water_inflow__discharge"]
# This array holds, for each node, the sum of sqrt(slope) x face width
# for each link/face.
self._grad_width_sum = grid.zeros("node")
# This array holds the prefactor in the algebraic equation that we
# will find a solution for.
self._alpha = grid.zeros("node")
# Instantiate flow router
self._flow_accum = FlowAccumulator(
grid,
"topographic__elevation",
flow_director="MFD",
partition_method="square_root_of_slope",
)
# Flag to let us know whether this is our first iteration
self._first_iteration = True
@property
def runoff_rate(self):
"""Runoff rate.
Parameters
----------
runoff_rate : float, optional (defaults to 1 mm/hr)
Precipitation rate, mm/hr. The value provide is divided by
3600000.0.
Returns
-------
The current value of the runoff rate.
"""
return self._runoff_rate
@runoff_rate.setter
def runoff_rate(self, new_rate):
assert new_rate > 0
self._runoff_rate = new_rate / 3600000.0 # convert to m/s
@property
def vel_coef(self):
"""Velocity coefficient."""
return self._vel_coef
@property
def depth(self):
"""The depth of water at each node."""
return self._depth
[docs]
def run_one_step(self, dt):
"""Calculate water flow for a time period `dt`."""
# If it's our first iteration, or if the topography may be changing,
# do flow routing and calculate square root of slopes at links
if self._changing_topo or self._first_iteration:
# Calculate the ground-surface slope
self._slope[self._grid.active_links] = self._grid.calc_grad_at_link(
self._elev
)[self._grid.active_links]
# Take square root of slope magnitude for use in velocity eqn
self._sqrt_slope = np.sqrt(np.abs(self._slope))
# Re-route flow, which gives us the downstream-to-upstream
# ordering
self._flow_accum.run_one_step()
self._nodes_ordered = self._grid.at_node["flow__upstream_node_order"]
self._flow_lnks = self._grid.at_node["flow__link_to_receiver_node"]
# (Re)calculate, for each node, sum of sqrt(gradient) x width
self._grad_width_sum[:] = 0.0
for i in range(self._flow_lnks.shape[1]):
self._grad_width_sum[:] += (
self._sqrt_slope[self._flow_lnks[:, i]]
* self._grid.length_of_face[
self._grid.face_at_link[self._flow_lnks[:, i]]
]
)
# Calculate values of alpha, which is defined as
#
# $\alpha = \frac{\Sigma W S^{1/2} \Delta t}{A C_r}$
cores = self._grid.core_nodes
self._alpha[cores] = (
self._vel_coef
* self._grad_width_sum[cores]
* dt
/ (self._grid.area_of_cell[self._grid.cell_at_node[cores]])
)
# Zero out inflow discharge
self._disch_in[:] = 0.0
# Upstream-to-downstream loop
for i in range(len(self._nodes_ordered) - 1, -1, -1):
n = self._nodes_ordered[i]
if self._grid.status_at_node[n] == 0:
# Solve for new water depth
aa = self._alpha[n]
cc = self._depth[n]
ee = (dt * self._runoff_rate) + (
dt
* self._disch_in[n]
/ self._grid.area_of_cell[self._grid.cell_at_node[n]]
)
self._depth[n] = newton(
water_fn,
self._depth[n],
args=(aa, self._weight, cc, self._depth_exp, ee),
)
# Calc outflow
Heff = self._weight * self._depth[n] + (1.0 - self._weight) * cc
outflow = (
self._vel_coef * (Heff**self._depth_exp) * self._grad_width_sum[n]
) # this is manning/chezy/darcy
# Send flow downstream. Here we take total inflow discharge
# and partition it among the node's neighbors. For this, we use
# the flow director's "proportions" array, which contains, for
# each node, the proportion of flow that heads out toward each
# of its N neighbors. The proportion is zero if the neighbor is
# uphill; otherwise, it is S^1/2 / sum(S^1/2). If for example
# we have a raster grid, there will be four neighbors and four
# proportions, some of which may be zero and some between 0 and
# 1.
self._disch_in[self._grid.adjacent_nodes_at_node[n]] += (
outflow * self._flow_accum.flow_director._proportions[n]
)
# TODO: the above is enough to implement the solution for flow
# depth, but it does not provide any information about flow
# velocity or discharge on links. This could be added as an
# optional method, perhaps done just before output.
if __name__ == "__main__":
import doctest
doctest.testmod()