**2. The governing equations and basic ocean dynamics**

In order to fully appreciate the role and importance of numerical ocean models, it is helpful to first understand some of the basic dynamics of the ocean circulation. Mathematically, investigating the ocean circulation can be considered as solving an initial boundary value problem described by the Navier Stokes equations. These form a set of nonlinear, partial differential equations which describes the motion of any Newtonian fluid. The core of this is the three dimensional equation for the conservation of momentum which is essentially an expression of Newton's second law of motion. The two fundamental forces that must be considered are the pressure gradient force and gravity. For geophysical fluids, rotation of the Earth is also important and therefore Coriolis force must also be added to the equations. To complete the description of the motion equations for mass conservation (continuity) and for the conservation of internal energy must also be added. The latter can be expressed in terms of density or in terms of temperature and salinity. To make these equations more tractable and directly applicable to the ocean circulation, various simplifications and approximations are applied. These simplifications are usually based on a scale analysis of the various terms in the equations. The two most common approximations are: (1) the vertical extent or depth of the fluid layer is much smaller than the horizontal scale of motion, and (2) the Boussinesq approximation in which the density variations are assumed to be small compared to the mean value and are therefore neglected except in the buoyancy term of the equation. As a result of the first approximation, the vertical component of the conservation of momentum can be reduced to a diagnostic equation for hydrostatic balance (i.e., the vertical component of the pressure gradient force exactly balances gravity or the weight of the fluid). The second approximation, which is roughly equivalent to assuming that seawater is incompressible, means that mass continuity can be reduced to a diagnostic equation for the conservation of volume (i.e., three dimensional nondivergence). The final set of the governing equations (usually referred to as the primitive equations) in Cartesian coordinates (x, y, z), includes seven equations as follows:

*Horizontal momentum* 

112 Hydrodynamics – Natural Water Bodies

primarily a descriptive science based on sparse and scattered observations. The quantitative aspects of physical and dynamical oceanography saw a major breakthrough with the publication of Henry Stommel's seminal work on the North Atlantic circulation (Stommel, 1948). With a simple mathematical model of the wind driven circulation he was able to elegantly explain the phenomenon of westward intensification (i.e., the formation of strong western boundary currents such as the Gulf Stream) as a result of the meridional variation

The idea of using numerical models to further expand the understanding of the intricacies and complexities of the ocean circulation was introduced nearly twenty years later in the pioneering work of Bryan & Cox (1967). As with Stommel's research, they too investigated the circulation of the North Atlantic Ocean which at the time was the most highly observed ocean basin. The purpose of the model was to solve an initial value problem based on a simplified version of the Navier-Stokes equations. Through their model they were able to study the interaction between the wind driven and the thermohaline components of the circulation. Their work drew heavily from the experience of numerical weather prediction which took nearly thirty years to develop the capability of producing skillful forecasts beginning with Richardson's (1922) original concept but unsuccessful attempt and continuing to Charney et al. (1950) producing the first successful 24 hr forecast. As computational capabilities have increased exponentially over the past thirty years, so too has ocean modeling developed from a tool for simplified and focused process studies to fully operational forecasting systems. In this sense, the distinction between process studies (or simulations) and a forecasting system can be explained as follows. In the former, the goal is to understand the physical basis of the process without regard to reproducing specific details at any particular instant in time. In the latter, attention is focused on being able to produce the most accurate simulation of a particular realization of the flow at a specific time. The development of models for process studies and simulations was a necessary step in the development of forecasting systems. Furthermore, the useful range of a forecast, which is closely related to the limit of predictability, is limited by the chaotic behavior of the fluid flow. One the other hand, longer term simulations for the projection of future climate change is perhaps the most common example today of a process study. In both modern process studies and forecasting systems, the initial focus of model development has been the circulation, but today major progress has been made in developing components for simulating and predicting the fundamental biogeochemical processes of the oceanic

The goal of this chapter is to present an overview of modern ocean modeling as a tool for basic research as well as for operational forecasting. Considering the rapid developments and extensive experience of the Mediterranean oceanographic research community from recent years, we will use the Mediterranean as the prototype to explain and demonstrate

In order to fully appreciate the role and importance of numerical ocean models, it is helpful to first understand some of the basic dynamics of the ocean circulation. Mathematically, investigating the ocean circulation can be considered as solving an initial boundary value problem described by the Navier Stokes equations. These form a set of nonlinear, partial differential equations which describes the motion of any Newtonian fluid. The core of this is

of the Coriolis force.

ecosystem as well.

these capabilities and successes in ocean modeling.

**2. The governing equations and basic ocean dynamics** 

