Modelling the ocean has become an essential component of coastal and navigational hazard prevention (beach erosion, pollutant transport, tidal or storm surge, ice drift, wave height). Moreover, the ocean being a large component of the global climate, to model its circulation is essential to obtaining a better understanding of the dramatic climatic changes which either have occurred or might occur. Finally, some theoretical studies require the use of models in order to understand fundamental physical processes which involve nonlinear dynamics and/or complex geometries, and which are beyond analytical approaches.

The first attempts to model the ocean circulation were made in the middle of the 20th century, after the work of Ekman Ekman05 who recognized the importance of the wind as the major source of mechanical forcing. Sverdrup sverdrup47 derived a simple law that relates the ocean currents to the curl of the wind friction. Stommel stom48 and Munk munk50 derived analytical models of the wind induced ocean circulation in closed rectangular basins using simplified dissipative laws. Both the Stommel and Munk models afford simple explanation of the westward intensification of oceanic currents, such as observed for the Gulf Stream or the Kuroshiwo. An important threshold in computer performances was reached in the late sixties, and this allowed for the first full prognostic three-dimensional studies of the ocean circulation [Bryan and Cox, 1967,Bryan, 1969]. These models were driven by mechanical forcing (the winds and a bottom drag) and by fluxes of salt and heat exchanged with the atmosphere. They could take into account the complexity of the geometry and the nonlinear nature of the oceanic currents. Ideally, these models should be able to fill the gaps in the data and give reasonable estimates of the ocean circulation. However, because of their inherent complexity, the poor knowledge of numerous physical processes and problems with the definition of coastlines, straits and sea-mounts, they drift easily from any reasonable state if no restoring terms are added to the equations for temperature and sanility. Thus, prognostic three-dimensional models may sometimes look like expensive interpolators and yield no very different results than simpler inverse models do. Nonetheless, they have produced useful estimates of the role of the oceans in the thermal global budget and the importance of the so-called conveyor belt.

Our main concern is the representation of irregular domains in numerical ocean models, and their influence on the dynamics of the currents. Models, so far, have only crudely represented these irregular boundaries, either in the vertical (the topography) or in the horizontal (the coastlines). Our objective in this thesis is to evaluate the accuracy of conventional numerical methods in the presence of irregular coastlines and to introduce more accurate alternatives.

Furthermore, we suspect that irregular coastlines have important but sometimes under-estimated influences of the dynamics of the currents flowing along them. The energetics of these currents are controlled by the transfer of energy to smaller scales by nonlinear interactions. These interactions are likely to take place along the western part of the oceanic basins where the currents are the strongest. In particular, as the geostrophic balance (the main assumption governing the dynamics of the currents) in these regions breaks down, we hint at possible interactions between the geostrophic and ageostrophic modes. We therefore focus on models which are simple enough to represent geostrophic-ageostrophic interactions on a large range of scales. A shallow water model seems appropriate for these goals. The dynamics are only two dimensional which allows for very high resolutions in the horizontal directions, as opposed to more complex three-dimensional models where the same level of resolution would be too expensive. In one particular context, we introduce a quasi-geostrophic (QG) model for comparison. This particular type of model represents only the geostrophic motions. Previous studies focusing on nonlinear interactions were usually done in very idealized and regular domains using QG models. Hence, we hope to make an interesting contribution by conducting relatively simple experiments involving idealized but still irregular coastlines and somewhat more complex dynamics compared to QG models.

The problem of accurately representing the geometry of the domain in ocean
models is divided in two sub-problems: representing the bottom
boundary (the topography) and representing the lateral boundaries (the
coastline). Topographical features (sills or sea-mounts) are essential to
the mixing of waters of different properties, origins and depths, and,
therefore, their influence extends to the largest scales. Coastlines
partially enclose the oceanic basins. Their presence is essential to the
comprehension of the oceanic currents (such as the westward
intensification of currents). It was early realized that a crude
vertical representation of the topography (** z**- or geopotential vertical
coordinate) could be detrimental to an accurate modelling of the ocean
circulation. In particular, waters of different properties tend to mix
over sills with dramatic consequences for the global ocean circulation if
the vertical discretization is too crude. To remedy this problem, vertical
terrain following coordinates were proposed, despite various known
limitations. However, the horizontal discretization has not received the
same level of scrutinity. Most of the modern oceanic models still crudely
represent coastlines. A crude horizontal discretization has
several consequences. Straits may be under-resolved and the associated
exchange of water modified: The strait of Gilbratar controls the
Mediterranean salt input into the Atlantic; the Bering strait controls the
freshwater input between the Arctic and the Pacific and the Indonesian
archipelago is notorious for being the location of the so-called return
flow of the vast thermo-haline conveyor belt which circles the globe.
A crude horizontal representation has also retardation effects for fast
oceanic modes (Kelvin waves) which propagate along coastlines
[Schwab and Beletsky, 1998]. For the wind-driven circulation, little is known about
the influence of crude horizontal representations.

