## A NUMERICAL TECHNIQUE FOR COUPLING THE MOMENTUM AND THE CONTINUITY EQUATIONS FOR SEMI-IMPLICIT 3D OCEAN MODELS

## Ali Shahabi

Department of Civil and Environmental Engineering Old Dominion University, Norfolk, VA ashah010@odu.edu

## Reza Ghiassi

Faculty of Civil Engineering, University College of Engineering University of Tehran, Tehran, Iran rghiassi@ut.ac.ir

August 5, 2024

## ABSTRACT

Semi-implicit methods are powerful and efficient tools for the three-dimensional modeling of coastal and oceanic processes. A semi-implicit finite difference method for 3D hydrostatic primitive equations is presented in this paper. The governing equations are time-discretized with an implicit treatment of barotropic pressure gradient and vertical viscosity. The discretized momentum equations along a water column are coupled with the depth-integrated continuity equation of the column to construct a linear system for free-surface elevations. The novelty of this work lies in formulating an efficient method with the order of complexity O(N) for coupling the momentum and the continuity equations. In this method, the horizontal velocity components are expressed in terms of neighbor free-surface elevations by some simple recursive formulas and then are substituted in the integrated continuity equation.

K eywords free-surface flow · shallow water · finite difference · semi-implicit · hydrostatic pressure · three-dimensional

## 1 Introduction

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 threedimensional 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 depthintegrated 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.

## 2 Governing Equations

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:

$$\frac { \partial u } { \partial t } + u \frac { \partial u } { \partial x } + v \frac { \partial u } { \partial y } + w \frac { \partial u } { \partial z } \ = - 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 ) + f v,$$

$$\frac { \partial v } { \partial t } + u \frac { \partial v } { \partial x } + v \frac { \partial v } { \partial y } + w \frac { \partial v } { \partial z } \ = - g \frac { \partial \xi } { \partial y } + \nu _ { h } \left ( \frac { \partial ^ { 2 } v } { \partial x ^ { 2 } } + \frac { \partial ^ { \bar { 2 } } v } { \partial y ^ { 2 } } \right ) + \frac { \partial } { \partial z } \left ( \nu _ { z } \frac { \partial v } { \partial z } \right ) - f u,$$

$$\frac { \partial u } { \partial x } + \frac { \partial v } { \partial y } + \frac { \partial w } { \partial z } \ = 0.$$

where u v , , and w are the velocity components in the x y , , and z directions respectively, t is the time, ξ is the free-surface elevation measured from the still water level, g is the gravity constant, f is the Coriolis parameter, and ν h , ν 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,$$

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.$$

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 \, d z + \frac { \partial } { \partial y } \int _ { - h } ^ { \xi } v \, d z = 0.$$

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

$$\nu \, \frac { \partial u } { \partial z } \Big | _ { \xi } = \tau _ { w x }, \quad \nu \, \frac { \partial v } { \partial z } \Big | _ { \xi } = \tau _ { w y },$$

where τ wx and τ 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 \, \frac { \partial u } { \partial z } \Big | _ { - h } = \tau _ { b x }, \quad \nu \, \frac { \partial v } { \partial z } \Big | _ { - h } = \tau _ { b y },$$

where τ b is the bed shear stress. The bed shear stress can be represented in the form of a quadratic friction law giving:

$$\tau _ { b x } = \frac { g \sqrt { u ^ { 2 } + v ^ { 2 } } } { C ^ { 2 } } u, \quad \tau _ { b y } = \frac { g \sqrt { u ^ { 2 } + v ^ { 2 } } } { C ^ { 2 } } v,$$

where C is Chezy's roughness coefficient.

## 3 Description of Proposed Algorithm

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. ∆ x , ∆ y , and ∆ 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 ∆ 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 ± 1 / 2 , j, k ) , ( i, j ± 1 / 2 , k ) and ( i, j, k ± 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.

Figure 1: Computational mesh and notations

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00990v1-figure-fdfe9e3c-220a-4523-a225-6f7bf3e5c123.jpg)

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:

$$\frac { u _ { i - \frac { 1 } { 2 }, j, k } ^ { n + 1 } - u _ { i - \frac { 1 } { 2 }, j, k } ^ { n } } { \Delta t } = F u _ { 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 - \frac { 1 } { 2 } } - u _ { i - \frac { 1 } { 2 }, j, k } ^ { n + 1 } } { \Delta z _ { i - \frac { 1 } { 2 }, j, k + \frac { 1 } { 2 } } } } { \Delta z _ { i - \frac { 1 } { 2 }, j, k } } }$$

$$\Delta \tau & \quad \cdot \, \frac { \cdot \, 2 \nu \vdots ^ { \prime \prime } } { 2 \nu \vdots ^ { \prime } } \quad \Delta x & \quad \Delta z _ { i - \frac { 1 } { 2 }, j, k } } \\ \frac { v _ { i, j - \frac { 1 } { 2 }, k } ^ { n + 1 } - v _ { i, j - \frac { 1 } { 2 }, k } ^ { n } } { \Delta t } & = F v _ { i, j - \frac { 1 } { 2 }, k } ^ { n } - g \frac { \zeta _ { i, j } ^ { n + 1 } - \zeta _ { i, j - 1 } ^ { n + 1 } } { \Delta x } + \frac { V _ { i, j - \frac { 1 } { 2 }, k - 1 } - \frac { v _ { i, j - \frac { 1 } { 2 }, k - 1 } ^ { n + 1 } } { 2 \alpha _ { i, j - \frac { 1 } { 2 }, k - \frac { 1 } { 2 } } } - V _ { i, j - \frac { 1 } { 2 }, k + \frac { 1 } { 2 } } ^ { n + 1 } \frac { 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 } }, \\ \dots \dots \, \wedge \, \tau \, \colon \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \delta \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dclots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \ dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \end{cases}$$

