TAMU Coastal Oceanography Modeling Group
2. Methodology
-
- A. Primitive Equations
-
- B. The Smagorinsky Scheme, The Mellor and Yamada Level-2.5
Turbulence Closure Model, and The Multidimensional Positive Definite
Advection Transport Algorithm
-
- C. Vertically Integrated Vorticity Balance Equation
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}