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