In order to study the influence of the choice of the numerical method, we propose to test several of them, and investigate which one best handles irregular coastlines. We therefore propose to test different staggerings of the finite difference (FD) method and several finite element (FE) formulations against a spectral element method. The test-cases we choose are very idealized in order to focus only on the dynamical aspects of two-dimensional flows (no physical parameterizations except for constant dissipative coefficients) and range from simple linear and non-linear test-cases in square domains, to linear and non-linear test-cases in smoothly irregular domains. The finite difference models range from the conventional Arakawa C-grid (preferentially used for regional studies -- e.g. Bleck and Boundra, 1981; Blumderg and Mellor, 1983), to the conventional Arakawa B-grid (preferentially used for global studies as in Bryan-Cox derived models -- e.g. Bryan, 1969; Cox, 1984), to the unconventional A-grid (Dietrich et al., 1993).

We choose to use the FE method because it has the decisive
characteristic of representing boundaries more efficiently than the more
conventional FD method and because it presents variable
resolution capabilities. Oceanographers, especially the tidal
community use FE models
to represent tidal interactions and resonances which occur at different
scales, from ocean basin to coastal inlets [Connor and Wang, 1974,Lynch and Gray, 1979,Walters and Cheng, 1979]. Some of the modern FE
models are derived from the earlier models formulated by Lynch and Werner lynch87 and
Le Provost *et al.* lepro94. Others used the QG
approximation and proposed a finite element formulation of the vorticity
equation for the general ocean circulation [Fix, 1975,Dumas *et al. , *1982,Myers, 1995]; results
were very encouraging. Unfortunately, no general circulation model based
on the primitive equations (explained below) has been proposed based on
the FE method and we try to determine the reasons for this
relative failure. On the other hand, the spectral element method (an
extension of the finite element method) has been used with a relative
success by Iskandarani et al. iskan95. Their method is based
on quadrangular elements; instead, we favor the use of triangular elements
which offer increased geometrical flexibility. Specifically, we propose to
test a spectral element model based on this discretization technique.
The apparent advantage of the spectral element method lies in the accepted
advantage of spectral methods (the accuracy and the fast convergence with
increasing resolution for regular problems) and the flexibility of an
irregular grid. However, as with the spectral method, there is always the
possibility that Gibbs oscillations appear when the fields being
approximated are too irregular or under-resolved (the classical example is
the step-function). We try to solve this problem by use of an adaptive
method which increases the resolution (the number of triangles) in regions
where the largest errors in the solution are observed (to be defined
later). Finally, since finite element methods are potentially more costly
than conventional finite difference methods due to the need for more
matrix inversions, the spectral element method may be a good alternative
because its enhanced accuracy (compared to finite elements) is not
severely offset by an excessive cost. In order to verify this statement,
we give an accuracy-to-cost function for all models.

