The model used in this study is called ECOM-si,
a semi-implicit variant
of the three-dimensional
Estuary, Coastal
and Ocean Model (ECOM) described by Blumberg and Mellor (1987).
ECOM-si was selected because it
includes a free surface, nonlinear advective terms, coupled
density and velocity fields, river runoff, heating and cooling
of the sea surface, a 2.5 level turbulence closure scheme
to represent vertical mixing (Mellor and Yamada, 1982),
and is designed to easily allow ``realistic'' simulations.
In addition, the combination of orthogonal curvilinear coordinates
in the horizontal plane and sigma-coordinates in the vertical
dimension allows grid refinement in regions of interest
without sacrificing the well-known characteristics of Cartesian grid
schemes.

The basic equations are expressed in a sigma coordinate system

where is the bottom topography and is the
surface elevation. The basic governing equations are presented
here in Cartesian coordinates to facilitate discussion. The
equations as expressed in curvilinear
coordinates may be found in Blumberg and Mellor (1987).

The continuity equation is

the x momentum equation is

and the y momentum equation is

where is the surface elevation,
u and v are the x and y components of velocity,
D is the total water depth ,
is the transformed vertical velocity (normal to
sigma surfaces),
is the vertical eddy viscosity,
is the water density,
is a reference water density,
and and are the horizontal viscous terms defined by

and

where is the horizontal viscosity. The parameterizations of
and are discussed in chapters 2.3 and 2.4,
respectively. The model also solves prognostically for temperature,
salt, concentration (used here to model sewage effluent),
turbulence kinetic energy and turbulence macroscale.

ECOM-si differs from the Blumberg and Mellor (1987) ECOM model in that
it uses a semi-implicit scheme for calculating
the free surface, therefore avoiding the gravity wave CFL condition
required by explicit schemes (eg Casulli, 1990).
This has the advantage that larger
time steps may be taken (on the order of minutes, rather than tens of
seconds). A potential disadvantage of implicit schemes is that they
more readily damp free wave motions, but in strongly forced and damped
shallow regions such as Massachusetts Bay, the effect is small. This
was determined by halving the time step and observing negligible differences
in simulation results. Another disadvantage is that because the
calculation of surface elevation requires solving a large matrix
equation at each time step and efficient solution of this equation
requires positive definiteness; boundary conditions for elevation
must be formulated in matrix form and must not destroy the positive
definiteness of the matrix. We use a combination of clamped
and gravity wave radiation conditions on the open boundary, made
possible by the implementation of the partially-clamped formulation
of Blumberg and Kantha (1985) discussed in chapter 2.8.