where ∆ 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:

$$\frac { u _ { i - \frac { 1 } { 2 }, j, 1 } ^ { n + 1 } - u _ { i - \frac { 1 } { 2 }, j, 1 } ^ { n } } { \Delta t } = F u _ { i - \frac { 1 } { 2 }, j, 1 } ^ { n } - g \frac { \zeta _ { i, j } ^ { n + 1 } - \zeta _ { i - 1, j } ^ { n + 1 } } { \Delta x } + \frac { \tau _ { w x } - V _ { i - \frac { 1 } { 2 }, j, j + \frac { 1 } { 2 } } ^ { n + 1 } - u _ { i - \frac { 1 } { 2 }, j, 1 } ^ { n + 1 } } { \Delta z _ { i - \frac { 1 } { 2 }, j, 1 } ^ { n + 1 } } } { \Delta z _ { i - \frac { 1 } { 2 }, j, 1 } ^ { n + 1 } },$$

$$\frac { v _ { i, j - \frac { 1 } { 2 }, 1 } ^ { n + 1 } - v _ { i, j - \frac { 1 } { 2 }, 1 } ^ { n } } { \Delta t } = F v _ { i, j - \frac { 1 } { 2 }, 1 } ^ { n } - g \frac { \zeta _ { i, j } ^ { n + 1 } - \zeta _ { i, j - 1 } ^ { n + 1 } } { \Delta y } + \frac { \tau _ { w y } - V _ { i, j - \frac { 1 } { 2 }, 1 + \frac { 3 } { 2 } } ^ { n + 1 } - v _ { i, j - \frac { 1 } { 2 }, 2 } } { \Delta z _ { i, j - \frac { 1 } { 2 }, 1 } },$$

and for the bottom layer cells, the equations become:

$$\frac { u _ { i - \frac { 1 } { 2 }, j, k } ^ { n + 1 } - u _ { i - \frac { 1 } { 2 }, j, k } ^ { n } } { \Delta t } = F u _ { 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 - \frac { 1 } { 2 } } - u _ { i - \frac { 1 } { 2 }, j, k } ^ { n + 1 } } { \Delta z _ { i - \frac { 1 } { 2 }, j, k } } - \tau _ { b x } } } { \Delta z _ { i - \frac { 1 } { 2 }, j, k } },$$

2

$$\frac { v _ { i, j - \frac { 1 } { 2 }, n k } ^ { n + 1 } - v _ { i, j - \frac { 1 } { 2 }, n k } ^ { n } } { \Delta t } = V _ { i, j - \frac { 1 } { 2 }, n k } ^ { n } - g \frac { \zeta _ { i, j } ^ { n + 1 } - \zeta _ { i, j - 1 } ^ { n + 1 } } { \Delta y } + \frac { V _ { i, j - \frac { 1 } { 2 }, n k - \frac { 1 } { 2 } } ^ { n + 1 } - v _ { i, j - \frac { 1 } { 2 }, n k - \frac { 1 } { 2 } } ^ { n + 1 } - \tau _ { b y } ^ { n + 1 } } { \Delta z _ { i, j - \frac { 1 } { 2 }, n k } }.$$

where the velocities in the time step n are used to calculate τ bx and τ 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:

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

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

where

$$A ^ { n } = \begin{bmatrix} b _ { 1 } & c _ { 1 } \\ a _ { 2 } & b _ { 2 } & c _ { 2 } \\ & \ddots & \ddots & \ddots \\ & & a _ { n k - 1 } & b _ { n k - 1 } & c _ { n k - 1 } \\ & & a _ { n k } & b _ { n k } & b _ { n k } \end{bmatrix}, \ \ U ^ { n + 1 } = \begin{bmatrix} u _ { 1 } ^ { n + 1 } \\ u _ { 2 } ^ { n + 1 } \\ \vdots \\ u _ { n k - 1 } ^ { n + 1 } \\ u _ { n k } ^ { n + 1 } \end{bmatrix}, \ \ V ^ { n + 1 } = \begin{bmatrix} v _ { 1 } ^ { n + 1 } \\ v _ { 2 } ^ { n + 1 } \\ \vdots \\ v _ { n k - 1 } ^ { n + 1 } \\ v _ { n k } ^ { n + 1 } \end{bmatrix}, \ \ D ^ { n + 1 } = \begin{bmatrix} d _ { 1 } \\ d _ { 2 } \\ \vdots \\ d _ { n k - 1 } \\ d _ { n k } \end{bmatrix}.$$

wherein the coefficients are:

$$s \, & \text{are} \colon & \quad a _ { k } = - \frac { v _ { z _ { k } - \frac { 1 } { 2 } } } { \Delta z _ { k - \frac { 1 } { 2 } } } \Delta z _ { k } }, \\ & \quad c _ { k } = \frac { v _ { z _ { k + \frac { 1 } { 2 } } } } { \Delta z _ { k + \frac { 1 } { 2 } } } \Delta z _ { k } }, \\ & \quad b _ { k } = 1 - ( a _ { k } + c _ { k } ), \\ & \quad d _ { i - \frac { 1 } { 2 }, j, k } = d _ { i - \frac { 1 } { 2 }, j, k } ^ { \prime } - g \Delta t \frac { \zeta ^ { n + 1 } - \zeta ^ { n + 1 } _ { i - 1, j } } { \Delta x }, \\ & \quad d _ { i, j - \frac { 1 } { 2 }, k } = d _ { i, j - \frac { 1 } { 2 }, k } ^ { \prime } - g \Delta t \frac { \zeta ^ { n + 1 } _ { i, j } - \zeta ^ { n + 1 } _ { i, j - 1 } } { \Delta y }, \\ & \quad d _ { i - \frac { 1 } { 2 }, j, k } = u _ { i - \frac { 1 } { 2 }, j, k } ^ { n } + \Delta t F u _ { i - \frac { 1 } { 2 }, j, k } ^ { n }, \\ & \quad d _ { i, j - \frac { 1 } { 2 }, k } = v _ { i, j - \frac { 1 } { 2 }, k } ^ { n } + \Delta t F v _ { i, j - \frac { 1 } { 2 }, k } ^ { n }. \\ \text{mic boundary conditions yields to} \colon$$

Now applying the dynamic boundary conditions yields to:

$$\text{ne dynamic boundary conditions yields to:} \\ a _ { 1 } = c _ { n k } = 0, \\ d ^ { \prime } _ { i - \frac { 1 } { 2 }, j, 1 } = u ^ { n } _ { i - \frac { 1 } { 2 }, j, 1 } + F \Delta t u ^ { n } _ { i - \frac { 1 } { 2 }, j, 1 } + \frac { \Delta t \tau _ { w x } } { \Delta z _ { i - \frac { 1 } { 2 }, j, 1 } }, \\ d ^ { \prime } _ { i - \frac { 1 } { 2 }, j, n k } = u ^ { n } _ { i - \frac { 1 } { 2 }, j, n k } + F \Delta t u ^ { n } _ { i - \frac { 1 } { 2 }, j, n k } - \frac { \Delta t \tau _ { b x } } { \Delta z _ { i - \frac { 1 } { 2 }, j, n k } }, \\ d ^ { \prime } _ { i, j - \frac { 1 } { 2 }, 1 } = v ^ { n } _ { i, j - \frac { 1 } { 2 }, 1 } + F \Delta t v ^ { n } _ { i, j - \frac { 1 } { 2 }, 1 } + \frac { \Delta t \tau _ { w y } } { \Delta z _ { i, j - \frac { 1 } { 2 }, 1 } }, \\ d ^ { \prime } _ { i, j - \frac { 1 } { 2 }, n k } = v ^ { n } _ { i, j - \frac { 1 } { 2 }, n k } + F \Delta t v ^ { n } _ { i, j - \frac { 1 } { 2 }, n k } + \frac { \Delta t \tau _ { b y } } { \Delta z _ { i, j - \frac { 1 } { 2 }, n k } }. \\ \text{lgorithm for solving Equs. (16-17) consists of three steps. First, a forward sweep eliminates the lower
coefficients matrix: Second, the rioht.hand side entries are partitioned into the known and unknown.$$

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 .

## 3.1 Step 1: Forward Sweep

In the first step, the lower diagonal is eliminated by a forward sweep. This step is similar to the well-known double-sweep algorithm [1]. By applying the forward sweep, the Eqs. (16) and (17) become as follows:

$$& \mathbf E _ { i - \frac { 1 } { 2 }, j } ^ { n } \mathbf U _ { i - \frac { 1 } { 2 }, j } ^ { n + 1 } = \mathbf F _ { i - \frac { 1 } { 2 }, j } ^ { n + 1 }, \\ & \mathbf E _ { i, j - \frac { 1 } { 2 } } ^ { n } V _ { i, j - \frac { 1 } { 2 } } ^ { n + 1 } = \mathbf F _ { i, j - \frac { 1 } { 2 } } ^ { n + 1 },$$

where wherein the entries are:

n

E

=



1













-

e

.

.

n

1

.

-

e

n

k

.

.

.





f



.

.

.

1

























n

+1

1

.

.

.

n

+1

k

.

.















1

f

.

n

+1

nk

$$- c _ { 1 }$$

$$e _ { k } ^ { n } & = \frac { - c _ { k } } { a _ { k } e _ { k - 1 } ^ { n } + b _ { k } }, \quad e _ { 0 } ^ { n } = 0, & ( 2 2 ) \\ f _ { k } ^ { n + 1 } & = \frac { d _ { k } ^ { n + 1 } - a _ { k } f _ { k - 1 } ^ { n + 1 } } { a _ { k } e _ { k - 1 } ^ { n } + b _ { k } }, \quad f _ { 0 } ^ { n + 1 } = 0. & ( 2 3 )$$



,

F

n

+1

=

f

.

## 3.2 Step 2: Decomposing F n +1

In the second step, the right-hand side matrix entries, i.e. Eq. (23), are partitioned into two parts in a way that each part can be defined recursively.

Lemma 1. The Eq. (23) can be defined as follows for the x and y directions respectively:

$$f _ { i - \frac { 1 } { 2 }, j, k } ^ { n + 1 } = \psi _ { i - \frac { 1 } { 2 }, j, k } ^ { ( 1 ) } + \psi _ { i - \frac { 1 } { 2 }, j, k } ^ { ( 2 ) } \left ( \zeta _ { i, j } ^ { n + 1 } - \zeta _ { i - 1, j } ^ { n + 1 } \right ),$$

$$f _ { i, j - \frac { 1 } { 2 }, k } ^ { n + 1 } = \psi _ { i, j - \frac { 1 } { 2 }, k } ^ { ( 1 ) } + \psi _ { i, j - \frac { 1 } { 2 }, k } ^ { ( 2 ) } \left ( \zeta _ { i, j } ^ { n + 1 } - \zeta _ { i, j - 1 } ^ { n + 1 } \right ),$$

where the recursive coefficient ψ (1) for both x and y directions is:

$$\psi _ { k } ^ { ( 1 ) } = \frac { d _ { k } ^ { \prime } - a _ { k } \psi _ { k - 1 } ^ { ( 1 ) } } { a _ { k } e _ { k - 1 } ^ { n } + b _ { k } }, \quad \psi _ { 0 } ^ { ( 1 ) } = 0,$$

and ψ (2) for x and y directions are:

$$\psi _ { i - \frac { 1 } { 2 }, j, k } ^ { ( 2 ) } = \frac { - g \frac { \Delta t } { \Delta x } - a _ { k } \psi _ { i - \frac { 1 } { 2 }, j, k - 1 } ^ { ( 2 ) } } { a _ { k } e _ { k - 1 } ^ { n } + b _ { k } }, \quad \psi _ { i - \frac { 1 } { 2 }, j, 0 } ^ { ( 2 ) } = 0,$$

$$\psi _ { i, j - \frac { 1 } { 2 }, k } ^ { ( 2 ) } = \frac { - g \frac { \Delta t } { \Delta y } - a _ { k } \psi _ { i, j - \frac { 1 } { 2 }, k - 1 } ^ { ( 2 ) } } { a _ { k } e _ { k - 1 } ^ { n } + b _ { k } }, \quad \psi _ { i, j - \frac { 1 } { 2 }, 0 } ^ { ( 2 ) } = 0.$$

Proof. Using mathematical induction, the Eq. (24) is proved and likewise the Eq. (25) can be proved. In the basis step, it is proved that the Eq. (24) holds for k = 1 . From the Eq. (23), f 1 is:

$$& f _ { 1 } ^ { n + 1 } = \frac { d _ { 1 } - a _ { 1 } f _ { 0 } ^ { n + 1 } } { a _ { 1 } e _ { 0 } + b _ { 1 } } = \frac { \left ( d _ { 1 } ^ { \prime } - g \frac { \Delta t } { \Delta x } ( \zeta _ { i, j } ^ { n + 1 } - \zeta _ { i - 1, j } ^ { n + 1 } ) \right ) - 0 } { 0 + b _ { 1 } } = \frac { d _ { 1 } ^ { \prime } } { b _ { 1 } } - \frac { g \Delta t } { b _ { 1 } \Delta x } \left ( \zeta _ { i, j } ^ { n + 1 } - \zeta _ { i - 1, j } ^ { n + 1 } \right ), \\ & \Rightarrow f _ { 1 } ^ { n + 1 } = \psi _ { 1 } ^ { ( 1 ) } + \psi _ { 1 } ^ { ( 2 ) } ( \zeta _ { i, j } ^ { n + 1 } - \zeta _ { i - 1, j } ^ { n + 1 } ).$$

In the inductive step, it is proved if the Eq. (24) holds for k , then it holds for k +1 . We have:

$$f _ { k + 1 } ^ { n + 1 } & = \frac { d _ { k + 1 } - a _ { k + 1 } f _ { k } ^ { n + 1 } } { a _ { k + 1 } e _ { k } ^ { n } + b _ { k + 1 } } \quad ( * ) \\ f _ { k } ^ { n + 1 } & = \psi _ { k } ^ { ( 1 ) } + \psi _ { k } ^ { ( 2 ) } ( \zeta _ { i, j } ^ { n + 1 } - \zeta _ { i - 1, j } ^ { n + 1 } ) \quad ( * * ).$$

Substituting (**) into (*) gives:

$$& f _ { k + 1 } ^ { n + 1 } = \frac { \left ( d _ { k + 1 } ^ { \prime } - g \frac { \Delta t } { \Delta x } ( \zeta _ { i, j } ^ { n + 1 } - \zeta _ { i - 1, j } ^ { n + 1 } ) \right ) - a _ { k + 1 } \left ( \psi _ { k } ^ { ( 1 ) } + \psi _ { k } ^ { ( 2 ) } ( \zeta _ { i, j } ^ { n + 1 } - \zeta _ { i - 1, j } ^ { n + 1 } ) \right ) } { a _ { k + 1 } e _ { k } ^ { n } + b _ { k + 1 } } \\ & = \left ( \frac { d _ { k + 1 } ^ { \prime } - a _ { k + 1 } \psi _ { k } ^ { ( 1 ) } } { a _ { k + 1 } e _ { k } ^ { n } + b _ { k + 1 } } \right ) + \left ( \frac { - g \frac { \Delta t } { \Delta x } - a _ { k + 1 } \psi _ { k } ^ { ( 2 ) } } { a _ { k + 1 } e _ { k } ^ { n } + b _ { k + 1 } } \right ) ( \zeta _ { i, j } ^ { n + 1 } - \zeta _ { i - 1, j } ^ { n + 1 } ) \\ & \Rightarrow f _ { k + 1 } ^ { n + 1 } = \psi _ { k + 1 } ^ { ( 1 ) } + \psi _ { k + 1 } ^ { ( 2 ) } ( \zeta _ { i, j } ^ { n + 1 } - \zeta _ { i - 1, j } ^ { n + 1 } ).$$

## 3.3 Step 3: Backward Sweep

In the third step, a backward sweep produces the horizontal velocities in terms of the free-surface elevations.

Lemma 2. The horizontal velocity components can be expressed as follows:

$$u _ { i - \frac { 1 } { 2 }, j, k } ^ { n + 1 } = \Omega _ { i - \frac { 1 } { 2 }, j, k } ^ { ( 1 ) } + \Omega _ { i - \frac { 1 } { 2 }, j, k } ^ { ( 2 ) } ( \zeta _ { i, j } ^ { n + 1 } - \zeta _ { i - 1, j } ^ { n + 1 } ),$$

$$v _ { i, j - \frac { 1 } { 2 }, k } ^ { n + 1 } & = \Omega _ { i, j - \frac { 1 } { 2 }, k } ^ { ( 1 ) } + \Omega _ { i, j - \frac { 1 } { 2 }, k } ^ { ( 2 ) } ( \zeta _ { i, j } ^ { n + 1 } - \zeta _ { i, j - 1 } ^ { n + 1 } ),$$

where the coefficients Ω (1) and Ω (2) are defined as:

$$\Omega _ { k } ^ { ( 1 ) } = \psi _ { k } ^ { ( 1 ) } + e _ { k } ^ { n } \Omega _ { k + 1 } ^ { ( 1 ) }, \quad \Omega _ { n k + 1 } ^ { ( 1 ) } = 0,$$

$$\Omega _ { k } ^ { ( 2 ) } = \psi _ { k } ^ { ( 2 ) } + e _ { k } ^ { n } \Omega _ { k + 1 } ^ { ( 2 ) }, \quad \Omega _ { n k + 1 } ^ { ( 2 ) } = 0.$$

Proof. A backward induction is used to prove Eqs. (27) and (28). In the basis step, it is proved that the Eq. (27) holds for k = n k . The Eq. (20) for k = n k is:

$$& u _ { n k } ^ { n + 1 } = f _ { n k } ^ { n + 1 } \\ & = \psi _ { n k } ^ { ( 1 ) } + \psi _ { n k } ^ { ( 2 ) } ( \zeta _ { i, j } ^ { n + 1 } - \zeta _ { i - 1, j } ^ { n + 1 } ) \\ & \Rightarrow u _ { n k } ^ { n + 1 } = \Omega _ { n k } ^ { ( 1 ) } + \Omega _ { n k } ^ { ( 2 ) } ( \zeta _ { i, j } ^ { n + 1 } - \zeta _ { i - 1, j } ^ { n + 1 } ).$$

In the inductive step, it is shown if the Eq. (27) holds for k +1 , it holds for k . the k -th row of the Eq. (20) is:

$$u _ { k } ^ { n + 1 } - e _ { k } u _ { k } ^ { n + 1 } = f _ { k } ^ { n + 1 } \ \ ( * ).$$

and Eq. (27) holds for k +1 is:

$$u _ { k + 1 } ^ { n + 1 } = \Omega _ { k + 1 } ^ { ( 1 ) } + \Omega _ { k + 1 } ^ { ( 2 ) } ( \zeta _ { i, j } ^ { n + 1 } - \zeta _ { i - 1, j } ^ { n + 1 } ) \ \ ( * * ).$$

Substituting (**) into (*) yields:

$$& u _ { k } ^ { n + 1 } - e _ { k } \left ( \Omega _ { k + 1 } ^ { ( 1 ) } + \Omega _ { k + 1 } ^ { ( 2 ) } ( \zeta _ { i, j } ^ { n + 1 } - \zeta _ { i - 1, j } ^ { n + 1 } ) \right ) = \psi _ { k } ^ { ( 1 ) } + \psi _ { k } ^ { ( 2 ) } ( \zeta _ { i, j } ^ { n + 1 } - \zeta _ { i - 1, j } ^ { n + 1 } ) \\ & \Rightarrow u _ { k } ^ { n + 1 } = \left ( \psi _ { k } ^ { ( 1 ) } + e _ { k } ^ { n } \Omega _ { k + 1 } ^ { ( 1 ) } \right ) + \left ( \psi _ { k } ^ { ( 2 ) } + e _ { k } ^ { n } \Omega _ { k + 1 } ^ { ( 2 ) } \right ) ( \zeta _ { i, j } ^ { n + 1 } - \zeta _ { i - 1, j } ^ { n + 1 } ) \\ & \Rightarrow u _ { k } ^ { n + 1 } = \Omega _ { k } ^ { ( 1 ) } + \Omega _ { k } ^ { ( 2 ) } ( \zeta _ { i, j } ^ { n + 1 } - \zeta _ { i - 1, j } ^ { n + 1 } ).$$

The Eq. (28) can be also proved similarly. Finally, the horizontal velocity components are determined in terms of free-surface elevations as follows (matrix form):

$$\mathsf U _ { i - \frac { 1 } { 2 }, j } ^ { n + 1 } & = \mathsf \Omega _ { i - \frac { 1 } { 2 }, j } ^ { ( 1 ) } + \mathsf \Omega _ { i - \frac { 1 } { 2 }, j } ^ { ( 2 ) } ( \zeta _ { i, j } ^ { n + 1 } - \zeta _ { i - 1, j } ^ { n + 1 } ).$$

$$V _ { i, j - \frac { 1 } { 2 } } ^ { n + 1 } & = \Omega _ { i, j - \frac { 1 } { 2 } } ^ { ( 1 ) } + \Omega _ { i, j - \frac { 1 } { 2 } } ^ { ( 2 ) } ( \zeta _ { i, j } ^ { n + 1 } - \zeta _ { i, j - 1 } ^ { n + 1 } ).$$

In order to determine ξ n +1 , the Eqs. (30) and (31) are substituted into the depth-integrated continuity equation. The discretized form of the Eq. (6) in a matrix notation takes the form:

$$\frac { \zeta _ { i, j } ^ { n + 1 } - \zeta _ { i, j } ^ { n } } { \Delta t } + \frac { ( \Delta z ) _ { i + \frac { 1 } { 2 }, j } ^ { \top } U _ { i + \frac { 1 } { 2 }, j } ^ { n + 1 } - ( \Delta z ) _ { i - \frac { 1 } { 2 }, j } ^ { \top } U _ { i - \frac { 1 } { 2 }, j } ^ { n + 1 } } { \Delta x } + \frac { ( \Delta z ) _ { i, j + \frac { 1 } { 2 } } ^ { \top } U _ { i, j + \frac { 1 } { 2 } } ^ { n + 1 } - ( \Delta z ) _ { i, j - \frac { 1 } { 2 } } ^ { \top } U _ { i, j - \frac { 1 } { 2 } } ^ { n + 1 } } { \Delta y } = 0.$$

where A is a five-diagonal matrix and B is the right-hand side vector.

Substitution of the Eqs. (30) and (31) into the Eq. (32) yields:

$$\zeta _ { i, j } ^ { n + 1 } + \frac { \Delta t } { \Delta x } \left \{ \left [ ( \Delta z ) _ { i + \frac { 1 } { 2 }, j } ^ { \top } } ^ { \top } \Omega _ { i + \frac { 1 } { 2 }, j } ^ { ( 2 ) } \right ] ( \zeta _ { i + 1, j } ^ { n + 1 } - \zeta _ { i, j } ^ { n + 1 } ) - \left [ ( \Delta z ) _ { i - \frac { 1 } { 2 }, j } ^ { \top } } ^ { \top } \Omega _ { i - \frac { 1 } { 2 }, j } ^ { ( 2 ) } ) \right ] ( \zeta _ { i, j } ^ { n + 1 } - \zeta _ { i - 1, j } ^ { n + 1 } ) \right \} \\ + \frac { \Delta t } { \Delta y } \left \{ \left [ ( \Delta z ) _ { i, j + \frac { 1 } { 2 } } ^ { \top } } ^ { \top } \Omega _ { i, j + \frac { 1 } { 2 } } ^ { ( 2 ) } \right ] ( \zeta _ { i, j + 1 } ^ { n + 1 } - \zeta _ { i, j } ^ { n + 1 } ) - \left [ ( \Delta z ) _ { i, j - \frac { 1 } { 2 } } ^ { \top } } ^ { \top } \Omega _ { i, j - \frac { 1 } { 2 } } ^ { ( 2 ) } \right ] ( \zeta _ { i, j } ^ { n + 1 } - \zeta _ { i, j - 1 } ^ { n + 1 } ) \right \} \\ = \zeta _ { i, j } ^ { n } - \frac { \Delta t } { \Delta x } \left \{ \left [ ( \Delta z ) _ { i + \frac { 1 } { 2 }, j } ^ { \top } } ^ { \top } \Omega _ { i + \frac { 1 } { 2 }, j } ^ { ( 1 ) } \right ] - \left [ ( \Delta z ) _ { i - \frac { 1 } { 2 }, j } ^ { \top } } ^ { \top } \Omega _ { i - \frac { 1 } { 2 }, j } ^ { ( 1 ) } \right ] \right \} \\ - \frac { \Delta t } { \Delta y } \left \{ \left [ ( \Delta z ) _ { i, j + \frac { 1 } { 2 } } ^ { \top } } ^ { \top } \Omega _ { i, j + \frac { 1 } { 2 } } ^ { ( 1 ) } \right ] - \left [ ( \Delta z ) _ { i, j - \frac { 1 } { 2 } } ^ { \top } } ^ { \top } \Omega _ { i, j - \frac { 1 } { 2 } } ^ { ( 1 ) } \right ] \right \}, \\ \text{which is a symmetric positive definite and diagonally dominant five-diagonal system.  This system can be solved}$$

which is a symmetric positive definite and diagonally dominant five-diagonal system. This system can be solved efficiently by a conjugate gradient method [4].

In summary, the proposed algorithm to calculate ξ n +1 takes the following simple steps:

- 1. Calculate a , b k k , c k and d ′ k from the Eqs. (21) and (22).
- 2. Calculate e k from the Eq. (25).
- 3. Calculate ψ (1) k and ψ (2) k from the Eqs. (29).
- 4. Calculate Ω (1) k and Ω (2) k from the Eqs. (32).
- 5. Solve the five-diagonal linear system (36).

## 3.4 Calculating Velocity Field

Once the new free-surface elevations have been determined, the new horizontal velocity components are directly determined by substituting the new free-surface elevations in the Eqs. (27) and (28). Eventually, the Eq. (3) is directly applicable to calculate the vertical velocity components. This equation is discretized as follows:

$$\frac { u _ { i + \frac { 1 } { 2 }, j, k } ^ { n + 1 } - u _ { i - \frac { 1 } { 2 }, j, k } ^ { n + 1 } } { \Delta x } + \frac { v _ { i, j + \frac { 1 } { 2 }, k } ^ { n + 1 } - v _ { i, j - \frac { 1 } { 2 }, k } ^ { n + 1 } } { \Delta y } + \frac { w _ { i, j, k + \frac { 1 } { 2 } } ^ { n + 1 } - w _ { i, j, k - \frac { 1 } { 2 } } ^ { n + 1 } } { \Delta z _ { i, k } } = 0.$$

Using the boundary condition (4) that implies w n +1 n k = 0 and starting from the bottom layer cells, the vertical velocity components are determined.

## 4 Method Validation and Efficiency

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

## 4.1 Standing Wave in a Rectangular Basin

In the first test case, a standing wave in a rectangular basin was simulated. The analytical solution for a small amplitude standing wave is given by the following relations [15]:

$$\xi = A \cos ( k _ { x } x ) \cos ( k _ { y } y ) \cos ( \omega t ),$$

$$( u, v ) = A \frac { g } { \sigma } \frac { \cosh k \left ( h + ( - z ) \right ) } { \cosh k h } \sin ( \sigma t ) \left ( ( k _ { x } \sin ( k _ { x } x ) \cos ( k _ { y } y ), k _ { y } \sin ( k _ { x } x ) \cos ( k _ { y } y ) ) \right ),$$

$$w = A \frac { g k } { \sigma } \frac { \sinh k \left ( h + ( - z ) \right ) } { \cosh k h } \cos ( k _ { x } x ) \cos ( k _ { y } y ) \sin ( \sigma t ),$$

where A is the wave amplitude, ( k , k x y ) = ( π L x , π L y ) which in this test L x and L y are taken as the basin dimensions, k = √ k 2 x + k 2 y , h is the mean water depth, z is the vertical position, and σ = k √ gh . Figure 2 shows the initial condition for the free-surface elevation. The analytical solution is derived in an ideal condition that the fluid is inviscid and the bottom friction is zero. These conditions have been also applied in the model. The wave amplitude and the

Figure 2: The initial condition of the free-surface in the 3D standing wave test case

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00990v1-figure-52d55ee0-f5cf-4e8c-b55e-cbc51f5ee4ca.jpg)

basin dimensions have been taken as 0.1m and 500 × 500 × 10 m. The computational space has been discretized with the cells of dimensions ∆ x = 10 m, ∆ = 10 y m, ∆ = 1 z m, and the time step was 0.05s.

Figure 3(a) shows the numerical prediction of the free-surface elevations compared with the analytical solution at t = 6 T . As shown in this figure, wave damping after six oscillations is small and acceptable. In figure 3(b), numerical and analytical solutions of the free-surface elevation are presented at x = 0 , y = 0 . Figure 3(c) shows the predicted and analytical values of the horizontal velocity u at x = 250 m, y = 0 . It can be seen that the simulated free-surface and horizontal velocity agree well with the analytical solution.

## 4.2 Wind-Driven Circulatory Flow

The second test case was a wind-driven circulatory flow in a closed rectangular basin. In this test, the advection and the horizontal diffusion were neglected and a linearized bottom friction was used. The analytical solution for a constant vertical eddy viscosity and a linearized bottom friction is [12]:

$$u = \frac { 1 } { 6 \nu _ { z } } g \frac { \partial \xi } { \partial \tilde { x } } \left ( 3 z ^ { 2 } - H ^ { 2 } \right ) + \frac { \tau _ { w } } { 2 \rho \nu _ { z } } \left ( H + 2 z \right ),$$

