Development and technical paper 18 Sep 2020
Development and technical paper  18 Sep 2020
Development of a semiLagrangian advection scheme for the NEMO ocean model (3.1)
 Environment and Climate Change Canada, 2121 route Transcanadienne ouest, Dorval, Québec, H9P 1J3, Canada
 Environment and Climate Change Canada, 2121 route Transcanadienne ouest, Dorval, Québec, H9P 1J3, Canada
Correspondence: Christopher Subich (christopher.subich@canada.ca)
Hide author detailsCorrespondence: Christopher Subich (christopher.subich@canada.ca)
As resolutions of ocean circulation models increase, the advective Courant number – the ratio between the distance travelled by a fluid parcel in one time step and the grid size – becomes the most stringent factor limiting model time steps. Some atmospheric models have escaped this limit by using an implicit or semiimplicit semiLagrangian formulation of advection, which calculates materially conserved fluid properties along trajectories which follow the fluid motion and end at prescribed grid points. Unfortunately, this formulation is not straightforward in ocean contexts, where the irregular, interior boundaries imposed by the shore and bottom orography are incompatible with traditional trajectory calculations.
This work describes the adaptation of the semiLagrangian method as an advection module for an operational ocean model. We solve the difficulties of the ocean's internal boundaries by calculating parcel trajectories using a timeexponential formulation, which ensures that all parcel trajectories remain inside the ocean domain despite strong accelerations near the boundary. Additionally, we derive this method in a way that is compatible with the leapfrog timestepping scheme used in the NEMOOPA (Nucleus for European Modelling of the Ocean, Océan Parallélisé) ocean model, and we present simulation results for a simplified test case of flow past a model island and for 10year free runs of the global ocean on the quarterdegree ORCA025 grid.
Recent work by Smith et al. (2018) has shown that over the medium term (up to 7 d), a coupled forecasting system involving ocean, ice, and atmospheric models can significantly improve forecasting skill over forecasts that extend initial ocean and ice conditions over the atmospheric forecast period. While this is an exciting development for the future of numerical weather prediction, coupling adds a new dimension to the computational cost. Developing a deployable forecast system, especially with regional or ensemble components, requires exploiting every reasonable opportunity for optimization. One straightforward optimization is to maximize the admissible time step of the ocean component, and we intend to improve the ocean time step limit by implementing a semiLagrangian advection module into the popular NEMOOPA (Nucleus for European Modelling of the Ocean, Océan Parallélisé; Madec, 2008, version 3.1) model, used in this coupled system. This module is intended as a dropin replacement for the model's other advection modules, and in particular it does not interfere with NEMO's timestepping algorithm (leapfrog).
1.1 Time step constraints in the ocean
A numerical model with an explicit timemarching scheme must generally limit its time step to satisfy a Courant–Friedrichs–Lewy (CFL) condition: information must not propagate more than a discretizationdefined maximum number of cells in a single step, leading to a maximum stable Courant number. For systems such as the Euler equations (for the atmosphere) or hydrostatic equations (as implemented by NEMOOPA), the information propagation speeds are controlled by the admissible wave modes of the systems, which become characteristic curves.
In the atmosphere, the most restrictive wave mode is that corresponding to sound waves. These waves are fast compared to atmospheric motions, and in response, atmospheric models generally treat sound waves either implicitly or through subcycling, especially in the most restrictive vertical direction. The secondmost stringent restriction comes from simple advection by winds in the upper atmosphere. At the Canadian Meteorological Centre (CMC), the atmospheric forecasting system (and atmospheric component of the coupled forecasting system) uses the GEM (Global Environmental Multiscale; Girard et al., 2014) model, which addresses this time step restriction through a semiLagrangian treatment of advection (Robert, 1982).
In the ocean, the Boussinesq assumption eliminates sound waves, but the model is left with the problem of surface gravity waves. Here, NEMOOPA takes a similar approach to that used by atmospheric models for sound waves, by either treating the surface pressure gradient in a timeimplicit manner (with a linearized free surface, used in this work) or by subcycling. The ocean lacks any direct equivalent to the atmosphere's strong upperair winds, and so advection by the background velocity and internal gravity wave modes compete as the next most limiting factor for the maximum stable time step. Lemarié et al. (2015) finds that the Courant number associated with vertical advection is more limiting than that associated with internal (baroclinic) gravity waves at resolutions of $\mathrm{1}/\mathrm{2}{}^{\circ}$, and the Courant number associated with horizontal advection catches up with that of gravity waves at resolutions of $\mathrm{1}/\mathrm{4}{}^{\circ}$ and finer.
1.1.1 Grid stretching
In order to cover the entire ocean in a single, continuous domain, global NEMOOPA model configurations typically use grids based on the ORCA “tripolar” grid (Madec and Imbard, 1996; Murray, 1996). This grid is defined in the Northern Hemisphere by an elliptical coordinate system, where the latitudelike coordinate is defined by ellipses with a shared pair of foci and the longitudelike coordinate is defined by the hyperbolas orthogonal to these ellipses. These coordinates match continuously at the Equator to lines of latitude and longitude in a Mercator projection. By placing the foci of the ellipses on land, the grid contains no singularities in the ocean domain.
Unfortunately, this placement causes an abundance of small grid cells in the north polar region, especially in the Canadian Arctic Archipelago. Figure 1 depicts this situation at a nominal $\mathrm{1}/\mathrm{4}{}^{\circ}$ resolution: the grid point spacing of 25–30 km near the Equator falls to 3–4 km in the archipelago. The areas in Fig. 1 with the narrowest grid spacing are also shallow seas, with depths of 200 m or less and nontidal currents of 15–30 cm s^{−1}. This grid stretching is of particular concern when adapted to regional models such as that of Lemieux et al. (2016), which refine this grid while retaining its tripolar structure to use conforming boundary conditions.
The coordinate system is also stretched in the vertical direction. Using the zlevel grid option of the NEMOOPA model, layers near the surface are spaced much more closely together than layers nearer the ocean bottom, in order to provide adequate resolution of the mixing layer. This stretching enhances the impact of vertical advection on the vertical Courant number, even if verticalcurrent magnitudes are low in absolute terms; Lemarié et al. (2015) notes that vertical advection provides a tighter bound on the time step than horizontal advection.
SemiLagrangian advection alleviates both vertical and horizontal Courant number restrictions by tracing fluid parcels in a Lagrangian, fluidfollowing coordinate system. This coordinate system is defined so that the end of each time step the fluid parcels arrive on the prescribed computational grid, and the properties of the fluid parcels (in the ocean setting, temperature, salinity, and horizontal velocity) at the end of the time step (on the computational grid) are found by interpolating the previousstep, gridded values to the origin point of each parcel's trajectory. This method provides an implicit treatment of advection, allowing time steps with advective Courant numbers greater than those usually permitted by explicit, Eulerianform models.
In this work, we describe the initial implementation of a semiLagrangian advection routine in NEMOOPA, based on the configuration of Smith et al. (2018). This configuration uses a linear free surface where the vertical coordinate does not move in time, but we believe that the described method can be generalized.
1.2 Existing work
In the atmosphere, semiLagrangian advection is a standard technique (Robert, 1982) for the implicit treatment of advection, but especially at large scales the effects of topography are relatively gentle. In particular, trajectory calculations can proceed under the assumption that the fluid parcel does not experience strong boundaryrelated acceleration. In the ocean domain this assumption is strongly violated, particularly for zlevel vertical grids where the bathymetry changes abruptly at lateral cell boundaries.
Some attempts have been made previously to incorporate semiLagrangian advection into the ocean context. The work of Casulli and Cheng (1992), which is used as part of the ELCOM lake and estuary model (Hodges and Dallimore, 2006), calculates parcel trajectories via a substepping approach, where fluid parcel trajectories are integrated via an explicit Euler method over many short steps per model time step. The twodimensional, unstructured shallow water model of Walters et al. (2007) takes a similar approach, where it also must take at least one substep per element boundary traversed by a fluid parcel.
In this work, we overcome this difficulty with an iterative trajectory calculation which reduces in the limit to an implicit trapezoidal rule. In addition, we also derive the semiLagrangian advection scheme in a form which calculates effective advective tendencies, such that the advection routine can operate as a “plugin” scheme for models which traditionally use Eulerian fluxes. We apply this to the NEMOOPA model, and we believe this algorithm may be useful when applied to other ocean models with a structured grid.
In exchange, however, the semiLagrangian advection formulation departs from NEMO's finitevolume interpretation of its tracer and velocity components. By tracing infinitesimal fluid parcels, semiLagrangian advection treats gridpoint values analogously to a finitedifference method, and as a consequence the scheme does not naturally offer conservation guarantees. This is not a primary concern for the short to mediumterm forecasting applications that form the direct target for this work, but extensions of the semiLagrangian scheme to ensure conservation (Lauritzen, 2007) may be needed before the technique is applicable to longerterm climate simulations.
Additionally, Leclair and Madec (2011) has developed an “arbitrary Lagrangian–Eulerian” vertical coordinate scheme, implemented in recent versions of NEMO. This scheme splits vertical motions into fast (high temporal frequency) and slow motions, and the former are treated by comoving vertical coordinate surfaces with a regridding step. This coordinate system reduces spurious diapycnal mixing caused by the highfrequency vertical motions, and its Lagrangian treatment of these motions relaxes the corresponding stability restriction.
1.3 Organization
We first introduce the time discretization of the semiLagrangian scheme in Sect. 2, in order to develop a formulation that remains compatible with the common leapfrog scheme. In Sect. 3, we begin to spatially discretize the semiLagrangian scheme by specifying the horizontal and vertical interpolation operators, and in Sect. 4 we complete the discretization by defining the trajectory calculations. We present preliminary numerical examples in Sect. 5, demonstrating the stability of the advection scheme.
The first requirement of a semiLagrangian advection scheme for the NEMOOPA model is that it be consistent with the model's overall timestepping approach: the advection scheme is but one component of the full model.
In version 3 of NEMOOPA, nondiffusive, nondamping processes such as advection are implemented via the leapfrog scheme (Mesinger and Arakawa, 1976), where at each time step a field f receives its new value at f^{A} (f “after” the time step) based on its value at the previous time step and forcing terms, which are all evaluated on the reference grid x_{ref}. This gives a schematic of
where f^{A} is the field calculated at time t_{0}+Δt, f^{B} is the field evaluated at the known prior time t_{0}−Δt (“before”), f^{N} is the field at the provided time t_{0} (“now”), and F is the forcing operator. The forcing operator includes advective processes at the now time level, but diffusive, damping, and hydrostatic pressure terms might be evaluated at either the before or after time levels.
This is an Eulerian approach to fluid motion, where tracer and momentum values are tracked along the fixed reference grid at all times, and fluid flows through this grid.
2.1 SemiLagrangian advection
In contrast, the Lagrangian advection schemes consider the fluid parcel to be the fundamental unit of discretization. In this perspective, if f is a property of a fluid parcel that is conserved along a trajectory^{1}, it satisfies the continuous equations:
where $\frac{\mathrm{D}}{\mathrm{D}t}={\partial}_{t}+\mathit{u}\cdot \mathrm{\nabla}$ is the material derivative, and F_{L} (Lagrangian righthand side) contains all the same forcing terms as F except those arising from tracer and momentum flux, which are included inside the material derivative.
Ordinarily, Eq. (2) is discretized so that F_{L} is evaluated following the Lagrangian particles in the moving coordinate frame x(t), satisfying the trajectory equation:
From an Eulerian point of view, Eq. (3) is a trivial identity based on the definition of the material derivative, but from the Lagrangian point of view, Eq. (3) must be solved to define x over time.
One technique for solving Eqs. (2) and (3) is the twotimelevel implicit semiLagrangian method, used in the GEM atmospheric model (Girard et al., 2014), among others. Here, the F_{L} terms are evaluated with a trapezoidal rule, discretizing Eqs. (2) and (3) as
The trajectory Eq. (4b) acts to implicitly define the paths of the traced fluid parcels, where each location on x_{ref} is associated with a corresponding departurepoint location x^{D}. Over the single time step, fluid parcels depart from x^{D} (which in general is not aligned with the grid) and arrive on the reference grid.
This offgrid, departure point evaluation of u and F_{L} is fundamental to Lagrangian and semiLagrangian methods, and f^{N}(x^{D}) (${\mathbf{F}}_{L}^{N}\left({\mathit{x}}^{\mathrm{D}}\right)$) can be written more simply as f^{D} (${\mathbf{F}}_{L}^{\mathrm{D}}$) for “departurepoint f (F)”. Neither the timeimplicit evaluations (generally) nor the offgrid evaluations (of nonadvective forcing) are compatible with the core structure of NEMOOPA, which considers advection to be just one of many independent operators influencing the F term of Eq. (1).
2.2 Reconciliation
Implementing semiLagrangian advection in NEMOOPA requires adopting as much of the framework of Eq. (1) as possible, without changing the evaluation of nonadvective forcing terms. Effectively, the semiLagrangian advection routine must ultimately supply a time trend that, from the perspective of the leapfrog timestep algorithm, is indistinguishable from a conventional fluxform advection operator.
To effect this, consider Eq. (2) without forcing terms (F_{L}=0). The function f is preserved following the flow, so this gives the simply written
This is approximated by taking one time step of Eq. (1) (with only advective forcing F_{adv}), but the latter involves integrating over the whole interval from t_{0}−Δt to t_{0}+Δt. Thus, we should identify f^{D} (and the departure points generally) not with the now time level in the leapfrog scheme but rather with the before time level. Doing so and equating Eqs. (5) and (1) gives
Equation (7) is prescriptive, and it gives the necessary trend for the leapfrog algorithm. Evaluating it requires f only at the already known, before time level and calculation of the departure points x^{D}. This calculation is further simplified by basing the departure points on the timecentered velocities u^{N}, and the exact algorithm for this calculation will be discussed in more detail in Sect. 4.
2.3 Effects of the Asselin filter
To prevent decoupling of odd and even time steps (damping the computational mode), NEMOOPA is typically configured to use the Asselin time filter (Asselin, 1972), which adds a small time damping proportional to $\frac{{\partial}^{\mathrm{2}}}{\partial {t}^{\mathrm{2}}}f$. Using the notation of Shchepetkin and McWilliams (2005) adapted to Eq. (1), the filter extends the timemarching scheme to the sequence:
Equation (8a) is the direct equivalent of Eq. (1), creating a provisional after value f^{A*}. Equation (8b) applies the filter (with a strength parameter ϵ) with this value and the previous step's provisional field to define a final now field, and finally Eqs. (8c) and (8d) are “bookkeeping” steps to shift field labels to become ready for the next time step. The forcing operator F^{N*} is evaluated based on the provisionally defined fields.
In applying this filter with the semiLagrangian forcing, Eq. (7) is oblivious to the presence of the filter or the difference between f^{N} and f^{N*}. Substituting Eq. (7) into Eq. (8a) and applying Eqs. (8c) and (8d) to Eqs. (8a) and (8b) gives the update equation:
In the case of onedimensional advection by a constant velocity u_{0}, the trajectory calculation is trivial, and
Since Eq. (9) is linear, we can also take its Fourier decomposition in space and consider only a single, arbitrary wave mode, giving $f=\widehat{f}\left(t\right)\mathrm{exp}\mathbf{i}kx$ for a timevarying coefficient $\widehat{f}$. Applying this to Eq. (9) casts the update in a matrix form as
The time stability of this filter is then governed by the eigenvalues of this matrix. Using the shorthand $\mathit{\omega}=k{u}_{\mathrm{0}}\mathrm{\Delta}t$, these eigenvalues are
and to leading order in ϵ, these eigenvalues have squared magnitudes of
signifying stability ($\left\mathit{\lambda}\right\le \mathrm{1}$) for all values of ω and thus all Courant numbers.
To perform the offgrid interpolations in Eq. (7) to find f^{B}(x^{D}), this method fits a cubic polynomial to the underlying function. If the single departure point ${\mathit{x}}_{\mathrm{d}}=({x}_{\mathrm{d}},{y}_{\mathrm{d}},{z}_{\mathrm{d}})$ lies within^{2} $({x}_{a},{x}_{a+\mathrm{1}})\times ({y}_{b},{y}_{b+\mathrm{1}})\times ({z}_{c},{z}_{c+\mathrm{1}})$ for integer values of a, b, and c coinciding with gridpoint locations, the full interpolation stencil consists of the gridindex cube $i\in [a\mathrm{1},a+\mathrm{2}]$, $j\in [b\mathrm{1},b+\mathrm{2}]$, and $k\in [c\mathrm{1},c+\mathrm{2}]$.
This grid cube contains up to 64 grid points where f(x) might be defined (subject to boundary conditions), and building a complete interpolation stencil would be cumbersome and inefficient. Instead, the interpolation procedure takes advantage of the tensorproduct nature of the grid to separate interpolation along each dimension:
To effect the onedimensional interpolations in Algorithm 1, we make use of the cubic Hermite polynomials (Hildebrand, 1974). On the interval $\mathrm{0}\le \mathit{\chi}\le \mathrm{1}$, these polynomials are
and a function f(χ) defined on this interval is interpolated via
Here, we prefer to use the cubic Hermite polynomials over simple Lagrange polynomial interpolation because the former choice allows greater freedom (via Eq. 15) in implementation. If f^{′} is approximated by a fourpoint finitedifference stencil, then Eq. (15) reduces to Lagrange interpolation. However, we can also make other choices for f^{′} to impose desirable properties: restricting f^{′} to have the same sign as the discrete difference imposes a type of slope limiting, and calculating f^{′} through a threepoint stencil provides for continuous derivatives. These approaches are discussed in more detail in the following sections.
Interpolation using the above algorithm involves appropriately defining the interval to be scaled to [0,1] and approximating f^{′} at the endpoints. Because of the high aspect ratio of oceanic flows and the special character of vertical motion in a stratified ocean, these approximations differ between the horizontal and vertical interpolations.
3.1 Horizontal interpolation
In the horizontal, the interpolation in Eq. (15) can be directly conducted in
gridindex space. Even when the underlying grid is mapped to the sphere, such as in the ORCA global
grid (Madec, 2008, chap. 16), the grid generally transitions smoothly and slowly from point to
point^{3}. The physical trajectory departure location x_{d} (and y_{d})
can be translated into a fractional grid index offset by dividing by the appropriate grid scale factor, available inside the NEMOOPA source code as one of e[12][tuv]
.
Achieving thirdorder accuracy inside Eq. (15) is possible, but doing so requires an equally accurate estimate for f^{′}. Unfortunately, interpolating successively in one dimension using the above algorithm does not allow for precomputation of these derivatives: after the vertical interpolation step, all of the function values need to be taken off grid, so any precomputed derivatives would themselves require interpolation. Instead, sufficiently accurate estimates of the derivative are available by applying a finitedifference formula to the function values themselves.
3.1.1 Derivative estimates
For notational simplicity, begin with the last step of the above algorithm where we have $f({x}_{\mathrm{d}},{y}_{j},{z}_{\mathrm{d}})$ and would like to estimate $f({x}_{\mathrm{d}},{y}_{\mathrm{d}},{z}_{\mathrm{d}})$. If y_{d} lies between y_{0} and y_{1}, then the fourpoint interpolation stencil implies that we have computed $f({x}_{\mathrm{d}},{y}_{j},{z}_{\mathrm{d}})$ for $j=\mathrm{1},\mathrm{0},\mathrm{1},\mathrm{2}$. To emphasize that this is now a onedimensional interpolation problem, let $g\left(j\right)=f({x}_{\mathrm{d}},{y}_{j},{z}_{\mathrm{d}})$, such that $f({x}_{\mathrm{d}},{y}_{\mathrm{d}},{z}_{\mathrm{d}})=g\left({j}^{\prime}\right)$ for some ${j}^{\prime}\in [\mathrm{0},\mathrm{1}]$. In this domain, g^{′}(0) and g^{′}(1) can be approximated by the finite differences:
which then substitute for the appropriate derivatives in Eq. (15).
These finite differences are exact expressions for the first derivative for polynomials up to third order in j, and their use essentially converts Eq. (15) to interpolation via Lagrange polynomials. The Hermite polynomial form, however, allows for an easier imposition of boundary conditions.
3.1.2 Boundary conditions
On the NEMOOPA zlevel grid, the lateral boundaries coincide with u and v points (velocity points), which are spaced halfway between t points (tracer points). Tracer points that lie inside the land region are masked (𝚝𝚖𝚊𝚜𝚔=0) as are velocity points that are at the edge of or within the masked region. This arrangement is illustrated for a sample region in Fig. 2.
The physical interpretation of the boundary varies with respect to the field being interpolated. For tracers, lateral boundaries imply noflux conditions for the purposes of advection, which in turn implies a zero derivative at the boundary. The normal velocity (u with respect to a boundary along the first grid dimension, v with respect to a grid boundary along the second) is obviously constrained to zero by geometry to give a Dirichlet boundary condition, whereas the tangential velocity can be set as freeslip, noslip, or some combination via a namelist entry. In the subsequent, we assume that velocity has a freeslip boundary condition, with boundary friction left for future work.
If a boundary occurs in the left portion of the interpolation stencil, there are a total of seven cases (see Algorithm 2).
For boundaries that occur in the right portion of the interpolation stencil, the values taken for ghost points are given symmetrically.
The combination of “the grid point is at the boundary” and “the boundary is of Neumanntype” is missing from Algorithm 2. This construction is forbidden by the grid structure of NEMOOPA, where tangential velocity is located one halfcell away from a boundary.
For twodimensional interpolation, Algorithm 2 applies independently to each dimension. When interpolating along x, the points $f({x}_{\mathrm{d}},{y}_{j},{z}_{\mathrm{d}})$ will each individually be either in the ocean domain and valid or in the land domain and masked, which provides the values necessary to compute $f({x}_{\mathrm{d}},{y}_{\mathrm{d}},{z}_{\mathrm{d}})$. This offgrid point itself must lie in water, which imposes a strong requirement on the trajectory calculations to be discussed in Sect. 4.
3.1.3 Slope limiting
As a final step, once values for the function and its derivative at the interval endpoints are specified, the derivative values are limited to help prevent new maxima in the interpolated function. In particular, if g(0) is a local minimum (maximum) among itself, g(−1), and g(1), then g^{′}(0) is set to zero if the above procedure finds that it would be negative (positive). A similar procedure applies symmetrically for g^{′}(1) if g(1) is a relative extremum.
This limiting is milder than methods derived from Bermejo and Staniforth (1992), which would strictly preserve positivity for any j^{′}, but it effectively limits excursions when j^{′} is close to 0 or 1. Without such limiting, numerical testing showed that semiLagrangian advection of temperature and salinity could cause weak instabilities near the coastline, where a locally extreme temperature or salinity could become “trapped” near the coast and slowly amplified.
3.2 Vertical interpolation
Vertical motion in the NEMOOPA model differs from horizontal motion in a number of respects:

Vertical gradients of temperature and density are much stronger than typical horizontal gradients, especially near the surface.

Typical vertical grids used with NEMOOPA are strongly stretched, with a higher resolution near the surface and a lower resolution in the deep ocean.

Vertical flow is often oscillatory, where vertical motion is driven by barotropic and baroclinic waves.
The horizontal interpolation described in Sect. 3.1 is thirdorder accurate; with the provided onesided formulas for calculating the endpoint derivatives it reduces to a fourpoint (cubic) Lagrangian interpolation process. However, the smooth field implied by this interpolation process is only C^{0} continuous: f(x_{j}−ϵ) “sees” ${f}_{j}^{\prime}$ calculated from f(x_{j−2}) to f(x_{j+1}), whereas f(x_{j}+ϵ) sees ${f}_{j}^{\prime}$ from f(x_{j−1}) to f(x_{j+2}).
We do not find this to be a practical concern for horizontal interpolation, since horizontal currents in most of the ocean tend to be dominated by relatively steady quasigeostrophic motions. In the vertical, however, we found that even lowamplitude oscillations caused by highfrequency gravity waves would cause the temperature and salinity fields to drift. The mechanism is that a fluid parcel displaced upwards by ϵ in one time step and downwards by ϵ in the next time step would see an effective diffusion proportional to the jump between the upward and downwardlooking vertical derivatives.
To maintain global accuracy, we impose C^{1} continuity in the vertical direction through an alternative treatment of the vertical derivative. Instead of applying Eqs. (16), we treat the physical depth (rather than grid index) as the relevant coordinate and construct a centered estimate of the derivative.
For a function f(z_{n}) defined at the z_{n} levels, define $\mathrm{\Delta}{f}^{+}=f\left({z}_{n+\mathrm{1}}\right)f\left({z}_{n}\right)$, $\mathrm{\Delta}{f}^{}=f\left({z}_{n}\right)f\left({z}_{n\mathrm{1}}\right)$, $\mathrm{\Delta}{z}^{+}={z}_{n+\mathrm{1}}{z}_{n}$, and $\mathrm{\Delta}{z}^{}={z}_{n}{z}_{n\mathrm{1}}$. These differences combine to give the estimated derivative:
which is accurate to O(Δz^{2}) for the derivative and accurately reproduces quadratic functions of z. In the limiting case of a constant Δz (equispaced vertical levels), this formula reduces to the classic centered difference.
Because vertical interpolation comes first in Algorithm 1, Eq. (17) need be evaluated only at grid points, and in fact it may be precomputed for the entire grid for a given function and time step. This is a key advantage of placing vertical interpolation first in the interpolation sequence, and it avoids duplication of work.
Whereas interpolation near the horizontal boundaries is complicated by the many combinations of grid staggering and physical boundary conditions, interpolation near the vertical boundaries is much simpler. On the NEMO grid, tracers and horizontal velocities lie along the same vertical level, and these levels are staggered one halfcell away from the boundaries. Likewise, the natural vertical boundary condition for both tracer and horizontal velocity fields is a noflux boundary condition; NEMOOPA models boundarylayer friction in another module. Interpolation near the boundaries then proceeds in two steps.
The first step is to define f_{z} at the top and bottom points in the water column, for which the
central difference formula of Eq. (17) is not directly valid. Here, we approximate the
physical noflux condition through a fictitious ghost point such that $\mathrm{\Delta}{f}^{}=\mathrm{0}$ at the top
boundary and $\mathrm{\Delta}{f}^{+}=\mathrm{0}$ at the bottom boundary, with the respective Δz matching the
layer thickness (e3t
).
The second step is to define how Eq. (15) applies to the interval between the grid level and the physical boundary. Here, the noflux boundary conditions reduce to even symmetry, and the derivative at the ghost points is the negative of the vertical derivative calculated for the inboundary point. Near the free surface, if the interpolation point is above the level of the free surface (above z=0), then it is clamped to the surface itself. Near the ocean bottom, if the interpolation point is below the level of the ocean bottom (below $z={z}_{max}$) then the point is masked and is treated as an “inside the boundary” point for the purposes of horizontal interpolation above.
3.2.1 Treatment of partial cells
Over most of the domain, this interpolation works well. Although there is no guarantee of positivity in the derivative formulation of Eq. (17), overshoots and the consequent generation of spurious maxima are limited. For the tests presented in Sect. 5, there was no need for slope limiting for vertical interpolation over most of the domain.
One exception to this rule is at the bottom boundary. Here, vertical levels are spaced far apart, but to better represent the ocean bottom the zlevel grid of NEMOOPA uses a partial cell configuration (Madec, 2008, sec. 5.9). For water columns where the bottommost cell is much deeper than its neighbours, a local (small) upwelling can cause an overshoot of temperature or salinity that spuriously increases the local density but does not diminish the upwelling. Over time, the maximaincreasing trend can accumulate and cause some points at the bottom boundary to reach implausibly cold temperatures (below $\mathrm{10}{}^{\circ}\phantom{\rule{0.125em}{0ex}}\mathrm{C}$, for example) or high salinities. In the absence of explicit horizontal diffusion (which would mix this maximum into more dynamically active regions), these spurious maxima do not generally corrupt the flow, although they obviously would corrupt wholeocean (or wholelevel) statistics such as average or extreme temperatures.
Near these boundary cells, vertical limiting is implemented in the simplest possible way: the interpolation of Eq. (15) is replaced with a constant, such that f(z)=f(z_{k}) over the interval from z_{k} downwards to the physical boundary.
Implementing this limiting over the whole bottom level is possible, but that is far stronger than necessary and leads to erroneous diffusion along gentle slopes. When the bottom layer is composed of partial cells of varying thickness, even interpolation along a horizontal plane (that is, without changing physical depth) requires vertical interpolation in grid space to find that constant level in adjoining columns. Imposing vertical limiting along the whole bottom level effects undesired horizontal diffusion, even though the problem solved by limiting is observed when adjoining cells have large relative thickness variations.
As a compromise between these two errors, we only apply the described limiting to vertical interpolation for cells at the bottom boundary which have a layer thickness greater than 1.75 times that of their “thinnest” neighbour.
This exact threshold is empirical, and other grids might require a retuning of this parameter. Ideally, the grid generation itself would avoid abrupt transitions in celllayer thicknesses, but adding such a restriction would make this advection scheme useless as a dropin replacement for the standard advection routines of NEMOOPA.
3.3 A numerical example
As a simple numerical example, consider the case of a tracer being advected in a rectangular, twodimensional domain by an internal wave and a background current. This tracer satisfies the advection equation
for some prescribed velocity field (u,w).
If this tracer field $\mathit{\sigma}(x,z,t)$ would be a function of z alone, i.e. $\overline{\mathit{\sigma}}\left(z\right)$, if not for the wave motion, then its motion is analytically given by
where η(x,z) is the isopycnal displacement, u_{0} is the xdirected background current (uniform in z), and c is the phase speed of the wave. Following Turkington et al. (1991), a streamfunction defined as
gives velocities
which are exact solutions of Eq. (18) for $\mathit{\sigma}(x,z,t)$.
To give an internal wave that respects noflux conditions at the top and bottom of the domain, we set
where k and m are horizontal and vertical wavenumbers respectively, and A is the wave amplitude. For a domain of size L_{x} in the horizontal (periodic) and L_{z} in the vertical, $k=\mathrm{2}\mathit{\pi}/{L}_{x}$ and $m=\mathit{\pi}/{L}_{z}$ give the lowest internal wave mode, used here.
In dimensional units, we take the model domain to be a channel L_{x}=1 km long and L_{z}=100 m deep with a background current of ${u}_{\mathrm{0}}=\mathrm{1}\phantom{\rule{0.125em}{0ex}}\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{1}}$, and we set $c=N/\sqrt{{k}^{\mathrm{2}}+{m}^{\mathrm{2}}}$ based on a mean buoyancy frequency of $N=\mathrm{0.03}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{1}}$, which corresponds to a 1 % density change from the surface to the bottom of the channel. With a wave amplitude of A=10 m, the maximum waveinduced current is about 10 % of u_{0}, and the phase speed is $c\approx \mathrm{0.94}\phantom{\rule{0.125em}{0ex}}\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{1}}$.
In order to represent the pycnocline found in many ocean waters, we choose^{4} $\overline{\mathit{\sigma}}\left(z\right)=\mathrm{tanh}(\frac{\mathrm{1}}{\mathrm{2}}\frac{\mathrm{1}}{\mathrm{10}}z{L}_{y}^{\mathrm{1}})$. The domain is discretized by N_{x}×N_{y} points, defined as
where $i=\mathrm{1},\mathrm{2},\mathrm{\cdots},{N}_{x}$, $j=\mathrm{1},\mathrm{2},\mathrm{\cdots},{N}_{z}$, and
This implements a stretched vertical coordinate that increases the vertical resolution in the vicinity of the pycnocline.
3.3.1 SemiLagrangian advection
In integrating this system with semiLagrangian advection, the leapfrog method reduces to an Euler method of twice the time step because there is no external forcing. The timediscrete equation is
where (x_{d(ij)},z_{d(ij)}) is the departure point of the trajectory that arrives at the grid point (x_{i},z_{j}), and the offgrid evaluation of σ proceeds via the interpolation processes described earlier without slope limiting.
The departure points are given by the trapezoidal rule^{5} with a timecentered evaluation of velocity:
where the velocities are evaluated exactly via Eq. (21). The overall system Eq. (25) is solved via simple iteration, with an initial guess given by setting $(x,z{)}_{\mathrm{d}\left(ij\right)}=(x,z{)}_{ij}$.
This algorithm is stable for large time steps, so we tested this system for time steps corresponding to Courant numbers of 0.2 and 2.1, with spatial grid resolutions between 40×4 and 2560×256. The final integration time was chosen to be ${t}_{\text{fin}}=\mathrm{5}{L}_{x}/{u}_{\mathrm{0}}$, which allowed the wave to propagate through the domain several times. Since the exact solution is analytically known, we recorded the maximum error experienced over the integration, and error convergence rates are shown in Fig. 3.
3.3.2 Fluxform advection
As a control, we also integrate this system in flux form, i.e. ${\mathit{\sigma}}_{t}\mathrm{\nabla}\cdot \left(\mathit{u}\mathit{\sigma}\right)=\mathrm{0}$, via centered differences, with σ evaluated at the midpoints between grid cells via a simple average, matching the central difference tracer advection scheme in NEMOOPA.
The velocity field given by Eq. (20) is divergence free, so this form of the equation is pointwise equivalent to Eq. (18). However, this no longer holds after discretization. In order to eliminate the divergence error, the velocity field is defined by creating the streamfunction at the staggered points $({x}_{i+\mathrm{1}/\mathrm{2}},{z}_{j+\mathrm{1}/\mathrm{2}})$ and defining discrete velocities u and w via the discrete equivalents to Eq. (21). With this modification, the discrete fluxform operator is equivalent to a discrete advection equation.
After leapfrog discretization in time, the discretized equation is
where $\mathrm{\Delta}{z}_{j}={z}_{j+\mathrm{1}/\mathrm{2}}{z}_{j\mathrm{1}/\mathrm{2}}=\frac{\mathrm{1}}{\mathrm{2}}({z}_{j+\mathrm{1}}{z}_{j\mathrm{1}})$. For the first time step, a single Euler step is taken of size Δt with timecentered velocities ($t=\mathrm{\Delta}t/\mathrm{2}$).
As usual, this leapfrog timestepping algorithm is only stable to a maximum Courant number of 1. With this staggered grid and vertical grid stretching, the Courant number can be defined by
For the modeone internal wave with background current used in this section, the maximum Courant number is reached at the top and bottom of the domain (where w=0), so $max\left(C\right)=max\left(u\right)/\mathrm{\Delta}x$.
We present results for Eq. (26) at a maximum Courant number of 0.2, chosen to give a “small time step” for later comparison with semiLagrangian results. The results are insensitive to the time step within the stable range, with less than 5 % change in maximumnorm error over the range $\mathrm{0.2}\le max\left(C\right)\le \mathrm{0.99}$.
3.3.3 Results
The error over time of this test case is shown in Fig. 3. As expected, each method achieves secondorder convergence. For the Eulerian advection control case, this is governed by its twopoint central difference scheme. For the semiLagrangian cases, the dominant contribution to the error field comes from the lowerorder vertical interpolation. While the semiLagrangian method has a higher order of accuracy for horizontal motion, here the problem is constructed such that horizontal and vertical motions are of equal importance.
As is often observed with semiLagrangian methods, the overall error of the scheme is somewhat lower for the highCFL case than for the lowCFL case. The interpolation used to evaluate σ off the grid introduces error with each interpolation, and the overall contribution of this error necessarily scales in proportion to the number of interpolations and inversely with the time step size.
Overall, this simplified test case supports the conclusion that the semiLagrangian treatment of advection is a viable replacement for fluxform advection. The semiLagrangian method achieves similar (for lowCFL flows) or better (for highCFL flows) error, and it remains stable for CFL values substantially larger than unity.
With the mechanism for evaluating the f^{B}(x^{D}) term in Eq. (7) established
in Sect. 3, the remaining half of the semiLagrangian advection algorithm is the
estimation of the x^{D} departure points. This corresponds to the positions at the before
time level (t_{0}−Δt) of those fluid parcels that will arrive on the reference grid
x_{ref} the after time level (t_{0}+Δt). One such upstream location
exists for each valid grid location, so in general x_{d} needs to be estimated for each
t
, u
, and v
point on the NEMOOPA grid to provide for (respectively) the
tracer and velocity advective forcings.
In general, calculation of the departure points is an implicit and nonlinear problem, requiring knowledge of the flow velocity at every subgrid place and time between the before and after time levels, before the flow at the latter has been computed. To make this problem tractable, we make a series of simplifying assumptions.
The first such assumption is to freeze the flow, such that trajectories are computed based on strictly the now velocities (that is, u, v, and w at the intermediate time level). This is consistent with the underlying leapfrog timestepping algorithm and the other advection schemes in NEMOOPA, where most fluxes are computed instantaneously with respect to the same now velocities. In physical terms, this constrains fluid parcels to follow paths based on estimated, instantaneous streamlines. In exchange, this decouples the trajectory computation from the after velocities and makes the process time explicit, which eliminates what would otherwise be a need to iterate the entire timestepping process.
4.1 Exponential integration
Ordinarily, the next assumption in the trajectory calculation is to approximate the particle paths, either by a straight line or by a lowdegree polynomial. In this case, the Lagrangian equation
is integrated with an approximate quadrature. Using the trapezoidal rule gives the approximation
where ${\mathit{x}}_{\mathrm{a}}=\mathit{x}({t}_{\mathrm{0}}+\mathrm{\Delta}t)$ is the ongrid “arrival” point of the fluid parcel. This approximation is secondorder in time, and it results in an iterative method where v(x_{d}) is interpolated, leading to a revised estimate of x_{d}.
Unfortunately, this approximation is not suitable for trajectory calculations in the general ocean because it does not appropriately handle flow near a solid boundary. Consider the case of twodimensional flow in the positive halfplane, where fluid velocities are prescribed as $(u,v)=(x,y)$. This forms an analytic continuation of flow near a boundary along the line x=0.
Now, apply equation Eq. (29) to the fluid parcel that arrives at ${\mathit{x}}_{\mathrm{a}}=(\mathrm{1},\mathrm{0})$. Along this streamline, v=0 by inspection so this equation reduces to one dimension and has the solution
For small values of Δt, this solution is reasonable. For Δt>1, however, this solution leads an unphysical trajectory, where the departure point is found to lie in the left halfplane (and thus lies inside the boundary).
The failure here is a specific example of Eq. (29) failing the Lipschitz trajectorycrossing criterion (Smolarkiewicz and Pudykiewicz, 1992), which requires ${u}_{x}\mathrm{\Delta}t<C\approx \mathrm{1}$. The trajectory implied by Eq. (30) crosses the trajectories of fluid parcels that arrive at ${\mathit{x}}_{\mathrm{a}}=(\mathrm{1}\pm \mathit{\u03f5},\mathrm{0})$, and the resulting advection loses its physical meaning.
This trajectorycrossing criterion is a physical limit for solutions which develop discontinuous shocks, such as those that can arise in simulations of the nondispersive, nonlinear shallow water equations. However, these shocks are not typical of threedimensional hydrostatic flows in the ocean, and they are certainly not universally seen at solid boundaries. The true trajectories of fluid parcels, if evaluated exactly, do not cross (and do not have origins inside the land domain), so a better approach is to directly integrate Eq. (28) without approximating the time derivative. Here, this onedimensional system reduces to the ordinary differential equation
with the boundary condition $x({t}_{\mathrm{0}}+\mathrm{\Delta}t)={x}_{a}=\mathrm{1}$. The solution to this equation is obviously of the form x(t)=Cexp (t) for some constant C, and applying the boundary condition gives $x\left(t\right)=\mathrm{exp}(t({t}_{\mathrm{0}}+\mathrm{\Delta}t\left)\right)$ and a departure point of ${x}_{\mathrm{d}}=\mathrm{exp}(\mathrm{2}\mathrm{\Delta}t)$.
This solution is very well behaved, lying exclusively in the right halfplane and asymptotically approaching the wall at x=0 as Δt→∞. This approach works when that of Eq. (29) fails because the direct integration properly captures the exponentialintime path of the fluid parcel.
A generalization of this approach forms the basis for trajectory calculation in this work. Since the solution of Eq. (28) is not analytically possible with an arbitrary velocity field, we exactly solve Eq. (28) based on an approximate, linearly varying velocity field. This is similar to an approach discussed by Walters et al. (2007), where within a single, twodimensional finiteelement cell the linear velocity form is exactly given by the underlying discretization rather than an approximation.
Assume that an arbitrary fluid parcel arrives at x_{a} and that we know the velocity there (v_{a}) and at another point v(x_{c})=v_{c}. We know that the fluid parcel must arrive at x_{a} travelling in the direction of ${\widehat{v}}_{a}={\mathit{v}}_{a}/{\mathit{\alpha}}_{a}$, with ${\mathit{\alpha}}_{a}=\Vert {\mathit{v}}_{a}\Vert $. Then v_{c} can be written in terms of this direction as ${\mathit{v}}_{c}={\mathit{\alpha}}_{c}{\widehat{v}}_{a}+{\mathit{\beta}}_{c}{\widehat{n}}_{a}$, for scalar α_{c} and β_{c} and some ${\widehat{n}}_{a}$ normal to ${\widehat{v}}_{a}$.
This forms a twodimensional system spanned by vectors ${\widehat{v}}_{a}$ and ${\widehat{n}}_{a}$. If we additionally make the assumption that v(x) varies linearly in this plane, we can construct a simplified, twodimensional coordinate system to solve Eq. (28). Here, the origin of the coordinate system corresponds to x_{a}, and the rotated coordinates $\widehat{x}$ and $\widehat{y}$ align with ${\widehat{v}}_{a}$ and ${\widehat{n}}_{a}$ respectively. This implies that x_{c} projects onto the point $({\mathit{x}}_{c}\cdot {\widehat{v}}_{a},{\mathit{x}}_{c}\cdot {\widehat{n}}_{a})=({x}_{c},{y}_{c})$. The linearly interpolated velocities lie strictly in this plane, so the equations of motion for a fluid parcel are
subject to the boundary condition that $x({t}_{\mathrm{0}}+\mathrm{\Delta}t)=y({t}_{\mathrm{0}}+\mathrm{\Delta}t)=\mathrm{0}$. Equation (32a) can be solved first, and applying the boundary condition $x({t}_{\mathrm{0}}+\mathrm{\Delta}t)=\mathrm{0}$ gives
Applying this to Eq. (32b) along with its boundary condition $y({t}_{\mathrm{0}}+\mathrm{\Delta}t)=\mathrm{0}$ gives
When the alongtrajectory acceleration is small ($\left\right({\mathit{\alpha}}_{c}{\mathit{\alpha}}_{a})\mathrm{\Delta}t/{x}_{c}\ll \mathrm{1}$), Eq. (33) reduces to a trapezoidal rule with secondorder accuracy in time.
4.1.1 Trajectory iteration
Evaluating Eq. (33) at $t={t}_{\mathrm{0}}\mathrm{\Delta}t$ and reprojecting the coordinates to the grid forms the basis of an iterative algorithm for trajectories.
This algorithm is ideally suited to cases that look like flow away from a stagnation point, where a fluid parcel is accelerating as it reaches the grid point at t_{0}+Δt. In those cases, the $({\mathit{\alpha}}_{c}{\mathit{\alpha}}_{a})/{x}_{c}$ terms will be positive, and the exponential terms will limit the size of the trajectory for finite Δt. In the opposite case, however, the exponential terms will tend to lengthen the trajectory. For large Δt or large deceleration, this effectively demands that Eqs. (32)–(33) extrapolate beyond the velocity sample at x_{c}, a potential source of instability.
To remedy this, a limiter is added to step 3 of Algorithm 3, whereby x(t_{0}−2Δt) is constrained to the greater^{6} of that from Eq. (33a) and $\mathrm{2}\mathrm{\Delta}tmax({\mathit{\alpha}}_{a},{\mathit{\alpha}}_{c})$. When limiting is necessary, it effectively reduces the time step used for the trajectory iteration, so for consistency a revised Δt^{′} is computed by inverting Eq. (33a) with the limited ${x}_{c}^{\prime}$, which is then used to evaluate Eq. (33b).
4.2 Underrelaxation and land boundaries
While the construction of Algorithm 3 guarantees that trajectories cannot converge to an outofboundary point, there are no guarantees that the algorithm remains in boundary during the iteration process or that the iteration converges. The problem of a divergent or oscillatory iteration is more likely when the underlying velocity field does not resemble the linearly approximated velocity field integrated by Eq. (32), as then each iteration might result in very different approximations.
Addressing the latter point first, this work heuristically applies underrelaxation when Algorithm 3 is slow to converge. After 10 local iterations, step 3 is replaced by ${\mathit{x}}_{c}\leftarrow \frac{\mathrm{1}}{\mathrm{2}}({\mathit{x}}_{c}+{\mathit{x}}_{c}^{\prime})$; after 20 iterations the righthand side becomes $\frac{\mathrm{1}}{\mathrm{4}}(\mathrm{3}{\mathit{x}}_{c}+{\mathit{x}}_{c}^{\prime})$, and after 30 iterations the righthand side becomes $\frac{\mathrm{1}}{\mathrm{8}}(\mathrm{7}{\mathit{x}}_{c}+{\mathit{x}}_{c}^{\prime})$. At 40 iterations, the trajectory is truncated by ending the iteration with the first indomain point returned from the process; this ensures some sort of advection even if the iterative process enters a limit cycle.
This underrelaxation also addresses the possibility that x_{c} might lie outside of the ocean domain. If x_{c} is masked, then there is no valid velocity to provide via offgrid interpolation, so instead of evaluating Eq. (33), ${\mathit{x}}_{c}^{\prime}$ is set to x_{a} in step 3 of Algorithm 3. This combines with the underrelaxation after 10 iterations to reduce the trajectory length until an inboundary point is found, whereupon iteration resumes normally.
These values for iteration counts and underrelaxation parameters are conservatively specified. In the numerical tests discussed in this work, the vast majority of trajectories converge after one or two iterations, without needing to resort to underrelaxation or trajectory truncation.
4.3 Velocity interpolation
The trajectory algorithm requires the offgrid interpolation of velocities at each iteration. In principle, these velocities can be interpolated using the interpolation process of Sect. 3. Doing so would be ideal for ultimate consistency with the final offgrid interpolation, but this process is also computationally expensive. In practice, it is more efficient to evaluate the offgrid velocity field in step 2 of Algorithm 3 using trilinear interpolation; doing so causes little change in the numerical test cases in this work.
Trilinear interpolation proceeds with the same order of operations as Algorithm 1: velocities are first interpolated in depth to the (x,y) corners of the grid box at the offgrid level, then along the x direction, and finally along the y direction. Each individual interpolation respects the relevant boundary condition, so for example the u velocity is considered to reflect symmetrically around a boundary in y and z but is constrained to zero at a boundary in x.
One complication of linear interpolation, however, is that the velocity points are staggered by half a cell with respect to the physical boundary. In two dimensions, if the tracer point T(0,0) (to use gridcell coordinates for the tracer grid denoted T) is a land point but T(0,1), T(1,0), and T(1,1) are all ocean points, then uvelocity point U(0,0) (denoting the uvelocity grid as U), halfway between T(0,0) and T(1,0), lies along the boundary. The boundary continues to U(0,0.5), whereupon U(0,0.5)–U(0,1) lies inside the ocean. This violates a basic assumption of linear interpolation, that the velocity should vary smoothly (and approximately linearly) within the u cell.
This causes two problems for trajectory computation. The first problem is that after repeated onedimensional interpolation, the boundary condition is no longer necessarily respected by the interpolated velocity, which can result in a trajectory iteration that “pushes” the departure point into the wall, causing nonconvergence. The second problem is that while the interpolation process guarantees continuity of the interpolated field at the cell corners, the boundary conditions can cause large discontinuities along the cell edges, again resulting in a convergence failure. In the above example, the interpolated velocity at $U(+\mathit{\u03f5},\mathrm{0.6})$ would be influenced by both U(0,1) and the zero velocity at the physical boundary of U(0,0), but the interpolated velocity at $U(\mathit{\u03f5},\mathrm{0.6})$ would be influenced by U(0,1) and its reflection at a ghost point. These problems are illustrated in Fig. 4a.
The solution to both of these problems is to blend the linearly interpolated function with a corner singularity solution. A bilinear function is a solution to Laplace's equation (∇^{2}f=0), so it is reasonable to consider corner solutions that are also solutions to Laplace's equation.
Without loss of generality, consider a grid cell defined by $(x,y)\in [\mathrm{0},\mathrm{1}{]}^{\mathrm{2}}$, such that there is a solid boundary along $(x=\mathrm{0},y<\mathrm{0.5})$ as depicted in Fig. 4. Treating the boundary as an infinite halfplane, with $f(\mathrm{0},y)=\mathrm{0}$ for y<0.5 and $f(\mathrm{0},y)=f(\mathrm{0},\mathrm{1})$ for y>0.5, the “corner” solution to Laplace's equation is
while bilinear interpolation would give
These two solutions are blended together, with Eq. (34) taking precedence along the solid boundary (x=0 and $\mathrm{0}\le y\le \mathrm{0.5}$) and Eq. (35) taking precedence along the x=1 and y=1 boundaries of the cell. This gives
where $\mathit{\sigma}=max(\mathrm{1}x,\mathrm{2}(y\mathrm{0.5}\left)\right)$.
The blended function exactly respects the solid boundary condition, and the discontinuity at the cell edges is significantly reduced. Blended functions for other configurations of the solid wall are given by applying the appropriate reflections and rotations to Eq. (36).
5.1 Flow past an island
To demonstrate the impacts of semiLagrangian advection on a simple test case with a lengthened time step, we first present the quasitwodimensional test case of isothermal flow past an interposed island.
This test case consists of a $\mathrm{280}\times \mathrm{70}\times \mathrm{3}$ point grid, with grid resolution $\mathrm{\Delta}x=\mathrm{\Delta}y=\mathrm{5}\phantom{\rule{0.125em}{0ex}}\mathrm{m}$ and Δz=10 m. A 50 m×50 m region (10 points ×10 points) is masked as land in the middle of the domain. The inflow boundary condition is set to $u=\mathrm{0.03}\phantom{\rule{0.125em}{0ex}}\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{1}}$, v=0; this was also imposed throughout the domain as an initial condition. The reference frame was also irrotational, with a Coriolis parameter of 0.
Relevant namelist parameters are given in Table 1, with parameters that differ between the control and semiLagrangian runs highlighted. The control run used fluxform velocity advection^{7} via the QUICKEST scheme (Leonard, 1979, 1991), whereas the semiLagrangian run used semiLagrangian advection of momentum in flux form as described in Sects. 3 and 4. To emphasize the dynamical differences between the advection schemes, both test cases were run with no explicit horizontal diffusion of momentum. Vertical mixing terms, largely irrelevant for this quasitwodimensional case, were set consistently with the ORCA025 simulations in Sect. 5.2.
Both series of runs used the implicit freesurface formulation (enabled with the compiletime key
key_dynspg_flt
), which damped the large initial surface gravity waves caused by the
imposition of the blocking island on the steadystate flow.
After the initial gravitywave adjustment, this test case quickly develops a set of recirculating vortices in the lee of the island. Over time these vortices grow in extent and would begin detaching to form a vortex street, but this does not happen before the 8000 s end of the simulation. Although there is no explicit horizontal diffusion of momentum in these runs, the flow regime is much more laminar than would be implied by the physical Reynolds number of 1.5×10^{6}, based on the freestream velocity, edgelength of the island, and molecular viscosity of water.
In moving around the box, the flow locally accelerates to a maximum steady velocity of about 0.05 m s^{−1}, and this maximum velocity is reached in the vicinity of the leadingedge corners of the box. The exact value of this maximum depends on both the simulation time and the time step, but our expected pattern holds: the control simulation is stable with a time step of 64 s, which corresponded to a maximum steady Courant number, i.e. $max\left({u}_{\text{steady}}\right)/\mathrm{\Delta}x$, above 0.6 (and a maximum transient CFL of 0.95), but it is unstable with a time step of 80 s.
SemiLagrangian advection maintains stability for much longer time steps. Figure 5 shows the freesurface height and flow streamlines for Δt between 5 and 160 s, and only the semiLagrangian method remains stable for 80 and 160 s time steps. For both advection schemes, the longer time step is associated with a more diffuse flow pattern, with lengthening (and less intense) recirculating vortices in the lee of the island.
This effect is stronger with semiLagrangian advection than with Eulerian advection. We attribute this to the nature of the flow at the leading edge of the island. Here, the dominant flow balance is cyclostrophic, where the pressure gradient at these corners balances the local vorticity. The operator splitting method used here treats the advective terms in a frame following the flow, but it can only apply the pressure force at the destination cell. This results in an inconsistency that grows with Δt, related to the forces in equation Eq. (2) being available only at the endpoint of the Lagrangian trajectory – an O(Δt) approximation to the integral.
This inconsistency is most evident in the 160 s time step case (bottom panel in Fig. 5b), where the maximum steady Courant number of 1.6 implies that fluid parcels are advected by about three grid cells over the 2Δt leapfrog step. There, the lowest pressure region at the leading edge of the flow has moved slightly further downstream.
In the full ocean, the geostrophic effect predominates, with a leadingorder balance between the pressure gradient and the Coriolis force (planetary vorticity), so we expect this issue to be less pronounced.
5.2 Global forced runs
To evaluate semiLagrangian advection in a more realistic forecasting setting, we conducted a preliminary series of free runs of the NEMOOPA model. The runs consisted of the following:

a control run, based on the configuration of Environment and Climate Change Canada's $\mathrm{1}/\mathrm{4}{}^{\circ}$ nominalresolution Global Ice/Ocean Prediction System (GIOPS) (Smith et al., 2016) with a 10 min model time step;^{8} tracers were advected with the model's tracer variance dissipation scheme (Lévy et al., 2001), and momentums were advected in vector form with the model's energy and enstrophy conserving scheme^{9} (Arakawa and Lamb, 1981);

a “semiLagrangian tracer” run, where momentum was advected as in the control scheme and the semiLagrangian advection described in this work was used for advection of salinity and temperature; additionally, this run disabled horizontal diffusion of salinity and temperature; and

a “semiLagrangian momentum and tracer” run, where momentum as well is advected with the semiLagrangian scheme; the configuration was otherwise the same as the semiLagrangian tracer run, save for a 15 min model time step.
The runs were all initialized at 1 October 2001 on the ORCA025 grid. The ocean was at rest, and temperature and salinity were given by the 2011 World Ocean Atlas climatology (Locarnini et al., 2013; Zweng et al., 2013). Atmospheric forcing was provided at 1 h intervals from Environment and Climate Change Canada's $\mathrm{1}/\mathrm{4}{}^{\circ}$ global atmospheric reforecast, and sea ice was modeled via coupling with version 4.0 of the CICE model (Hunke and Dukowicz, 1997), with dynamically active (moving) ice. Selected namelist parameters are provided in Table 2.
As with Sect. 5.1, the test cases used NEMOOPA's linear freesurface with a timeimplicit solver, and tidal forcing was not present in these configurations. In a typical time step, the vast majority of semiLagrangian trajectories converged in one iteration (mean 1.004 over the semiLagrangian tracer run). A very small minority of cells required an extended number of iterations or underrelaxation as described in Sect. 4.2, but this did not affect the overall trajectorycalculation performance because convergence was measured (and iterations limited) on a percell basis.
Each run continued through late 2009. For reasons of space efficiency, we recorded the twodimensional sea surface height, temperature, and salinity fields for each model day, and we preserved every fifth dailymean, threedimensional output of temperature, salinity, and horizontal ocean velocity.
For short and mediumterm forecasts, the operational coupled forecasting systems at CMC are constrained by observations and periodic reinitialization. With a focus on this forecasting horizon, the objectives of these long freeruns were as follows:

to provide a test of model stability with semiLagrangian advection, in terms of both avoiding crashes and providing plausible ocean fields;

to check for any largescale conservation errors, which might be difficult to correct given the sparsity of observation data for the deep ocean; and

to note any qualitative improvement or deterioration in the effective resolution of the model.
This first goal of model stability was met in part by the successful completion of these runs. Use of semiLagrangian advection for both tracers and momentum allowed us to increase the effective time step from 10 min (with typical maximum Courant number^{10} of 0.2, found in the vertical direction) to 15 min (Courant number 0.3). Further increases led to instability and model crashes not from the advection component, but from the ice model. In this version of the model, the ocean–ice stress is coupled in a timeexplicit way between the water and ice components. Concurrent work towards a timeimplicit coupling has given encouraging preliminary results on further time step increases.
The use of semiLagrangian advection also gives global flows qualitatively similar to the control run, and average transports in the Atlantic overturning circulation and Circumpolar current are comparable between the control and semiLagrangian runs (Fig. 6). The semiLagrangian runs appear to result in a slightly weaker overturning circulation and a slightly stronger circumpolar current than the control run, but these results may not be robust to retuned physical parameterizations. Using semiLagrangian advection for the velocity components results in a significant decrease to overall ocean kinetic energy (Fig. 7), both during and after the spinup period.
The cause of this energy disparity is under investigation, but we believe the most likely cause is the application of slope limiting to the u and v fields independently. Future work will focus on taking a more nuanced approach to filtering, but this effect may not be very significant in a shorterterm forecast setting with frequent reinitializations from an analysis.
The second goal of global conservation was met. Although semiLagrangian advection does not guarantee conservation of tracers, the impact on the global balance of temperature and salinity was small. Figure 8 shows the evolution of oceanaverage temperature and salinity over time in these runs, and the effect of nonconservation attributable to the semiLagrangian advection of tracers is comparable to the magnitude of uncertainty in the global balance of atmospheric forcing – the imbalance seen in the control run. Each case saw an overall temperature drift of about 0.04 K over the simulated period, with the semiLagrangian cases having a slight warming trend against the control run's slight cooling trend, and all three runs showed a very small increase in oceanaverage salinity, by about 0.01 PSU.
The temperature change vs. depth over the simulated period is shown in Fig. 9. Both the control and semiLagrangian runs showed a warming trend in the surface layers, but the semiLagrangian runs showed temperature stability in fluid layers below 1000 m depth, whereas the control run showed a cooling trend in these waters.
Despite the energy shortfall with semiLagrangian advection of momentum, we see tentative signs that the method increases the model's effective resolution. Figure 10 shows one particular sea surface temperature realization, from the 31 December 2005 of each test case, along with the magnitude of the temperature gradient. The largescale flows are similar between the control and semiLagrangian runs (and most similar between the control and semiLagrangian tracer run), but the semiLagrangian runs have noticeably stronger gradients in the sea surface temperature, in patterns that resemble smallerscale eddies.
This work has derived a semiLagrangian advection scheme for the NEMOOPA model. After advecting a tracer or momentum field along estimated fluid parcel trajectories, it calculates a time trend to provide to the remainder of the model; in this way the semiLagrangian scheme serves as a dropin replacement for other tracer and (fluxform) momentum schemes.
The development of this advection module relied on several new or newly applied algorithms that might be relevant to other ocean models or other domains:

The “semiLagrangian trend” form of equation Eq. (7) might be useful in other models when researchers wish to implement semiLagrangian advection after the fact, without disrupting the calculation of other forcing terms.

The Hermite interpolation form in Sect. 3, especially combined with the C^{1}continuous estimate of the vertical derivative in Sect. 3.2 might find application in other domains where, as in the ocean, one dimension (the vertical) is more oscillatory than others.

The exponential integration of trajectories in Eq. (4) may be useful in other applications that feature strong accelerations over trajectories. In particular, it forbids trajectorycrossing in onedimensional flows, and here that property ensures that trajectories remain inside the ocean domain.

The “corner solution” treatment of velocity for trajectory calculations near corners might find use in other applications with staggered velocity components.
Overall, we find that the semiLagrangian method is effective at extending the realizable time step in the NEMOOPA model. In the simple domain of Sect. 5.1, this method resulted in a stable simulation with advective Courant numbers in excess of 1. Although we only extended the time step from 10 to 15 min for the semiLagrangian momentum run in Sect. 5.2, this limitation was imposed by the ice model. Disabling ice dynamics allowed us to increase the time step to 30 min, but this would have made the results incomparable with those of the control and semiLagrangian tracer runs. Preliminary work with the CICE sea model and implicit coupling of the ice–ocean stress seems to allow us to relax the icerelated time step restriction.
6.1 Performance and implementation
In spite of this increased time step, the semiLagrangian method by itself does not yet improve overall computational performance. The semiLagrangian momentum and tracer run of Sect. 5.2 took approximately 1 h of computational time per 5 d of simulated time, using 128 Intel Xeon E5530 processors at 2.4 GHz. With a 10 min time step, the semiLagrangian tracer run took approximately 50 min for the same 5 d of simulated time, whereas the control run took just 30 min. We expect these results to improve with further numerical optimization work. In particular, we did not take great care to ensure that loops were vectorized where possible, and it is much more difficult for compilers to automatically vectorize the pointbypoint semiLagrangian computations compared to volume flux calculations in the traditional advection schemes.
About onethird of the additional computational cost comes from trajectory iterations, and the remainder comes from the cubic interpolation. This suggests that the relative cost of semiLagrangian advection will be lower than presented here if trajectories can be reused for multiple tracer species (such as biogeochemical constituents). Additionally, it suggests that a further optimization may be to reuse tracer trajectories for momentum advection, at least away from the boundaries where interpolating (staggering) trajectories might be reasonable. It seems unlikely that optimization will reduce the pertimestep penalty to the 20 % value seen by Ritchie et al. (1995) for an atmospheric model – owing to the lack of threedimensional implicit equations and expensive physical parameterizations elsewhere in NEMO0PA – but we are hopeful that semiLagrangian advection will nonetheless improve overall system performance.
The parallel (MPI) implementation of this algorithm was straightforward. With the relatively modest
increase in Courant number for the cases in this work, we simply needed to increase the
interprocessor lateral halo (parameters jpreci
and jprecj
) to three points, which
was sufficient to allow a fluid parcel arriving at a processor's edge to apply the full
interpolating stencil for the Courant numbers reached in the presented simulations. This increase
in halo size was small compared to the processor tile size of about 50×260 grid points for
the runs in Sect. 5.2. Extending this to support very large horizontal Courant
numbers, however (if another solution could be found to stabilize baroclinic waves), would require
either prohibitively large halo sizes or additional interprocessor communication to track fluid
parcels that cross MPI boundaries.
6.2 Qualitative comments on results
Although the semiLagrangian method does not guarantee tracer conservation, we see no evidence that its implementation here leads to a degradation relevant in a weekly to seasonal forecast setting. In particular, the deepwater temperature stability shown in Fig. 9 is an encouraging sign that semiLagrangian advection will preserve the deepwater structure that is weakly constrained by data. Even small imbalances, however, might become significant over the decadetocentury timescales of climate simulations. Further work will be necessary to characterize this method before we can safely recommend semiLagrangian advection in such settings.
For the test cases in Sect. 5.2, semiLagrangian advection of tracers appears to slightly increase the effective resolution of the model. However, this effect is much more mixed when momentum is also advected with the semiLagrangian method, in part because the underlying currents differ. Both of these differences will be the subject of future study, with the specific intention of assessing these effects in the setting of shorterterm forecasts. We speculate that the overall loss of kinetic energy with semiLagrangian advection of momentum is attributable to the use of the slope limiter: limiting each component of velocity separately may be causing unrealistic diffusion of smallerscale structures in the presence of background vorticity. We hope to address this issue with more selective limiting.
6.3 Future development
Finally, the development in this paper implicitly assumes that the coordinate system is static with time. This is not the case in NEMOOPA when using its nonlinear freesurface option, which necessarily implies timevarying vertical levels. Adapting the semiLagrangian method to this more general coordinate system will be a focus of future work, which will be required to apply this advection scheme to higherresolution domains that require tidepermitting simulations.
Additionally, future versions of NEMO intend to move to a thirdorder Runge–Kutta timestepping algorithm (Shu and Osher, 1997), which constructs a full time step as a linear combination of forward Euler steps. We expect that the semiLagrangian “advective trend” of Eq. (7) can be adapted to this framework in a straightforward manner by basing the calculated trend on the currentstep values of tracers and velocities, but the adaptation may require care to preserve the higherorder temporal accuracy of the overall scheme.
The modified NEMO (CeCILL license, version 2.0) code along with scripts and data used in this paper are available at https://doi.org/10.5281/zenodo.3905100 (Subich et al., 2020). The modified CICE 4.0 used in Sect. 5.2 is not redistributable under a free license, but it has been made available for the topical editors and anonymous reviewers.
All authors were responsible for the concept. CS, PP, and FD were responsible for the necessary software development. CS contributed to the original draft preparation, and all authors contributed to review and editing.
The authors declare that they have no conflict of interest.
The authors would like to acknowledge Francois Roy of Environment Canada, who provided considerable technical support, especially with the forcing runs of Sect. 5.2, and Jerome Chanut of Mercator Ocean, who provided helpful comments on a draft of this paper. The authors also are grateful for the peer review provided by Florian Lemarie and Mike Bell, which greatly improved this paper from its original submission.
The authors also owe a debt of inspiration to the semiLagrangian dynamical core of MC2 (Thomas et al., 1998).
This paper was edited by Simone Marras and reviewed by Florian Lemarie and Mike Bell.
Arakawa, A. and Lamb, V. R.: A Potential Enstrophy and Energy Conserving Scheme for the Shallow Water Equations, Mon. Weather Rev., 109, 18–36, https://doi.org/10.1175/15200493(1981)109<0018:APEAEC>2.0.CO;2, 1981. a
Asselin, R.: Frequency Filter for Time Integrations, Mon. Weather Rev., 100, 487–490, https://doi.org/10.1175/15200493(1972)100<0487:FFFTI>2.3.CO;2, 1972. a
Bermejo, R. and Staniforth, A.: The Conversion of SemiLagrangian Advection Schemes to QuasiMonotone Schemes, Mon. Weather Rev., 120, 2622–2632, https://doi.org/10.1175/15200493(1992)120<2622:TCOSLA>2.0.CO;2, 1992. a
Casulli, V. and Cheng, R. T.: Semiimplicit finite difference methods for threedimensional shallow water flow, Int. J. Numer. Methods Fluid., 15, 629–648, https://doi.org/10.1002/fld.1650150602, 1992. a
Ducousso, N., Le Sommer, J., Molines, J. M., and Bell, M.: Impact of the “Symmetric Instability of the Computational Kind” at mesoscale and submesoscalepermitting resolutions, Ocean Model., 120, 18–26, https://doi.org/10.1016/j.ocemod.2017.10.006, 2017. a
Ferry, N., Parent, L., Masina, S., Storto, A., Zuo, H., and Balmaseda, M.: PRODUCT USER MANUAL For Global Ocean Reanalysis Products GLOBALREANALYSISPHYS001009 GLOBALREANALYSISPHYS001011 GLOBALREANALYSISPHYS001017, Tech. Rep. 2.9, Copernicus Marine Service, available at: https://resources.marine.copernicus.eu/documents/PUM/CMEMSGLOPUM001009011017.pdf, (last access: 12 June 2020), 2016. a
Girard, C., Plante, A., Desgagné, M., McTaggartCowan, R., Côté, J., Charron, M., Gravel, S., Lee, V., Patoine, A., Qaddouri, A., Roch, M., Spacek, L., Tanguay, M., Vaillancourt, P. A., and Zadra, A.: Staggered vertical discretization of the Canadian Environmental Multiscale (GEM) model using a coordinate of the loghydrostaticpressure type, Mon. Weather Rev., 142, 1183–1196, 2014. a, b
Hildebrand, F.: Introduction to Numerical Analysis: Second Edition, McGrawHill Inc., 1974. a
Hodges, B. R. and Dallimore, C.: Estuary, Lake and Coastal Ocean Model: ELCOM v2.2 science manual, Tech. rep., Center for Water Research, University of Western Australia, 2006. a
Hollingsworth, A., Kållberg, P., Renner, V., and Burridge, D.: An internal symmetric computational instability, Q. J. Roy. Meteor. Soc., 109, 417–428, 1983. a
Hunke, E. C. and Dukowicz, J. K.: An elastic–viscous–plastic model for sea ice dynamics, J. Phys. Oceanogr., 27, 1849–1867, https://doi.org/10.1175/15200485(1997)027<1849:AEVPMF>2.0.CO;2, 1997. a
Lauritzen, P. H.: A Stability Analysis of FiniteVolume Advection Schemes Permitting Long Time Steps, Mon. Weather Rev., 135, 2658–2673, https://doi.org/10.1175/MWR3425.1, 2007. a
Leclair, M. and Madec, G.: $\stackrel{\mathrm{\u0303}}{z}$Coordinate, an Arbitrary Lagrangian–Eulerian coordinate separating high and low frequency motions, Ocean Model., 37, 139–152, https://doi.org/10.1016/j.ocemod.2011.02.001, 2011. a
Lemarié, F., Debreu, L., Madec, G., Demange, J., Molines, J. M., and Honnorat, M.: Stability constraints for oceanic numerical models: implications for the formulation of time and space discretizations, Ocean Model., 92, 124–148, https://doi.org/10.1016/j.ocemod.2015.06.006, 2015. a, b
Lemieux, J.F., Beaudoin, C., Dupont, F., Roy, F., Smith, G. C., Shlyaeva, A., Buehner, M., Caya, A., Chen, J., Carrieres, T., Pogson, L., DeRepentigny, P., Plante, A., Pestieau, P., Pellerin, P., Ritchie, H., Garric, G., and Ferry, N.: The Regional Ice Prediction System (RIPS): verification of forecast sea ice concentration, Q. J. Roy. Meteor. Soc., 142, 632–643, https://doi.org/10.1002/qj.2526, 2016. a
Leonard, B. P.: A stable and accurate convective modelling procedure based on quadratic upstream interpolation, Comput. Methods Appl. Mech. Eng., 19, 59–98, https://doi.org/10.1016/00457825(79)900343, 1979. a
Leonard, B. P.: The ULTIMATE conservative difference scheme applied to unsteady onedimensional advection, Comput. Methods Appl. Mech. Eng., 88, 17–74, https://doi.org/10.1016/00457825(91)90232U, 1991. a
Lévy, M., Estublier, A., and Madec, G.: Choice of an advection scheme for biogeochemical models, Geophys. Res. Lett., 28, 3725–3728, https://doi.org/10.1029/2001GL012947, 2001. a
Locarnini, R. A., Mishonov, A. V., Antonov, J. I., Boyer, T. P., Garcia, H. E., Baranova, O. K., Zweng, M. M., Paver, C. R., Reagan, J. R., Johnson, D. R., Hamilton, M., and Seido, D.: World Ocean Atlas 2013, Volume 1: Temperature, 2013. a
Madec, G.: NEMO ocean engine, Note du Pôle de modélisation, Institut PierreSimon Laplace (IPSL), France, No. 27, 2008. a, b, c
Madec, G. and Imbard, M.: A global ocean mesh to overcome the North Pole singularity, Clim. Dynam., 12, 381–388, https://doi.org/10.1007/BF00211684, 1996. a
Mesinger, F. and Arakawa, A.: Numerical methods used in atmospheric models, Global Atmospheric Research Programme (GARP), World Meteorological Organization, 1976. a
Murray, R. J.: Explicit Generation of Orthogonal Grids for Ocean Models, J. Comput. Phys., 126, 251–273, https://doi.org/10.1006/jcph.1996.0136, 1996. a
Ritchie, H., Temperton, C., Simmons, A., Hortal, M., Davies, T., Dent, D., and Hamrud, M.: Implementation of the SemiLagrangian Method in a HighResolution Version of the ECMWF Forecast Model, Mon. Weather Rev., 123, 489–514, https://doi.org/10.1175/15200493(1995)123<0489:IOTSLM>2.0.CO;2, 1995. a
Robert, A.: A semiLagrangian and semiimplicit numerical integration scheme for the primitive meteorological equations, J. Meteorol. Soc. Jpn. Ser. II, 60, 319–325, 1982. a, b
Roy, F., Chevallier, M., Smith, G. C., Dupont, F., Garric, G., Lemieux, J.F., Lu, Y., and Davidson, F.: Arctic sea ice and freshwater sensitivity to the treatment of the atmosphereiceocean surface layer, J. Geophys. Res.Oceans, 120, 4392–4417, https://doi.org/10.1002/2014JC010677, 2015. a
Shchepetkin, A. F. and McWilliams, J. C.: The regional oceanic modeling system (ROMS): a splitexplicit, freesurface, topographyfollowingcoordinate oceanic model, Ocean Model., 9, 347–404, https://doi.org/10.1016/j.ocemod.2004.08.002, 2005. a
Shu, C.W. and Osher, S.: Efficient Implementation of Essentially Nonoscillatory ShockCapturing Schemes, II, in: Upwind and HighResolution Schemes, edited by: Hussaini, M. Y., van Leer, B., and Van Rosendale, J., Springer, Berlin, Heidelberg, 328–374, https://doi.org/10.1007/9783642605437_14, 1997. a
Smith, G., Liu, Y., Colan, D. S., Reszka, M., Roy, F., Dupont, F., Lemieux, J.F., Beaudoin, C., Lellouche, J.M., Garric, G., Testut, C.E., Pellerin, P., Ritchie, H., Lu, Y., and Davidson, F.: The CMC Global Ice Ocean Prediction System v2.1, Tech. rep., Canadian Centre for Meteorological and Environmental Prediction, available at: http://collaboration.cmc.ec.gc.ca/cmc/cmoi/product_guide/docs/tech_notes/technote_giops210_e.pdf, (last access: 15 September 2020), 2016. a
Smith, G. C., Bélanger, J.M., Roy, F., Pellerin, P., Ritchie, H., Onu, K., Roch, M., Zadra, A., Colan, D. S., Winter, B., Fontecilla, J.S., and Deacu, D.: Impact of Coupling with an Ice–Ocean Model on Global MediumRange NWP Forecast Skill, Mon. Weather Rev., 146, 1157–1180, https://doi.org/10.1175/MWRD170157.1, 2018. a, b
Smolarkiewicz, P. K. and Pudykiewicz, J. A.: A Class of SemiLagrangian Approximations for Fluids, J. Atmos. Sci., 49, 2082–2096, https://doi.org/10.1175/15200469(1992)049<2082:ACOSLA>2.0.CO;2, 1992. a
Subich, C., Pellerin, P., Smith, G., and Dupont, F.: Code repository for “Development of a semiLagrangian advection scheme for the NEMO ocean model (3.1)”, a GMD discussion paper, Zenodo, https://doi.org/10.5281/zenodo.3905100, 2020. a
Thomas, S. J., Girard, C., Benoit, R., Desgagné, M., and Pellerin, P.: A new adiabatic kernel for the MC2 model, Atmos.Ocean, 36, 241–270, https://doi.org/10.1080/07055900.1998.9649613, 1998. a
Turkington, B., Eydeland, A., and Wang, S.: A Computational Method for Solitary Internal Waves in a Continuously Stratified Fluid, Stud. Appl. Math., 85, 93–127, https://doi.org/10.1002/sapm199185293, 1991. a
Walters, R. A., Lane, E. M., and Henry, R. F.: SemiLagrangian methods for a finite element coastal ocean model, Ocean Model., 19, 112–124, https://doi.org/10.1016/j.ocemod.2007.06.008, 2007. a, b
Zweng, M., Reagan, J., Antonov, J., Locarnini, R., Mishonov, A., Boyer, T., Garcia, H., Baranova, O., Johnson, D., Seidov, D., and Biddle, M.: World Ocean Atlas Volume 2: Salinity, 2013. a
This is true for temperature, salinity, and momentum provided the ocean is treated as an incompressible fluid. This assumption is satisfied by NEMOOPA's adoption of the Boussinesq approximation.
If x_{d} lies along an edge or corner of this interval, then at least one of the resulting interpolations will be trivial. In that case, the choice of which neighbouring interval x_{d} lies “within” is arbitrary.
This is not necessarily the case, however, for grids that have manually specified, nonsmooth regions of enhanced resolution. In such cases a more nuanced treatment of interpolation would be advisable.
Since this section tests advection alone, the scaling of σ is not dynamically relevant. In fact, the wave structure of Eq. (22) corresponds to an exact internal mode of the incompressible Navier–Stokes equations for a linear stratification.
The trajectory calculation scheme of Sect. 4 could be used instead, but since the overall trajectory lengths are small compared to the length scales of the velocity field (L_{x} and L_{z}), that method would give equivalent results to the trapezoidal rule.
Since the rotated x axis is aligned with the fluid velocity at x_{a}, x_{c} is generally negative in the rotated frame.
This choice of velocity advection provided the best results for the control run, of the advection models supported in NEMO version 3.1.
This time step is shorter than other commonlyused ORCA025 configurations, such as in the ocean reanalysis of Ferry et al. (2016). This shorter time step is required to stabilize the coupling of ocean–ice stress with the CICE model, where following Roy et al. (2015) the ice–ocean drag coefficient is larger than typically considered. We chose to maintain this configuration and coupling approach to provide for the cleanest likeforlike comparisons against the operational configuration
For compatibility with the operational model, as run in this work, the scheme did not include the “fix” for the Hollingsworth instability (Hollingsworth et al., 1983) reported in Ducousso et al. (2017). This instability is more prominent at higher resolutions, and we do not believe it meaningfully impacted the results as presented in this section.
Defining the Courant numbers in each direction as $max\left(\rightu\left\right)/\mathtt{\text{e1u}}$, $max\left(\rightv\left\right)/\mathtt{\text{e2v}}$, and $max\left(\rightw\left\right)/\mathtt{\text{e3w}}$.