This model describes the surface dynamics of an enclosed water-body following an initial perturbation. It is an implementation of the "bathtub" model described in The Dynamics of Coastal Models by Clifford J. Hearn (Cambridge).
The model is 1-dimensional and therefore has only one spatial dimension in which flow needs to be resolved (i.e. $x$ rather than $x$, $y$, and $z$). We don't have to deal with the Coriolis force nor, since flow is depth-averaged, variations with depth.
In the model, water moves because of two forces: pressure and friction. The pressure force occurs because of horizontal variations in the water depth which produce a horizontal pressure gradient causing water to flow horizontally. The friction force acts to slow down the flow in response to friction at the seabed.
When the water moves, the water surface changes. This is because water flow removes water volume from one place and adds it to another (we assume water is not created or destroyed, but is "conserved". It is also not compressed). Therefore, there is a feedback between flow and water level: the flow alters the water surface, and the resulting new water surface and its' associated pressure gradient modify the flow. The model tracks these mutually dependent processes through time.
Water flow is described by the depth-averaged momentum equation:
$$ \frac{\partial u}{\partial t} = -g\frac{\partial \eta}{\partial x} - \frac{u}{\tau H}$$where $u$ is the depth-averaged current velocity (m3 m-2 s-1 or m s-1), $g$ is the acceleration due to gravity, $\eta$ (greek letter eta) is the water level above some datum (m), $H$ is the total water depth (m), and $\tau$ (greek letter tau) is a linear frictional damping time (i.e. the reciprocal of a friction coefficient). The first term on the right-hand side is the force due to the pressure gradient and the second term represents friction. The $H$ quotient in the friction term effectively averages the bed friction over the entire water column.
The response of the water surface is described by the depth-averaged volume continuity equation:
$$ \frac{\partial \eta}{\partial t} = -\frac{\partial Hu}{\partial x}$$This equation basically conserves water volume. It states that the rate of water level change at a given point in space (left-hand side) is equal to the volume flux at that point (right-hand side), i.e. the gradient of the depth-average current velocity multiplied by the water depth.
First, we'll start off by importing the required numerical and plotting libraries. numpy allows us to operate on vectors and matrices. matplotlib gives us plotting functionality.
In [102]:
import matplotlib.pyplot as plt
import numpy as np
# plot in this page rather than in separate window
%matplotlib inline
Now, let's set up the numerical scheme. We'll divide the space domain of the model into a number of discrete cells (in the $x$ direction), and we'll divide time up into discrete steps. Eventually, we'll iterate through each timestep and for each cell in the spatial domain we'll solve the depth-averaged equations, determining how the water moves in that time step and how the water surface responds.
We can adjust the numerical parameters as desired, but we need to ensure that the timestep is sufficiently small relative to the distance step so that water does not flow through a whole cell during a single timestep. Such a case would cause the model to become unstable. To ensure this, we set $dt$ as a function of $dx$. This means that is we decrease $dx$, $dt$ automatically decreases too. See this Wikipedia page for more information on this constraint.
In [130]:
N = 25 # number of cells in the x direction
dx = 1.0/N # distance between cells, i.e. size of distance step
dt = 0.05*dx # size of time step
Next, we must specify some model parameters. This model differs from that of Hearn by explicitly including the acceleration due to gravity in the pressure term of the depth-averaged momentum equation (the Hearn implementation omits this by converting to non-dimensional units). So we set the value for $g$,
In [131]:
g = 9.8
and a value for the frictional damping time ($\tau$),
In [132]:
tau = 0.05
Next we have to specify how the water level is described. We call the water depth at any given location in the model domain, confusingly, the water height, i.e. the height of the water surface from the bed. So the variation in water height across the model domain at any point in time defines the water surface at that time.
Instead of using a single variable for water height, we split the total water height up in to two components: (1) the height to a fixed water level datum; and (2) the height of the water surface relative to this datum. The reason for doing this is that the water level needs to be measured relative to a horizontal plane if differences in water height across space (and therefore the pressure gradient) are to be resolved. If we simply used the total water depth to measure differences in water level across space we would be implicitly defining the seabed as the reference plane. But the seabed is scarcely completely flat and in the case of a non-flat seabed it would not be clear whether changes in water depth were due to variations in the seabed or variations in the water surface. We're going to assume a flat seabed in this model, but it is good practice to used an independent, fixed, horizontal datum against which to define both the water surface and the seabed.
We're calling the height to the datum $h$ and the water level relative to the datum $\eta$. Since we're just assuming uniform bathymetry, we'll set $h$ to a constant value:
In [1]:
h = 10
We now define $\eta$, the water-level relative to $h$. Since this is a variable which varies across the model domain, we must hold a value for it for every cell in the $x$ direction. Therefore we'll initialise a vector of length $N$ (the size of the domain) to hold a water-level value for each spatial cell as the model propagates through time:
In [134]:
eta = np.zeros(N)
We used the numpy .zeros() function which initialises a vector containing entirely zeros values. If we leave the eta variable like this we won't have much of an interesting model - we've effectively defined a flat water surface ($h + 0$ for all $x$). What we want is to start the model with the water suface in some sort of "disturbed", non-equlibrium, state. So let's alter the values in the eta vector.
In [161]:
for x in range(0,N):
eta[x] = 2.0 - (4.0/25.0)*x
We iterated through all the values in eta using Python's for loop statement. In the loop, we simply defined the water level as a linear function of $x$ with a y-intercept of $+2$ and a gradient of $4/N$. This basically gives use a slope of 4 across the model domain, meaning that the water level varies from $+2$ to $-2$.
Let's create a plot of this initial condition to check it.
In [162]:
fig = plt.figure()
wlevel = fig.add_subplot(111, xlim=(1, N-1), ylim=(0, np.max(h + eta)))
wlevel.plot(range(N), h + eta, lw=2, color='k')
wlevel.grid()
plt.xlabel('distance (x)')
plt.ylabel('Water height (H)')
plt.show()
That looks okay. Notice this is the entire water height (or depth), since we added eta to the constant datum height h.
The last bit of set up we need to do is to create some containers for each of the other variables which we need to solve for during each time step. The are the speed ($u$) and the total water height ($H$) which both appear in each of the depth-averaged equation we are using. So, as before we'll use the numpy .zeros() function to create placeholder vectors for these values. In contrast to the $\eta$ variable, for these variables we are going to include an extra value - i.e. one more than the number of spatial cells - so that we can handle the flow at the model boudary in an appropriate way. We cannot have water flowing through the boundary. So we create two $N+1$ sized vectors for $u$ and $H$:
In [99]:
u = np.zeros(N+1)
H = np.zeros(N+1)
In [ ]:
In [100]:
def time_step():
# Update local water heights based on last elevations
for i in range(1,N):
H[i] = h + eta[i]
# Handle boundary water heights
H[N] = H[N-1]
H[0] = h
# Calculate local current velocities
for i in range(1,N):
u[i] = u[i] + dt * (-g*((eta[i] - eta[i-1])/dx) - u[i]/(tau * H[i]))
# Calculate new local elevations (forward differences)
for i in range(0,N):
eta[i] = eta[i] - dt * ((H[i+1] * u[i+1] - H[i]*u[i])/dx)
if eta[i] < -h:
eta[i] = -h
In [101]:
fig = plt.figure()
# Limit x-axis to hide stationary points
ax = fig.add_subplot(111, xlim=(1, N), ylim=(8, 12))
ax.grid()
for _ in range(1):
time_step()
ax.plot(range(N+1), H, lw=2, color='k')
ax.text(N+1, H[-1], 't = 1', horizontalalignment='left')
for _ in range(9):
time_step()
ax.plot(range(N+1), H, lw=1, color='k')
ax.text(N+1, H[-1], 't = 10', horizontalalignment='left')
for _ in range(10):
time_step()
ax.plot(range(N+1), H, '-', lw=1, dashes=[1, 1], color='k')
ax.text(N+1, H[-1], 't = 20', horizontalalignment='left')
for _ in range(30):
time_step()
ax.plot(range(N+1), H, '-', lw=1, dashes=[2, 2], color='k')
ax.text(N+1, H[-1], 't = 50', horizontalalignment='left')
for _ in range(450):
time_step()
ax.plot(range(N+1), H, '-', lw=1, dashes=[3, 3], color='k')
ax.text(N+1, H[-1], 't = 500', horizontalalignment='left')
for _ in range(1000):
time_step()
ax.plot(range(N+1), H, '-', lw=1, dashes=[4, 4], color='k')
ax.text(N+1, H[-1], 't = 1500', horizontalalignment='left')
plt.show()
In [232]:
import matplotlib.animation as animation
In [ ]: