Modeling strategy

Following the observations of Chereskin (1983), Hibiya (1986, 1988) succesfully used a numerical method based on the direct solution of the 2D Euler equations to study the flow over Stellwagen Bank. The same idea was followed by Lamb in a series of papers on the generation of ISWs in the Georges Bank area. Solving the Euler equations in place of approximate equations (such as the KdV or eKdV, refs...) eliminates the (severe) restriction imposed by the several approximations involved in deriving simplified models, at the cost of forgoing analytical solutions and requiring rather expensive computations. Even though 3D effects both at large scales, such as the finiteness of the Bank, the curvature of the wave fronts as well as at small scales (turbulence) are likely important, a numerical solution of the 3D Euler or Navier-Stokes equations with sufficient resolution is out of question, at least for the imediate future. On the other hand, previous work by Hibiya and Lamb have shown that valuable insight can be obtained by looking at the 2D dynamics, provided, as is the case here, that the curvature of both the Bank and the shoaling are small compared to the horizontal distances involved (about 50 Km in the present case).

Hence, to model the shoaling of ISWs in Massachusetts Bay we consider a 2D channel with a variable bottom (see fig. 6). The section between $x_1=0$ Km and $x_2=45$ Km realistically models the bottom topography observed along the track connecting the moorings. East of $x_2$ the bottom is flattened out. We assume that the flow is governed by the 2D Euler equations for conservation of mass and momentum, using the Bousinnesq approximation, and we fix a frame of reference $(x,z)$ with $x$ positive in the offshore direction, $z$ positive upward in the vertical direction and the surface located at $z=0$, and the flow occupies the region $0\geq z\geq D(x),0\geq x\geq L$ (see Figure 6). The channel is open on the eastern side (right), while is it closed on the western side. To account for rotation, we allow the existence of a velocity $v$ normal to the plane $xz$, but we restrict all quantities to depend on $x$ and $z$ alone. We write the total density as the sum of a steady part plus a perturbation term $\rho=\rho_0(1+\overline\rho(x,z)+\rho'(x,z,t)),$ and we solve for $\rho'$. The flow is bounded by a rigid, porous lid at the surface and by the topography at the bottom, with zero normal flow at the bottom boundary and a prescribed velocity at the top. The latter trick removes the extra complication of having to deal with a time dependent domain without substancially alter the physics.

If we introduce a streamfunction $\psi$, such that

$$(u,w)=\left({\partial\psi\over\partial z},-{\partial\psi\over\partial x}\right),$$ (1)

and let $\omega=-\nabla^2\psi$ be the vorticity in the direction normal to the plane of the flow, the governing equations in the f-plane approximation can be written as

$${\partial \rho'\over\partial t}-J(\psi,\rho') = J(\psi,\overline\rho) ,\cr {\partial \omega\over\partial t}-J(\psi,\omega)+f{\partial v\over\partial z}$$ (2)

where $J(\cdot,\cdot)$ denotes the Jacobian operator, $V$ is the geostrophic jet associated to horizontal mean density gradients.

The flow is forced by prescribing the total (irrotational) transport as a function of space and time. In this formulation, this amounts to setting $\psi=0$ at the bottom and at the left boundary, $\psi=Q_t(t)x/L,$ at the surface and $\psi=Q_t(t)z/(D(L))$ at the right boundary, where $Q_t(t)$ is the flux entering the domain from the right. Notice that this choice of boundary conditions implies an irrotational flow which is consistent with a surface moving at speed $d\eta/dt=Q_t(t)/L$. In this way, the barotropic tide goes from being maximum at the eastern side to zero as we approach the shore.