Circulation modelling

An ocean model is a numerical and computational tool used to understand and predict ocean variables (Griffies, 2006), providing a discrete solution of the geophysical fluid dynamic equations. It represents a rigorous way of linking the ocean state parameters through mathematical equations representing the physics that governs the oceans.

In the next subsections, we will introduce the different components of an OGCM, that is part of the OOFS (steps 1 and 2 as in Figure 4.1), focusing on mathematical equations, numerical methods, and spatial discretization techniques. A list of available numerical ocean models is provided in Table 5.1 in Section 5.4.3. Data assimilation methods used in OOFSs are instead presented in Section 5.5.

The Navier-Stokes equations represent the fundamental laws of fluid dynamics; they are based on conservation of momentum, conservation of mass, and an equation of state.

Oceans are also represented by the following equations (although with some significant simplifications as explained in Madec et al., 2022):

• Spherical Earth approximation: the geopotential surfaces are assumed to be oblate spheroids that follow the Earth’s bulge, and are approximated by spheres which gravity is locally vertical (parallel to the Earth’s radius) and independent from latitude;
• Thin-shell approximation: the ocean depth is neglected compared to the Earth’s radius;
• Turbulent closure hypothesis: the turbulent fluxes - which represent the effect of small-scale processes on the large scale - are expressed in terms of large scale features;
• Boussinesq hypothesis: density variations are neglected, except in their contribution to buoyancy force:
  
• Hydrostatic hypothesis: the vertical momentum equation is reduced to a balance between the vertical pressure gradient and the buoyancy force (this removes convective processes from the initial Navier-Stokes equations and so convective processes must be parameterized instead):
  
• Incompressibility hypothesis: the 3D divergence of the velocity vector U is assumed to be zero:
  
• Neglect of additional Coriolis terms: the Coriolis terms that vary with the cosine of latitude are neglected.

Because the gravitational force dominates in the equations of large-scale motions, it is useful to choose an orthogonal set of unit vectors (i,j,k) linked to the Earth such that k is the local upward vector and (i,j) are 2 vectors orthogonal to k. Let us define additionally: U the vector velocity, T the potential temperature, S the salinity, _ρ the insitu density. The vector invariant form of the primitive equations in the (i,j,k) vector system provides the following equations: 

• The momentum balance:
  
• The heat and salt conservation equations:
  

where ∇ is the generalised derivative vector operator in (i,- j,k) directions, t is the time, z is the vertical coordinate, ρ is the in-situ density given by Eq. 5.1, _ρ0 is the reference density, p is the pressure, f=2Ω∙k is the Coriolis acceleration (where Ω is the Earth’s angular velocity vector) and g is the gravitational acceleration. D^u, D^T and D^s are the parameterizations of small-scale physics for momentum, temperature and salinity, while F^u, F^T and F^s are surface forcing terms. 

OGCMs are able to resolve the mesoscale in some regions but not in others; additionally, once applied for climate research, they cannot entirely reproduce the rich mesoscale eddy activity we observe in reality. For this reason, mixing associated with sub-grid scale turbulence needs to be parameterized. 

A common problem an ocean modeller is facing when he/ she deals with primitive equations is the numerical discretization in space and time. As described in Hallberg (2013), numerical ocean models need to represent the effects of mesoscale eddies, which are the typical horizontal scales of less than 100 km and timescales in the order of a month. When defining the spatial grid for the numerical integration of the primitive equations, it is important to account for the ratio of a model’s grid spacing to the deformation radius, defined as:

where c_g is the first-mode internal gravity wave speed, f is again the Coriolis parameter, and β is its meridional gradient (Chelton et al., 1998). 

Figure 5.2 shows the ocean model resolution required for the baroclinic deformation radius to be twice the grid spacing, based on an eddy-permitting ocean model after one year of spin-up from climatology (Hallberg, 2013).

The next step towards the setup of a numerical model is the discretization phase, which involves the spatial discretization and the equation discretization.

The spatial discretization consists in defining a grid or mesh that would represent the space continuum with a finite number of points where the numerical values of the physical variables must be determined. In Section 5.4.2.1-2, basic concepts for dealing with horizontal grids and vertical discretization will be introduced. Once the mesh is defined, we move to the final step related to the primitive equations discretization by using numerical methods, which consist in transforming the mathematical model into an algebraic, nonlinear system of equations for the mesh-related unknown quantities. The concepts on the basis of the time stepping are treated in Section 5.4.2.3. With the definition of the time-dependent numerical formulation, we finally select the discretization method to use for the equations, described in Section 5.4.2.4.

