A Numerical Technique for Coupling the Momentum and the Continuity Equations for Semi-Implicit 3D Ocean Models

Numerical modeling of free-surface flows has been widely used to study the hydrodynamics and the transport processes of oceans, coastal waters, and rivers. Many of the geophysical flow models are based on Hydrostatic Primitive Equations (HPE). The HPEs relax the assumption of hydrostatic pressure, which reduces the vertical momentum equation to the hydrostatic pressure distribution. Large-scale phenomena in the free-surface flows are accurately described by the HPEs [14]. Several three-dimensional hydrostatic models have been developed and applied successfully to various practical problems. Different numerical methods, horizontal and vertical grid systems, and time discretization approaches have been employed to develop models.

In geophysical flow modelling, a reasonable balance between accuracy, stability, and computational efficiency is generally achieved when time discretization of the barotropic pressure gradient and the vertical eddy viscosity in momentum equations are implicit. This is due to the severe stability limitations that arise from these terms [21]. Different solution methods have been presented yet for the HPEs.

Some models (e.g., POM [2], ROMS [17], and FVCOM [9]) utilize mode-splitting technique [13], in which the pressure gradient terms are not often treated implicitly. However, in order to retain computational efficiency, the free-surface elevations are determined externally from the depth-integrated equations (external mode) and the velocity components are calculated from the main three-dimensional equations (internal mode). Therefore, a large time-step can be used for the internal mode because it is no longer limited by CFL stability conditions. Nevertheless, models using this technique have errors associated with using two systems of equations. Recent efforts in this field have led to the development of more schemes to solve the three-dimensional hydrostatic systems.

A free-surface correction method developed by Chen presents a predictor-corrector method for the HPEs [10]. In this method, first an intermediate free-surface is obtained with the explicit and implicit treatment of the pressure gradient and the vertical viscosity terms. Then the intermediate free-surface is corrected by solving a five-diagonal linear system for free-surface changes over the time step.

A semi-implicit finite difference method has been first presented by Casulli and Cheng which solves the three-dimensional HPEs directly with an implicit treatment of the pressure gradient and the vertical viscosity [4]. The semi-implicit algorithm is accurate, stable, and straightforward. It calculates the free-surface elevations and the horizontal velocity components simultaneously. This method has been also used to develop quasi-hydrostatic [6], non-hydrostatic [11], unstructured grid [8], and isopycnal coordinate [5] models. In recent years, semi-implicit methods have been used in several studies. The semi-implicit models such as TRIM-3D [4], UnTRIM [8], ECOM-si [3], SUSTAN [11], and ELCIRC [22] have enjoyed great success in simulating geophysical flows.

In the semi-implicit methods, the solution approach is to couple the horizontal momentum equations with the depth-integrated continuity equations in a way that a linear system is constructed for the free-surface. For this purpose, the horizontal velocity components should be derived in terms of neighbor water elevations. By substituting them in the depth-integrated continuity equation, a five-diagonal (or sparse in case of unstructured grid) linear system for free-surface is achieved. In the discretized momentum equation, each horizontal velocity component is related with both the neighbor water levels and the neighbor horizontal velocities (along the column). A basic way to solve such a system for the horizontal velocity is to arrange all the momentum equations along a water column in a tridiagonal linear system and inverting the coefficient matrix. Doing so is however computationally expensive and unjustified in practice.

In this paper, an algorithm is presented to construct a linear system for free-surface elevations. A three-dimensional shallow water model is developed to illustrate the details and validate the algorithm. The governing equations are discretized by finite difference method. In the proposed method, the horizontal velocity components are determined in terms of the free-surface elevations by some simple recursive formulas. The derived formulas are then substituted in the depth-integrated continuity equation to calculate the free-surface elevations.

