TAMU Coastal Oceanography Modeling Group

1. Numerical Model Formulation
A. Orthogonal Grid Generation
B. Model Parameters and Boundary Conditions
C. Data Sources
D. Lagrangian Particle Trajectories

A. Orthogonal Grid Generation

The present model employs an orthogonal curvilinear coordinate system in the horizontal. The advantage of the orthogonal curvilinear coordinate system is that it flexibly allows ones to control the grid resolution throughtout the model domain, where finer grid resolution prevails over frontal regions with strong gradients of velocity, and coarser grid resolution in others without increasing the computational cost. For application to the Texas-Louisiana shelf, the model uses a 172 x 66 curvilinear orthogonal grid in the horizontal (Fig. 4). \begin{figure} \centering \psfig{figure=figure_4.ps,height=9.0in,bbllx=1.25in,bblly=+0.25in,bburx=10in,bbury=10in} \addcontentsline{lof}{figure}{\protect\numberline{4}{The 172 $\times$ 66 curvilinear orthogonal model grid of the Texas-Louisiana shelf. The mesh resolution varies from 3.8 km on the inner shelf to 19.7 km near the open boundary.}} \addtocounter{figure}{4} \addtocounter{page}{0} \end{figure} The mesh resolution is non-uniform horizontally, finer on the inner shelf and gradually coarser near the open boundary, varying spatially from 3.8 to 19.7 km. The grid resolution over the inner shlef is considerably uniform with average spacing of about 6 km which is sufficient to resolve scales of $O(R_{i})\sim$ 20 - 50 km. To avoid unphysical waves caused by sharp change of grid resolution, a 10-grid-point transition zone with grid spacing increasing linearly is deployed near the open boundary (Chen, 1992). Also, the model consists of 14 levels with irregular spacing in the vertical, which is established by log distributions at the top and bottom and a linear distribution between (Mellor, 1996). Therefore, the distribution of the vertical $\sigma$ coordinate includes finer resolution near the surface and the bottom in order to resolve the surface wind-mixing layer and the bottom boundary layer (Table.~\ref{tab:sigma1}).

Table. 1. The vertical sigma coordinate distribution
Level sigma Level sigma
1 0.000 8 -0.500
2 -0.025 9 -0.600
3 -0.050 10 -0.700
4 -0.100 11 -0.800
5 -0.200 12 -0.900
6 -0.300 13 -0.950
7 -0.400 14 -1.000