Modelling the ocean is very challenging due to the coexistence of many
physical processes at various spatial and time scales, from the lowest
scales (salt intrusion and viscous boundary layers of a few centimeters),
to surface waves induced by wind or wave breakings, tides, geostrophic
eddies, to the general ocean circulation. Since all these processes can
interact with each other, it is virtually impossible to reproduce and
isolate with a high degree of realism any of these processes. Most often,
approximations and parameterizations are used to represent the small scale
(or sub-grid) phenomena and limit the explicit motions of the model to the
scales of interest. In this hierarchy of approximations, the barotropic
QG model represents the leading largest scale and lowest
frequency approximation. There is no vertical structure and no horizontal
divergent motions such as gravity waves. Then, somewhat more complex is
the shallow water model. The variables are the horizontal velocity,
**( u,v)** and the elevation of the free-surface, .
It allows for
divergent motions but still does not permit vertical structure. Layered
models are extended versions of shallow water models and allow for crude
vertical (baroclinic) variations. For more realistic vertical structures,
the so-called primitive equations are used. They are based in the
incompressible Navier-Stokes equations and use the Boussinesq and
hydrostatic assumptions. Further improvement can be gained by using an
non-hydrostatic model which can represent the small scale convection.
However, the limitation imposed by computer performance fixes the length
scales and the physical processes which can be
explicitly resolved. Some features,
such as synoptic eddies, are very difficult to resolve in global
circulation models. These eddies are of the order of 10 to 100 km and are
relatively small compared to the basin scale (10,000 km). Nonetheless,
some authors [Holland and Lin, 1975,Treguier, 1992] stress the importance of
representing explicitly the role of the eddies in the transfer of energy
between the different scales and their positive influence on the (more
realistic) mean fields. These processes can not be perfectly mimicked by
the alternative strategy of using eddy-parameterizations and, therefore,
this strategy is argued to be flawed [Lesieur, 1997]. Moreover, these
parameterizations use coefficients difficult to adjust to real
observations when these coefficients are not simply ``cosmetic''.
Therefore, a good general ocean model should be eddy-resolving. However,
since the required resolution for a general circulation model of the ocean
is too high (10 km at mid-latitudes), models should have, at least,
variable resolution capabilities, in the sense that they should
have capabilities
to follow and resolve isolated eddies or westward boundary currents, while the
rest of the domain is discretized at a coarser resolution.

It may be important to represent other physical processes. The Topex-Poseidon
satellite mission, for instance, renewed interest in the surface ocean
large scale activity: tides, Kelvin and Rossby waves and synoptic eddies
whose signature were measurable on the surface elevation field of the
ocean. Hence, a good general or regional ocean model should also have a
moving free-surface which allows for fast barotropic gravity waves. This
was also pointed out by studies of the vertical eddy-viscosity over the
rough topography of the Atlantic ridge [Polzin *et al. , *1997]. The large values
of the eddy-viscosity observed over the ridge, probably induced by
external tides, imply that the general circulation must interact with the
tides, i.e. the QG physics
interacts with the large scale gravity waves (ageostrophic horizontally
divergent
dynamics). From a theoretical point of view, it seems also more and more
necessary to include ageostrophic motions in numerical models, even in the
extra-tropics. The difficulty comes from explaining the cascade of energy
down to the molecular viscosity scale where the energy can be dissipated.
Indeed, the two dimensional (and QG) dynamics tend to cascade the energy
up to the Rhines arrest scale (50 to 200 km) in typical basins and not
down. Therefore, there is no clear mechanism that cascades down and
dissipates the energy in QG dynamics. This mechanism may
come from the non-linear interactions between the geostrophic and
ageostrophic modes. This may be visible from spectral analysis
[Stammer, 1997] which show no particular cut-off frequencies or
wave-numbers separating geostrophic and ageostrophic modes
^{1.1}.

The presence of irregular coastlines may also be important for the interactions of the geostrophic and ageostrophic modes. First, because it provides a forcing source at various wave numbers and, moreover, because the westward side of ocean basins is the location where the geostrophic approximation is most likely to break down, i.e., where the transfer of energy is most likely to occur. These preceding arguments imply that general circulation models should allow for the interactions between the geostrophic and ageostrophic motions. The simplest system that allows for such interactions is the shallow water equation system.