The conservation laws of mass and momentum that describe the fluid motion are known as Reynolds-Averaged Navier-Stokes equations (RANS). The three-dimensional hydrostatic primitive equations (HPE) are derived from the Navier-Stokes equations by involving the Coriolis force and applying the hydrostatic approximation. Provided that the density is constant, a conservative form of the equations is expressed as:

$\displaystyle\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}+w\frac{\partial u}{\partial z}$

$\displaystyle=-g\frac{\partial\xi}{\partial x}+\nu_{h}\left(\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial y^{2}}\right)+\frac{\partial}{\partial z}\left(\nu_{z}\frac{\partial u}{\partial z}\right)+fv,$

(1)

$\displaystyle\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}+w\frac{\partial v}{\partial z}$

$\displaystyle=-g\frac{\partial\xi}{\partial y}+\nu_{h}\left(\frac{\partial^{2}v}{\partial x^{2}}+\frac{\partial^{2}v}{\partial y^{2}}\right)+\frac{\partial}{\partial z}\left(\nu_{z}\frac{\partial v}{\partial z}\right)-fu,$

(2)

$\displaystyle\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}$

$\displaystyle=0.$

(3)

where $u$, $v$, and $w$ are the velocity components in the $x$, $y$, and $z$ directions respectively, $t$ is the time, $\xi$ is the free-surface elevation measured from the still water level, $g$ is the gravity constant, $f$ is the Coriolis parameter, and $\nu_{h}$, $\nu_{z}$ are the horizontal and vertical viscosity coefficients. Because of the small effects of the horizontal viscosity terms in shallow water environments, it is justified to use a constant coefficient for the horizontal viscosity [21]. The vertical eddy viscosity coefficient is calculated using the Prandtl’s mixing-length theory [18].

The mathematical expression for the kinematic boundary condition at the impermeable bottom is:

$$u\frac{\partial h}{\partial x}+v\frac{\partial h}{\partial y}+w=0,$$

(4)

where $h$ is the water depth measured from the still water level. At the moving free surface, the kinematic boundary condition is:

$$\frac{\partial\xi}{\partial t}+u\frac{\partial\xi}{\partial x}+v\frac{\partial\xi}{\partial y}=w.$$

(5)

The depth-integrated continuity equation is derived by using the kinematic boundary conditions and integrating the continuity Eq. (3) over the water depth giving:

$$\frac{\partial\xi}{\partial t}+\frac{\partial}{\partial x}\int_{-h}^{\xi}u\,dz+\frac{\partial}{\partial y}\int_{-h}^{\xi}v\,dz=0.$$

(6)

The dynamic boundary condition at the free surface is determined by wind shear stress giving:

$$\nu\left.\frac{\partial u}{\partial z}\right|_{\xi}=\tau_{wx},\quad\nu\left.\frac{\partial v}{\partial z}\right|_{\xi}=\tau_{wy},$$

(7)

where $\tau_{wx}$ and $\tau_{wy}$ are the prescribed wind stress components [19]. Similarly, at the impermeable bottom boundary, the dynamic boundary condition is specified by the bed shear stress and is expressed as:

$$\nu\left.\frac{\partial u}{\partial z}\right|_{-h}=\tau_{bx},\quad\nu\left.\frac{\partial v}{\partial z}\right|_{-h}=\tau_{by},$$

(8)

where $\tau_{b}$ is the bed shear stress. The bed shear stress can be represented in the form of a quadratic friction law giving:

$$\tau_{bx}=\frac{g\sqrt{u^{2}+v^{2}}}{C^{2}}u,\quad\tau_{by}=\frac{g\sqrt{u^{2}+v^{2}}}{C^{2}}v,$$

(9)

where $C$ is Chezy’s roughness coefficient.