$$\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} + w \frac{\partial u}{\partial z} - fv = -\frac{1}{\rho\_0} \frac{\partial p}{\partial x} + DIFF(u) \tag{1}$$

$$\frac{\partial v}{\partial t} + u \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y} + w \frac{\partial v}{\partial x} + fu = -\frac{1}{\rho\_0} \frac{\partial p}{\partial y} + DIFF(v) \tag{2}$$

Where u, v, w are the velocity components in the x, y, z directions, t is time, ρ is the density (the subscript 0 indicates the mean value), � = ������ is the Coriolis parameter (Ω is the rotation rate of the Earth and φ is the latitude), p is the pressure, and *DIFF*(ψ) is the diffusion given by ����(�) <sup>=</sup> � �� ��� �� ��� + � �� ��� �� ��� + � �� ��� �� �� �, where *Ah* and *Az* are the horizontal and vertical diffusion coefficients, respectively;

*Vertical momentum* (hydrostatic equation)

$$\frac{\partial p}{\partial x} = -\rho g \tag{3}$$

Where g is gravity;

*Mass continuity* 

$$
\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} = \mathbf{0} \tag{4}
$$

**Conservation of internal energy** (can be written in terms of density or temperature and salinity)

$$\frac{\partial T}{\partial t} + \mu \frac{\partial T}{\partial x} + \nu \frac{\partial T}{\partial y} + \mathcal{w} \frac{\partial T}{\partial z} = DIFF(T) \tag{5}$$

Numerical Modeling of the Ocean Circulation:

North America.

**3. Numerical ocean modeling** 

(2006).

From Process Studies to Operational Forecasting – The Mediterranean Example 115

component of diffusion. The solution is the so called Ekman spiral for the surface layer in which the current magnitude decays with depth and the current vector rotates clockwise with depth in the northern hemisphere. Another example is the investigation of vertical convective mixing and its role in the deepening of the surface mixed layer through the use of the one-dimensional version of Eq. (5) and (6) in which the local time derivative is

