**3. Numerical ocean modeling**

114 Hydrodynamics – Natural Water Bodies

డ௬ ݓ డௌ

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.,

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

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

డ௭ ൌ ܨܨܫܦሺܵሻ (6)

ߩ ൌ ߩሺܵǡ ܶǡ ሻ (7).

డ௫ ݒ డௌ

డௌ డ௧ ݑ డௌ

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

*Equation of state* 

Orlanski, 1976).

for temporal changes.

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.

Numerical Modeling of the Ocean Circulation:

invert a large tridiagonal matrix at every time step.

treated computationally the same as sub grid scale processes.

**4. The Mediterranean Sea – a laboratory ocean basin** 

energy conservation.

From Process Studies to Operational Forecasting – The Mediterranean Example 117

velocity components are shifted one half of a grid point in their respective directions relative to the mass variables (Arakawa & Lamb, 1977). This arrangement ensures that the numerical equations will preserve certain integral properties of the continuous equations such as

In the time integration of the equations, the computational stability of the differencing scheme must also be considered. This means that the time step and spatial grid spacing must be chosen in such a way as to properly resolve the motion of the fastest moving waves that can be simulated by the model (usually the free surface or external mode gravity wave). Most models use some type of centered explicit or split explicit scheme that is also second order accurate. In explicit schemes all variables can be advanced to the forward time step based on the values known at the present and/or backward time step. In split explicit schemes, the terms in the momentum equations that are identified with the fastest moving waves are integrated separately with a shorter time step. This is done to improve the numerical efficiency and execution time of the model. An alternative is to use a fully implicit time scheme in which advancing the model to the forward step involves simultaneously knowing the values of the variables at the forward, present, and backward time steps. This has the advantage of the scheme being absolutely stable (i.e., no numerical amplification) but the disadvantage of being computationally cumbersome and slow due to the need to

The final point to consider in the construction of an ocean model is how to account for the unresolved scales of motion. The grid spacing or resolution of a model limits the explicitly resolved processes. Strictly speaking, from a mathematical perspective the shortest length scale that can be explicitly resolved by a model is 2*Δx*. Practically however, scales shorter than ~4-6 *Δx* are often misrepresented due to numerical damping or phase speed errors which may be an inherent characteristic of certain differencing schemes. However there are also important processes which may occur on scales smaller than the grid spacing which affect the larger scale flow. The classic example of this is vertical convective mixing forced by the wind or by induced by static instability when the water is cooled from above. This mixing will transport various properties of the water but is accomplished by small scale turbulent eddies which are not usually explicitly resolved. Such processes are accounted for by adding sub grid scale parameterizations, often in the form of a diffusion term but with a diffusion coefficient that is several orders of magnitude larger than the value for molecular diffusion. The quasi-empirical method for computing these eddy diffusion coefficients is usually referred to as a turbulence close scheme (e.g., Mellor & Yamada, 1982; Pacanowski & Philander, 1981). An analogous term is usually included for horizontal mixing. Finally, ecosystem models, which are becoming more common components of ocean models, require the addition of advection-diffusion equations, similar to Eqs. (5) and (6), for the relevant biogeochemical variables in addition to all of the biogeochemical processes which are

Since the early 20th century the Mediterranean Sea was known to be a concentration basin where excess evaporation drives a basin wide thermohaline cell in which less saline water enters from the Atlantic Ocean through the Strait of Gibraltar (Nielsen, 1912). This surface water becomes more saline and denser. It sinks to a depth of ~ 250 m and then returns to the strait where it is carried by a subsurface outflow back to the Atlantic. During the past 25-30

We must always critically examine the results to be sure the model is doing what it should. With this in mind we present a very brief survey of the most commonly used numerical methods used in ocean modeling today.

The governing equations presented in the previous section form a set of time dependent, hyperbolic partial differential equations. Numerical methods for solving such equations have been developed and have appeared in the mathematical literature over many years. As noted in the introduction, the approach to constructing numerical ocean models and the choice of particular methods has benefited greatly from and closely followed the development of atmospheric models and numerical weather prediction, which preceded ocean modeling by 10-15 years. The most common method used today in ocean models is finite differencing. In recent years finite elements have also become popular. For various reasons, spectral methods have not been widely used, perhaps due to the associated difficulties of dealing with irregular boundaries (i.e., coastlines). Without loss of generality, we will use the finite difference method to illustrate the fundamental principles of a typical approach to ocean modeling. A detailed presentation of ocean modeling methodology and applications can be found in the recent books of Kampf (2009, 2010).

The first step in developing a model is to define the domain of interest and divide it into a discrete set of grid points in space. The goal of the model is to approximate the continuous equations with a compatible set of algebraic equations which are solved at the grid points. This requires accurate methods for representing the first and second order spatial derivatives based on the values of the relevant variables at the grid points. If we consider a dependent variable *y* as a function of the spatial coordinate *x*, say *y(x)*, then the gradient or first derivative �� �� ⁄ can be approximated from a Taylor series expansion around the *i*-th grid point *xi* as

$$\frac{\partial \boldsymbol{y}}{\partial \boldsymbol{x}} \approx \frac{\boldsymbol{y}\_{l+1} - \boldsymbol{y}\_{l}}{\Delta \boldsymbol{x}} + \mathcal{O}(\Delta \boldsymbol{x}) \tag{8}$$

or

$$\frac{\partial y}{\partial x} \approx \frac{y\_t - y\_{t-1}}{\Delta x} + \mathcal{O}(\Delta x) \tag{9}$$

which are referred to as the forward and backward differencing schemes, respectively, and where *Δx* is the grid spacing. These schemes are first order accurate as indicated by the truncation error *O(Δx)*. By averaging these two schemes we obtain the more accurate centered differencing scheme

