Methodology

2. Methodology

The numerical model used in this study is a modified version of the Princeton Ocean Model (POM), a three-dimensional, primitive-equation, estuarine and coastal ocean circulation model, which was initially developed by Blumberg and Mellor (1983, 1987). The finite difference model uses a centered scheme in space, a leap-frog scheme in time, and an energy-enstrophy conserving Arakawa C-grid (Arakawa and Lamb, 1977). Incorporating a $\sigma$-coordinate (terrain-following) system in the vertical and a curvilinear coordinate system in the horizontal, the model is able to represent significant variability of realistic bottom topography and irregular coastal geometry. The model has a free surface boundary so that tidally driven currents and the propagation of surface gravity waves can be simulated.

Physically, turbulent mixing or horizontal diffusion is described as a function of velocity shear and a mixing or diffusion coefficient. Although there is some uncertainty in chosing diffusion coefficients, the model predicts the realistic parameterization of the vertical mixing according to the level-2.5 turbulence closure model of Mellor and Yamada (1974, 1982), which is based on the kinetic energy and the length scale of the turbulence. For additional details, see Appendix A. The horizontal diffusivity is calculated according to the Smagorinsky (1963) scheme, which depends on the velocity gradients and the horizontal grid spacing. In addition, the model has been modified by including a flux-corrected transport (FCT) scheme as used in the model ECOM-si (Blumberg, 1994). The FCT submodel, known as the Multidimensional Positive Definite Advection Transport Algorithm (MPDATA) (Smolarkiewicz, 1984; Smolarkiewicz and Clark, 1986; Smolarkiewicz and Grabowski, 1990), has been used for simulating the transport of salinity and temperature, and is second-order accurate in both time and space and also positive definite. For additional details, see Appendix B. The governing equations and numerical techniques are described in detail by Blumberg and Mellor (1987), Blumberg (1994) and Mellor (1996). A brief description of the numerical ocean model is given here.

A. Primitive Equations

The fundamental equations of the ocean motions are governed by the incompressible Navier-Stokes equations on a rotating spherical earth with two assumptions: hydrostatic assumption and Boussinesq approximation (see, e.g., Gill, 1982). Together with the salinity and temperature equations, the system is sometimes referred to as the primitive equations. In an orthogonal Cartesian coordinate system, in which $x, y$ and $z$ increase respectively eastward, westward and upward, the governing equations can be written as \begin{eqnarray} \frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} + w\frac{\partial u}{\partial z} - fv & = & -\frac{1}{\rho_{o}}\frac{\partial p}{\partial x} + \frac{\partial} {\partial z}\left(K_{m}\frac{\partial u}{\partial z}\right) + F_{x} \label {eq:mom1} \\ \nonumber \\ \frac{\partial v}{\partial t} + u\frac{\partial v}{\partial x} + v\frac{\partial v}{\partial y} + w\frac{\partial v}{\partial z} + fu & = & -\frac{1}{\rho_{o}}\frac{\partial p}{\partial y} + \frac{\partial} {\partial z}\left(K_{m}\frac{\partial v}{\partial z}\right) + F_{y} \label {eq:mom2} \\ \nonumber \\ \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} & = & 0 \label {eq:cont} \\ \nonumber \\ \rho g & = & -\frac{\partial p}{\partial z} \label {eq:hyd} \\ \nonumber \\ \frac{\partial\theta}{\partial t} + u\frac{\partial\theta}{\partial x} + v\frac{\partial\theta}{\partial y} + w\frac{\partial\theta}{\partial z} & = & \frac{\partial}{\partial z}\left(K_{h}\frac{\partial\theta} {\partial z}\right) + F_{\theta} \label {eq:tem} \\ \nonumber \\ \frac{\partial S}{\partial t} + u\frac{\partial S}{\partial x} + v\frac{\partial S}{\partial y} + w\frac{\partial S}{\partial z} & = & \frac{\partial}{\partial z}\left(K_{h}\frac{\partial S} {\partial z}\right) + F_{S} \label {eq:sal} \\ \nonumber \\ \rho & = & \rho\left(\theta,S\right) \end{eqnarray} where $u, v$ and $w$ are the $x, y$ and $z$ components of mean velocity, respectively; $\rho_{o}$ is the reference or mean density, $\rho$ the perturbation density, $g$ the gravitational acceleration, $p$ the pressure, $S$ the salinity, $\theta$ the potential temperature, $K_{m}$ the vertical eddy viscosity, and $K_{h}$ the vertical eddy viscosity for heat and salt. The potential density, $\rho_{total} = \rho_{o} + \rho$, is calculated using the potential temperature and salinity according to an equation of state (Fofonoff, 1962), in which $\rho$ is independent of pressure.

The hydrostatic assumption is based on an elementary scale analysis for large scale ocean flows, in which the horizontal scales are large in comparison with the ocean depth. As a result, the vertical momuntum equation can be simplified to the virtually hydrostatic relation (\ref{eq:hyd}). The Boussinesq approximation leads to the result that the density variations can be neglected in the horizontal momuntum equations but must be included in the hydrostatic equation, where buoyancy effects are important. Equation (\ref{eq:cont}) is the equation of continuity, simplified from the mass conservation equation by assuming that the fluid is incompressible and isentropic so that \begin{eqnarray} \frac{D\rho}{Dt} = 0 \end{eqnarray} where the material derivative is \begin{eqnarray} \frac{D(\ )}{Dt} = \frac{\partial(\ )}{\partial t} + u\frac{\partial(\ )}{\partial x} + v\frac{\partial(\ )}{\partial y} + w\frac{\partial(\ )}{\partial z} \ ; \end{eqnarray} i.e., that the velocity is nondivergent (see, e.g., Gill, 1982).

$F_{x}, F_{y}, F_{\theta}$ and $F_{S}$ represent the horizontal viscosity and diffusion terms, which can be written as \begin{eqnarray} F_{x} & = & \frac{\partial}{\partial x}\left(\tau_{xx}\right) + \frac{\partial}{\partial y}\left(\tau_{xy}\right) \label{eq:f1} \\ \nonumber \\ F_{y} & = & \frac{\partial}{\partial x}\left(\tau_{xy}\right) + \frac{\partial}{\partial y}\left(\tau_{yy}\right) \end{eqnarray} where \begin{eqnarray} \tau_{xx} & = & 2A_{m}\frac{\partial u}{\partial x} \\ \nonumber \\ \tau_{xy} & = & \tau_{yx} \ = \ A_{m}\left(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial y}\right) \\ \nonumber \\ \tau_{yy} & = & 2A_{m}\frac{\partial v}{\partial y} \end{eqnarray} Also, \begin{eqnarray} F_{\theta} & = & \frac{\partial}{\partial x}\left(A_{h}\frac{\partial \theta}{\partial x} \right) + \frac{\partial}{\partial y}\left(A_{h}\frac{\partial \theta}{\partial y}\right) \\ \nonumber \\ F_{S} & = & \frac{\partial}{\partial x}\left(A_{h}\frac{\partial S}{\partial x} \right) + \frac{\partial}{\partial y}\left(A_{h}\frac{\partial S}{\partial y}\right) \label{eq:f2} \end{eqnarray} where $A_{m}$ and $A_{h}$ are the horizontal diffusivities and are functions of $(x, y, t)$, obtained by the Smagorinsky's formula.

In the absence of surface and bottom heat fluxes, the vertical boundary conditions can be expressed as follows:

(i) at the free surface, $z = \eta (x,y,t)$ \begin{eqnarray} w = \frac{\partial \eta}{\partial t} + u\frac{\partial \eta}{\partial x} + v\frac{\partial \eta}{\partial y} \\ \nonumber \\ K_{m}\left(\frac{\partial u}{\partial z}, \frac{\partial v}{\partial z}\right) = (\frac{\tau_{sx}}{\rho_{o}}, \frac{\tau_{sy}}{\rho_{o}}) \\ \nonumber \\ K_{h}\left(\frac{\partial \theta}{\partial z}, \frac{\partial S}{\partial z}\right) = (0, 0) \end{eqnarray}