In numerical methods, we can use:

• Structured grids
• Unstructured grids

A mesh is structured when the grid cells have the same number of sides and the same number of neighbouring cells. Typically, in ocean models three kinds of grids may be used (Figure 5.3): the Arakawa-A grid, the Arakawa-B grid and the Arakawa C-grid. In the Arakawa-A grid (Figure 5.3A), all variables are evaluated at the same location. Then, the B and C grids have been developed respectively for coarse and fine resolution models. In the Arakawa-B grid (Figure 5.3B) both u (Northwards current component, in orange) and v (Eastwards current component, in green), for example, are evaluated at the same point and the velocity points are situated at the point that is equidistant from the four nearest elevation points (Elevation, in blue). In the Arakawa-C grid (Figure 5.3C), the u points lie east and west of elevation points, while the vpoints lie north and south of the elevation points.

Unstructured grids (Figure 5.4C) allow one to tile a domain using more general geometrical shapes (most commonly triangles) that are pieced together to optimally fit details of the geometry. They are extremely attractive for ocean modelling, especially for coastal models, in which the high-quality representation of geometrical features of a given domain is essential, and from the numerical point of view they may reach a significant level of complexity (Griffies et al., 2000). 

Besides their ability to better represent coastlines, unstructured grid approaches also offer the possibility to smoothly increase the resolution over a region of interest or depending on physical parameters (Sein et al., 2017). This is also possible with structured curvilinear grids (for example, see the BLUElink Australian prediction model grid in Brassington et al., 2005, and Figure 5.4A), though with likely more constraints on the grid deformation properties. However, in any of the two cases, numerical stability is dictated by the smallest grid element, which substantially increases the computational problem. An additional difficulty is that sub-grid parameterizations have to be valid throughout the domain, whatever the grid size and eddy resolution regime are (Hallberg, 2013). In the structured grid case, block structured refinement techniques enable to circumvent some of the aforementioned difficulties by allowing a stepwise change (over a given grid patch) of the space and time resolutions (by integer factors, Figure 5.5B). Parameterizations and numerical schemes can also be changed accordingly. Grid exchanges can either be “one-way” if finer grids only receive information at their dynamical boundaries from the outer grid, or “two-way” if they also feed information back to the underlying mesh. In the latter case, data transferred at each model time step allows for a nearly seamless transition at the interface and possibly guarantees perfect conservation of prognostic quantities (Debreu et al., 2012). 

Several libraries do facilitate the implementation of block structured refinement. Among them, the AGRIF library (Debreu et al., 2008) has been successfully used in HYCOM, MARS, NEMO and ROMS models. It is noteworthy that refinement techniques can eventually be adaptive, hence refinement regions can move over the course of the model integration (Blayo and Debreu, 1999). Resolution is in that case increased only where needed, depending on a local numerical or physical criterion, to save computing resources. The use of AMR techniques in realistic ocean models is nevertheless still poorly documented.

The problem of vertical discretization is connected to physical processes that the modeler wants to resolve and it must address questions related to: a) the representation of pressure gradients; b) the representation of sub-grid scale processes; c) the need to concentrate the resolution in a specific region (e.g. the shelf, the coastal areas, etc.); and d) the comparison with observations. Griffies et al. (2000) distinguished among three traditional approaches (Figure 5.5):

• Depth/geopotential vertical coordinates; 
• Terrain-following;
• Potential density (isopycnic) vertical coordinates. 

Geopotential (z-) coordinates (Figure 5.5A) have been largely used in ocean and atmospheric models because of their simplicity and straightforward nature for parameterizing the surface boundary layer. On the contrary, they are not able to adequately represent the effect of topography on the large-scale ocean models. Terrain-following coordinate systems (Figure 5.5B) are used especially in coastal applications, where bottom boundary layers and topography need to be well resolved. As z-coordinates, they suffer from spurious diapycnal mixing due to problems with numerical advection. In isopycnic vertical coordinates (Figure 5.5C), the potential density is referred to a given pressure. This system basically divides the water column into distinct homogeneous layers, which thicknesses can vary from place to place and from one time step to the next. This choice of coordinate works well for modelling tracer transport, which tends to be along surfaces of constant density. While both layered and isopycnal models use density as the vertical coordinate, there are subtle differences between the two types. Griffies et. al. (2000) and Chassignet et al. (2006), provide a discussion on the advantages and disadvantages of each vertical coordinate system.

Once the model is set from the spatial point of view and discretization in horizontal and vertical is defined, the time step for the computation needs to be considered as well. In the numerical schemes used to integrate the primitive equations, the time step must be small enough to guarantee computational stability. The Courant-Friedrichs-Lewy criterion (CFL) is the stability condition that states that the velocity c at which the information is propagating at times the time step ∆t must be less than the horizontal grid spacing ∆x:

Where C is the Courant number and Cmax depends on the specific used scheme: explicit schemes allow to advance the solution to the next time level, one spatial grid point at a time, and are quite simple to implement (Kantha and Clayson, 2000); in an implicit time-stepping scheme, the solution at the next time level must be derived for all grid points simultaneously. These schemes are computationally more intensive, but are unconditionally stable, thus permitting larger time steps to be taken than would otherwise be required.

Three families of methods are available for discretizing the space derivatives that enters in the primitive equations:

• Finite difference Method (FDM);
• Finite Volume Method (FVM);
• Finite Elements Method (FEM). 

Here we provide an introduction to each method but for more detailed explanation refer to Hirsch, 2007. 

The FDM is based on the properties of the Taylor expansions: it corresponds to an estimation of a derivative by the ratio of two differences according to the theoretical definition of the derivative, like the following:

If we remove the limit in Eq. 5.9, we obtain a finite difference: additionally, if ∆x is “small” but finite, the expression on the RHS of Eq. 5.9 is an approximation of the exact value of ux. Since ∆x is finite, an error is introduced, called truncation error, which goes to zero for ∆x tending to zero. The power of ∆x with which this error tends to zero, is caller order of accuracy of the difference approximation and can be obtained by a Taylor series of u(x+∆x) around point x (Eqq. 5.10 and 5.11):

Equation 5.11 shows that:

• The RHS of Eq. 5.9 is an approximation of the first derivative u_x in the point x;
• The remaining terms in the RHS represent the error associated with this formula. 

If we restrict the truncation error to its dominant term, that is the lower power of ∆x, we see that this approximation for u(x) goes to zero like the first power of ∆x and is said to be the first order in ∆x:

Where O(∆x) is the truncation error. 

The FVM is a numerical technique by which the integral formulation of the conservation laws is discretized directly in the physical space. It is based on cell-averaged values, which makes this method totally different from FDM and FEM where the main numerical quantities are the local function values at the mesh points. For each cell, a local finite volume, also called control volume, is associated to each mesh point and applies the integral conservation law to this local volume. For this reason, the FVM is considered a conservative method. The essential property of this formulation is the presence of the surface integral and the fact that the time variation of a generic variable u inside the volume only depends on the surface values of the fluxes.

The FVM requires:

• The subdivision of the mesh, obtained from the space discretization, into finite small volumes, one control volume being associated to each mesh point;
• The application of the integral conservation law to each of these finite volumes. 

The FEM originates from the field of structural analysis and it has two common points with the FVM:

• The space discretization is considered a set of volumes or cells, called elements;
• It requires an integral formulation as a starting point that can be considered as a generalisation of the FVM. 

The FEM requires:

• Discretization of the spatial domain into a set of elements of arbitrary shapes;
• In each element, a parametric representation of the unknown variables, based on families for interpolating or shape functions, associated to each element or cell is defined.

Such nice properties of the FEM as conservation of energy, that is common for all variational methods of solving differential equations, treatment of boundary conditions, and flexibility of irregular meshes have made them quite attractive, since they are also well suited to parallel computing. For this reason, it is considered as an interesting alternative to FDM commonly used in ocean modelling (Danilov et al., 2004).

In Table 5.1, are summarised some of the most used ocean models that integrate numerically the primitive equations for a wide range of spatial domains, from global ocean to coastal scales.

Table 5.1. List of available ocean models used from global to coastal scales.

The need to resolve the small scales of ocean circulation in coastal seas, aswell as the impracticability to run models at sufficiently high resolution and detailed physics at global scales, led to the development of downscaling approaches for both the direct modelling and the data assimilation problems.

Two families of modelling approaches can be distinguished: (1) models running at global scales with mesh refinement in the coastal areas of interest; and (2) one-way or two-way nesting of coastal models into regional or global ones. In practice, the first one is achieved by setting variable-mesh grids, such as unstructured or curvilinear structured grids (as discussed in 5.4.2.1). To our knowledge, only 2D (i.e. barotropic) unstructured models dedicated to storm surges and/ or tides modelling, such as the tidal atlas FES2014 (Lyard et al., 2021), are running over the global ocean and satisfy the resolution requirements in shallow waters. In the second approach, the large-scale global (or regional) model, i.e. the ‘parent’ model, provides open-boundary conditions to the coastal (‘child’) model; in case of two-way nesting, both models are coupled and the child model returns an estimate of the ocean state at its boundary, which is used in turn to force the parent simulation. General resolution issues for both approaches and practical considerations are discussed in Greenberg et al. (2007).

However, nesting methods do not just consist in reproducing the large-scale solution with more details. Indeed, the child model may represent different processes from those solved by the parent model (e.g. tides, surface gravity waves, etc.) or may rely on different parameterizations or parameters. Besides, due to the strong nonlinearity of the ocean flow, the internal variability of the child model may decouple from that of the parent, leading to divergent solutions (Katavouta and Thompson, 2016). Figure 5.6 shows an example of spectral nudging in the Gulf of Maine as in Katavouta and Thompson (2016). The spatial domain is given in Figure 5.6-top: the black box represents the bounding box of the regional model GoMSS (NEMO, 1/36° horizontal resolution), which is nested into the HYCOM+NCODA global 1/12° analysis system. GoMSS+ is the regional configuration with spectral nudging where temperature and salinity variables are directly updated. By adopting such a nesting approach, the regional configuration significantly improves the quality of the solution as shown in Figure 5.6-bottom: it represents the sea surface temperature snapshots for 22 July 2012 based on satellite (“Obs”), the global system (“Global system”), the regional system (“GoMSS”), and that implementing the spectral nudging (“GoMSS+”). The GoMSS+ exhibits an improved version of the coastal sea surface temperature representation, which is typical for a higher resolution model that takes into account coastal processes (e.g., tides). At the same time, it is able to capture the warm slope water and cold shelf waters as shown in the observations, which are well represented in the global model thanks to data assimilation. For further details, please refer to Katavouta and Thompson (2016).

In 2007, the GODAE Coastal and Shelf seas Working Group (De Mey et al., 2007) noted that: “It is becoming increasingly clear that specifying the offshore boundary conditions of coastal models by using forecasts from a hydrodynamical large-scale ocean model has the potential (1) to provide better local estimates by adding value to GODAE products, (2) to extend predictability on shelves, and (3) to enhance the representativeness of local observations.” Despite considerable efforts since 2007 on both coastal modelling capabilities and nesting methods, downscaling still raises obvious numerical and physical issues. In the following paragraphs an attempt has been made, but not exhaustively, to present the various difficulties that arise and the solutions found in the literature to address them. The coastal ocean is subject to both local (e.g. atmosphere, river mouths) and remote forcings (e.g. astronomical potential, coastal waveguide, wind fetch, biogeochemical connectivity). Therefore, the boundaries of a coastal model, which also intercept strong bathymetry gradients, play a critical role. In addition, solving primitive equations on a limited area domain with OBC does not lead to a unique physically realistic solution. Consequently, a variety of ad hoc methods to set-up practical OBC have been developed with a dependence upon flow dynamics, model resolution, types of information at the open boundaries, etc., as reviewed by Blayo and Debreu (2005). A simple view of the OBC issues consists in viewing the problem because of inconsistencies between the parent and child models which, as mentioned previously, arise due to different physics of the model, to different forcing (e.g. atmospheric, runoff, bathymetry), and to truncated information at the open boundary. The last refers to the fact that the parent information is provided as discrete fields in space and time (e.g. daily or hourly averages); high-frequency motions are therefore filtered out or aliased.