A finite difference method is used to discretize the governing equations and the boundary conditions on a staggered grid layout in a Cartesian coordinate system. The computational domain is discretized by $n_{i}$, $n_{j}$, and $n_{k}$ cells in the $x$, $y$, and $z$ directions, respectively. $\Delta x$, $\Delta y$, and $\Delta z$ are the cell sizes in $x$, $y$, and $z$ directions. In the $x$ and $y$ directions, a uniform cell size is used, while in the vertical direction $\Delta z$ varies with $k$. The height of top layer cells also varies with time and horizontal location to simulate the free-surface. Each cell is denoted by $(i,j,k)$ at center and the cell faces are denoted by $(i\pm 1/2,j,k)$, $(i,j\pm 1/2,k)$ and $(i,j,k\pm 1/2)$, as shown in Figure 1. The scalar quantities are defined at the cell centers, while the velocity components are located at the faces.

Discretization of the governing equations is based on the semi-implicit method [4] in which the barotropic pressure gradient and the vertical viscosity terms are discretized implicitly. However, the convective, the Coriolis, and the horizontal viscosity terms are finite differenced explicitly. A central difference scheme is used for spatial discretization. Notably, the nonlinear advection terms are solved based on the method described in [16] The general discretized form of the Eqs. (1) and (2) can be expressed as:

$\displaystyle\frac{u_{i-\frac{1}{2},j,k}^{n+1}-u_{i-\frac{1}{2},j,k}^{n}}{\Delta t}=Fu_{i-\frac{1}{2},j,k}^{n}-g\frac{\zeta_{i,j}^{n+1}-\zeta_{i-1,j}^{n+1}}{\Delta x}+\frac{V_{i-\frac{1}{2},j,k-\frac{1}{2}}\frac{u_{i-\frac{1}{2},j,k-1}^{n+1}-u_{i-\frac{1}{2},j,k}^{n+1}}{\Delta z_{i-\frac{1}{2},j,k-\frac{1}{2}}}-V_{i-\frac{1}{2},j,k+\frac{1}{2}}\frac{u_{i-\frac{1}{2},j,k}^{n+1}-u_{i-\frac{1}{2},j,k+1}^{n+1}}{\Delta z_{i-\frac{1}{2},j,k+\frac{1}{2}}}}{\Delta z_{i-\frac{1}{2},j,k}}$

(10)

$\displaystyle\frac{v_{i,j-\frac{1}{2},k}^{n+1}-v_{i,j-\frac{1}{2},k}^{n}}{\Delta t}=Fv_{i,j-\frac{1}{2},k}^{n}-g\frac{\zeta_{i,j}^{n+1}-\zeta_{i,j-1}^{n+1}}{\Delta y}+\frac{V_{i,j-\frac{1}{2},k-\frac{1}{2}}\frac{v_{i,j-\frac{1}{2},k-1}^{n+1}-v_{i,j-\frac{1}{2},k}^{n+1}}{\Delta z_{i,j-\frac{1}{2},k-\frac{1}{2}}}-V_{i,j-\frac{1}{2},k+\frac{1}{2}}\frac{v_{i,j-\frac{1}{2},k}^{n+1}-v_{i,j-\frac{1}{2},k+1}^{n+1}}{\Delta z_{i,j-\frac{1}{2},k+\frac{1}{2}}}}{\Delta z_{i,j-\frac{1}{2},k}},$

(11)

where $\Delta t$ is the time step and $F$ is an explicit operator including the explicit terms. There are several ways to discretize the explicit terms. We have used a second-order upwind explicit method for the convection terms and a central explicit method for the horizontal viscosity in the model.

By applying the dynamic boundary conditions, the discretized momentum equations for top layer cells take the form:

$\displaystyle\frac{u_{i-\frac{1}{2},j,1}^{n+1}-u_{i-\frac{1}{2},j,1}^{n}}{\Delta t}=Fu_{i-\frac{1}{2},j,1}^{n}-g\frac{\zeta_{i,j}^{n+1}-\zeta_{i-1,j}^{n+1}}{\Delta x}+\frac{\tau_{wx}-V_{i-\frac{1}{2},j+\frac{1}{2}}\frac{u_{i-\frac{1}{2},j,1}^{n+1}-u_{i-\frac{1}{2},j,2}^{n+1}}{\Delta z_{i-\frac{1}{2},j+\frac{1}{2}}}}{\Delta z_{i-\frac{1}{2},j,1}},$

