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*:

where $u$ is the depth-averaged current velocity (m^{3} 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*:

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.

`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

```
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

```
In [134]:
```eta = np.zeros(N)

`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()

```
```

`eta`

to the constant datum height `h`

.

`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 [ ]:
```

$$u_{x,t} = u_{x,t-1} + dt (-g\frac{\eta_{x,t-1} - \eta_{x-1,t-1}}{dx} - \frac{u_{x,t-1}}{\tau \times H_{x,t-1}})$$$$\eta_{x,t} = \eta_{x,t-1} - dt\frac{H_{x+1,t-1} \times u_{x+1,t-1} - H_{x,t-1} \times u_{x,t-1}}{dx}$$

```
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 [ ]:
```