(ii) at the bottom, $z = - H (x,y)$ \begin{eqnarray} w = -u\frac{\partial H}{\partial x} - v\frac{\partial H}{\partial y} \\ \nonumber \\ K_{m}\left(\frac{\partial u}{\partial z}, \frac{\partial v}{\partial z}\right) = (\frac{\tau_{bx}}{\rho_{o}}, \frac{\tau_{by}}{\rho_{o}}) \\ \nonumber \\ K_{h}\left(\frac{\partial \theta}{\partial z}, \frac{\partial S}{\partial z}\right) = (0, 0) \end{eqnarray} where ($\tau_{sx}, \tau_{bx}$) and ($\tau_{sy}, \tau_{by}$) are the $x$ and $y$ components of the surface wind and the bottom stresses, respectively. The surface wind stress is calculated by \begin{eqnarray} \tau_{s} & = & \rho_{a} C_{D} W_{D}^{2} \end{eqnarray} where $\rho_{a}$ is the density of air, $C_{D} \simeq 1.4 \times 10^{-3}$ the drag coefficient and $W_{D}$ the wind speed in $m s^{-1}$. The bottom stress is determined according to the logarithmic ``law of the wall,'' given by \begin{eqnarray} K_{m} \left(\frac{\partial u}{\partial z}, \frac{\partial v}{\partial z}\right) & = & C_{z}\left[u^{2} + v^{2}\right]^{1/2} (u, v) \end{eqnarray} $C_{z}$ is the bottom drag coefficient, defined as \begin{eqnarray} C_{z} & = & max\left[\frac{\kappa^{2}}{[ln(z/z_{o})]^{2}}, \ 0.0025\right] \end{eqnarray} where $\kappa = 0.4$ is the von K\'{a}rm\'{a}n constant and $z_{o}$ is the roughness parameter.

The model employs a sigma $(\sigma)$ or terrain-following coordinate system in the vertical, in which local water depths are transformed into a non-dimensional, multi-layer system in which the vertical coordinate $\sigma$ ranges from $\sigma = 0$ at $z = \eta$ (at the surface) to $\sigma = -1$ at $z = -H$ (at the bottom) (Phillips, 1957; Blumberg and Mellor, 1987). In other words, the irregular bottom topography is mapped to a flat bottom domain based on the $\sigma$-coordinate transformation, which is given as

\begin{eqnarray} x^{*} = x, \ y^{*} = y, \ \sigma = \frac{z - \eta}{H + \eta}, \ t^{*} = t \ . \end{eqnarray} Let the total water depth $D = H + \eta$ and define a new vertical velocity as \begin{eqnarray} \omega = (\eta + H)\frac{D}{Dt^{*}}\left(\frac{z - \eta}{\eta + H}\right) \nonumber \end{eqnarray} or \begin{eqnarray} \omega & = & w - u\left(\sigma\frac{\partial D}{\partial x^{*}} + \frac{\partial \eta}{\partial x^{*}}\right) - v\left(\sigma\frac{\partial D}{\partial y^{*}} + \frac{\partial \eta}{\partial y^{*}}\right) - \left(\sigma\frac{\partial D}{\partial t^{*}} + \frac{\partial \eta}{\partial t^{*}}\right) \end{eqnarray} Physically, the new vertical velocity, $\omega$, is the velocity component normal to $\sigma$ surfaces and satisfies the vertical boundary conditions as \begin{eqnarray} \omega (x^{*}, y^{*}, 0, t^{*}) = 0, \ \omega (x^{*}, y^{*}, -1, t^{*}) = 0 \label{eq:wbc} \end{eqnarray}

In addition to the $\sigma$-coordinate system, a horizontal, orthogonal, curvilinear coordinate system is also used in the model, which enables variable grid resolutions and which can be applied to a shelf model with complex coastal geometry by finite difference methods. The expression or transformation of the governing equations in an orthogonal curvilinear coordinate is developed by Batchelor (1967) or Blumberg (1994) and is not repeated here.

B. The Smagorinsky Scheme, The Mellor and Yamada Level-2.5 Turbulence Closure Model, and The Multidimensional Positive Definit Advection Transport Algorithm

Mathematically the primitive equations are not closed because there are eleven unknown variables, $u, v, w, p, \rho, A_{h}, A_{m}, K_{h}, K_{m}, \theta$, and $S$, in seven equations. Two subset models, therefore, are introduced into the model: the Smagorinsky scheme (1963) and the level-2.5 turbulence closure model of Mellor and Yamada (1974, 1982). The horizontal diffusivities, $A_{m}$ and $A_{h}$, are determined according to the Smagorinsky scheme, which can be expressed as

\begin{eqnarray} A_{m} & = & A_{h} \ = \ C\Delta x\Delta y \left[\left(\frac{\partial u}{\partial x}\right)^{2} + \left(\frac{\partial v}{\partial y}\right)^{2} + \frac{1}{2}\left(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}\right)^{2}\right]^{1/2} \end{eqnarray} where the parameter $C$ is non-dimensional and typically in the range of 0.10 to 0.20 (Mellor, 1996). Here, spatially and temporally, $A_{m}$ and $A_{h}$ vary with the grid resolution and the local velocity gradients. They decrease as the grid resolution becomes finer and become small if the velocity gradients are small.

