The primitive equations describe the behaviour of a geophysical fluid at space and time scales larger than a few kilometres in the horizontal, a few meters in the vertical and a few minutes. They are usually solved at larger scales: the specified grid spacing and time step of the numerical model. The effects of smaller scale motions (coming from the advective terms in the Navier-Stokes equations) must be represented entirely in terms of large-scale patterns to close the equations. These effects appear in the equations as the divergence of turbulent fluxes ( fluxes associated with the mean correlation of small scale perturbations). Assuming a turbulent closure hypothesis is equivalent to choose a formulation for these fluxes. It is usually called the subgrid scale physics. It must be emphasized that this is the weakest part of the primitive equations, but also one of the most important for long-term simulations as small scale processes in fine balance the surface input of kinetic energy and heat.
The control exerted by gravity on the flow induces a strong anisotropy between the lateral and vertical motions. Therefore subgrid-scale physics D , and in (2.1a), (2.1d) and (2.1e) are divided into a lateral part D , and and a vertical part D, and . The formulation of these terms and their underlying physics are briefly discussed in the next two subsections.
The model resolution is always larger than the scale at which the major sources of vertical turbulence occur (shear instability, internal wave breaking...). Turbulent motions are thus never explicitly solved, even partially, but always parameterized. The vertical turbulent fluxes are assumed to depend linearly on the gradients of large-scale quantities (for example, the turbulent heat flux is given by , where is an eddy coefficient). This formulation is analogous to that of molecular diffusion and dissipation. This is quite clearly a necessary compromise: considering only the molecular viscosity acting on large scale severely underestimates the role of turbulent diffusion and dissipation, while an accurate consideration of the details of turbulent motions is simply impractical. The resulting vertical momentum and tracer diffusive operators are of second order:
Lateral turbulence can be roughly divided into a mesoscale turbulence associated with eddies (which can be solved explicitly if the resolution is sufficient since their underlying physics are included in the primitive equations), and a sub mesoscale turbulence which is never explicitly solved even partially, but always parameterized. The formulation of lateral eddy fluxes depends on whether the mesoscale is below or above the grid-spacing ( the model is eddy-resolving or not).
In non-eddy-resolving configurations, the closure is similar to that used for the vertical physics. The lateral turbulent fluxes are assumed to depend linearly on the lateral gradients of large-scale quantities. The resulting lateral diffusive and dissipative operators are of second order. Observations show that lateral mixing induced by mesoscale turbulence tends to be along isopycnal surfaces (or more precisely neutral surfaces McDougall ) rather than across them. As the slope of neutral surfaces is small in the ocean, a common approximation is to assume that the `lateral' direction is the horizontal, the lateral mixing is performed along geopotential surfaces. This leads to a geopotential second order operator for lateral subgrid scale physics. This assumption can be relaxed: the eddy-induced turbulent fluxes can be better approached by assuming that they depend linearly on the gradients of large-scale quantities computed along neutral surfaces. In such a case, the diffusive operator is an isoneutral second order operator and it has components in the three space directions. However, both horizontal and isoneutral operators have no effect on mean ( large scale) potential energy whereas potential energy is a main source of turbulence (through baroclinic instabilities). Gent and Mcwilliams  have proposed a parameterisation of mesoscale eddy-induced turbulence which associates an eddy-induced velocity to the isoneutral diffusion. Its mean effect is to reduce the mean potential energy of the ocean. This leads to a formulation of lateral subgrid-scale physics made up of an isoneutral second order operator and an eddy induced advective part. In all these lateral diffusive formulations, the specification of the lateral eddy coefficients remains the problematic point as there is no really satisfactory formulation of these coefficients as a function of large-scale features.
In eddy-resolving configurations, a second order operator can be used, but usually the more scale selective biharmonic operator is preferred as the grid-spacing is usually not small enough compared to the scale of the eddies. The role devoted to the subgrid-scale physics is to dissipate the energy that cascades toward the grid scale and thus to ensure the stability of the model while not interfering with the resolved mesoscale activity. Another approach is becoming more and more popular: instead of specifying explicitly a sub-grid scale term in the momentum and tracer time evolution equations, one uses a advective scheme which is diffusive enough to maintain the model stability. It must be emphasised that then, all the sub-grid scale physics is included in the formulation of the advection scheme.
All these parameterisations of subgrid scale physics have advantages and drawbacks. There are not all available in NEMO. In the -coordinate formulation, five options are offered for active tracers (temperature and salinity): second order geopotential operator, second order isoneutral operator, Gent and Mcwilliams  parameterisation, fourth order geopotential operator, and various slightly diffusive advection schemes. The same options are available for momentum, except Gent and Mcwilliams  parameterisation which only involves tracers. In the -coordinate formulation, additional options are offered for tracers: second order operator acting along surfaces, and for momentum: fourth order operator acting along surfaces (see §9).
The lateral Laplacian tracer diffusive operator is defined by (see Appendix B):
For iso-level diffusion, and are zero. reduces to the identity in the horizontal direction, no rotation is applied.
For geopotential diffusion, and are the slopes between the geopotential and computational surfaces: they are equal to and , respectively (see (2.22) ).
For isoneutral diffusion and are the slopes between the isoneutral and computational surfaces. Therefore, they are different quantities, but have similar expressions in - and -coordinates. In -coordinates:
The normal component of the eddy induced velocity is zero at all the boundaries. This can be achieved in a model by tapering either the eddy coefficient or the slopes to zero in the vicinity of the boundaries. The latter strategy is used in NEMO (cf. Chap. 9).
The lateral bilaplacian tracer diffusive operator is defined by:
The Laplacian momentum diffusive operator along - or -surfaces is found by applying (2.7e) to the horizontal velocity vector (see Appendix B):
Such a formulation ensures a complete separation between the vorticity and horizontal divergence fields (see Appendix C). Unfortunately, it is not available for geopotential diffusion in coordinates and for isoneutral diffusion in both - and -coordinates ( when a rotation is required). In these two cases, the and fields are considered as independent scalar fields, so that the diffusive operator is given by:
As for tracers, the fourth order momentum diffusive operator along or -surfaces is a re-entering second order operator (2.41) or (2.41) with the eddy viscosity coefficient correctly placed:
geopotential diffusion in -coordinate:
geopotential diffusion in -coordinate:
Gurvan Madec and the NEMO Team
NEMO European Consortium2017-02-17