$$\frac { \partial \xi } { \partial x } = \frac { 3 } { 2 } \frac { \tau _ { w } } { \rho g H } \frac { ( 2 \nu _ { z } + k H ) } { ( 3 \nu _ { z } + k H ) },$$

where H is the water depth and k is the linearized bottom friction coefficient. In this test, τ w was set to 0.1 N/m 2 in both numerical and analytical solutions. The other parameters were: H = 40 m, ν z = 0 03 . m /s, 2 k = 0 005 . m/s, and ρ = 1000 Kg/m . The basin dimensions were 2500×2500×40m, which were discretized into the cells with 3 ∆ x = 50 m, ∆ x = 50 m, ∆ z = 2 m, and the time step was 2s. A comparison of the horizontal velocity profile between the numerical model and the analytical solution is shown in Figure 4.

## 4.3 Efficiency of the Method

Assuming that the number of horizontal layers is N , any methods employed to couple the momentum and the continuity equations will finish at O N ( ) , O N ( 2 ) , etc. Generally, the smaller the order of complexity of the algorithm, the faster it will run. The standing wave test case in section 4.1 was repeated for a different number of horizontal layers to investigate the order of complexity of the method and all CPU times were recorded.

Figure 5(a) shows the simulation time for a different number of horizontal layers. As shown in this figure, the simulation time increases linearly and the method is of the least order of complexity O N ( ) .

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00990v1-figure-05da5927-096a-407b-b2f2-67086d8e7164.jpg)

Time (s)

Figure 3: Comparisons of simulated free-surface and horizontal velocity u with the analytical solution for the 3D linear standing wave test case: (a) free-surface elevations at y = 0 and t = 6 T ; (b) free-surface elevation at ( x, y ) = (0 , 0) ; (c) horizontal velocity u at ( x, y ) = (250 0) , .

For the purpose of comparison, a fast tridiagonal matrix inversion algorithm [20] with the order of complexity O N ( 2 33 . ) was also used to couple the equations in the test case. The result shows how important is the order of complexity of the employed coupling method in the computational efficiency of the model.

## 5 Summary

A semi-implicit model is developed to simulate three-dimensional HPEs. The pressure gradient terms and the vertical viscosity terms are discretized implicitly in the horizontal momentum equations, while the other terms are treated explicitly with respect to time. A method is presented to couple the horizontal momentum equation and the depthintegrated continuity equation which constructs a linear system to calculate the free-surface elevations. In the proposed algorithm, some simple recursive formulas have been formulated to express the horizontal velocity components in terms of the neighbor free-surface elevations which allows them to be directly substituted in the depth-integrated continuity equation.

The algorithm was validated by two idealized classical test cases: a uni-nodal standing wave in a rectangular basin and a wind-driven circulatory flow in a closed rectangular basin. It has been shown that the computation time increases linearly when the number of horizontal layers increases. It means that the order of complexity of the method is O N ( ) . Although the most expensive part of computations in each time step is to solve the five-diagonal linear system of the free-surface, it has been shown that the order of complexity the employed method to couple the equations considerably affects the computational efficiency of the model.

Figure 4: Comparison of the simulated horizontal velocity profile with the analytical solution at x = 1250 m, y = 1250 m for the wind-driven circulatory flow case

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00990v1-figure-281509ea-1999-47a0-9638-4d78a35201fa.jpg)

Figure 5: Computation time versus the number of horizontal layers

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00990v1-figure-93a40438-6617-4b21-a684-860b08f3cfa7.jpg)

## References