The vertical eddy viscosities, $K_{h}$ and $K_{m}$, are parameterized according to the level-2.5 (second-order) turbulence closure model of Mellor and Yamada (1982) as modified by Galperin et al. (1988), in which the vertical eddy viscosities are calculated based on the turbulent kinetic energy and the turbulent macroscale equations. The Mellor and Yamada level-2.5 turbulence closure model is derived by starting from the Reynolds stress and turbulent heat flux equations under the hypotheses of Rotta (1951), Kolmogorov(1941) and a near isotropic condition, where the Reynolds stress is generated due to the exchange of momentum in the turbulent mixing process. To make the turbulence equations closed, all empirical constants are obtained by assuming that turbulent heat production is primarily balanced by turbulent dissipation. As a result, the vertical eddy viscosities are determined according to the local Richardson number given as

\begin{eqnarray} R_{i} = - \frac{g}{\rho_{o}} \frac{\partial \rho}{\partial z} \left/ \left[\left(\frac{\partial u} {\partial z}\right)^{2} + \left(\frac{\partial v}{\partial z}\right)^{2}\right] \right . \end{eqnarray} A critical Richardson number, $R_{i} = 0.19$, is found, above that turbulence and mixing cease to exist (Mellor and Yamada, 1982). A detailed description of the level-2.5 turbulence closure model is given in the Appendix A.

In numerical modeling of ocean currents or geophyiscal phenomena, it is often necessary to solve the advection-diffusion transport equation for conservative tracers such as salinity or density, which are of course the positive-definite in nature. In fact, second- or high-order numerical schemes can produce unphysical oscillations appearing in regions of high-gradient density, such as ahead of and behind fronts, due to numerical dispersion. These dispersive ripples may cause positive-definte, conservative tracers to become negative near high-gradient regions. Hence, to ensure positivity of numerical solutions, the model incorporates a FCT-like, positive-definite submodel, the Multidimensional Positive Definite Advection Transport Algorithm (MPDATA) (Smolarkiewicz, 1984; Smolarkiewicz and Clark, 1986; Smolarkiewicz and Grabowski, 1990), to solve the transport equation for the density field.

The general concept behind the positive definite schemes is as follows: large numerical diffusion is first generated by a numerical truncation error to guarantee positivity and stability, and then a corrective term of anti-diffusion is introduced to balance out the numerical dissipation in the first step. The MPDATA scheme is developed in this manner and is based on an upwind scheme. In comparison with other positive-definite schemes, the MPDATA scheme is less time consuming because the advection transport equation is solved by the first-order upwind algorithm which is the most computationally efficient of any of the advective schemes. The accuracy of calculations can be increased by reapplying a number of corrective iterations, but that would considerably increase the time consumption of the scheme. The formulation of the MPDATA scheme is discussed in detail in the Appendix B.

C. Vertically Integrated Vorticity Balance Equation

The vorticity of a fluid is defined as the curl of the velocity field (see, e.g., Gill, 1982). The vorticity is a characteristic expression of the tendency for a fluid to rotate. For a three-dimensional case, the vertically integrated vorticity balance equation can be obtained by taking the curl of the vertically averaged momuntum equations (e.g., Ezer and Mellor, 1994). The vertically averaged equations in the $\sigma$ coordinate can be written as follows:

\begin{eqnarray} \frac{\partial \bar{u}D}{\partial t} + A_{x} - f\bar{v}D & = & \nonumber \\ & & \hspace{-1.5in}-Dg\frac{\partial \eta} {\partial x} - \frac{gD}{\rho_{o}}\int^{0}_{-1}\int^{0}_{\sigma}\left[D\frac{\partial \rho}{\partial x} - \frac{\partial D}{\partial x}\sigma^{'}\frac{\partial \rho}{\partial \sigma}\right]d\sigma^{'}d\sigma + \frac{1}{\rho_{o}}\frac{\partial}{\partial z}(\tau_{sx} - \tau_{bx}) \label{eq:vert1} \\ \nonumber \\ \frac{\partial \bar{v}D}{\partial t} + A_{y} + f\bar{u}D & = & \nonumber \\ & & \hspace{-1.5in}-Dg\frac{\partial \eta} {\partial y} - \frac{gD}{\rho_{o}}\int^{0}_{-1}\int^{0}_{\sigma}\left[D\frac{\partial \rho}{\partial y} - \frac{\partial D}{\partial y}\sigma^{'}\frac{\partial \rho}{\partial \sigma}\right]d\sigma^{'}d\sigma + \frac{1}{\rho_{o}}\frac{\partial}{\partial z}(\tau_{sy} - \tau_{by}) \label{eq:vert2} \end{eqnarray} where $A_{x}$ and $A_{y}$ represent the $x$ and $y$ components of horizontal advection and diffusion terms, respectively, and $\bar{u}, \bar{v}$ are the vertically average velocities. Based on Stokes' theorem, the vorticity can be evalutated by taking a line integral around each closed small element for equations (\ref{eq:vert1}) and (\ref{eq:vert2}) and then dividing by the area of the element.

For physical interpretation, the pressure gradient terms in equations (\ref{eq:vert1}) and (\ref{eq:vert2}) may be expressed as follows in a $z$-coordinate system.

\begin{eqnarray} \tilde{\Phi}_{x} & = & \int^{\eta}_{-H}\int^{\eta}_{z}\frac{\partial b}{\partial x}dz^{'}dz \label{eq:b1}\\ \nonumber \\ \tilde{\Phi}_{y} & = & \int^{\eta}_{-H}\int^{\eta}_{z}\frac{\partial b}{\partial y}dz^{'}dz \end{eqnarray} where $b = g\rho/\rho_{o}$ is the buoyancy. It can be shown that \begin{eqnarray} \tilde{\Phi}_{x} + Dg\frac{\partial \eta}{\partial x} & = & \frac{\partial\Phi}{\partial x} + D\frac{\partial P_{b}}{\partial x} \\ \nonumber \\ \tilde{\Phi}_{y} + Dg\frac{\partial \eta}{\partial y} & = & \frac{\partial\Phi}{\partial y} + D\frac{\partial P_{b}}{\partial y} \label{eq:b2} \end{eqnarray} where $\Phi = \int^{\eta}_{-H}zbdz$ is the potential energy and $P_{b} = g\eta + \int^{\eta}_{-H}bdz$ is the bottom pressure.

Now $\partial/\partial x(\ref{eq:vert2}) - \partial/\partial y(\ref{eq:vert1})$ gives the vertically integrated vorticity balance equation (using (\ref{eq:b1})$\sim$ (\ref{eq:b2})) \begin{eqnarray} \frac{\partial}{\partial t}\left(\frac{\partial\bar{v}D}{\partial x} - \frac{\partial\bar{u}D}{\partial y} \right) + \frac{\partial A_{y}}{\partial x} - \frac{\partial A_{x}}{\partial y} + \frac{\partial (f\bar{u}D)} {\partial x} + \frac{\partial (f\bar{v}D)}{\partial y} & = & \nonumber \\ \nonumber \\ & &\hspace{-4.0in}\frac{\partial P_{b}}{\partial x}\frac{\partial D}{\partial y} - \frac{\partial P_{b}} {\partial y}\frac{\partial D}{\partial x} + \frac{1}{\rho_{o}}\frac{\partial^{2}}{\partial x\partial z}( \tau_{sy} - \tau_{by}) - \frac{1}{\rho_{o}}\frac{\partial^{2}}{\partial y\partial z}(\tau_{sx} - \tau_{bx}) \ . \label{eq:vort} \end{eqnarray} The terms in the equation (\ref{eq:vort} are respectively referred to as the tendency term, the advection and diffusion terms, the Coriolis term, the bottom pressure torque term, and the surface and bottom stress curl terms. For a steady-state flow, the transport equation becomes $\partial(\bar{u}D)/\partial x + \partial(\bar{v}D)/\partial y = 0$, and thus the term containing the Coriolis parameter reduces to the beta term $\beta\bar{v}D$, where $\beta = df/dy$. Physically, the bottom pressure torque term contains a component of the JEBAR term that plays an important role in the potental vorticity balance, defined as \begin{eqnarray} \mbox{JEBAR} & = & J\left(\Phi, \frac{1}{D}\right) \end{eqnarray} where $J(A,B)$ is the Jacobian operator defined by \begin{eqnarray} J(A, B) & = & \frac{\partial A}{\partial x}\frac{\partial B}{\partial y} - \frac{\partial A}{\partial y}\frac{\partial B}{\partial x} \ . \end{eqnarray}