Note
This page was generated from a jupyter notebook.
The Implicit Kinematic Wave Overland Flow Component¶
Overview¶
This notebook demonstrates the KinwaveImplicitOverlandFlow Landlab component. The component implements a two-dimensional kinematic wave model of overland flow, using a digital elevation model or other source of topography as the surface over which water flows.
Theory¶
The kinematic wave equations are a simplified form of the 2D shallow-water equations in which energy slope is assumed to equal bed 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):
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 discharge depends on the local depth, bed-surface gradient \(\mathbf{S}=-\nabla\eta\) (this is the kinematic wave approximation; \(\eta\) is land surface height), and a roughness factor \(C_r\):
Reads may recognize this as a form of the Manning, Chezy, or Darcy-Weisbach equation. If \(\alpha = 5/3\) then we have the Manning equation, and \(C_r = n\) is “Manning’s n”. If \(\alpha = 3/2\) then we have the Chezy/Darcy-Weisbach equation, and \(C_r = 1/C = (f/8g)^{1/2}\) represents the Chezy roughness factor \(C\) and the equivalent Darcy-Weisbach factor \(f\).
Numerical solution¶
The solution method used by this component 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:
where \(H\) is water depth at a given grid node, \(t\) indicates time step number, \(\Delta t\) is time step duration, \(Q_{in}\) is total inflow discharge, \(Q_{out}\) is total outflow discharge, \(A\) is cell area, and \(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:
where \(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:
where \(N\) is the number of outflow faces (i.e., faces where the ground slopes downhill away from the cell’s node), and \(W_i\) is the width of face \(i\).
We use the depth at the cell’s node, so this simplifies to:
Notice that we know everything here except \(H'\). The reason we know \(Q_{out}\) is that it equals \(Q_{in}\) (which is either zero or we calculated it previously) plus \(RA\).
We define \(H\) in the above as a weighted sum of the “old” (time step \(t\)) and “new” (time step \(t+1\)) depth values:
If \(w=1\), the method is fully implicit. If \(w=0\), it is a simple forward explicit method.
When we combine these equations, we have an equation that includes the unknown \(H^{t+1}\) and a bunch of terms that are known. If \(w\ne 0\), it is a nonlinear equation in \(H^{t+1}\), and must be solved iteratively. We do this using a root-finding method in the scipy.optimize library.
In order to implement the algorithm, we must already know which of neighbors of each node are lower than the neighbor, and what the slopes between them are. We accomplish this using the FlowAccumulator and FlowDirectorMFD components. Running the FlowAccumulator also generates a sorted list (array) of nodes in drainage order.
The component¶
Import the needed libraries, then inspect the component’s docstring:
[ ]:
import matplotlib.pyplot as plt
import numpy as np
from tqdm.notebook import trange
from landlab import RasterModelGrid
from landlab.components.overland_flow import KinwaveImplicitOverlandFlow
from landlab.io import esri_ascii
[ ]:
print(KinwaveImplicitOverlandFlow.__doc__)
The docstring for the __init__ method will give us a list of parameters:
[ ]:
print(KinwaveImplicitOverlandFlow.__init__.__doc__)
Example 1: downpour on a plane¶
The first example tests that the component can reproduce the expected steady flow pattern on a sloping plane, with a gradient of \(S_p\). We’ll adopt the Manning equation. Once the system comes into equilibrium, the discharge should increase width distance down the plane according to \(q = Rx\). We can use this fact to solve for the corresponding water depth:
which implies
This is the analytical solution against which to test the model.
Pick the initial and run conditions
[ ]:
# Process parameters
n = 0.01 # roughness coefficient, (s/m^(1/3))
dep_exp = 5.0 / 3.0 # depth exponent
S = 0.01 # slope of plane
R = 72.0 # runoff rate, mm/hr
# Run-control parameters
run_time = 240.0 # duration of run, (s)
nrows = 5 # number of node rows
ncols = 11 # number of node columns
dx = 2.0 # node spacing, m
dt = 10.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")
# set topography
elev[grid.core_nodes] = S * (np.amax(grid.x_of_node) - grid.x_of_node[grid.core_nodes])
Plot topography, first in plan view…
[ ]:
grid.imshow(elev)
…then as a cross-section:
[ ]:
plt.plot(grid.x_of_node, elev)
plt.xlabel("Distance (m)")
plt.ylabel("Height (m)")
plt.grid(True)
[ ]:
# Instantiate the component
olflow = KinwaveImplicitOverlandFlow(
grid, runoff_rate=R, roughness=n, depth_exp=dep_exp
)
[ ]:
# 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
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
[ ]:
# Compare with analytical solution for depth
Rms = R / 3.6e6 # convert to m/s
nc = grid.number_of_node_columns
x = grid.x_of_node[grid.core_nodes][: nc - 2]
Hpred = (n * Rms * x / (S**0.5)) ** 0.6
plt.plot(x, Hpred, "r", label="Analytical")
plt.plot(
x,
grid.at_node["surface_water__depth"][grid.core_nodes][: nc - 2],
"b--",
label="Numerical",
)
plt.xlabel("Distance (m)")
plt.ylabel("Water depth (m)")
plt.grid(True)
plt.legend()
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.
[ ]:
# Process parameters
n = 0.1 # roughness coefficient, (s/m^(1/3))
dep_exp = 5.0 / 3.0 # depth exponent
R = 72.0 # runoff rate, mm/hr
# Run-control parameters
rain_duration = 240.0 # duration of rainfall, s
run_time = 480.0 # duration of run, s
dt = 10.0 # time-step size, s
dem_filename = "../hugo_site_filled.asc"
# Derived parameters
num_steps = int(run_time / dt)
# set up arrays to hold discharge and time
time_since_storm_start = np.arange(0.0, dt * (2 * num_steps + 1), dt)
discharge = np.zeros(2 * num_steps + 1)
[ ]:
# 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 surface_water_inflow__discharge at these 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).
[ ]:
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]}")
[ ]:
# Instantiate the component
olflow = KinwaveImplicitOverlandFlow(
grid, runoff_rate=R, roughness=n, depth_exp=dep_exp
)
[ ]:
discharge_field = grid.at_node["surface_water_inflow__discharge"]
for i in trange(num_steps):
olflow.run_one_step(dt)
discharge[i + 1] = np.sum(discharge_field[outlet_nodes])
[ ]:
plt.plot(time_since_storm_start[:num_steps], discharge[:num_steps])
plt.xlabel("Time (s)")
plt.ylabel("Discharge (cms)")
plt.grid(True)
[ ]:
grid.imshow(
grid.at_node["surface_water__depth"], cmap="Blues", colorbar_label="Water depth (m)"
)
Now turn down the rain and run it a bit longer…
[ ]:
olflow.runoff_rate = 1.0 # just 1 mm/hr
for i in trange(num_steps, 2 * num_steps):
olflow.run_one_step(dt)
discharge[i + 1] = np.sum(discharge_field[outlet_nodes])
[ ]:
plt.plot(time_since_storm_start, discharge)
plt.xlabel("Time (s)")
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 in the topography).