\begin{table}[h] \caption{The vertical $\sigma$ coordinate distribution}\label{tab:sigma1} \begin{center} \begin{tabular}[h]{ccccc} \hline\hline Level & $\sigma$ & & Level & $\sigma$ \\ \hline 1 & 0.000 & \hspace{0.5in} & 8 & -0.500 \\ 2 & -0.025 & \hspace{0.5in} & 9 & -0.600 \\ 3 & -0.050 & & 10 & -0.700 \\ 4 & -0.100 & & 11 & -0.800 \\ 5 & -0.200 & & 12 & -0.900 \\ 6 & -0.300 & & 13 & -0.950 \\ 7 & -0.400 & & 14 & -1.000 \\ \hline \end{tabular} \end{center} \end{table} The program for generating the curvilinear orthogonal grids was initially developed by Wilkins (1988) using the power conformal transformation technique of Ives and Zacharias (1987). The general conformal mapping method can be interpreted by introducing the concept of complex coordinates (see, e.g., Reid, 1954). Let we assume that \begin{eqnarray} Z & = & x + i \ y \label{eq:z1} \\ \omega & = & \xi_{1} + i \ \xi_{2} \end{eqnarray} and \begin{eqnarray} \omega & = & f(Z) \label{eq:z2} \end{eqnarray} where $\xi_{1}$ and $\xi_{2}$, corresponding respectively to $u$ and $v$, are the orthogonal computational coordinates in the mapped complex $\omega$-plane. If we differentiate $\omega$ with respect to $x$ and $y$, then we can obtain \begin{eqnarray} \frac{\partial^{2}x}{\partial \xi_{1}^{2}} + \frac{\partial^{2}y}{\partial \xi_{1}^{2}} & = & 0 \label{eq:l1} \\ \nonumber \\ \frac{\partial^{2}x}{\partial \xi_{2}^{2}} + \frac{\partial^{2}y}{\partial \xi_{2}^{2}} & = & 0 \label{eq:l2} \end{eqnarray} which are the Laplace equations preserving orthogonality. Following Ives and Zacharias (1987), equation (\ref{eq:z1}) can be expressed in the polar coordinates $r, \vartheta$ as follows: \begin{equation} Z = x + i \ y = r \ cos\vartheta + i \ r \ sin\vartheta = r \ e^{i \vartheta} \ . \end{equation} Subsequently, equation (\ref{eq:z2}) is replaced by the power transformation written as \begin{equation} \omega = Z^{n} = r^{n} e^{i n \vartheta} \end{equation} and let \begin{equation} n = \frac{\varphi}{\alpha} \end{equation} where $\varphi$ is the desired angle at the hinge point in the $\omega$ plane, and $\alpha$ is the included angle at the hinge point in the $Z$ plane. To map a closed contour to a rectangle, the desired angle $\varphi$ is set to $\pi$ for all boundary points, except for the four corner points where it is equal to $\pi/2$. The boundary points between the physical and computational planes are thus solved by the power $n$ tranformation technique. Once the boundary points are determined in both the physical and computational planes, the grid is filled in orthogonally by solving the Laplace equations (\ref{eq:l1}) and ({\ref{eq:l2}) for the boundary value problems of $x$ and $y$, using the fast Poisson solver SEPELI of Swartztrauber and Sweet (1975). Since the procedure of the conformal mapping method is based on Cartesian coordinates, a tranformation between the spherical coordinates (in longitude and latitude) and a physical plane is required. Several map projection methods have been developed, but up to date any projection method is area-conserving in some certain regions. Three conformal map projections are generally used for mapping different regions in meteorology: the Lambert conformal, the polar steorographic and the Mercator projection methods. The Lambert conformal projection method is adopted here because it is area-conserving for middle latitude regions. The description of three projection methods is given in detail in the Appendix C.

B. Model Parameters and Boundary Conditions

As mentioned previously, the model employs a free surface boundary so that the dynamics of ocean motions resolved by the primitive equations or the shallow water model contains fast propagating surface gravity waves and slow moving internal gravity waves as well. For computational efficiency, the model utilizes a mode splitting technique (Simons, 1974; Madala and Piacsek, 1977) in which the vertically integrated equations (external mode) are separated from the three dimensional momentum equations (internal mode). As a result, the surface elevation is solved by the transport equation. Based on the {\em Courant-Friedrichs-Levy} (CFL) theory of the computational stability condition, the time step of the external mode, vertically integrated transport equations is according to \begin{equation} \Delta t_{e} \le \frac{1}{C_{t}} \left(\frac{1}{\Delta x^{2}} + \frac{1}{\Delta y^{2}} \right)^{-1/2} \end{equation} where $C_{t} = 2(gH)^{1/2} + \bar{u}_{max}$; $\bar{u}_{max}$ is the expected maximum average velocity. For the internal mode, it allows much larger time steps because the fast propagating external mode effects have been removed. The time step criteria is given as \begin{equation} \Delta t_{i} \le \frac{1}{C_{T}} \left(\frac{1}{\Delta x^{2}} + \frac{1}{\Delta y^{2}} \right)^{-1/2} \end{equation} where $C_{T} = 2C + u_{max}$; $C$ is the maximum internal gravity wave phase speed ($\sim 2 \ m s^{-1}$) and $u_{max}$ is the maximum advective velocity. In the POM model, the time dependent terms are discretized by a leap-frog or center-time scheme. Unphysically computational mode wave solutions can be introduced because the solutions at odd and even time steps can be slowly diverging or time-splitting. Technically, the time splitting can be removed by averaging or smoothing the numerical solution to essentially suppress the computational mode. Here the numerical solution is smoothed at each time step based on the Asselin (1972) filter given as follows: \begin{equation} T_{s} = T + \frac{\alpha}{2} ( T^{n+1} - 2T^{n} + T^{n-1} ) \end{equation} where $T_{s}$ is the smoothed solution and the parameter $\alpha = 0.10$. At the side walls, the normal gradients of salinity and temperature are set to zero so that no advective and diffusive heat and salt fluxes cross these solid boundary. Also, the velocities normal to the land boundaries are set to zero, and the tangential velocities are zero at a half-grid point inside the walls according to the C-grid system. One of the difficulties in modeling coastal circulation is implementation of the open boundary conditions. Without predescribed boundary values, the boundary values are usually solved by Sommerfeld-type radiation conditions. The present model consists of one open boundary on the south from Rio Grande, Texas to Apalachicola, Florida and off the 2000 m isobath in the central portion as shown in Fig. 4. Unknown variables at the open boundary are salinity, temperature, sea surface elevation, and the tangential components of external and internal velocities to the boundary. For the temperature and salinity the open boundary conditions are basically Sommerfeld radiation conditions with an upwind scheme sense (Mellor, 1996), which can be expressed as \begin{equation} \frac{\partial \phi}{\partial t} - c\frac{\partial \phi}{\partial y} = 0 \label{eq:sommer} \end{equation} here $\phi$ represents $\theta$ and $S$, and $c = v_{b+1}$; $b$ is the boundary point. \begin{eqnarray} \phi^{n+1}_{b} & = & \left\{ \begin{array}{ll} \phi^{n-1}_{b} - \mu \ (\phi^{n-1}_{b} - \phi_{BDY}) \ , & \mbox{if \ $c \ge 0$} \vspace{0.15in} \\ \phi^{n-1}_{b} - \mu \ (\phi^{n-1}_{b+1} - \phi^{n-1}_{b}) \ , & \mbox{if \ $c < 0$} \end{array} \right. \label{eq:cond} \end{eqnarray} where $\phi_{BDY}$ is given boundary values, and $\mu = 2 c \Delta t_{i}/\Delta y$. Also, if $c < 0$, a vertcal mixing is added as follows: \begin{equation} \phi^{n+1}_{b,k} = \tilde\phi^{n+1}_{b,k} - \frac{w^{n}_{b+1,k} \Delta t_{i}}{\Delta z} (\phi^{n}_{b+1, k-1} - \phi^{n}_{b+1, k+1}) \end{equation} where $\tilde\phi$ is the solution of (\ref{eq:cond}). For the surface elevation and external velocity, an implicit gravity-wave radiation condition is applied at the open boundary (Chapman, 1985), which is based on equation (\ref{eq:sommer}) except $\phi$ represents $\eta$ and $\bar{v}$, and $c = (gH)^{1/2}$ is wave phase speed. Thus, \begin{equation} \phi^{n+1}_{b} = (\phi^{n}_{b} + \mu_{e} \phi^{n+1}_{b+1})/(1+\mu_{e}) \end{equation} where $\mu_{e} = c \Delta t_{e}/\Delta y$. The internal part of the velocity field is then obtained by using an implicit Orlanski (1976) radiation boundary condition at the boundary, where \begin{equation} C_{L} = \frac{\phi^{n-1}_{b+1} - \phi^{n+1}_{b+1}}{\phi^{n+1}_{b+1} + \phi^{n-1}_{b+1} - 2\phi^{n}_{b+2}} \end{equation} and \begin{eqnarray} \mu_{i} & = & \left\{ \begin{array}{ll} 1 \ , & \mbox{if $C_{L} \ge 1$} \\ C_{L} \ , & \mbox{if $0 < C_{L} < 1$} \\ 0 \ , & \mbox{if $C_{L} \le 0$} \ . \end{array} \right. \end{eqnarray} Thus, the boundary condition becomes \begin{equation} \phi^{n+1}_{b} = \left[ \phi^{n-1}_{b}(1 - \mu) + 2 \mu \phi^{n}_{b+1} \right]/(1 + \mu) \end{equation} In addition to the radiation boundary conditions, a hyperviscous ``sponge'' layer (Israeli and Orszag, 1981; Chen, 1992; Oey, 1993) is applied to the external and internal velocity fields near the open boundary in order to release energy out of the model domain and to reduce wave reflection from the open boundary. The present model includes 6 grid points for the sponge layer along the open boundary, constructed as Newtonian damping terms, $-\gamma u$ and $-\gamma v$, in the horizontal momentum and the vertically integrated momentum equations (Sarmiento and Bryan, 1982), where energy gradually vanishes toward the open end and is passed out through the boundary by applying the radiation conditions. The Newtonian damping coefficient $\gamma$ is the form suggested by Israeli and Orszag (1981) or Chapman (1985) \begin{equation} \gamma(j) = \left[ \frac{\gamma_{max}}{n}(b-j) \right]^{m} \ , \ j = b \pm 1,2,\ldots n \end{equation} where $\gamma_{max}$ is the maximum $\gamma$ at outer edge of the sponge layer, and $n$ is the number of grid points for the sponge layer. The appropriate values of $\gamma_{max}$ and $m$ are found by trial and error, in which $\gamma_{max} = 0.4$ and $m = 2$. Also, all parameters used in the model are listed in (Table.~\ref{tab:para1}).
Table. 2. Parameters used in the present model
external time step 24 seconds
internal time step 360 seconds
coefficient used in Smagorinsky's formula 0.12
background vertical viscosity 2 x 10-5 m2 s-1
bottom friction coefficient 0.0025
bottom roughness parameter 0.003 m
filter coefficient 0.1
Newtonian damping coefficient 0.4
\begin{table}[h] \caption{Parameters used in the present model}\label{tab:para1} \begin{center} \begin{tabular}{lcccl} \hline\hline external time step & & $\Delta t_{e}$ & & 24 seconds \\ internal time step & & $\Delta t_{i}$ & & 360 seconds \\ coefficient used in Smagorinsky's formula & & $C$ & & 0.12 \\ background vertical viscosity & & $\nu$ & & $2 \times 10^{-5} m^{2} s^{-1}$ \\ bottom friction coefficient & & $C_{z}$ & & 0.0025 \\ bottom roughness parameter & & $z_{o}$ & & 0.003 m \\ filter coefficient & & $\alpha$ & & 0.1 \\ Newtonian damping coefficient & & $\gamma_{max}$ & & 0.4 \\ \hline \end{tabular} \end{center} \end{table} It has to be mentioned here that the transport equation for solving the surface elevation does not include a damping term because of keeping the conservation of mass. Freshwater discharging from two major rivers onto the Texas-Louisiana continental shelf are considered in the present study: the Atchafalaya and the Mississippi rivers. The model river inflows are represented by point sources at the coast or the outermost interior grid point according to the Arakawa C-grid system, where there are one grid cell for the Atchafalaya river and six grid cells for the Mississippi river as shown in Fig. 1. The model implementation of river discharge is a river/dam and onshore intake/outfall discharge method as used in the model ECOM-si (Blumberg, 1994), in which the volume transport of river discharge is $Q_{r}(t)$ with unit $[L^{3}T^{-1}]$, distributing over a grid cell. Hence, the external mode of the river inflow velocity $V_{re}(t)$ is obtained according to \begin{eqnarray} V_{re}(t) = \frac{Q_{r}(t)}{ D \ \Delta L_{x \ or \ y}} \end{eqnarray} and the internal mode of the river inflow velocity $V_{ri}$ is \begin{eqnarray} V_{ri,k}(t) = \frac{Q_{r}(t) \times P_{k}}{\Delta D_{k} \ \Delta L_{x \ or \ y}} \end{eqnarray} where $\Delta L_{x \ or \ y}$ is the grid spacing in the direction perpendicular to the river inflow velocity, $P_{k}$ is the percentage of river discharge at each vertical $\sigma$-level, and $\Delta D_{k}$ is the distance between two $\sigma$-levels. Also, predescribed values of salinity and temperature for the river inflow are required. Once river forcing is established, potential energy becomes available to drive density currents.

C. Data Sources

The bottom topography is obtained from the Dynalysis of Princeton high resolution bathymetry database with $0.01^{\circ} \times 0.01^{\circ}$ resolution, provided by James Herring, and an inverse distance square formula is used to interpolate that onto each model grid point. Since the present study is focused on the inner continental shelf, the model depth can be held constant at 500 m for actual depths greater than 500 m and the minimum depth is set to 3 m over the shallow water region, for numerical efficiency. Also, the bottom topography is slightly smoothed by Mellor's (1996) method (Fig. 5), \begin{figure} \centering \psfig{figure=figure_5.ps,height=9.0in,bbllx=-2.50in,bblly=+1.50in,bburx=6.50in,bbury=10.50in,angle=180} \addcontentsline{lof}{figure}{\protect\numberline{5}{Smoothed bottom topography. Depth contours are in meters. The circles indicate the hydrographic stations of LATEX A cruise H06; crosses indicate the hydrographic stations of LATEX B cruise 4; stars indicate the wind observing stations.}} \addtocounter{figure}{5} \addtocounter{page}{0} \end{figure} in which the bottom slope of any two adjacent grid points is less than or equal to a parameter, in order to reduce the numerical errors caused by the vertical $\sigma$-coordinate transformation over steep bottom topography (Haney, 1991). Basically, only the effects of river discharge and wind forcing to the inner-shelf circulation are concerned in the present study. The volume transports of the Atchafalaya and Mississippi rivers are included according to the daily observational data of U.~S. Geological Survey (1994) and U.~S. Army Corps of Engineers (1994) (Fig. 2). The percentage of the Mississippi river discharge distributing to each pass is listed in Table.~\ref{tab:rivt} (U.~S. Army Corps of Engineers, 1984).
Table. 3. The percentage of the Mississippi river discharge for each pass
river mouth % of the total river discharge
Southwest Pass 33
South Pass 18
Pass a la Loutre 30
Grand Pass 5
Baptise Pass 4
Main Pass 10
\begin{table} \caption{The percentage of the Mississippi river discharge for each pass}\label{tab:rivt} \begin{center} \begin{tabular}{clcccc} \hline\hline &river mouth & & & & \% of the total river discharge\\ \hline &Southwest Pass & & & & 33 \\ &South Pass & & & & 18 \\ &Pass \`{a} la Loutre & & & & 30 \\ &Grand Pass & & & & 5 \\ &Baptise Pass & & & & 4 \\ &Main Pass & & & & 10 \\ \hline \end{tabular} \end{center} \end{table} Surface meteorological data is obtained from the meteorological buoys deployed by the LATEX program and is supplemented by other meteorological datasets during the summer-fall period of 1993 as shown in Fig. 5 (see, Wang, 1996), in which hourly wind speed and direction were measured at about 3 m above sea level. The analytical results show that wind stress is considerably uniform over the shelf during the summer and fall period (Wang, 1996). Therefore, the vector wind stress is assumed to be spatially uniform over the entire model domain at any given time, which is taken from the data of the observing station L17 at (268.04$^{\circ}$E, 29.19$^{\circ}$N) July and that of the station L22 at (266.04$^{\circ}$E, 28.35$^{\circ}$N) Aug.-Oct., 1993. The initial conditions including salinity and temperature are used according to two kinds of data sets: the monthly mean climatological salinity and temperature database (Levitus, 1982), with resolution $1^{\circ} \times 1^{\circ}$, and the high-quality hydrographic data obtained from cruise H06 (26 July - 7 Aug., 1993) of LATEX-A (Texas-Louisiana Shelf Circulation and Transport Processes) at Texas A\&M University and that obtained from cruise 4 (13 - 21 July, 1993) of LATEX-B (Mississippi River Plume Hydrography) at Louisiana State University. The hydrographic data were taken by a Sea-Bird SBE-9 conductivity-temperature-depth (CTD) meter, where the locations of the hydrographic stations are shown in Fig. 5. The initial conditions are thus generated by an objective bi-linear method, composed from the hydrographic data from LATEX cruises, augmented by climatological data in regions of the model domain not covered by the LATEX data set.

D. Lagrangian Particle Trajectories

It is well known that wave propagation transports energy but not particles. Water particles move in an approximately circular path but not closed, and there exists a small component of forward propagation, particularly in large-amplitude waves. The small net forward displacement of particle propagation is called the Stokes drift, in which the Stokes drift velocity, derived from a first order approximation to the Lagrangian velocity in the derivatives of the Eulerian velocity field, can be analytically expressed in terms of the difference between the Eulerian mean and the Lagrangian mean velocity (Longuet-Higgins, 1969; Zimmerman, 1979); \begin{eqnarray} u_{\mbox{\scriptsize Lagrange}} = u_{\mbox{\scriptsize Euler}} + u_{\mbox{\scriptsize Stokes}} \ . \end{eqnarray} The Lagrangian trajectory ${\bf y}({\bf x}_{o}, t)$ of a particle starting at location ${\bf x}_{o}$ can be written as \begin{eqnarray} {\bf y}({\bf x}_{o}, t) = {\bf x}_{o} + \int^{t}_{0}{\bf u}_{\ell}({\bf x}_{o}, t^{'}) dt^{'} \end{eqnarray} where ${\bf u}_{\ell}$ is the Lagrangian velocity. In the present study, the Lagrangian velocity is obtained from that at any given time, solved by the momentum equations. Many studies suggest that the physical process of water mixing can be revealed by studies of particle trajectories such as ``chaotic stirring'' or turbulent dispersion caused by tide-induced forcing over variable bottom topography (Zimmerman, 1976, 1986; Ridderinkhof, 1992), and the exchange or mixing rate of tidal flow through a narrow strait (Awaji, 1980; Awaji, 1982). For all experiments, a group of Lagrangian particles, distributed from surface to bottom, are released initially in the region near the Mississippi river mouth in order to study (1) the process of water exchange in the frontal zone and (2) the influence of variable bottom topography on the mixing process. The arrangement of marked particles is 9 particles in a grid cell over 6 $\times$ 6 grid points horiztontally, initially deployed at the positions of near the surface, the middle depth, and the bottom in the vertical. As evolution in time, the velocity fields are given by a weighted trilinear interpolation using the velocity fields around particles, and the trajectories of suspended particles are calculated by a bilinear interpolation method.