09-sci-underflow.qmd

# Underflow Model {#sec-underflow}

## Overview

The underflow model in ELCOM is based on the numerical approach proposed by @Brad99 with finite-volume spatial discretization, conservation equations advanced by a predictor-corrector algorithm and @Roe81 approximate Riemann solver for fluxes through cell faces. The conservation equations of @Brad99 have been modified to use salinity and temperature instead of turbidity, and terms have been added to model 1) effects of the earth's rotation, 2) gradients in ambient fluid salinity and 3) momentum of the background circulation. The entrainment relationship has been changed to the formulation of @Dall01, developed from a parameterisation of the turbulent kinetic energy equation. The general grid transformations of @Brad99 are replaced with a simple Cartesian-grid approach for compatibility with the Cartesian grid in ELCOM.

## Governing Equations

The underflow model governing equations are the vertically-integrated equations of volume, momentum and salinity based on derivations of @Elli59. Assuming unity values for the shape factors while including the ambient fluid properties, the conservation equations of volume, x momentum, y momentum and salinity are respectively given by

$$
\frac{\partial(H)}{\partial t}+\frac{\partial(UH)}{\partial x}+\frac{\partial(VH)}{\partial y}=E\sqrt{U^2+V^2}
$$ {#eq-uf-volume}

$$
\begin{aligned}
\frac{\partial(UH)}{\partial t}+\frac{\partial(U^2H)}{\partial x}+\frac{\partial(UVH)}{\partial y} &= fVH-C_DU\sqrt{U^2+V^2}- \frac{1}{2}\frac{\partial(g'H^2)}{\partial x} \\
&\quad -g'HS_x+U_{amb}E\sqrt{U^2+V^2}
\end{aligned}
$$ {#eq-uf-xmom}

$$
\begin{aligned}
\frac{\partial(VH)}{\partial t}+\frac{\partial(UVH)}{\partial x}+\frac{\partial(V^2H)}{\partial y} &= -fUH-C_DV\sqrt{U^2+V^2}- \frac{1}{2}\frac{\partial(g'H^2)}{\partial y} \\
&\quad -g'HS_y+V_{amb}E\sqrt{U^2+V^2}
\end{aligned}
$$ {#eq-uf-ymom}

$$
\frac{\partial(SH)}{\partial t}+\frac{\partial(SUH)}{\partial x}+\frac{\partial(SVH)}{\partial y}=S_{amb}E\sqrt{U^2+V^2}
$$ {#eq-uf-salinity}

where $H$ is the height of the underflow; $U$ and $V$ are vertically averaged velocities in the $x$ and $y$ directions; $S$ is the vertically averaged salinity; $g'$ is the reduced gravity term $= g(\rho-\rho_{amb})/\rho_{amb}$, $f$ is the Coriolis parameter; $C_D$ is the bottom drag coefficient; $E$ is the entrainment coefficient and $S_x$ and $S_y$ are the bottom slopes in the $x$ and $y$ directions. The subscripts $amb$ refer to the ambient value of the variable. Ambient values used by the underflow model are obtained from the salinity and velocity at the bottom boundary cells of the 3D model, and are used to model entrainment through the last terms in the momentum and salinity equations.

The entrainment coefficient ($E$) in the governing equations is given by the entrainment law of @Dall01. The equation has been altered to account for the ambient flow and is given by

$$
E=\frac{C_KC_D^{3/2}+C_s\left(\frac{\Delta u}{\sqrt{U^2+V^2}}\right)^{3/2}}{Ri+ 10\left(C_KC_D^{3/2}+C_s\left(\frac{\Delta u}{\sqrt{U^2+V^2}}\right)^{3/2}\right)}
$$ {#eq-entrainment}

where

$$
\Delta u =\sqrt{(U^2-U_{amb}^2)+(V^2-V_{amb}^2)}
$$ {#eq-delta-u}

and $Ri$ is the bulk Richardson number given by

$$
Ri=\frac{g'H}{U^2+V^2}
$$ {#eq-Ri}

The parameterisation constants $C_K$ and $C_S$ account for the bottom production and shear production terms respectively.

## Model Implementation

The model is implemented using a finite-volume method. The predictor-corrector scheme (explicit prediction and implicit correction) as outlined by @Brad99 is used to obtain second-order temporal accuracy. The predictor solution is obtained from the primitive form of the governing equations. Following @Brad99, the primitive forms of governing equations are used for computational efficiency and can be written as:

$$
\frac{\partial H}{\partial t}+U\frac{\partial H}{\partial x}+H\frac{\partial U}{\partial x}+V\frac{\partial H}{\partial y} +H\frac{\partial V}{\partial y}=E\sqrt{U^2+V^2}
$$ {#eq-prim-vol}

$$
\begin{aligned}
\frac{\partial U}{\partial t}+g'\frac{\partial H}{\partial x}+U\frac{\partial U}{\partial x}+\frac{1}{2}H\frac{\partial g'}{\partial x}
+V\frac{\partial U}{\partial y} &= \frac{1}{H}\left(U_{amb}E\sqrt{U^2+V^2}-C_DU\sqrt{U^2+V^2} \right) \\
&\quad +fVH-gS_x
\end{aligned}
$$ {#eq-prim-xmom}

$$
\begin{aligned}
\frac{\partial V}{\partial t}+g'\frac{\partial H}{\partial y}+V\frac{\partial V}{\partial y}+\frac{1}{2}H\frac{\partial g'}{\partial y}
+U\frac{\partial V}{\partial x} &= \frac{1}{H}\left(V_{amb}E\sqrt{U^2+V^2}-C_DV\sqrt{U^2+V^2} \right) \\
&\quad -fUH-gS_y
\end{aligned}
$$ {#eq-prim-ymom}

$$
\frac{\partial S}{\partial t}+U\frac{\partial S}{\partial x}+V\frac{\partial S}{\partial y}=\frac{1}{H}(S-S_{amb})E\sqrt{U^2+V^2}
$$ {#eq-prim-sal}

An approximate (predictor) solution for the conservation of volume at time $t+\Delta t/2$ in cell $i$, is given by

$$
H_{i,j}^{n+1/2}=H_{i,j}^{n}-\frac{\Delta t}{2}
\left[
U\frac{\overline{\Delta H_x}}{\Delta x}
+V\frac{\overline{\Delta H_y}}{\Delta y}
+H\frac{\overline{\Delta U_x + \Delta U_y}}{\Delta x}
-E\sqrt{U^2+V^2}
\right]_{i,j}^n
$$ {#eq-predictor}

where superscripts represent the time step and subscripts represent the center of a grid cell in the horizontal underflow model grid. Note that $n+1/2$ represents the computational time of the predictor step. Similar expressions (not shown for brevity) can be written for predictor conservation equations of salinity and momentum. Overbars in @eq-predictor denote average gradients at the cell center computed from gradients across the cell faces in the $x$ and $y$ directions. For example

$$
\overline{\Delta H_x}=\text{avg}(H_{i,j}-H_{i-1,j},H_{i+1,j}-H_{i,j})=\text{avg}(a,b)
$$ {#eq-avg-gradient}

where the averaging function used in the underflow model is the $\beta$ family of averages given by

$$
\overline{\Delta H_x}=
\begin{cases}
\text{sign}(a)\min(\max(|a|,|b|),\beta \min(|a|,|b|)) & ab>0 \\
0 & ab\le0
\end{cases}
$$ {#eq-superbee}

For all simulations, $\beta$ is equal to two, which corresponds to a Superbee averaging function.

The corrector solution is a coupled, implicit solution of a discrete form of the conservative governing equations using the predictor-step fluxes. For example, the conservation of volume equation can be written as

$$
\begin{aligned}
\Delta x\Delta y \frac{H^{n+1}-H^{n}}{dt} &= F_F^{n+1/2}\Delta y+F_B^{n+1/2}\Delta y+F_R^{n+1/2}\Delta x+F_L^{n+1/2}\Delta x \\
&\quad +E\sqrt{U^{(n+1/2)^2}+V^{(n+1/2)^2}}
\end{aligned}
$$ {#eq-corrector}

where $F_F$, $F_R$, $F_B$ and $F_L$ are the volume fluxes through the front, right, back and left cell faces respectively. Here the nomenclature is consistent with the ELCOM code, where the front is the cell face in the positive $x$ direction, back is in negative $x$, right is positive $y$, and left is negative $y$. The fluxes are defined as positive in the direction of increasing computational coordinates.

The fluxes at each cell face are calculated using the approximate Riemann solver developed by @Roe81. These fluxes are given by

$$
F_{F(i,j)}=F_{B(i+1,j)}=\frac{1}{2}\left(f_{F(i,j)}^{n+1/2} + f_{B(i+1,j)}^{n+1/2} - \hat{R}|\hat{\Lambda}||\Delta \hat{V}|\right)
$$ {#eq-fluxTerm}

where the $f_{F(i,j)}$ and $f_{B(i+1,j)}$ refer to the estimated fluxes at the front face of cell $(i,j)$ and the back face of $(i+1,j)$ respectively. The third term on the right hand side is the Riemann flux term. The estimated fluxes at the front of cell $(i,j)$ are determined from the extrapolation of the primitive variables at the centre of the cell to the front face of the cell. For example the height at the front face of cell $(i,j)$ is given by

$$
H_{F(i,j)}^{n+1/2}=H_{(i,j)}^{n+1/2}+\frac{1}{2}\overline{\Delta H_x}_{i,j}^n
$$ {#eq-H-front}

and the estimated volume flux at the front face of cell $(i,j)$ is given by

$$
f_{F(i,j)}^{n+1/2} =H_{F(i,j)}^{n+1/2}U_{F(i,j)}^{n+1/2}
$$ {#eq-f-front}

Similar expressions can be written for the other faces and fluxes.

In matrix form the Riemann flux terms are detailed further in @Dall01 and @Brad99.

Two boundary conditions are used in the underflow model, solid wall and inflow. Inflows are implemented using ghost cells located outside the computational domain. The specified flow rate, salinity and underflow height are placed within the ghost cell for each timestep. Within the ghost cell all gradient terms are set to 0. The fluxes through the interface between the ghost cell and the interior cell are calculated using the Riemann solver as described above. Solid wall boundaries are cell faces where the fluxes are set to zero.

## Front Tracking

Following @Brad99 a front-tracking scheme is used to eliminate artificial propagation of the density front resulting from discretization. Cells in the underflow model are defined as being 'underflow', 'ambient' or 'frontal'. 'Underflow' cells are cells that are filled by the underflow; 'ambient' cells contain ambient water without underflow water, while 'frontal' cells are the boundary cells between the underflow and the ambient fluid that contain both types of water. In a frontal cell, fluxes are only calculated for faces that are bounded by an underflow cell. The front propagation distance in the $x$ direction, $X_f$, from time $n$ to $n+1$ within a frontal cell is given by

$$
X_{f}^{n+1} =X_{f}^{n}+U_f^{n+1}\Delta t
$$ {#eq-front}

where $U_f$ is the front velocity. An analogous formula is used for the frontal distance in the $y$ direction. If $X_f^{n+1}$ is greater than the grid dimension ($\Delta x$) then the cell is changed to an underflow cell and any adjacent ambient cell becomes a new frontal cell. While the concept is straightforward in one dimension, in 2D a cell may be a frontal cell with respect to both the $x$ and $y$ coordinate directions. This is simply a bookkeeping problem for the model, that is handled by designating any cell with both ambient and underflow water as a 'frontal' cell then checking whether one or two faces are involved in frontal propagation. At the end of each timestep, any ambient cell that has frontal cells on both $x$ or both $y$ faces is also designated as an underflow cell. This simplifies the computation for two fronts that are propagating in opposite directions and meet in a single ambient grid cell.

As pointed out by @Brad99 there are a number of possibilities for the frontal velocities $U_f$ and $V_f$. We follow their approach wherein the frontal velocity is based on the computed values of $U$ and $V$ in the frontal cells. We add the effect of the ambient flow on front propagation by computing the frontal velocities as

$$
U_{f} =U+U_{amb}
$$ {#eq-front-vel}

This formulation provides a simple and robust first-order approximation for calculating the frontal propagation. Due to the explicit nature of the front-tracking scheme, the front must not propagate further than one grid cell in any timestep. If this condition is violated, the underflow model timestep is automatically halved and the underflow is recomputed. This acts recursively so that the underflow solution will always be stable. Note the timestep for the ambient 3D flow solver is not affected by the change in timestep of the underflow model.

## Coupling the Model to ELCOM

This section discusses the coupling of the 2D underflow model with the 3D solver ELCOM. As an underflow travels along the bottom boundary of a water body it has two effects on the ambient water. Firstly it entrains ambient water into itself and secondly it displaces ambient fluid away from the bottom boundary. The proposed coupling scheme attempts to model these effects via the use of inflows and outflows across the bottom boundary of the 3D domain. The fundamental idea behind the scheme is to insert the net change in the underflow volume during each timestep through the bottom boundary of the 3D domain over the appropriate cells while still tracking the salinity of the underflow with the 2D model. In effect the water that is inserted into the 3D domain is designed to act as a 'place holder' that will be replaced with the underflow water once the underflow reaches its insertion level. A schematic of the coupling procedure is shown for a homogeneous ambient water body in @fig-coupling. The three basic steps of the scheme are:

1. Removal of the entrained water from the 3D domain (@fig-coupling a)
2. Insertion of the net volume increase in each underflow cell into the corresponding ELCOM bottom boundary cell (@fig-coupling b)
3. Insertion of salinity once the plume has reached its level of neutral buoyancy (@fig-coupling c and d)

Step one and two are implemented every timestep, step three only occurs once the underflow reaches the insertion level.

![Profiles showing coupling of the 2D underflow model to the 3D solver. The net volume for each underflow cell is calculated for each timestep (a). This volume is inserted as an inflow across the bottom boundary of the 3D domain (b). Once the underflow reaches its insertion level the salinity within a volume at the bottom of each column in the 3D domain is replaced with the salinity of the underflow. (d) shows the representation of the underflow within the 3D domain after insertion.](images/coupling.pdf){#fig-coupling}

Both the entrainment and the net change in volume of an underflow cell (@fig-coupling a) are easily calculated from the underflow model. After each underflow timestep is completed the entrained volume is removed from the appropriate cell on the bottom boundary of the 3D domain. The net change in volume of the underflow cell is then inserted into the 3D domain as an inflow across the bottom boundary of the corresponding 3D column. So as not to introduce any unwanted currents in the 3D domain this inflow is given a salinity equal to the salinity at the bottom of the ELCOM column. The idea is to insert the underflow volume over the appropriate cells while having the minimal effect on the ambient motions. In practice, because the entrained and the insertion fluid both have a salinity given by the salinity at the bottom of the 3D domain, step 1 and 2 are combined to form a single inflow for each column. The volume of the inflow is given by the change in volume of the cell minus the volume entrained into the cell and is equal to the net volume flux into the cell given by the modeled Riemann fluxes. If the net volume flux into the underflow cell is negative (e.g. at the tail of the underflow) then the appropriate volume is removed from the ELCOM cell as an outflow.

It can easily be shown that the sum of the net volume fluxes into all underflow cells, and hence the total volume flux into ELCOM, is equal to the specified inflow flux at the side boundary. This ensures that conservation of volume is maintained in the 3D solver and that the modeled free surface height in ELCOM accurately represents the supplied inflows.

Once the underflow reaches its insertion level the salinity in the 3D domain is altered to reflect the salinity of the underflow. The criteria for insertion are defined as when any underflow body cell is less dense than the ambient ELCOM cell 'above it' or when the height of the underflow is greater than 5 times ELCOM's vertical grid size.

The idea behind the salinity insertion subroutine is to replace the salinity within a certain volume at the bottom of each column ($i,j$) in the 3D domain such that after the insertion the salinity in the 3D domain represents the height and salinity predicted by the underflow model at each cell ($i,j$). If the 3D domain has a uniform salinity before the underflow event occurs then the insertion is trivial. For this case both the volume input at the bottom boundary of the ELCOM domain and the entrained water into the underflow have a salinity equal to the uniform ambient salinity. It can be shown that for this scenario conservation of salinity is achieved by replacing the ambient salinity at the bottom of each column ($i,j$) up to a height equal to the underflow height at cell ($i,j$) with the underflow salinity ($S_{i,j}$). The replacement of salinity is achieved using the following formulae:

$$
S^{ELCOM} =
\begin{cases}
S & z_{top}<H \\
\frac{S(H-z_{bot})}{dz} +\frac{S^{ELCOM}(z_{top}-H)}{dz} & z_{bot}<H \le z_{top} \\
S^{ELCOM} & z_{bot}>H
\end{cases}
$$ {#eq-insertSal}

where $S^{ELCOM}$ is the salinity in the ELCOM cell, $dz$ is the vertical grid size and $z_{top}$ and $z_{bot}$ are the height above the bottom of the top and bottom faces of the ELCOM cell. The representation of the underflow salinity in the 3D solver after insertion is shown in @fig-coupling(d).

When the ambient salinity is not uniform the above scheme is not conservative for salinity. The problem is the total amount of salinity being replaced is not necessarily equal to the salinity that was input into the 3D domain before the time of insertion. There are two possible solutions to this; either the volume of water being replaced or the underflow salinity can be altered to arbitrarily balance salinity. In ELCOM the former was chosen for two reasons. Firstly, if the underflow salinity is changed there is the possibility that the new salinity will be unrealistically high. Secondly, because the height of the underflow in ELCOM must be an integer number of cells high a small change in the volume of water replaced in the ELCOM domain will have a small effect on the representation of the underflow within ELCOM. The height over which salinity is replaced in a column is given by

$$
H_{insert}=H+\frac{S_{removed}-S_{added}}{S\Delta x\Delta y}
$$ {#eq-h-insert}

where $S_{removed}$ is the total salinity replaced in the ELCOM column, $S_{added}$ is the salinity that has been added through the bottom of the ELCOM column and $H$ and $S$ are the predicted height and salinity of the underflow. Once $H_{insert}$ has been calculated the salinity within the 3D domain can be modified according to @eq-insertSal above with $H$ replaced by $H_{insert}$. Note that if $S_{removed} = S_{added}$, as in the ambient uniform case, @eq-h-insert predicts an insertion height equal to the height of the underflow.

The two major parts of the insertion routine can be summarised as:

1. At each time step a volume is inserted through the bottom boundary of each column in the 3D domain. This volume is equal to the net volume flux into the corresponding underflow cell. The salinity of the inserted water is given by the salinity in the bottom cell of the ELCOM column.
2. Once the underflow has reached its level of insertion then the salinity in a certain volume at the base of every column in the 3D domain is replaced by the salinity of the corresponding underflow cell. The volume of water replaced is chosen to ensure conservation of salinity.

It should be noted that once the insertion of salinity has taken place and the tracking of the underflow has essentially been handed over to ELCOM that numerical convective entrainment will still occur. Two scenarios lead to the insertion of salinity. The first is when the height of the underflow is large compared to the vertical grid size in which case the effect of numerical mixing is small compared to the volume of the flow. The second scenario is when the underflow reaches neutral buoyancy. For this scenario the effect of convective entrainment is reduced because as water flows over a step in the bathymetry the water below it will have a similar salinity.

## Plunge Point Analysis

In order to model the plunging region of an inflow event a simple first-order approximation was made to the underflow model to account for the error in the momentum equation. This modification deals with the symptoms of the momentum error rather than addressing the error directly. Once the underflow volume fluxes are calculated the underflow height for each underflow cell at timestep $n+1$ is first estimated by neglecting the entrainment into the flow such that

$$
\Delta x\Delta y \frac{H^{n+1}-H^{n}}{dt}=F_F^{n+1/2}\Delta y+F_B^{n+1/2}\Delta y+F_R^{n+1/2}\Delta x+F_L^{n+1/2}\Delta x
$$ {#eq-plunge-height}

The predicted underflow height for each cell is then compared to the corresponding depth of the reservoir in the 3D domain. If the underflow height is greater than the water depth then the volume fluxes are manipulated so that the $H^{n+1}$ is equal to the depth of the reservoir. More specifically the volume flux through the face in the direction of the underflow front is increased to provide the necessary change. Scalar fluxes are altered to reflect the change in volume flux. In order to ensure that no mixing takes place before the plunge point of the underflow the entrainment coefficient for these cells is set to zero. Rather than altering the momentum equation to deal with a plunging flow, the conservation of volume equation is utilised to give

$$
U^{n+1}=\frac{F_F^{n+1/2}+F_B^{n+1/2}}{H^{n+1/2}}
$$ {#eq-plunge-U}

$$
V^{n+1}=\frac{F_R^{n+1/2}+F_L^{n+1/2}}{H^{n+1/2}}
$$ {#eq-plunge-V}

A further complication can arise during the solution of the conservation equations. The predicted height at time $t = t+\Delta t$ is given by the conservation of volume equation

$$
\begin{aligned}
\Delta x\Delta y \frac{H^{n+1}-H^{n}}{dt} &= F_F^{n+1/2}\Delta y+F_B^{n+1/2}\Delta y+F_R^{n+1/2}\Delta x+F_L^{n+1/2}\Delta x \\
&\quad +E\sqrt{U^{(n+1/2)^2}+V^{(n+1/2)^2}}
\end{aligned}
$$ {#eq-plunge-full}

If the entrainment term is large it is still possible for $H^{n+1}$ to be greater than the depth of the water column. This is easily dealt with by modifying the entrainment coefficient, $E$, where the underflow height is greater than the water depth to ensure $H^{n+1}$ is equal to the depth.