- [1] Anderson J. (1995). Computational Fluid Dynamics (Vol. 206). New York: McGraw-Hill. doi:https://doi.org/10.1201/9781351124027.
- [2] Blumberg A. F., &amp; Mellor G. L. (1987). A description of a three-dimensional coastal ocean circulation model. Three-dimensional coastal ocean models, 4, 1-16. doi:https://doi.org/10.1029/CO004p0001.
- [3] Blumberg A. F. (1994). A primer for ECOM-si. Technical report of HydroQual 66.
- [4] Casulli V., &amp; Cheng R. T. (1992). Semi-implicit finite difference methods for three-dimensional shallow water flow. International Journal for Numerical Methods in Fluids, 15(6), 629-648. doi:https://doi.org/10.1002/fld.1650150602.
- [5] Casulli V. (1997). Numerical simulation of three-dimensional free surface flow in isopycnal coordinates. International Journal for Numerical Methods in Fluids, 25(6), 645-658. doi:https://doi.org/10.1002/(SICI)10970363(19970930)25:6%3C645::AID-FLD579%3E3.0.CO;2-L.
- [6] Casulli V., &amp; Stelling G. S. (1998). Numerical Simulation of 3D Quasi-Hydrostatic Free-Surface Flows. Journal of Hydraulic Engineering, 124(7), 678-686. doi:https://doi.org/10.1061/(ASCE)0733-9429(1998)124:7(678).

- [7] Casulli V. (1999). A Semi-Implicit Finite Difference Method for Non-Hydrostatic Free-Surface Flows. International Journal for Numerical Methods in Fluids, 30(May 1998), 425-440. doi:https://doi.org/10.1002/(SICI)10970363(19990630)30:4%3C425::AID-FLD847%3E3.0.CO;2-D.
- [8] Casulli V., &amp; Walters R. A. (2000). An unstructured grid three-dimensional model based on the shallow water equations. International Journal for Numerical Methods in Fluids, 32(3), 331-348. doi:https://doi.org/10.1002/(SICI)1097-0363(20000215)32:3%3C331::AID-FLD941%3E3.0.CO;2-C.
- [9] Chen C., Liu H., &amp; Beardsley R. C. (2003). An unstructured finite volume three-dimensional primitive equation ocean model: Application to coastal ocean and estuaries. Journal Atmospheric and Oceanic Technology, 20, 159-186. doi:https://doi.org/10.1175/1520-0426(2003)020%3C0159:AUGFVT%3E2.0.CO;2.
- [10] Chen X. (2003). A free-surface correction method for simulating shallow water flows. Journal of Computational Physics, 189(2), 557-578. doi:https://doi.org/10.1016/S0021-9991(03)00234-1.
- [11] Fringer O. B., Gerritsen M., &amp; Street R. L. (2006). An unstructured-grid finite-volume nonhydrostatic parallel coastal ocean simulator. Ocean Modelling, 14(3-4), 139-173. doi:https://doi.org/10.1016/j.ocemod.2006.03.006.
- [12] Huang, H. &amp; Spaulding, M. (1995). A three-dimensional numerical model of estuarine circulation and water quality induced by surface discharges. Estuarine, Coastal and Shelf Science, 40, 357-380.
- [13] Madala R. V., &amp; Piacseki S. A. (1977). A semi-implicit numerical model for baroclinic oceans. Journal of Computational Physics, 23(2), 167-178. doi:https://doi.org/10.1016/0021-9991(77)90119-X.
- [14] Marshall J., Hill C., Perelman L., &amp; Adcroft A. (1997). Hydrostatic, quasi-hydrostatic, and nonhydrostatic ocean modeling. Journal of Geophysical Research: Oceans, 102(C3), 5733-5752. doi:https://doi.org/10.1029/96JC02776.
- [15] Mei C. C. (1989). The applied dynamics of ocean surface waves. Singapore: World Scientific.
- [16] Shahabi, A. and Ghiassi, R. (2022). A robust second-order godunov-type method for Burgers' equation. International Journal of Applied and Computational Mathematics, 8(2), 82. Springer. doi:https://doi.org/10.1007/s40819021-01171-7.
- [17] Shchepetkin A. F., &amp; McWilliams J. C. (2005). The regional oceanic modeling system (ROMS): a splitexplicit, free-surface, topography-following-coordinate oceanic model. Ocean Modelling, 9(4), 347-404. doi:https://doi.org/10.1016/j.ocemod.2004.08.002.
- [18] Smagorinsky J. (1963). General Circulation Experiments With the Primitive Equations. Monthly Weather Review, 91(3), 99-164.
- [19] Taylor G. I. (1916). Skin friction of the wind on the earth's surface. Proc. R. Soc. Lond. A, 92(637), 196-199.
- [20] Usmani R. A. (1994). Inversion of a tridiagonal jacobi matrix. Linear Algebra and Its Applications, 212-213, 413-414.
- [21] Vreugdenhil C. B. (1994). Numerical methods for shallow-water flow (Vol. 13). Springer Science &amp; Business Media.
- [22] Zhang Y., Baptista A. M., &amp; Myers E. P. (2004). A cross-scale model for 3D baroclinic circulation in estuary-plume-shelf systems: I. Formulation and skill assessment. Continental Shelf Research, 24(18), 2187-2214. doi:https://doi.org/10.1016/j.csr.2004.07.021.
Ali Shahabi Reza Ghiassi
A Numerical Technique for Coupling the Momentum and the Continuity Equations for Semi-Implicit 3D Ocean Models
Semi-implicit methods are powerful and efficient tools for the
three-dimensional modeling of coastal and oceanic processes. A semi-implicit
finite difference method for 3D hydrostatic primitive equations is presented in
this paper. The governing equations are time-discretized with an implicit
treatment of barotropic pressure gradient and vertical viscosity. The
discretized momentum equations along a water column are coupled with the
depth-integrated continuity equation of the column to construct a linear system
for free-surface elevations. The novelty of this work lies in formulating an
efficient method with the order of complexity O(N) for coupling the momentum
and the continuity equations. In this method, the horizontal velocity
components are expressed in terms of neighbor free-surface elevations by some
simple recursive formulas and then are substituted in the integrated continuity
equation.