(12)

$\displaystyle\frac{v_{i,j-\frac{1}{2},1}^{n+1}-v_{i,j-\frac{1}{2},1}^{n}}{\Delta t}=Fv_{i,j-\frac{1}{2},1}^{n}-g\frac{\zeta_{i,j}^{n+1}-\zeta_{i,j-1}^{n+1}}{\Delta y}+\frac{\tau_{wy}-V_{i,j-\frac{1}{2},1+\frac{1}{2}}\frac{v_{i,j-\frac{1}{2},1}^{n+1}-v_{i,j-\frac{1}{2},2}^{n+1}}{\Delta z_{i,j-\frac{1}{2},1+\frac{1}{2}}}}{\Delta z_{i,j-\frac{1}{2},1}},$

(13)

and for the bottom layer cells, the equations become:

$\displaystyle\frac{u_{i-\frac{1}{2},j,k}^{n+1}-u_{i-\frac{1}{2},j,k}^{n}}{\Delta t}=Fu_{i-\frac{1}{2},j,k}^{n}-g\frac{\zeta_{i,j}^{n+1}-\zeta_{i-1,j}^{n+1}}{\Delta x}+\frac{V_{i-\frac{1}{2},j,k-\frac{1}{2}}\frac{u_{i-\frac{1}{2},j,k-1}^{n+1}-u_{i-\frac{1}{2},j,k}^{n+1}}{\Delta z_{i-\frac{1}{2},j,k-\frac{1}{2}}}-\tau_{bx}}{\Delta z_{i-\frac{1}{2},j,k}},$

(14)

$\displaystyle\frac{v_{i,j-\frac{1}{2},nk}^{n+1}-v_{i,j-\frac{1}{2},nk}^{n}}{\Delta t}=Fv_{i,j-\frac{1}{2},nk}^{n}-g\frac{\zeta_{i,j}^{n+1}-\zeta_{i,j-1}^{n+1}}{\Delta y}+\frac{V_{i,j-\frac{1}{2},nk-\frac{1}{2}}\frac{v_{i,j-\frac{1}{2},nk-1}^{n+1}-v_{i,j-\frac{1}{2},nk}^{n+1}}{\Delta z_{i,j-\frac{1}{2},nk-\frac{1}{2}}}-\tau_{by}}{\Delta z_{i,j-\frac{1}{2},nk}}.$

(15)

where the velocities in the time step $n$ are used to calculate $\tau_{bx}$ and $\tau_{by}$.

A detailed description of the algorithm is presented below. For brevity, the $i$ and $j$ indices are not shown if they are unnecessary. The Eqs. (12)-(15) can be written in the following form:

$\displaystyle\mathbf{A}_{i-\frac{1}{2},j}^{n}\mathbf{U}_{i-\frac{1}{2},j}^{n+1}=\mathbf{D}_{i-\frac{1}{2},j}^{n+1},$

(16)

$\displaystyle\mathbf{A}_{i,j-\frac{1}{2}}^{n}\mathbf{V}_{i,j-\frac{1}{2}}^{n+1}=\mathbf{D}_{i,j-\frac{1}{2}}^{n+1},$

(17)

where