The example of tides is particularly enlightening on these limitations. Even though the parent model resolves tides, forcing the child with the parent tidal waves (either barotropic or both barotropic and baroclinic) implies the availability of the large-scale forcing at very high frequency (a few minutes). In practice, especially for operational systems, this is very difficult to achieve as it requires huge storage capacities. Therefore, coastal models are usually forced by low-frequency dynamics and tidal constituents, both of which not necessarily stemming from the same parent models (tidal constituents are often chosen from accurate global tidal atlases). Herzfeld and Gillibrand (2015) noted that conditions for incoming tidal waves may be reflective for the low-frequency external data and propose OBC based on dual relaxation time scales. Furthermore, the difference of bathymetry and representation of the coastline between the parent and child models may lead to large inconsistencies between the tidal solutions in both models, with a risk of spurious patterns developing in the coastal domain close to the open boundaries (e.g. rim currents). Toublanc et al. (2018) proposed a simple approach that reduces such inconsistencies by pre-processing the tidal forcing using a 2D simulation with a dedicated 2D tidal model. At last, filtering out the high-frequency 3D incoming information by using for instance hourly or daily averages from the parent simulation, may lead to a loss of energy in the coastal domain, in particular because of the missing internal waves forcing, as recently evidenced by Mazloff et al. (2020).

Another difficulty in one-way nesting arises from the possibility that the child model develops an internal variability that diverges from the parent’s one. In many operational systems, global or large-scale solutions stem from a data assimilation system in which the mesoscale dynamics are constrained by satellite data (e.g. altimetry). If no data assimilation is performed in the coastal domain, the developing mesoscale (and a fortiori submesoscale) may deviate from reality leading to the undesirable case in which the parent solution is closer to observations at large-scale and mesoscale than the child. Sandery and Sakov (2017) report that even with data assimilation, increasing the resolution does not automatically improve the skill of the forecast, because of the inverse cascade of unconstrained submesoscale towards mesoscale. Methods such as spectral nudging are developed to ensure that the large-scale patterns, e.g. eddies or meandering jets that are accurately represented in the parent model, are maintained in the child; an example of such method can be found in Katavouta and Thompson (2016).

A last but not least issue concerns quantifying the errors in the child simulations due to the nesting process. The errors originate from the OBC scheme (numerical implementation and physical assumptions) and from the uncertainties on the parent forcing fields. In the latter case, the question is how the parent model errors are downscaled. Ensemble approaches can help to characterise and estimate the downscaling of parent errors, as for instance explored in Ghantous et al. (2020). 

Figure 5.7 shows an example of ensemble downscaling of a coastal ocean model (Symphonie model, 500 m resolution) for the south-east Bay of Biscay in an ensemble of a regional model (NEMO, 1/36°) (Ghantous et al., 2020). Figure 5.7a presents the regional domain, in particular the parent domain over the map, while the blue box is the domain of the child model. Figures 5.7b-d show the ensemble spread (standard-deviation) in sea surface height (SSH) in the domain of the child model for ensembles of 50 members. In particular, Figure 5.7b is the parent ensemble, generated by perturbing the wind in the parent domain; Figure 5.7c is the child ensemble, generated by perturbing the wind in the child domain; Figure 5.7d is the child ensemble generated by perturbing both the wind and the OBC conditions (the OBC perturbations stem from the parent ensemble). The numerical experiment reveals that, on average over the period of study, the spread in SSH is greatest where the mesoscale eddies are present (in the deeper area of the domain). It also reveals that the contribution from the OBC uncertainties is larger than the impact of local wind uncertainties. It is a valuable result for the next generation of ensemble data assimilation systems.

An example of nesting capacities of circulation modelling in short-term forecast is shown in Figure 5.8. This is the result of downscaling the Copernicus Marine Service Iberia-Biscay-Ireland – Monitoring and Forecasting Centre (IBI-MFC, 🔗1) product on a higher spatial grid; in the bottom panel it can be seen a detail of surface currents in the Gulf of Cadiz and Alborán Sea.

The downscaling approach is extremely powerful to allow the modeller to set up an OOFS at high resolution, and every OOFS may be used to build another OOFS in a seamless way. In Section 5.9 can be found an initial but exhaustive list of OOFSs’ providers from which the modeller may select to nest her/his new OOFS.