Note
This page was generated from a jupyter notebook.
2D Surface Water Flow component¶
Overview¶
This notebook demonstrate the usage of the river flow dynamics Landlab component. The component runs a semi-implicit, semi-Lagrangian finite-volume approximation to the depth-averaged 2D shallow-water equations of Casulli and Cheng (1992) and related work.
Theory¶
The depth-averaged 2D shallow-water equations are the simplification of the Navier-Stokes equations, which correspond to the balance of momentum and mass in the fluid. It is possible to simplify these equations by assuming a well-mixed water column and a small water depth to width ratio, where a vertical integration results in depth-averaged equations. These require boundary conditions at the top and bottom of the water column, which are provided by the wind stress and the Manning-Chezy formula, respectively:
where \(U\) is the water velocity in the \(x\)-direction, \(V\) is the water velocity in the \(y\)-direction, \(H\) is the water depth, \(\eta\) is the water surface elevation, \(Cz\) is the Chezy friction coefficient, and \(t\) is time. For the constants \(g\) is the gravity acceleration, \(\epsilon\) is the horizontal eddy viscosity, \(\mathbf{f}\) is the Coriolis parameter, \(\gamma_T\) is the wind stress coefficient, and \(U_a\) and \(V_a\) are the prescribed wind velocities.
Numerical representation¶
A semi-implicit, semi-Lagrangian, finite volume numerical approximation represents the depth averaged, 2D shallow-water equations described before. The water surface elevation, \(\eta\), is defined at the center of each computational volume (nodes). Water depth, \(H\), and velocity components, \(U\) and \(V\), are defined at the midpoint of volume faces (links). The finite volume structure provides a control volume representation that is inherently mass conservative.
The combination of a semi-implciit water surface elevation solution and a semi-Lagrangian representation of advection provides the advantages of a stable solution and of time steps that exceed the CFL criterion. In the semi-implicit process, \(\eta\) in the momentum equations, and the velocity divergence in the continuity equation, are treated implicitly. The advective terms in the momentum equations, are discretized explicitly. See the cited literature for more details.
The component¶
Import the needed libraries:
[ ]:
import matplotlib.pyplot as plt
import numpy as np
from IPython.display import clear_output
from tqdm import trange
from landlab import RasterModelGrid
from landlab.components import RiverFlowDynamics
from landlab.io import esri_ascii
Information about the component¶
Using the class name as argument for the help function returns descriptions of the various methods and parameters.
[ ]:
help(RiverFlowDynamics)
Examples¶
Example 1: Flow in a rectangular channel 6.0 m long¶
This first basic example illustrates water flowing through a rectangular channel 1.0 \(m\) wide and 6.0 \(m\) long. Our channel is made in concrete, so we choose a Manning’s roughness coefficient equal to 0.012 \(s/m^\frac{1}{3}\), and it has a slope of 0.01 \(m/m\).
We specify some basic parameters such as the grid resolution, time step duration, number of time steps, and the domain dimensions by specifying the number of columns and rows.
[ ]:
# Basic parameters
mannings_n = 0.012 # Manning's roughness coefficient, [s/m^(1/3)]
channel_slope = 0.01 # Channel slope [m/m]
# Simulation parameters
n_timesteps = 1000 # Number of timesteps
dt = 0.1 # Timestep duration, [s]
nrows = 20 # Number of node rows
ncols = 60 # Number of node cols
dx = 0.1 # Node spacing in the x-direction, [m]
dy = 0.1 # Node spacing in the y-direction, [m]
Create the grid:
[ ]:
# Create and set up the grid
grid = RasterModelGrid((nrows, ncols), xy_spacing=(dx, dy))
Create the elevation field and define the topography to represent our rectangular channel:
[ ]:
# The grid represents a basic rectangular channel with slope equal to 0.01 m/m
te = grid.add_field(
"topographic__elevation", 1.0 - channel_slope * grid.x_of_node, at="node"
)
te[grid.y_of_node > 1.5] = 2.5
te[grid.y_of_node < 0.5] = 2.5
We show a top view of the domain:
[ ]:
# Showing the topography
grid.imshow("topographic__elevation")
The channel is empty at the beginning of the simulation, so we create the fields for the water surface elevation, depth and velocity:
[ ]:
# We establish the initial conditions, which represent an empty channel
h = grid.add_zeros("surface_water__depth", at="node")
# Water velocity is zero in everywhere since there is no water yet
vel = grid.add_zeros("surface_water__velocity", at="link")
# Calculating the initial water surface elevation from water depth and topographic elevation
wse = grid.add_field("surface_water__elevation", te, at="node")
Then, we specify the nodes at which water is entering into the domain, and also the associated links. These are going to be the entry boundary conditions for water depth and velocity. In this case, water flows from left to right at 0.5 \(m\) depth, with a velocity of 0.45 \(m/s\):
[ ]:
# We set fixed boundary conditions, specifying the nodes and links in which the water is flowing into the grid
fixed_entry_nodes = np.array([300, 360, 420, 480, 540, 600, 660, 720, 780, 840, 900])
fixed_entry_links = grid.links_at_node[fixed_entry_nodes][:, 0]
# We set the fixed values in the entry nodes/links
entry_nodes_h_values = np.array([0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5])
entry_links_vel_values = np.array(
[0.45, 0.45, 0.45, 0.45, 0.45, 0.45, 0.45, 0.45, 0.45, 0.45, 0.45]
)
And now we show the boundary condition in the cross-section:
[ ]:
plt.plot(
grid.y_of_node[fixed_entry_nodes], entry_nodes_h_values + te[fixed_entry_nodes]
)
plt.plot(grid.y_of_node[grid.nodes_at_left_edge], te[grid.nodes_at_left_edge])
plt.title("Cross-section")
plt.xlabel("Distance [m]")
plt.ylabel("Elevation [m]")
plt.axis([0.25, 1.75, 0.75, 2.75])
plt.grid(True)
We construct our component by passing the arguments we defined previously:
[ ]:
# Finally, we run the model and let the water fill our channel
rfd = RiverFlowDynamics(
grid,
dt=dt,
mannings_n=mannings_n,
fixed_entry_nodes=fixed_entry_nodes,
fixed_entry_links=fixed_entry_links,
entry_nodes_h_values=entry_nodes_h_values,
entry_links_vel_values=entry_links_vel_values,
)
And finally, we run the simulation for 100 timesteps (10 seconds).
[ ]:
# Set the animation frequency to n_timesteps if you
# don't want to plot the water depth
# display_animation_freq = n_timesteps
display_animation_freq = 5
grid.imshow("surface_water__depth", output=True)
for timestep in trange(n_timesteps):
rfd.run_one_step()
if timestep % display_animation_freq == 0:
clear_output(wait=True) # This will clear the previous image
grid.imshow("surface_water__depth", output=True)
Exploring the water depth results at the latest time:
[ ]:
grid.imshow("surface_water__depth")
And the water surface elevation:
[ ]:
grid.imshow("surface_water__elevation")
Example 2: Surface water flowing over a DEM¶
On this case, we will import a digital elevation model (DEM) for a side-channel of the Kootenai River, Idaho, US.
[ ]:
# Getting the grid and some parameters
asc_file = "DEM-kootenai_37x50_1x1.asc"
with open(asc_file) as fp:
grid = esri_ascii.load(fp, name="topographic__elevation")
te = grid.at_node["topographic__elevation"]
Again, we specify some basic parameters such as the time step number and duration. For simplicity, we will keep our previous Manning’s coefficient. Notice that we already loaded all the required libraries.
[ ]:
# Basic parameters
mannings_n = 0.012 # Manning's roughness coefficient, [s/m^(1/3)]
# Simulation parameters
n_timesteps = 75 # Number of timesteps
dt = 1.0 # Timestep duration, [s]
Let’s see our new topography:
[ ]:
# Showing the topography
grid.imshow("topographic__elevation")
Our side-channel is empty at the beggining of the simulation, so we create the proper fields:
[ ]:
# We establish the initial conditions, which represent an empty channel
h = grid.add_zeros("surface_water__depth", at="node")
# Water velocity is zero in everywhere since there is no water yet
vel = grid.add_zeros("surface_water__velocity", at="link")
# Calculating the initial water surface elevation from water depth and topographic elevation
wse = grid.add_field("surface_water__elevation", te, at="node")
Then, we specify the nodes at which water is entering into the domain, and also the associated links. These are going to be our entry boundary conditions for water depth and velocity. On this case, water flows from right to left:
[ ]:
# We set fixed boundary conditions, specifying the nodes and links in which the water is flowing into the grid
fixed_entry_nodes = grid.nodes_at_right_edge
fixed_entry_links = grid.links_at_node[fixed_entry_nodes][:, 2]
# We set the fixed values in the entry nodes/links
entry_nodes_h_values = np.array(
[
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.04998779,
0.05999756,
0.03997803,
0.0,
0.0,
0.0,
0.05999756,
0.10998535,
0.12994385,
0.09997559,
0.15997314,
0.23999023,
0.30999756,
0.36999512,
0.45996094,
0.50994873,
0.54998779,
0.59997559,
0.63995361,
0.65997314,
0.65997314,
0.60998535,
0.5,
0.13995361,
0.0,
]
)
entry_links_vel_values = np.array(
[
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
-2.58638018,
-2.58638018,
-2.58638018,
0.0,
0.0,
0.0,
-2.58638018,
-2.58638018,
-2.58638018,
-2.58638018,
-2.58638018,
-2.58638018,
-2.58638018,
-2.58638018,
-2.58638018,
-2.58638018,
-2.58638018,
-2.58638018,
-2.58638018,
-2.58638018,
-2.58638018,
-2.58638018,
-2.58638018,
-2.58638018,
0.0,
]
)
Now we can plot our entry boundary condition in the cross-section:
[ ]:
plt.plot(
grid.y_of_node[fixed_entry_nodes], entry_nodes_h_values + te[fixed_entry_nodes]
)
plt.plot(grid.y_of_node[grid.nodes_at_right_edge], te[grid.nodes_at_right_edge])
plt.title("Entry cross-section")
plt.xlabel("Distance [m]")
plt.ylabel("Elevation [m]")
plt.grid(True)
Then we create the component by passing the arguments defined previously:
[ ]:
# Finally, we run the model and let the water fill our channel
rfd = RiverFlowDynamics(
grid,
dt=dt,
mannings_n=mannings_n,
fixed_entry_nodes=fixed_entry_nodes,
fixed_entry_links=fixed_entry_links,
entry_nodes_h_values=entry_nodes_h_values,
entry_links_vel_values=entry_links_vel_values,
)
And we run 75 time steps of 1 \(s\) duration (around 1 minute of computing time):
[ ]:
# Set the animation frequency to n_timesteps if you
# don't want to plot the water depth
# display_animation_freq = n_timesteps
display_animation_freq = 5
grid.imshow("surface_water__depth", output=True)
for timestep in trange(n_timesteps):
rfd.run_one_step()
if timestep % display_animation_freq == 0:
clear_output(wait=True) # This will clear the previous image
grid.imshow("surface_water__depth", output=True)
Finally, we can explore the results by plotting the resulting water depth:
[ ]:
grid.imshow("surface_water__depth")
### |
And that’s it! |
Nic |
e work completing
this tutorial.
You know now how
to use the
|