Also for physical reasons, even eddy-resolving models need an explicit
parameterization of dissipation. Although very crude, this is usually done
through an explicit eddy-viscosity (Laplacian operator) parameterization.
Such a parameterization requires an arbitrary choice for the dynamical
boundary condition at the walls (a problem which is exacerbated when
higher order dissipations are employed). We consider herein only two
boundary conditions. One is the free-slip boundary condition and
corresponds to fluids being free to slip along lateral boundaries. There
is some ambiguity as to the precise definition of free-slip. Pedlosky
(1987, p. 183) takes the point of view that it corresponds
to there being no viscous flux of tangential momentum across the boundary
(i.e.,
at the
wall). We take a less stringent definition by simply forcing the normal
derivative of the tangential velocity to be zero (for instance,
on a meridionally oriented wall). The latter choice is
the one generally found in the literature. On straight walls, the two
definitions are equivalent, and correspond to vertical vorticity vanishing
at the wall. On curved boundaries, either definition mentioned above
results in non-zero relative vorticity ^{1.2}. This bears some
dynamical consequences that we discuss below. The second boundary
condition is the so-called no-slip boundary condition which corresponds to
fluids that do not slip along walls (the tangential velocity is zero), and
leads to strong shear along walls. This boundary condition is considered
to be the ``real'' one because it is the one observed in laboratory
experiments at microscales. However, at the resolution used in modern
ocean models (1 to 100 km), it is not clear which, if either, boundary
condition is appropriate. The present trend in ocean modelling is to go
towards higher Reynolds number (smaller lateral viscosity) and higher
resolution along with no-slip boundary conditions. But the level of
resolution is still very far from sufficient to represent realistic
viscous boundary layers, although inertial boundary layers are definitely
becoming more realistic. On the other hand, authors experienced problems
with the free-slip boundary condition. In an idealized square basin and
barotropic ocean, it is observed that under single-gyre wind forcing and
decreasing eddy-viscosity, the oceanic currents tend to jump to
unrealisticly high values with no signs of transient eddies, whatever the
resolution of the model is. Hence, the free-slip boundary condition
prevents the transient activity which usually allows for reasonable mean
currents by transporting the excess of negative vorticity from the
interior, through the inertial layer to the walls, where it is dissipated
[Pedlosky, 1996]. Using a barotropic QG approximation, Dengg
dengg92 showed that free-slip flows tend not to separate from
a cape compared to the clear separation observed with no-slip flows,
whatever the value of the wind forcing. For all of these reasons, it
seems safer to use the no-slip boundary condition, even if it requires
unrealisticly large viscosity values in order to resolve the boundary
viscous sub-layer. Nonetheless, these studies fail to realize some
important issues. The absence of transients and separation, under the
free-slip boundary condition, is connected to the fact that most of those
models fail to produce relative vorticity at the walls. This is
because the relative vorticity at the boundary is specified to be zero for
reasons of simplicity. Therefore, those studies fail to note that, even
under free-slip, flows can produce relative vorticity simply
because of the coastline curvature. Hence, absence of transients in a
square basin is due to the idealized straight walls. If the walls were
curved (or, in more general sense, irregular, as they are in nature) ,
there is a chance that transients would appear and play the important role
of transporting excess of negative vorticity from the interior to the
walls. For the same reason, separation of the oceanic currents around a
cape can occur because the cape is round and can produce the necessary
vorticity required for separation. Of course, due to their
fractal nature, the curvature of the coastline depends on the sampling
resolution chosen to represent the coastline. Therefore, the knowledge
about coastline curvature is subjective. This, in itself, would be a good
reason for not considering the free-slip boundary condition for practical
ocean modelling. Nonetheless, we would like to revisit the debate between
free-slip and no-slip and investigate if, at least from a theoretical
point of view (when the curvature is known), use of free-slip is
permissible. For this reason, we discard all baroclinic processes and only
consider an idealized model of the ocean, the so-called shallow water
model, in presence of irregular boundaries and driven solely by wind
friction. The description of both geostrophic and ageostrophic motions in
this model allows for observations of the interactions between the two
kinds of motions which may be important, especially in presence of
irregular boundaries. In terms of physics and geometrical
representations, our study contrasts with and should be an improvement
over earlier theoretical studies based on the QG approximation and
rectangular basins.

The thesis is organized as follows. In Chapter 2, we present the different numerical methods and, in Chapter 3, we test them for simple test-cases in order to understand the effective truncation order and cost of these methods in presence of irregular domains. In Chapter 4, we further analyze the issue of discretization in FD models and, in particular, how it relates to vorticity budgets of the whole basin. In Chapter 5, we investigate the inertial run-away problem under free-slip boundary conditions in irregular domains. Conclusions are presented in Chapter 6.

**Get the PS or PDF version here
**