$$\frac{\partial y}{\partial x} \approx \frac{y\_{l+1} - y\_{l-1}}{2\Delta x} + \mathcal{O}(\Delta x^2) \tag{10}$$

which is second order accurate. Higher order schemes that are even more accurate can be formed by various weighted combinations of the respective Taylor series expansions, although most models use centered differencing. The second derivative is approximated by the three point stencil

$$\frac{\partial^2 y}{\partial x^2} \approx \frac{y\_{l+1} - 2y\_l + y\_{l-1}}{\Delta x^2} \tag{11}$$

The location of the dependent variables is a matter of choice. They can all be co-located at the grid points, or can be located on a staggered grid in which certain variables are shifted by half of a grid point. A commonly staggered grid is the Arakawa C-grid in which the

We must always critically examine the results to be sure the model is doing what it should. With this in mind we present a very brief survey of the most commonly used numerical

The governing equations presented in the previous section form a set of time dependent, hyperbolic partial differential equations. Numerical methods for solving such equations have been developed and have appeared in the mathematical literature over many years. As noted in the introduction, the approach to constructing numerical ocean models and the choice of particular methods has benefited greatly from and closely followed the development of atmospheric models and numerical weather prediction, which preceded ocean modeling by 10-15 years. The most common method used today in ocean models is finite differencing. In recent years finite elements have also become popular. For various reasons, spectral methods have not been widely used, perhaps due to the associated difficulties of dealing with irregular boundaries (i.e., coastlines). Without loss of generality, we will use the finite difference method to illustrate the fundamental principles of a typical approach to ocean modeling. A detailed presentation of ocean modeling methodology and

The first step in developing a model is to define the domain of interest and divide it into a discrete set of grid points in space. The goal of the model is to approximate the continuous equations with a compatible set of algebraic equations which are solved at the grid points. This requires accurate methods for representing the first and second order spatial derivatives based on the values of the relevant variables at the grid points. If we consider a dependent variable *y* as a function of the spatial coordinate *x*, say *y(x)*, then the gradient or first derivative �� �� ⁄ can be approximated from a Taylor series expansion around the *i*-th

which are referred to as the forward and backward differencing schemes, respectively, and where *Δx* is the grid spacing. These schemes are first order accurate as indicated by the truncation error *O(Δx)*. By averaging these two schemes we obtain the more accurate

which is second order accurate. Higher order schemes that are even more accurate can be formed by various weighted combinations of the respective Taylor series expansions, although most models use centered differencing. The second derivative is approximated by

��� <sup>≈</sup> �������������

The location of the dependent variables is a matter of choice. They can all be co-located at the grid points, or can be located on a staggered grid in which certain variables are shifted by half of a grid point. A commonly staggered grid is the Arakawa C-grid in which the

�� � �(��) (8)

�� � �(��) (9)

��� � �(���) (10)

��� (11).

applications can be found in the recent books of Kampf (2009, 2010).

�� �� <sup>≈</sup> �������

�� �� <sup>≈</sup> �������

��

�� <sup>≈</sup> ���������

���

methods used in ocean modeling today.

grid point *xi* as

centered differencing scheme

the three point stencil

or

velocity components are shifted one half of a grid point in their respective directions relative to the mass variables (Arakawa & Lamb, 1977). This arrangement ensures that the numerical equations will preserve certain integral properties of the continuous equations such as energy conservation.

In the time integration of the equations, the computational stability of the differencing scheme must also be considered. This means that the time step and spatial grid spacing must be chosen in such a way as to properly resolve the motion of the fastest moving waves that can be simulated by the model (usually the free surface or external mode gravity wave). Most models use some type of centered explicit or split explicit scheme that is also second order accurate. In explicit schemes all variables can be advanced to the forward time step based on the values known at the present and/or backward time step. In split explicit schemes, the terms in the momentum equations that are identified with the fastest moving waves are integrated separately with a shorter time step. This is done to improve the numerical efficiency and execution time of the model. An alternative is to use a fully implicit time scheme in which advancing the model to the forward step involves simultaneously knowing the values of the variables at the forward, present, and backward time steps. This has the advantage of the scheme being absolutely stable (i.e., no numerical amplification) but the disadvantage of being computationally cumbersome and slow due to the need to invert a large tridiagonal matrix at every time step.

The final point to consider in the construction of an ocean model is how to account for the unresolved scales of motion. The grid spacing or resolution of a model limits the explicitly resolved processes. Strictly speaking, from a mathematical perspective the shortest length scale that can be explicitly resolved by a model is 2*Δx*. Practically however, scales shorter than ~4-6 *Δx* are often misrepresented due to numerical damping or phase speed errors which may be an inherent characteristic of certain differencing schemes. However there are also important processes which may occur on scales smaller than the grid spacing which affect the larger scale flow. The classic example of this is vertical convective mixing forced by the wind or by induced by static instability when the water is cooled from above. This mixing will transport various properties of the water but is accomplished by small scale turbulent eddies which are not usually explicitly resolved. Such processes are accounted for by adding sub grid scale parameterizations, often in the form of a diffusion term but with a diffusion coefficient that is several orders of magnitude larger than the value for molecular diffusion. The quasi-empirical method for computing these eddy diffusion coefficients is usually referred to as a turbulence close scheme (e.g., Mellor & Yamada, 1982; Pacanowski & Philander, 1981). An analogous term is usually included for horizontal mixing. Finally, ecosystem models, which are becoming more common components of ocean models, require the addition of advection-diffusion equations, similar to Eqs. (5) and (6), for the relevant biogeochemical variables in addition to all of the biogeochemical processes which are treated computationally the same as sub grid scale processes.