$\displaystyle\mathbf{A}^{n}=\begin{bmatrix}b_{1}&c_{1}&&&\\ a_{2}&b_{2}&c_{2}&&\\ &\ddots&\ddots&\ddots&\\ &&a_{nk-1}&b_{nk-1}&c_{nk-1}\\ &&&a_{nk}&b_{nk}\end{bmatrix},\quad\mathbf{U}^{n+1}=\begin{bmatrix}u_{1}^{n+1}\\ u_{2}^{n+1}\\ \vdots\\ u_{nk-1}^{n+1}\\ u_{nk}^{n+1}\end{bmatrix},\quad\mathbf{V}^{n+1}=\begin{bmatrix}v_{1}^{n+1}\\ v_{2}^{n+1}\\ \vdots\\ v_{nk-1}^{n+1}\\ v_{nk}^{n+1}\end{bmatrix},\quad\mathbf{D}^{n+1}=\begin{bmatrix}d_{1}\\ d_{2}\\ \vdots\\ d_{nk-1}\\ d_{nk}\end{bmatrix}.$

wherein the coefficients are:

		$\displaystyle a_{k}=-\frac{v_{z_{k-\frac{1}{2}}}\Delta t}{\Delta z_{k-\frac{1}{2}}\Delta z_{k}},$		(18)
		$\displaystyle c_{k}=\frac{v_{z_{k+\frac{1}{2}}}\Delta t}{\Delta z_{k+\frac{1}{2}}\Delta z_{k}},$
		$\displaystyle b_{k}=1-(a_{k}+c_{k}),$
		$\displaystyle d_{i-\frac{1}{2},j,k}=d^{\prime}_{i-\frac{1}{2},j,k}-g\Delta t\frac{\zeta_{i,j}^{n+1}-\zeta_{i-1,j}^{n+1}}{\Delta x},$
		$\displaystyle d_{i,j-\frac{1}{2},k}=d^{\prime}_{i,j-\frac{1}{2},k}-g\Delta t\frac{\zeta_{i,j}^{n+1}-\zeta_{i,j-1}^{n+1}}{\Delta y},$
		$\displaystyle d^{\prime}_{i-\frac{1}{2},j,k}=u_{i-\frac{1}{2},j,k}^{n}+\Delta tFu_{i-\frac{1}{2},j,k}^{n},$
		$\displaystyle d^{\prime}_{i,j-\frac{1}{2},k}=v_{i,j-\frac{1}{2},k}^{n}+\Delta tFv_{i,j-\frac{1}{2},k}^{n}.$

Now applying the dynamic boundary conditions yields to:

		$\displaystyle a_{1}=c_{nk}=0,$		(19)
		$\displaystyle d^{\prime}_{i-\frac{1}{2},j,1}=u_{i-\frac{1}{2},j,1}^{n}+F\Delta tu_{i-\frac{1}{2},j,1}^{n}+\frac{\Delta t\tau_{wx}}{\Delta z_{i-\frac{1}{2},j,1}},$
		$\displaystyle d^{\prime}_{i-\frac{1}{2},j,nk}=u_{i-\frac{1}{2},j,nk}^{n}+F\Delta tu_{i-\frac{1}{2},j,nk}^{n}-\frac{\Delta t\tau_{bx}}{\Delta z_{i-\frac{1}{2},j,nk}},$
		$\displaystyle d^{\prime}_{i,j-\frac{1}{2},1}=v_{i,j-\frac{1}{2},1}^{n}+F\Delta tv_{i,j-\frac{1}{2},1}^{n}+\frac{\Delta t\tau_{wy}}{\Delta z_{i,j-\frac{1}{2},1}},$
		$\displaystyle d^{\prime}_{i,j-\frac{1}{2},nk}=v_{i,j-\frac{1}{2},nk}^{n}+F\Delta tv_{i,j-\frac{1}{2},nk}^{n}+\frac{\Delta t\tau_{by}}{\Delta z_{i,j-\frac{1}{2},nk}}.$

The proposed algorithm for solving Eqs. (16-17) consists of three steps. First, a forward sweep eliminates the lower diagonal of the coefficients matrix. Second, the right-hand side entries are partitioned into the known and unknown parts and each part is formulated recursively. Third, a backward sweep solves the system for $U^{n+1}$ and $V^{n+1}$.

Two classical test cases have been undertaken to validate the proposed algorithm and to verify the mass and momentum conservation of the developed model.