In contrast to the vertical column models, other processes in which the horizontal variations are important or of interest can be investigated using two dimensional, depth integrated versions of the equations. To study the wind driven gyres in the upper ocean, (Stommel (1948) and Munk (1950) both started with the geostrophic form of Eqs. (1) and (2) with the addition of a frictional drag term as an alternative to horizontal diffusion. They took advantage of the non-divergence of the geostrophic flow to recast the equations as a single equation for vorticity. Through the solutions of the equations they were able to explain the underlying dynamics of the observed circulation in the North Atlantic Ocean. The general anticylonic (clockwise) gyre was driven by the curl of the wind stress (i.e., change in the direction from the easterly Trade Winds in the tropics and subtropics to the Westerlies in the mid-latitudes). The appearance of the intense western boundary current (i.e., the Gulf Stream) was a result of intensification due the accumulation of anticylonic relative vorticity by the wind stress and anticylonic planetary vorticity due to meridional variations of the Coriolis parameter but bounded by the damping effect of friction with the east coast of

The various examples presented above are meant to demonstrate some of the basic and salient features of ocean processes which have been investigated over the past 100 years through the use of various simplified versions of the governing equations for geophysical fluid dynamics. It does not even scratch the surface of the vast body of scientific literature in this rapidly expanding and exciting field of study. An in depth survey of these processes can be found in many of the excellent modern textbooks published in recent years such as Vallis

As noted in the introduction, the rapid development of computer technology over the past few decades has encouraged the massive development and advancement of numerical ocean models since the original effort of Bryan & Cox (1967). The main advantage of numerical modeling as compared to simplified process studies is that the numerical models are based on the more complete form of the governing equations presented in the previous section. This allows us to investigate more complex flow regimes and processes than in the past. In fact some of the simplifications such as hydrostatic balance in the governing equations are also being removed in recent models, thereby restoring the full time dependent equation for the vertical component of velocity. This is driven by the interest in and capabilities to investigate and simulate smaller scale processes. A model has the potential to fill in the many gaps left by limited in situ observations, subject of course to the computational and mathematical limitations of any model. One disadvantage of using more and more complex models is that it becomes more difficult to isolate and understand specific dynamical processes and thereby we develop the tendency to use a model as a black box. Even when running the most complex models we must never lose sight of what exactly the model is doing. Successful completion of a simulation does not guarantee proper results.

balanced by the vertical component diffusion (e.g., Martin, 1985).

$$\frac{\partial \mathcal{S}}{\partial t} + \boldsymbol{u} \frac{\partial \mathcal{S}}{\partial \boldsymbol{x}} + \boldsymbol{\nu} \frac{\partial \mathcal{S}}{\partial \boldsymbol{y}} + \boldsymbol{w} \frac{\partial \mathcal{S}}{\partial \boldsymbol{z}} = DIFF(\mathcal{S}) \tag{6}$$

Where *T* and *S* are the temperature and salinity, respectively; *Equation of state* 

$$
\rho = \rho(\mathbb{S}, T, \mathfrak{p}) \tag{7}
$$

Details of the derivation of the governing equations can be found in any text book on geophysical fluid dynamics such as Cushman-Roisin (1994). In order to solve the equations it is necessary to specify appropriate spatial boundary conditions and the top, the bottom and the sides of the domain as well as initial conditions. There is no general formulation of the boundary conditions since they depend upon the particular problem being addressed. Examples of boundary conditions at the top include wind stress for the momentum equations, or heat and mass fluxes for the internal energy equations. The bottom boundary conditions usually consist of frictional drag and no vertical mass flux. Lateral boundary conditions may be as simple as no flow at the coastline or some type of wave radiation condition at an open lateral boundary which allows waves to escape with no reflection (e.g., Orlanski, 1976).

The equations as they appear above describe a wide range of atmospheric and oceanic motions (except sound waves which are filtered out by the Boussinesq approximation). To study particular phenomena or processes, they can be further simplified, usually through additional scale analysis which leads to neglecting other terms. In some cases analytical solutions can be found, but in most cases numerical approaches are necessary. A very powerful and widely used simplification of Eq. (1) and (2) is geostrophic flow in which the local time derivative, the nonlinear advections terms, and diffusion are neglected. The remaining leading order terms, which roughly balance each other, are the Coriolis force and the horizontal pressure gradient force (last term on the left hand side and first term of the right hand side, respectively). The immediate implication is that the currents must flow parallel to the isobars (lines of constant pressure) rather than from high pressure to low pressure zones as in non-rotating fluids. This also means that the currents can be diagnosed directly from the pressure or mass field. The practical importance of this is that it is much easier and cheaper to measure the mass field variables (i.e., temperature, salinity, and pressure) than to measure the motion field (currents). Consequently the vast majority of physical oceanographic measurements consists of the three dimensional distribution of temperature and salinity. Combining the geostrophic equations with the hydrostatic equation allows us to compute the vertical shear of the currents due to horizontal pressure or density gradients. The two main weaknesses of the geostrophic approximation are that it breaks down in tropical areas where the Coriolis force is very weak, and it does not allow for temporal changes.

Another common method for simplifying the equations is to reduce their spatial dimensionality. For example, the primary external forcing of the ocean originates in the atmosphere and is applied from above (winds and heat flux). Consequently, the vertical gradients of the primary dependent variables in Eq. (1), (2), (5), and (6) are much larger than the horizontal gradients. It is therefore quite common to study the importance of this stratification through the use of one-dimensional water column models. A classic example is the study of the wind forced surface boundary layer by Ekman (1905) in which Eqs. (1) and (2) are reduced to steady state equations balancing the Coriolis force with the vertical component of diffusion. The solution is the so called Ekman spiral for the surface layer in which the current magnitude decays with depth and the current vector rotates clockwise with depth in the northern hemisphere. Another example is the investigation of vertical convective mixing and its role in the deepening of the surface mixed layer through the use of the one-dimensional version of Eq. (5) and (6) in which the local time derivative is balanced by the vertical component diffusion (e.g., Martin, 1985).

In contrast to the vertical column models, other processes in which the horizontal variations are important or of interest can be investigated using two dimensional, depth integrated versions of the equations. To study the wind driven gyres in the upper ocean, (Stommel (1948) and Munk (1950) both started with the geostrophic form of Eqs. (1) and (2) with the addition of a frictional drag term as an alternative to horizontal diffusion. They took advantage of the non-divergence of the geostrophic flow to recast the equations as a single equation for vorticity. Through the solutions of the equations they were able to explain the underlying dynamics of the observed circulation in the North Atlantic Ocean. The general anticylonic (clockwise) gyre was driven by the curl of the wind stress (i.e., change in the direction from the easterly Trade Winds in the tropics and subtropics to the Westerlies in the mid-latitudes). The appearance of the intense western boundary current (i.e., the Gulf Stream) was a result of intensification due the accumulation of anticylonic relative vorticity by the wind stress and anticylonic planetary vorticity due to meridional variations of the Coriolis parameter but bounded by the damping effect of friction with the east coast of North America.

The various examples presented above are meant to demonstrate some of the basic and salient features of ocean processes which have been investigated over the past 100 years through the use of various simplified versions of the governing equations for geophysical fluid dynamics. It does not even scratch the surface of the vast body of scientific literature in this rapidly expanding and exciting field of study. An in depth survey of these processes can be found in many of the excellent modern textbooks published in recent years such as Vallis (2006).
