A Primer on the Vertical Lagrangian‐Remap Method in Ocean Models Based on Finite Volume Generalized Vertical Coordinates

This paper provides a primer on the mathematical, physical, and numerical foundations of ocean models that are formulated using finite volume generalized vertical coordinate equations and that use the vertical Lagrangian‐remap method to evolve the ocean state. We consider the mathematical structure of the governing ocean equations in both their strong formulation (partial differential equations) and weak formulation (finite volume integral equations), thus enabling an understanding of their physical content and providing a physical‐mathematical framework to develop numerical algorithms. A connection is made between the Lagrangian‐remap method and the ocean equations as written using finite volume generalized vertical budgets. Thought experiments are offered to exemplify the mechanics of the vertical Lagrangian‐remap method and to compare with other methods used for ocean model algorithms.


Introduction
A finite volume generalized vertical coordinate (GVC) formulation of the ocean equations, when combined with the vertical Lagrangian-remap method, enables a wide suite of vertical coordinate options. It also naturally allows for an exactly conserving realization of wetting and drying by following the numerical algorithms required for vanishing layers in isopycnal models. Since its inception for numerical ocean models by Bleck (2002), the vertical Lagrangian-remap method has realized scientific and engineering maturity in the ocean modeling community. In particular, numerical algorithms based on the vertical Lagrangian-remap method form the dynamical core for community ocean models such as HYCOM (Bleck, 2002) and MOM6 (Adcroft et al., 2019). These models have been configured in a variety of applications related to process physics, realistic operational prediction, and comprehensive earth system science. Use of the vertical Lagrangian-remap method for ocean models is motivated by the numerical accuracy and flexibility it offers relative to alternatives. In response to this promise, we here offer a primer on ocean models that use the vertical Lagrangian-remap method as framed within finite volume GVC equations.

An Outline of the Lagrangian-Remap Method
As depicted in Figure 1, the Lagrangian-remap method has three main steps: • Lagrangian step: fluid equations are solved in a Lagrangian reference frame whereby the grid mesh follows the fluid flow; • regrid/regularization step: the grid mesh is moved to a specified target location; and • remap step: the fluid state is conservatively remapped onto the target grid.
The regrid step generally specifies a movement of the grid relative to the fluid thus complementing the Lagrangian step whereby the grid moves with the fluid. The combined regrid-remap steps can be numerically identical to advection in an Eulerian method (Dukowicz & Baumgardner, 2000). When using GVCs, important differences from a fixed grid approach arise in practice due to the movement of the grid with the flow in the Lagrangian step, thus facilitating the reduction of cross-grid flow and the corresponding spurious calculation of cross-grid transport that occurs with traditional approaches (e.g., Griffies, Pacanowski, & Hallberg, 2000).
As realized in ocean models, the vertical Lagrangian-remap method considers only the vertical direction for the Lagrangian step, with fixed horizontal positions for the mesh. In the ocean modeling community the vertical Lagrangian-remap method is sometimes considered a particular realization of the Arbitrary Lagrangian-Eulerian (ALE) method (Donea et al., 2004;Hirt et al., 1974). Indeed, ALE involves three steps very much aligned with those outlined above for Lagrangian remapping. However, ALE generally prescribes the grid motion, and that motion is not typically Lagrangian. Furthermore, Lagrangian methods date back to the origins of computational fluid mechanics (von Neumann & Richtmyer, 1950) whereas ALE evolved as a particular means to address limitations of early Lagrangian methods (Noh, 1964). We here follow terminology used in the computational physics community by referring to the method utilized by Bleck (2002) as vertical Lagrangian remapping rather than ALE.

Focus for This Paper and Presentation Style
Finite volume-based ocean equations written in GVCs provide the theoretical basis for many ocean models developed for applications ranging from the coast to the globe. These equations are key to a formulation of ocean model dynamical cores based on the vertical Lagrangian-remap method. Our central goal for this paper is to accessibly present the physical, mathematical, and numerical basis for such models. That is, we aim to articulate concepts and methods forming the foundation of ocean models, thus allowing the methods to be understood sufficiently well that they can be scrutinized by a broad range of readers and in turn improved upon by future innovations. This paper thus supports a community understanding of modern numerical ocean models, and by doing so we hope to enhance the utility, transparency, integrity, and maturity of the science going into ocean models as well as the science emerging from their simulations.
To approach these goals we present the material in the form of a primer or tutorial. The reader is expected to have a firm understanding of ocean fluid mechanics yet not to have extensive experience with ocean model methods nor numerical algorithms. Correspondingly, we synthesize intellectual threads that underlie ocean model dynamical cores, with particular emphasis on finite volume formulations, GVCs, and vertical

Finite Volume Equations for Moving Fluid Regions
One can formulate the equations of geophysical fluid mechanics by developing partial differential equations describing the motion of tiny ("infinitesimal") elements of seawater, and we summarize the physical content of that description in Appendix A. However, differential equations do not provide our desired starting point for a finite volume formulation of discrete ocean model equations. We instead formulate integral budget equations for arbitrary closed moving fluid regions of finite extent, with our presentation inspired by similar treatments given by Donea et al. (2004) and Ringler (2011). The resulting finite volume budgets for scalars (e.g., seawater mass, tracer mass, and potential enthalpy) and vectors (e.g., linear momentum and vorticity) offer a weak formulation of the ocean fluid mechanical equations, whereas the corresponding partial differential equations provide a strong formulation. The terminology refers to the need to make more assumptions about smoothness of the flow state when working with the partial differential equations (i.e., "strong assumptions"), whereas the finite volume integral equations generally allow for nonsmooth flow features (i.e., "weak assumptions"). For many purposes of geophysical fluid mechanics, one can seamlessly move between the weak and strong formulations via the Leibniz-Reynolds transport theorem as detailed in These equations constitute the weak form of the scalar budget equations for a finite region of fluid, with their connection to the strong form provided in Appendix C. The left-hand side of Equation 1 expresses the time tendency for region-integrated mass of a scalar whose mass concentration (mass of scalar per mass of seawater) is C. Examples include salinity, S, Conservative Temperature times the constant heat capacity, c p Θ, and biogeochemical trace species. Likewise, the left-hand side of Equation 2 is the time tendency for the total seawater mass in the region, with the in situ density of seawater. The right-hand side to these equations is the surface integral for the transport of these scalar properties through the region boundary. The first transport term arises from the difference between the fluid velocity, v, and the velocity of the boundary, v (b) , with (v − v (b) ) ·n d referred to as the dia-surface transport of seawater mass (dimensions of mass per time) across a surface (see section D0.4 for more details), and C (v−v (b) )·n d the corresponding dia-surface transport of tracer mass. The dia-surface transport is projected onto the outward normal at the boundary so that the scalar budgets are unaffected by transport tangential to the boundary. The second transport term in the tracer Equation 1 arises from the SGS flux across the region boundary as parameterized by J. Notably, this flux generally incorporates surface and bottom boundary fluxes via Neumann boundary conditions placed on J.
There is no subgrid flux in the seawater mass Equation 2 since the mass of a fluid element moves according to the (barycentric) seawater velocity (see section A0.4 for more details). That is, the seawater fluid velocity embodies the mass weighted motion of matter constituents within an element of seawater. Correspondingly, setting the region boundary velocity so that its normal component matches that of the seawater velocity, v (b) ·n = v ·n ⇐⇒ Lagrangian region, yields the scalar budgets for a Lagrangian region, (v), with the corresponding budget equation referred to as the Reynolds transport theorem The budget Equation 5 means that the seawater mass contained in a Lagrangian region is fixed in time, whereas Equation 4 means that the tracer and potential enthalpy content of a Lagrangian region are altered only via fluxes crossing the region boundary as embodied by J.

Budget Equation for Linear Momentum
The linear momentum of a fluid region, ∫  v dV, is modified through body forces integrated over the region, plus contact forces integrated over the region boundary d dt where f body is the body force per mass acting at a point in the fluid and is the contact force per area acting on the boundary to the region. We are concerned with body forces from the Coriolis force and the effective gravitational force (central gravity plus planetary centrifugal) with the angular velocity of the rotating earth and Φ the geopotential. Contact forces arise from pressure, p (an isotropic compressive force), kinetic stresses from momentum transport, and friction due to molecular viscosity, Reynolds stresses, and boundary stresses (e.g., from winds and bottom drag) with I the identity tensor, T the frictional stress tensor (a second-order symmetric tensor), and ⊗ denoting an outer product (see Equation A12 for more details).
Gauss's divergence theorem supports a dual accounting of the net effects from contact forces acting on a region. Namely, they can be formulated as the area integral of the contact forces acting along the region boundaries or as the volume integral of the divergence of stresses within the region volume. Mathematically, this repackaging is given by which in turn allows for the momentum budget to take the form d dt in which we introduced the frictional force per volume As for scalar fields we can specialize to a Lagrangian region in which the linear momentum is not affected by the kinetic stress d dt

Finite Volume Ocean Model Grid Cell Budgets
We here specialize the finite volume scalar and linear momentum budgets from section 2 to develop budget equations for numerical ocean model grid cells such as those shown in Figure 2. The grid cells have top and bottom interfaces defined by surfaces of GVCs denoted by s = s(x, , z, t), with discrete values labeled s k ± 1/2 and with k an integer increasing downward. Along with interior transport of seawater between grid cells, we must also account for the transport of matter, enthalpy, tracers, and momentum across the surface and bottom boundaries of the ocean domain.
In these finite volume budgets, by virtue of the weak formulation, the grid is sufficiently described by the volume boundary positions alone, and there are no properties associated with a point at the center of a cell. Thus, in finite volume models, there is a general emphasis on reconstruction methods for calculating transports across surfaces rather than on finite difference methods for gradients between points.

Framing the Budgets
In the ocean models of concern here, the lateral grid cell boundaries have a velocity, v (b) , that vanishes. Consequently, the corresponding dia-surface transport through the cell boundary (also called "cell face") arises just from the fluid velocity, v. If a grid cell boundary abuts the solid-earth bottom, then the no-normal flow condition, prevents any resolved fluid flow from crossing that boundary. Hence, the only contribution to budgets from the bottom arises from SGS processes and/or boundary fluxes such as bottom frictional stresses and geothermal heating.
Lateral boundaries for interior grid cells are vertical so that the corresponding dia-surface mass transport arises from the horizontal seawater velocity, u = (u, v, 0). For cell boundaries defined by a constant s-surface, we make use of the GVC dia-surface transport detailed in section D0.4. In that treatment, we define w ( . as the dia-surface volume transport of seawater crossing a constant s-surface, per unit horizontal area. The dia-surface transport w ( . s) serves as a "vertical velocity" in that multiplication by horizontal area element yields the instantaneous transport of seawater crossing the undulating s-surface. Following details in section D0.5, we introduce a corresponding dia-surface SGS transport to account for unresolved processes such as dia-surface diffusion that cause tracers and momentum to cross an s-surface.
At the surface ocean boundary, dia-surface mass transport arises from precipitation, evaporation, and surface imposed river runoff (some models impose river runoff at depth as well). Dia-surface transport of tracer and momentum across the ocean surface arises from air-sea and ice-sea boundary fluxes, including turbulent and radiative heating (with shortwave radiation penetrating below the surface layer), turbulent stresses from winds, transport of enthalpy and tracers within the matter crossing the boundary (i.e., precipitation, evaporation, and runoff carry enthalpy, momentum, and tracer content), and other boundary tracer fluxes.

Scalar Budgets
For a grid cell that does not abut the ocean bottom, the finite volume budget for a scalar, given in general by Equation 1, takes on the following form for an ocean model grid cell d dt In this equation we introduced the advective plus subgrid flux components with the flux C w ( . s) + J ( . s) the dia-surface flux that crosses the top and bottom interfaces to the grid cell. We also introduced the difference operators used to compute the difference of fluxes across the cells where i, j, k are discrete cell indices and their half-values denote quantities estimated on cell faces. The swapped signs for the vertical index, k, arise since k increases downward whereas the geopotential coordinate, z, increases upward (see Figure 2). Generalization to a grid cell that abuts the ocean bottom follows

Journal of Advances in Modeling Earth Systems
10.1029/2019MS001954 by adjusting the corresponding boundary flux, with numerical models accounting for boundaries through the use of land-sea masking. The accuracy of a numerical realization of the budget (14) is a function of the accuracy used to estimate the flux components along the grid cell faces.

Cell Averaged Fields
The finite volume budgets developed thus far consist of continuous fields integrated over a cell volume or cell face. A corresponding finite volume numerical model provides a discrete representation of the grid cell averages, with fluxes computed at the cell faces estimated through methods that utilize the cell averaged values. In this subsection we introduce grid cell averaged fields with a focus on issues related to the vertical direction, though note that similar issues arise for the horizontal.
We introduce cell averaged fields through the following cell volume, seawater mass, tracer mass, and linear momentum Now assume the horizontal area of the grid cell, is static so that the grid cell volume is given by in which case The cell averaged tracer content is modified by the convergence of fluxes integrated over the area of the cell faces. We thus introduce the area averaged flux components The horizontal flux components are thickness weighted, with the thickness weighting reflecting the general case of a grid cell with a distinct thickness on each of its lateral cell faces (e.g., Figure 2). In a discrete model, the area averaged fluxes at cell faces are estimated via algorithms that make use of the cell averaged fields. There are many algorithms for this estimation with details outside of our scope.

Budget for Cell Averaged Scalars
Making use of the averages defined in section 3.3 in the finite volume scalar budget (14) renders the discrete cell averaged equation ] .
The s-subscript on the left-hand side reminds us that the partial time derivative is computed while holding the horizontal position fixed yet while vertically following a grid cell defined by surfaces of constant s. That is, the time derivative is a mix of Eulerian (fixed horizontal positions) and Lagrangian (moving vertical positions). Likewise, the horizontal derivatives are computed with s held fixed since they are computed layer-wise. We further explore the mathematics of such GVC derivative operators in Appendix D. The corresponding seawater mass budget follows by setting the tracer concentration to a spatial constant, in which case the SGS fluxes all vanish so that

Budget for Cell Averaged Linear Momentum
Following from the discussion in section 2.2, we write the grid cell budget for linear momentum in the form d dt where we chose to write the pressure and friction contributions in their volume integral form whereas advection is written as an area integral of the kinetic stresses acting on the grid cell boundaries. Now assume the geopotential takes on the form with g a constant gravitational acceleration, and assume a hydrostatic balance whereby the vertical pressure force is balanced by the weight of fluid in the cell For this case the horizontal momentum budget for a hydrostatic fluid is d dt with ∇ z p the horizontal pressure gradient and F horiz the horizontal SGS force per volume. In accord with the hydrostatic balance, we made the traditional approximation in Equation 27 by keeping just the local vertical component to the planetary rotation with the latitude and f the Coriolis parameter. We also wrote the velocity vector as where u =x u +ŷ v is the horizontal velocity.
Now introduce the grid cell averaged horizontal pressure gradient and friction as well as the area averaged advective fluxes along the cell faces

Journal of Advances in Modeling Earth Systems
Substituting into the horizontal momentum Equation 27 for a hydrostatic fluid renders the grid cell budget ] . (31)

The Vertical Lagrangian-Remap Method
The finite volume GVC ocean equations developed in section 3 can be time stepped using a variety of numerical algorithms. We are interested here in three broad categories based on (i) those used in geopotential models (quasi-Eulerian); (ii) those inspired by isopycnal models (vertical Lagrangian) (Adcroft & Hallberg, 2006;Bleck, 2002); and (iii) those used in vertical ALE models as proposed by Leclair and Madec (2011) and Petersen et al. (2015). For example, the MITgcm of Adcroft and Campin (2004), the ROMS code of Shchepetkin and McWilliams (2005), the NEMO code of Madec (2008), and the MOM5 code of Griffies (2012) are formulated as GVC ocean models that use quasi-Eulerian numerical methods whereby the dia-surface transport is diagnosed based on mass continuity (see section 5.2). In contrast, vertical Lagrangian-remap methods as used by HYCOM (Bleck, 2002) and MOM6 (Adcroft et al., 2019) have direct control over the motion of coordinate surfaces, as does the vertical ALE method of Leclair and Madec (2011) and Petersen et al. (2015) though using distinct approaches. In this section we offer a conceptual description of the vertical Lagrangian-remap method, and then in section 5 we compare numerical realizations of quasi-Eulerian, vertical ALE, and vertical Lagrangian remapping.

Essence of ALE and Vertical Lagrangian Remapping
The finite volume budgets from section 3 form the theoretical basis for the ALE method as reviewed by Donea et al. (2004), as well as Lagrangian-remap methods. Both methods consider a grid mesh that moves as the fluid evolves (see Figure 1), with Lagrangian methods allowing the grid to follow the flow whereas ALE typically prescribes a distinction between fluid and grid motion. In essence, the methods are differentiated by specifications of the velocity, v (b) , here measuring the velocity of the grid lines. A static grid (v (b) = 0) represents a trivial limit of ALE in which the grid is Eulerian. A moving grid is Lagrangian by setting v (b) ·n = v·n, whereas ALE grids are typically prescribed otherwise.

Concerning the Choice of Grids
A key constraint on the utility of grids is that they remain nonsingular in order to uniquely estimate the fluid state: The surfaces of constant GVC should never go vertical, should never cross over other coordinate surfaces, nor should they overturn. Furthermore, we generally wish to utilize grid surfaces in a "smart" manner that allows for accurate estimates of fluid properties (e.g., density, tracers, and velocity). These concerns are central to the problem of designing grids for realistic geophysical fluid motion (e.g., boundary currents, jets, fluid instabilities, and turbulent motions) with inhomogeneous properties (e.g., vast range of density structure between the tropics vs. high latitudes). One approach that follows the vertical ALE method was proposed by Leclair and Madec (2011) and Petersen et al. (2015). Their method (detailed in section 5) retains a suitable gridding by damping their moving grid back to a low-pass filtered state that remains nonsingular. For this approach the first arrow in Figure 1 is not a Lagrangian step and the second arrow (regridding) is absent. s) = 0, as well as evolution of the ocean state. The Lagrangian step is assumed to be induced here by an adiabatic gravity wave so that the ocean state experiences no diffusive mixing. Values from above and below the view have advected in vertically. The third panel shows the regrid step, whereby we identify a target s-grid (solid horizontal lines) according to prescribed target geopotentials (which are here identical to the initial geopotentials shown in the first panel). The fourth panel shows the remap step, whereby the ocean state, originally known at the old grid, has now been estimated at the vertical position of the target grid via an interpolation algorithm. Under a perfect remapping the ocean state does not evolve. Rather, remapping adjusts the location of the s-coordinate surfaces and offers a new estimate of the ocean state on the new grid. However, due to truncation errors arising from finite resolution, the representation of the state on the target grid is not a perfect transformation between the two grids. The final panel shows the new state as represented on the new grid, now ready to move forward another time step.

Illustrating the Method
We illustrate the vertical Lagrangian-remap method in Figure 3, which is a special case of the three-dimensional Lagrangian-remap method depicted in Figure 1. Here, the first step is Lagrangian just in the vertical. With a freely moving grid, the grid must be periodically reinitialized or regularized onto a target grid (this is the regrid step) to ensure it remains nonsingular and/or continues to offer an accurate representation of the fluid flow (e.g., by not losing too many grid layers as can occur in isopycnal models). A regrid step necessitates a corresponding remap step whereby the fluid state is estimated on the new grid. The remap updates the discrete representation of the ocean state by transforming (via interpolation) that representation onto the target grid. In principle it does not change the ocean state, whereas in practice there are changes due to interpolation errors. In section 4.3, we discuss the remaining steps to the algorithm illustrated in Figure 3. Dukowicz and Baumgardner (2000) applied the notions of Lagrangian remapping to horizontal transport in sea ice models. They emphasized that remapping is operationally identical to advective transport. Furthermore, the method does not require regridding/remapping every time step. Instead, the regrid/remap step can occur after a few Lagrangian steps, particularly if the flow is laminar and/or the time steps are small.

Example Ocean and Sea Ice Realizations
Vertical Lagrangian remapping is a natural extension of methods used for isopycnal layered ocean models (see Bleck, 1998, for a review). Correspondingly, the numerical realization of the vertical Lagrangian-remap method by Bleck (2002) facilitates a remapping of the ocean state over an arbitrary number of vertically displaced grid cells. This capability makes use of numerical technology developed to allow coordinate layers to vanish or inflate, which is an essential feature of isopycnal coordinate models used for realistic simulations. In turn, these algorithms lead to a natural realization of conservative wetting and drying with corresponding applications in studies of estuaries and moving ice-shelf grounding lines (e.g., Goldberg et al., 2012aGoldberg et al., , 2012b. Finally, we note that Leclair and Madec (2012b) and Petersen et al. (2015) present an implementation of vertical ALE within two ocean models. Their method is described in section 5.

The Hydrostatic Ocean Equations
In presenting the vertical Lagrangian-remap method, we work with the finite volume grid cell budgets for a hydrostatic fluid using GVCs. We derived these budget equations in section 3, where the equations integrate the cell averaged fields and we summarize the equations in the form

Journal of Advances in Modeling Earth Systems
Note that we reduced notational clutter by dropping the "cell" script on fields that are defined as cell averages (section 3.3), and yet we maintain the overbar on fluxes to emphasize that they are estimated on cell faces using the grid cell averaged fields. We also introduced a shorthand notation for the horizontal divergence computed in constant s-layers The dia-surface volume flux, w ( . s) (red term in the prognostic Equations 32a-32d), measures the flux of seawater that crosses the moving GVC surface (see section D0.4 for details). Hence, if the normal component to the s-surface velocity matches that of seawater fluid elements, then there is zero seawater that flows across the s-coordinate surface. Time stepping the ocean equations with w ( . s) = 0 thus constitutes a vertical Lagrangian update to the ocean state since the s-coordinate surfaces follow the flow when w ( . s) = 0. In this manner we see how the GVC equations are ideally suited for a vertical Lagrangian time stepping of the ocean equations. We expand on this point next.

The Vertical Lagrangian Step
The first step of the vertical Lagrangian-remap method (second panel in Figure 3) results from having the generalized vertical coordinate, s, follow the vertical motion of fluid elements so that w ( . s) = 0. A surface maintaining w ( . s) = 0 means that no net seawater mass crosses the surface. If the fluid is homogeneous, then w ( . s) = 0 also means the surface is material; that is, no matter crosses it. However, we are concerned with real fluids containing multiple constituents and irreversible mixing. Mixing of fluid constituents acts to exchange matter between fluid elements while fixing the mass (or volume for an incompressible Boussinesq fluid) of each fluid element. Hence, the vertical Lagrangian step allows the s-surface to follow the motion of a seawater fluid element (the barycentric velocity of section A0.4) even as matter and momentum are irreversibly exchanged across the surface. Operationally, the Lagrangian step is realized by time stepping the fluid state with w ( . s) = 0 in the ocean Equations 32a-32f. To allow matter (tracers) and momentum to irreversibly cross the s-surfaces, we retain diffusion and friction in the tracer and momentum equations. This is also the step where boundary fluxes are incorporated to the fluid.

The Vertical Regrid Step
As emphasized above, the Lagrangian step allows for both reversible and irreversible changes to the ocean state, all while having the s-coordinate surfaces vertically follow fluid elements. If continued indefinitely, the s-surfaces will generally drift far from their initial positions. If the grid were moving in three dimensions, it would become corrugated and filamentary, especially in the presence of turbulent motion. However, so long as the horizontal time stepping of transport is stable to thus ensure layer thicknesses remain positive, then the vertical Lagrangian method is more forgiving because all that can happen is the vertical resolution becomes heterogeneous. To ensure an accurate representation of the discrete model equations, it is necessary to include a second step that we call the vertical regrid step (also referred to as coordinate reinitialization or regularization), with this step depicted in the third panel of Figure 3.
The vertical regrid step occurs by laying down a new s-grid according to a prescribed regridding method associated with a target coordinate layout. There is great freedom in choosing this target grid. For example, we may choose a target grid according to surfaces of constant geopotential placed at prescribed vertical positions, as illustrated in the third panel of Figure 3. Any sensible alternative is available as well, such as prescribed potential density classes or a mixture of geopotentials, isopycnals, and terrain-following surfaces as in Bleck (2002). Indeed, as noted by Bleck (2002), there is no need to provide an analytical expression for the target grid; one only needs an algorithmic prescription. Furthermore, the vertical resolution of the target grid can differ between time steps. In this manner the vertical Lagrangian-remap method can be considered "coordinate free."

The Vertical Remap Step
The ocean state is known on the grid that moves with the flow under the Lagrangian step. Upon regridding, the ocean state must be estimated on the target grid via interpolation (or extrapolation if needed near boundaries). This estimation is known as the vertical remap step whereby we map the ocean state from its Lagrangian-displaced grid onto the target grid. Every regrid step requires a remap step. The remap step is depicted in the fourth panel of Figure 3. Under a perfect realization of remapping, the ocean state does not evolve since no fluid elements move nor does any trace matter. The position of the s-grid changes during regridding, thus necessitating a new representation of the ocean state on the new grid. Furthermore, as seen in section 5.4, vertical remapping corresponds operationally to dia-surface advection.
In practice, there is spurious evolution during the remapping step since a discrete remapping incurs discretization error. We see such errors in the final panel of Figure 3 whereby the fluid state has been fully remapped onto the new grid in preparation for the next time step. As they degrade the integrity of the simulated ocean state, such numerical errors also affect the mechanical energy by modifying the background potential energy, much like errors in traditional advection schemes (Ilicak et al., 2012). Reducing these numerical errors in vertical remapping, while maintaining scalar conservation (i.e., integrated scalar properties remain unchanged by the remapping algorithm), is a primary concern for vertical remapping methods such as those developed by White and Adcroft (2008). We provide more comments on this point in section 6.

Some Examples
We illustrate the vertical Lagrangian-remap method by considering some examples for a Boussinesq fluid with a linear equation of state where buoyancy is proportional to Conservative Temperature. Although trivial oceanographically, these examples help to conceptually establish how the method works in practice. 4.6.1. Regridding/Remapping to Geopotentials As already noted, Figure 3 illustrates the vertical Lagrangian-remap method for target geopotential vertical coordinates where the s-coordinate surfaces are initially horizontal. Furthermore, we assume that they are initially aligned with horizontal isopycnals. The Lagrangian step allows for the coordinate surfaces to evolve while maintaining fixed volume within the coordinate layers. It also allows for mixing to act on the tracer and velocity fields. The figure illustrates the Lagrangian step induced by a gravity wave in the absence of mixing. The wave causes the s-coordinate surfaces to undulate along with the isopycnals. Since the s-surfaces are displaced by the wave, the regridding step reinitializes the s-coordinate surfaces to their target geopotentials. The remapping step estimates the ocean state at the target positions thus completing the time stepping of the ocean state. Note that if the target grid was isopycnal, then there would be no regrid/remap step since the s-coordinate surfaces remain fixed to their original isopycnal surfaces in the absence of mixing or boundary buoyancy forcing.

Regridding/Remapping to Isopycnals
Now consider the case of a wave in the presence of vertical diffusion and assume an isopycnal target grid. As before, we assume that the Lagrangian step is induced by a wave that causes the s-coordinate surfaces to undulate just like the second panel of Figure 3. Additionally, vertical mixing causes the isopycnals to spread. Yet the coordinate surfaces do not feel diffusion since they evolve with w ( . s) = 0. Hence, the isopycnals and s-coordinate surfaces separate due to mixing in the fluid. The regrid step returns the s-coordinate surfaces to their target isopycnals, and the remap step estimates the state within the new isopycnal positions (see section 6.1 for more on the challenges of finding robust positions for isopycnals).

Zero Fluid Particle Motion With Nonzero Mixing
To further emphasize the essence of the method, drop the wave motion yet retain vertical fluid mixing that acts on vertically stratified horizontal isopycnals. There is no motion of fluid elements in this case, and yet vertical mixing causes the isopycnals to move vertically. The Lagrangian step retains s-isolines aligned with their initial geopotentials since there is no fluid motion. In the case of a geopotential target grid, there is no need to regrid or to remap since there is no fluid motion. All that happens is that density spreads relative to the static geopotentials. The complement holds for an isopycnal target grid whereby vertical mixing causes the s-coordinates to lose their initial potential density, thus requiring a regrid step to return the s-grid to its target isopycnal grid. A regrid step in turn requires a remap step to estimate the ocean state on the target grid.
Exposing a bit of the maths, we set w = 0 in Equation D18 (equation relating vertical fluid motion to motion of the grid and dia-surface transfer), so that the vertical position, z, of a potential density surface, , evolves according to is the buoyancy relative to a reference temperature Θ 0 , the constant thermal expansion coefficient, and N 2 = b∕ z the squared buoyancy frequency. Vertical diffusion causes a nonzero w ( . ) so that the isopycnal surfaces move during the Lagrangian step (which includes diapycnal diffusion in the tracer equation). In contrast, the s-surfaces remain fixed during the Lagrangian step since there is no change in the s-layer thicknesses. It is only during the regridding step that the s-surfaces are brought back to their target -surfaces. They are brought back by moving the s-surfaces vertically by an amount determined by ∫ Δt 0 w ( . ) dt, with Δt the time step between regridding.

Summarizing Some Conceptual Points
When allowing the grid to move as in the vertical Lagrangian-remap method, motion of fluid elements and thus the transport of fluid properties are measured relative to the moving grid. This transport is the dia-surface transport detailed mathematically in section D0.4. If the grid lines are fixed, as in an Eulerian method, then advective transport arises from fluid moving across the fixed grid lines. For the vertical Lagrangian-remap method the grid lines move vertically with the fluid during the Lagrangian step. Vertical advective transport arises only during the remap step as part of estimating the fluid state on the new grid. If there is no need to vertically remap, then there is no vertical advection (e.g., adiabatic flow with an isopycnal target grid). In short, an Eulerian method watches fluid move through fixed grids whereas the vertical Lagrangian method first follows the flow and then (if needed) vertically regrids relative to the fluid and remaps onto the target grid.
We emphasize that during the Lagrangian step the fluid experiences the full suite of mixing and boundary fluxes that lead to irreversible evolution of the ocean state. So although the vertical position of the grid lines moves with fluid particles during the Lagrangian step (as if there was no mixing), the fluid state evolves in the presence of mixing and boundary fluxes. Hence, the vertical Lagrangian step is not generally a reversible step.

Comparison of Numerical Algorithms
In this section we expose certain algorithm details required to realize the vertical Lagrangian-remap method in a numerical model. In the process we compare this algorithm to a quasi-Eulerian algorithm and a vertical ALE algorithm. The quasi-Eulerian algorithm (see Chapter 12 of Griffies, 2004 for a summary) follows the general strategy of time stepping used in traditional geopotential coordinate models, with "quasi" used since the algorithm is applicable to some GVCs so that grid cells can have moving top and bottom interfaces. The vertical ALE algorithm we feature is based on that of Leclair and Madec (2011) and Petersen et al. (2015).

Focusing on the Thickness and Tracer Equations
Our aim is to sketch elements of the three distinct numerical algorithms and to compare them side by side. Our primary concern with questions of numerical stability is to note the need to resolve cross-grid vertical motions in time according to a vertical CFL condition with the quasi-Eulerian and vertical ALE methods, but not with the vertical Lagrangian-remapping method. Additionally, we are not concerned with numerical accuracy, noting that while the chosen forward Euler scheme is only first order accurate in time overall and might be insufficient for certain applications, this limitation can be compensated within the spatial discretization. Table 1 Summary of the Three Numerical Algorithms Detailed in Sections 5.2,5.3,and 5.4 Quasi-Eulerian: Thickness tendency determined by an analytic vertical coordinate Horz advection thickness update Vert advection thickness update Vert advection tracer update Vertical ALE: General thickness tendency with a general vertical coordinate Horz advection thickness update Vert advection thickness update Vert advection tracer update Vertical Lagrangian remap: Coordinate free Layer motion via convergence of horz advection Remap tracer using dia-surface transport Note. For brevity we consider the Boussinesq equations for layer thickness and thickness-weighted tracer concentration. The quasi-Eulerian (QE) algorithm is a generalization of the traditional Eulerian algorithms originating with Bryan (1969), with the vertical grid defined analytically. The vertical ALE (V-ALE) algorithm allows for more general grid and does not require an analytic expression for the vertical grid. One example is the algorithm of Leclair and Madec (2011) and Petersen et al. (2015). The vertical Lagrangian-remap (V-LR) algorithm is based on the method detailed in section 4 and is generally coordinate free. Each algorithm has been written using operator splitting between vertical and horizontal to better compare.
To focus on the essential similarities and differences between the algorithms, we emphasize the distinct treatment by the algorithms of w ( . s) , the dia-surface velocity, and the grid velocity, w grid , with w grid determining how the layer thicknesses evolve and hence how the coordinate interfaces move. For this purpose it is sufficient to consider just the prognostic equations for layer thickness and thickness-weighted tracer concentration in a hydrostatic Boussinesq ocean; to ignore all SGS terms and boundary fluxes; and to make use of a forward Euler time stepping scheme from time step (n) to time step (n + 1). Further terms such as those from the momentum equation add complexity important for the full ocean equations yet without providing new conceptual pieces to the algorithm discussion.
Hence, we consider the thickness equation and thickness-weighted tracer equation in the continuous form where C is a generic tracer concentration (e.g., salinity, Conservative Temperature, and biogeochemical tracer), and the updated tracer concentration is obtained by The thickness-weighted velocity is updated as per the momentum equation, but again, its consideration offers no new conceptual element to the following discussion of numerical algorithms. 10.1029/2019MS001954

Quasi-Eulerian (QE): Vertical Coordinate Determines h/ t
The quasi-Eulerian method is characterized by a thickness tendency that is directly related to the ocean free surface tendency (for Boussinesq fluids) or bottom pressure tendency (for non-Boussinesq fluids), with this relation determined by the vertical coordinate choice. For example, the rescaled geopotential coordinate (Adcroft & Campin, 2004;Stacey et al., 1995) has its layer thickness and corresponding time tendency In these equations, z = (x, , t) is the ocean free surface and z = −H(x, ) is the bottom (see Figure 2). Likewise, for terrain-following coordinates, with s = z * ∕H, the thickness tendency is h∕ t = ds ( ∕ t). The free surface time tendency for a Boussinesq fluid is given by with U = ∫ −H u dz the depth integrated horizontal flow and Q m the surface mass flux from precipitation, evaporation, runoff, etc. Consequently, if we know − ∇ ·U+Q m / 0 , then we can update both the free surface (Equation 40) and the interior layer thicknesses ( h/ t). We can then diagnose the dia-surface flow by vertically integrating starting from the boundary condition w ( . s) = 0 at the ocean bottom. Note that this bottom kinematic boundary condition holds for arbitrary topography since there is no flow normal to the solid-earth bottom.
To bring the algorithm arising from the quasi-Eulerian method in-line with the two subsequent algorithms, we measure the motion of the grid surfaces via The quantity, w grid , offers a convenient means to identify the distinct approaches used to determine motion of the grid. In the quasi-Eulerian algorithm detailed in this subsection, the grid motion is determined by motion of the vertical coordinate surfaces. This grid motion is more general for the vertical ALE and vertical Lagrangian-remap methods. We emphasize that this grid motion is distinct from motion of seawater crossing the grid surfaces, which is measured by w ( . s) .
Some quasi-Eulerian models perform an incremental update for the tracer concentration by operator splitting the horizontal and vertical advection into the two steps (e.g., Easter, 1993). We present the quasi-Eulerian algorithm using this split since it allows for a clean comparison with the vertical ALE and vertical Lagrangian-remap algorithms discussed in sections 5.3 and 5.4. The following summarizes the quasi-Eulerian algorithm , and likewise for the tracer equation. This sanity check verifies that the operator split algorithm offers a consistent discretization of the thickness and tracer equations. Furthermore, note that setting tracer concentration to a constant reduces the thickness-weighted tracer updates to those for the thickness, with this reduction required to ensure compatibility between updates of the layer thickness and updates of the tracer concentration.

Journal of Advances in Modeling Earth Systems
The key takeaway points for the quasi-Eulerian algorithm are the following.
• It uses a diagnostic calculation for w ( . s) via the continuity Equation 44b. • The time tendency of the layer thickness is known through the free surface Equation 40. More generally, it is known through an analytic specification of the vertical coordinate. • All vertical motion, including gravity waves, must be resolved for advection. Hence, an explicit time realization of vertical advective transport requires the time step to respect the vertical advective CFL condition This condition concerns time-explicit schemes and becomes rate limiting when cell thicknesses become small and/or when dia-surface motion becomes large. Consequently, there have been a variety of algorithms developed to alleviate this time step constraint, such as the Courant number-dependent implicit scheme of Shchepetkin (2015); the flux-form semi-Lagrangian approach of Lin and Rood (1996); and the cell integrated semi-Lagrangian (CISL) scheme of Nair and Machenhauer (2002). As seen in sections 5.3 and 5.4, vertical ALE and vertical Lagrangian remapping address the time step constraint in a complementary manner by aiming to reduce the cross-coordinate flow.

Vertical ALE (V-ALE): Arbitrary h/ t
The quasi-Eulerian method can be generalized to allow the vertical grid to be unconstrained in its time tendency so that it is not necessarily tied to the free surface motion. In this manner, the grid cell thicknesses are no longer defined by an analytic vertical coordinate. Instead, they are prescribed by an algorithm as per a vertical ALE method.
To motivate our realization of vertical ALE, return to the z * coordinate of Equation 37. Observe that if the flow arises solely from free surface motion (i.e., barotropic flow), then there is zero dia-surface transport so that, in this special case, z * is a Lagrangian coordinate. Thez coordinate of Leclair and Madec (2011) takes this idea one step further by absorbing relatively fast motion from internal layer undulations (e.g., isopycnal layers with gravity waves) into the layer motion. They do so by prescribing a method for evolution of the cell thicknesses that preferentially captures the faster motion within the grid and allows the slower motion to be handled much like a traditional z * algorithm.
One means to understand the Leclair and Madec (2011) method is to split the velocity for a fluid element into a slow or low-pass filtered velocity, v (L) , and a fast or high-pass filtered velocity, v (H) , so that Consequently, the dia-surface velocity in Equation (D14) takes the form where v (s) is the velocity of a point on the s-coordinate surface. As before, a purely Lagrangian approach is realized by a grid that maintains w ( . s) = 0, meaning that the grid follows fluid elements,n · v (s) =n · v. An intermediate approach considers v (L) as a low-pass filtered version of the full velocity so that it captures the relatively slow motion. The complement velocity, v (H) = v − v (L) , captures the fast motion such as from internal gravity waves. To minimize numerically induced spurious dia-surface transport associated with the fast motion, Leclair and Madec (2011) prescribe a grid that matches the fast motion so that Hence, dia-surface transport arises just from slower motions presumably simpler to accurately realize using a z or z * grid This method is realized algorithmically by introducing a target layer thickness, h target , determined by an algorithm such as that detailed by Leclair and Madec (2011) and Petersen et al. (2015). This target thickness determines the vertical grid motion via The remainder of this vertical ALE algorithm is identical to quasi-Eulerian given by Equations 44a-44f. The key takeaway points for the vertical ALE algorithm are the following.
• Any chosen fluid motion can be absorbed by the vertical grid. Correspondingly, the grid does not need to be defined analytically. • An explicit time realization of vertical advective transport requires the time step to respect the vertical advective CFL condition |w ( . s) | Δt∕h < 1. However, the grid motion is generally chosen to render |w ( . s) | smaller than with the quasi-Eulerian algorithm, thus allowing for a larger Δt and/or a smaller layer thicknesses. It is only that portion of the motion that is not captured by the grid that must be resolved for tracer and velocity advection. • The vertical ALE algorithm allows for a Lagrangian grid by choosing so that w ( . s) = 0, such as for an adiabatic isopycnal model. • As shown by Leclair and Madec (2011) and Petersen et al. (2015), absorbing fast motion into the grid helps to reduce spurious cross-grid transport in the presence of linear gravity waves, in which case there is a reduced need to perform vertical advection. However, the method does not generally help reduce spurious transport associated with balanced motions such as geostrophic turbulence (e.g., Griffies, Pacanowski, & Hallberg, 2000;Ilicak et al., 2012).

Vertical Lagrangian Remapping (V-LR)
As noted by Adcroft and Hallberg (2006), an algorithmic realization of the vertical Lagrangian-remap method can be written as an operator split algorithm whereby the first group of steps is Lagrangian (here just incorporating horizontal transport) and the second and third groups of steps involve vertical regridding and remapping. We detailed this algorithm in section 4, and it takes on the following form for the present discrete equation set h † = h (n) + Δt Δ s w grid horz advection thickness update, [hC] † = [hC] (n) − Δt ∇ s · [hC u] (n) horz advection tracer update, [hC] (n+1) = [hC] † − Δt Δ s (w ( . s) C † ) remap tracer using dia-surface transport.
Equations 52a-52c constitute the Lagrangian portion of the algorithm, corresponding to the second panel ("Lagrangian") in Figure 3. In particular, Equation 52b  Equation 52d is the regrid step, corresponding to the third panel in Figure 3. This is the key step in the algorithm with the new grid defined by the new thickness h (n + 1) . The new thickness is prescribed by the target values for the vertical grid, h (n+1) = h target . The prescribed target grid thicknesses are then used to diagnose the dia-surface velocity according to This step, and the remaining step (52f) for tracers, constitutes the remapping portion of the algorithm as depicted in the fourth panel of Figure 3. For example, if the prescribed coordinate surfaces are geopotentials, then w ( . s) = w and h target = h (n) , in which case the remap step reduces to Cartesian vertical advection. Furthermore, notice the relative lateness of the diagnosis of w ( . s) for the vertical Lagrangian-remapping algorithm as compared to the quasi-Eulerian and vertical ALE algorithms. The reason is that during the first half of vertical Lagrangian remapping, the equations are time stepped with w ( . s) = 0. It is only for developing the remapping step that we need to diagnose w ( . s) according to Equation 53.
This realization of the vertical Lagrangian-remap method allows for the use of numerical algorithms that remove the vertical CFL constraint. Conceptually, this freedom is afforded since vertical regrid step rearranges the grid without transmitting physical signals. If one designs the remapping algorithm to allow for remapping the ocean state over more than a single grid cell in the vertical, then there are no corresponding vertical CFL constraints for either tracer or velocity. This property is shared by other Lagrangian transport schemes (e.g., Lin, 2004). It is naturally facilitated by the vertical Lagrangian-remap method, but it is not a necessary component to an algorithmic implementation, nor is it an exclusive feature of vertical Lagrangian remapping.
As for any generalized vertical coordinate layer model, it is essential to maintain positive definite layer thicknesses. Positive layer thicknesses can be ensured for time-explicit methods if the horizontal (within layer) CFL constraint is maintained by the horizontal transport. By doing so, the layer thicknesses can properly respond to the horizontal divergence of mass within a layer. This horizontal CFL constraint is shared with the quasi-Eulerian and vertical ALE methods.
The key takeaway points for algorithms based on the vertical Lagrangian-remap method are the following.
• If the numerical remapping algorithm allows for remapping the ocean state over more than a single grid cell in the vertical, as in MOM6 and HYCOM, then there are no corresponding vertical CFL constraints for either tracer or velocity. Positive layer thicknesses are ensured for time-explicit methods if the horizontal (within layer) CFL constraint is maintained by the horizontal transport. • Equations 52a-52f can be time stepped and subcycled in three groups, each with distinct time steps.
1. The Lagrangian dynamics in Equations 52a-52c form the first group, and it typically requires the smallest time step sufficient to resolve gravity waves and inertial oscillations. 2. The Lagrangian tracer update (52c) is the second, which is constrained by horizontal advection and diffusion. 3. The regrid/remap equations 52d-52f form the third group. This update is constrained by drift of the coordinate surfaces and is thus a function of the vertical coordinates and flow regime. The tracer updates are less constrained than the Lagrangian dynamics, thus allowing for a longer tracer time step than dynamics time step while still satisfying the CFL criteria. This feature is of particular use when incorporating many dozens of biogeochemical tracers found in comprehensive earth system and ecosystem models. Bleck (2002) pioneered the vertical Lagrangian-remap method for ocean models and realized it algorithmically within the HYCOM code. It has taken the subsequent years for the broader modeling community to catch up to Bleck's insights and developments and to improve on certain of his original algorithmic details. Given the broad and growing community of ocean modelers making use of HYCOM and MOM6, we suspect that the vertical Lagrangian-remap method will see more use and further algorithmic improvements. Furthermore, we hope that this primer will help to remove some mysteries surrounding the method among the community of model developers, practitioners, and interested spectators. We here close the main part of this paper by presenting some challenges and opportunities forming the focus of ongoing research.

Regridding, Remapping, and Spurious Mixing
The numerical integrity of a model using a vertical Lagrangian-remapping dynamical core relies on two key properties of the simulation: (A) the amount of flow crossing surfaces of constant generalized vertical coordinate and (B) the accuracy of vertical remapping. When either are done poorly, the simulation degrades through spurious mixing (Griffies, Pacanowski, & Hallberg, 2000), with the degradation especially pertinent for climate simulations where errors can accumulate over time. As suggested by Griffies, Böning, et al. (2000), the use of GVCs has emerged as a promising method to reduce such errors. The vertical Lagrangian-remap method is a suitable framework to build hybrid coordinate models to address the spurious mixing problem while maintaining a dynamical core flexible enough to simulate the wide variety of flow regimes and water mass structures found in the World Ocean. We do note that a practical realization remains a work in progress, and we here touch upon some of the issues.

Vertical Coordinate Choices and Regridding
We generally aim to minimize the amount of flow crossing vertical coordinate surfaces for the purpose of reducing the amount of spurious cross-surface transport. Doing so requires a "smart" choice to the vertical coordinate that reduces the need for regridding and the associated remapping. For example, if the physical flow has relatively small levels of diapcynal mixing, such as for stably stratified flows in the ocean interior, then an isopycnal coordinate is preferred over a geopotential coordinate, thus allowing for the regrid/remap step to be largely absent. The real ocean, however, is not so simple. Rather, a simulation of the global ocean encounters a complex suite of stratification regimes that span surface and bottom boundary layers, interior ocean regimes of stable stratification, and low-latitude and high-latitude processes. Simulating these multiple regimes is made even more complicated by the nonlinear equation of state for seawater Hallberg, 2005;IOC et al., 2010;Sun et al., 1999).
Acknowledging the complex nature of the ocean's vertical density stratification, Bleck (2002) designed his "HYCOM" vertical coordinate as a regional hybrid of pressure in the upper ocean, where strong mixing leads to weak if not negative vertical density stratification; potential density referenced to 2,000 dbar ( 2,000 ) in the ocean interior, where density stratification is stable; and terrain following near the bottom to represent the bottom boundary layer. Another approach was proposed by Hofmeister et al. (2010), who introduced a prognostic equation for the vertical coordinate that resulted in a vertical coordinate with nontrivial vertical stratification while ensuring lateral coordinate surface smoothness. Other algorithms have been proposed for atmospheric models by Thuburn (2018a, 2018b). At this stage of research, the choice of a hybrid vertical coordinate for realistic ocean climate simulations remains largely an art rather than a science.

Vertical Remapping and Spurious Mixing
As noted in section 4.5, in principle the vertical remap step does not alter the ocean state. However, in practice numerical errors create spurious evolution. Reducing numerical errors requires both a suitable vertical gridding and an accurate remapping algorithm. The remapping scheme from Bleck (2018a) was rather low order and thus led to a nontrivial amount of spurious mixing (Megann et al., 2010). However, White and Adcroft (2008) and White et al. (2009) developed a higher-order conservative scheme that greatly ameliorated the problems with Bleck (2002).
In general, remapping errors degrade buoyancy gradients that in turn affect the mechanical energy, most notably the potential energy. Ilicak et al. (2012) measured changes in potential energy and isolated that portion of the change from spurious mixing. Gibson et al. (2017) went one step further using MOM6 by partitioning the spurious mixing into that arising from lateral processes (e.g., along coordinate advection) and that portion arising from vertical remapping. This ability to separately attribute causes for spurious 10.1029/2019MS001954 mixing is natural in models based on vertical Lagrangian remapping given the operator split between lateral processes and vertical remapping.
Although we are not ready to conclude that the spurious mixing problem has been solved, we suggest that the vertical Lagrangian-remap method provides a robust and flexible framework for enabling "smart" vertical coordinates and highly accurate interpolation schemes, thus allowing for the reduction of spurious mixing to physically acceptable levels. We thus conjecture that with further advances in vertical coordinate construction, the vertical Lagrangian-remap method will facilitate the removal of spurious mixing from the suite of numerical biases that plague realistic ocean climate models. Adcroft et al. (2019) offer evidence to support this optimistic conjecture.

Conservative Wetting-Drying
The key algorithmic feature required for conservative wetting-drying, as utilized in HYCOM and MOM6, is to allow the vertical remapping (or vertical advection) to span an arbitrary number of vertical coordinate layers. In so doing the vertical CFL constraint is removed, thus allowing for vanishing/inflating grid layers and the conservation of scalar properties to machine precision. HYCOM and MOM6 arrived at this technology through implementing numerical algorithms historically developed for vanishing and inflating layers in isopycnal models. Alternatively, one can make use of general vertical remapping/advection algorithms within models based on quasi-Eulerian or vertical ALE methods. Machine precision is required for century to millennial scale climate and earth system applications, in which conservation of mass, potential enthalpy, salt, and biogeochemical tracers is essential. Such precision is ensured by a conservative remapping/advection scheme that allows for vanishing and inflating grid layers.
One application of conservative wetting and drying algorithms concerns the dynamical modification of the land-sea boundary through vanishing and reappearing vertical coordinate layers. Models with such capabilities enable mechanistic studies of interactive and coupled ocean/ice-shelf processes whereby the ice-shelf grounding line evolves (Goldberg et al., 2012a(Goldberg et al., , 2012b. This system is critical to robustly simulate in global climate models aiming to address questions concerning sea level rise from melting land ice. It also facilitates the study of the transition between glacial and interglacial periods whereby massive changes to the land-sea boundaries occurred. Finally, conservative wetting-drying algorithms support the study of tidal fluctuations in shallow seas and coastal estuaries, thus enabling a wealth of applications for coastal oceanography.

Nonhydrostatic Ocean Dynamical Cores
As noted by Hallberg (2006, 2012b), the vertical Lagrangian-remap method has only been realized for ocean circulation models making the hydrostatic approximation. However, there is nothing fundamental that precludes using the method for nonhydrostatic modeling. Supporting evidence can be found in the nonhydrostatic applications reviewed by Donea et al. (2004); through its use for the nonhydrostatic atmospheric modeling by Chen et al. (2013) and Kavćić and Thuburn (2018a); by the nonhydrostatic isopycnal ocean model for simulating nonlinear internal solitary gravity waves by Vitousek and Fringer (2014); and by the layered nonhydrostatic approach of Popinet (2020). The key condition is that one must enable regridding at a sufficiently fine time resolution to maintain monotonic vertical grid layering and thus to prevent grid singularities. We thus conjecture that it is only a matter of time before we see a vertical Lagrangian-remap-based nonhydrostatic ocean dynamical core for realistic global and regional applications.

Appendix A: Differential Equations for Fluid Elements
In this appendix we present the fluid mechanical equations of the ocean as expressed by partial differential equations for an infinitesimal fluid element. These equations represent the strong formulation of the ocean equations, with this discussion complementing the finite volume weak formulation from section 2. Here, we conceive of a seawater fluid element as a constant mass region with length scale at or smaller than L micro ≈ 10 −3 m. This microscale is where molecular viscous forces play an important role in the momentum equation. That is, the Reynolds number based on molecular viscosity is order unity at the microscale. Hence, the equations described in this section are formally appropriate for motions extending down to the microscale. However, we expose certain features of the subgrid scale (SGS) operators required for coarse-grained equations used for most simulations. Although relying on textbooks (e.g., Olbers et al., 2012;Vallis, 2017) for derivation of these equations, we here describe their physical content with an eye toward numerical models.

A1. Partial Differential Equations for the Ocean Fluid
A seawater fluid element maintains a constant mass, yet its matter, momentum, and enthalpy are exchanged with other fluid elements through both reversible and irreversible processes. Partial differential equations describing the evolution of fluid elements are based on Newton's laws of motion, thermodynamics, and mass conservation. Applying Newton's second law to a fluid element leads to the equation for linear momentum, v V, where V is the volume of the fluid element, its mass density, and v its velocity. The linear momentum for seawater moving on a rotating and graviting earth is affected by body forces (gravity, centrifugal, and Coriolis) and contact forces (pressure and friction). Conservation of mass for a fluid element, V, leads to the continuity equation, and mass conservation for trace constituents as exemplified by the mass concentration/fraction of salt, S, leads to the salinity equation. The evolution of potential enthalpy per mass, c p Θ, occurs under the assumption of local thermodynamic equilibrium and is realized by a prognostic equation for Conservative Temperature, Θ, with c p the constant specific heat capacity (McDougall, 2003). Finally, we make use of an empirically determined equation of state for seawater mass density, , with IOC et al. (2010) providing the state-of-the-science version (see Roquet et al., 2015, for an efficient numerical implementation).
In mathematical form, the coupled nonlinear partial differential equations for the ocean can be written The Lagrangian time derivative, D/Dt, is related to the Eulerian (or laboratory frame) time derivative via with v·∇ the advection operator and v the barycentric (center of mass) velocity for the fluid element.

A2. Linear Momentum Equation
The linear momentum of a fluid element is its velocity times its mass, v ( V), and its time change is given by D(v V)∕Dt = ( V) Dv∕Dt, which follows since we are interested in the dynamics of constant mass fluid elements so that D( V)∕Dt = 0 (which leads to the continuity Equation A2 discussed in section A3). Through Newton's law of motion, we equate time changes in the linear momentum to forces acting on the fluid element. There are two types of forces acting on a continuous media: external/body forces and internal/contact forces. Body forces act throughout the volume of a fluid element whereas contact forces arise from the mechanical interactions between fluid elements and so act on fluid element boundaries.

A2.1. Body Forces From Coriolis and Geopotential (Gravity + Centrifugal)
Observing fluid motion in the rotating planetary reference frame requires body forces from Coriolis and the geopotential. The Coriolis force per volume is written −2 ∧ v, with the angular velocity of the planet directed along the rotational axis through the north pole. We assume that is constant in time, which is a very accurate assumption for the time scales of order 10 4 years or less commonly considered in ocean climate studies. For example, adding/removing ice sheets changes | | by roughly one part in 10 5 .
The effective gravitational force per volume is given by − ∇Φ, where the geopotential Φ combines accelerations from planetary gravitational (gravitational acceleration between fluid element and the earth) and planetary centrifugal (centrifugal acceleration due to motion on the rotating planet). In many applications we set Φ = gz so that − ∇Φ = − gẑ. In these expressions, z is the geopotential vertical position (z = 0 defines the geoid as discussed in Gregory et al., 2019),ẑ is the local vertical unit vector, and g is a constant effective gravitational acceleration commonly taken to be roughly g = 9.8ms −2 . The geopotential becomes a space-and time-dependent field when simulating ocean tides (e.g., Arbic et al., 2018) and when rotational and deformation properties of the planet change due to mass redistributions such as that occurring from ice sheet melt (e.g., Gregory et al., 2019).

A2.2. Contact Forces From Pressure, Viscosity, and Motion
Three contact forces concern us. The first arises from pressure, p, which acts in a compressive manner thus leading to the force wheren  is the oriented surface area on the boundary of the fluid region, withn the outward normal. The second expression exposed the Cartesian tensor indices a = 1, 2, 3. We assume local thermodynamic equilibrium so that the mechanical pressure, p, in a moving fluid is the same as the thermodynamic pressure used in the equation of state (e.g., section 4.5 of Kundu et al., 2012). The exception to this assumption occurs for Boussinesq fluids, in which we only need the mechanical pressure since the equation of state is evaluated with pressure approximated by p eos (z) = − 0 gz, with 0 the constant Boussinesq reference density (Young, 2010).
The second contact force arises from viscous and/or SGS stresses. These stresses can be organized into a second-order stress tensor, T, so that the corresponding contact force is where the second expression exposes the Cartesian tensor indices and makes use of the implied summation convention where repeated indices are summed over their range. We assume that there are no internal sources of angular momentum so that T is a symmetric tensor, with ≈ 10 −6 m 2 s −1 the kinematic molecular viscosity and S the rate of strain tensor (dimensions T −1 ). More general forms appear for coarse-grained equations with the molecular viscosity converted to a viscosity tensor (e.g., see Smagorinsky, 1993, as well as section 17.4 of Griffies, 2004). In general, T is known as a frictional stress tensor with constraints on the viscosity tensor such that T is a sign-definite sink of kinetic energy.
The third contact force arises from the advective transport of linear momentum, which acts as a mechanical momentum flux. This nonlinear term is referred to as self-advection, and it allows for turbulent motion and the turbulence cascade. In the context of contact forces we refer to this term as a kinetic stress (see section 13.5.2 of Thorne & Blandford, 2017), whose form is revealed when moving to the laboratory or Eulerian reference frame. In this case the material time derivative of velocity takes on the Eulerian flux-form where v⊗v is the outer product of the velocity vector with components The kinetic stress is given by − v⊗v, and the corresponding contact force is When coarse-grained, the kinetic stress leads to SGS turbulent Reynolds stresses.

10.1029/2019MS001954
Contact forces are in local mechanical equilibrium in accord with Newton's third law (action/reaction law), for example, section 4.3 of Griffies (2004) and section 2.3 of Olbers et al. (2012). Hence, the net contact force acting on a finite fluid region is due only to the contact forces acting on the region boundary. This property of contact forces leads to the notion of a form stress, which arises from the horizontal component of the pressure force acting on a sloped surface such as an interior fluid interface (interfacial form stress), the ocean surface (atmosphere/ocean form stress), or the ocean bottom (topographic form stress). We have more to say about form stress in section B0.1.

A2.3. Eulerian Flux-Form Momentum Equation
Bringing all forces together leads to the Eulerian flux-form momentum equation with I the 3 × 3 identity tensor. By definition, body forces appear as sources of momentum whereas the contact forces arise from the divergence of stresses. This distinction has implications when formulating the momentum budget for finite volumes in section 2.2. As a relation between geometric objects (vectors and tensors), the momentum budget (A14) is independent of coordinate representation.

A3. Mass Conservation for Seawater and Trace Constituents
As noted above, we formulate the ocean equations for fluid elements that maintain a constant mass. The Lagrangian expression for mass conservation is D( V)∕Dt = 0, which takes on its Eulerian flux-form as the continuity equation The corresponding Eulerian flux-form equation for trace matter, such as salt, is given by the salinity equation The flux J(S) parameterizes SGS motions typically in the form of an advection and/or diffusion process. The simplest form arises from molecular diffusion with S ≈ 1.2 × 10 −9 m 2 s −1 the molecular diffusivity for salt (e.g., page 56 of Olbers et al., 2012). Note that we ignore the cross-diffusion process (e.g., section 2.5 of Olbers et al., 2012).

A4. Barycentric Velocity
Seawater is a multicomponent fluid that is well approximated as a binary mixture of freshwater and salt. Following DeGroot and Mazur (1984) (see their Chapter II) and Olbers et al. (2012) (see their section 2.2), we interpret the velocity, v, as the barycentric velocity. An important consequence of this kinematic interpretation is that the SGS matter flux, J(S), vanishes when the matter concentration is spatially constant This kinematic property is trivially realized for molecular diffusion as given by Equation A17. Correspondingly, the continuity Equation A15 has no SGS operator even after coarse graining; that is, there is no diffusion of seawater mass, and hence, there are no mass sources/sinks in the ocean interior. This property means that a diffusive flux of salt crossing the boundary of a fluid element balances an oppositely directed diffusive flux of freshwater, thus leaving the fluid element with a constant mass but with a nonconstant salt and freshwater content. Nurser and Griffies (2019) explore the associated consequences for surface boundary fluxes of water and salt.
We maintain the kinematic property (A18) for the cell averaged finite volume (weak formulation) of the ocean equations as described in section 3.3. For the continuum equations one can make use of density-weighted averages as proposed by Hesselberg (1926) and

A5. Conservative Temperature Equation
As noted by McDougall (2003), the most suitable scalar field for representing ocean heat transport is given by the potential enthalpy divided by the specific heat capacity, thus defining the Conservative Temperature, Θ. The Eulerian flux-form equation for Θ takes on the same form as that for salinity in Equation (A16): Likewise, the mathematical form of the subgrid flux of Conservative Temperature, J(Θ), is generally the same as that for salinity in Equation (A16) (McDougall, 2003), though with possibly distinct dianeutral diffusivities in the presence of double-diffusion (e.g., Schmitt, 1994), in which the molecular diffusive flux of Conservative Temperature occurs with Θ ≈ 1.4 × 10 −7 m 2 s −1 ≈ 100 S (e.g., page 55 of Olbers et al., 2012). Additionally, temperature experiences a contribution from penetrative solar radiation in the upper ocean.

A6. Evolving the Ocean State
The ocean state is determined by the three components to the velocity, v = ux + vŷ + wẑ; the in situ density, ; the salinity, S; the Conservative Temperature, Θ; and the pressure, p. Evolution of these seven fields is governed by the six prognostic equations A1-A4 along with the equation of state (A5). The subgrid fluxes and stress tensor are determined by the resolved model fields according to constitutive relations (e.g., the Newtonian relation between stress and strain rate) or parameterizations/closures. Furthermore, boundary fluxes of momentum, matter, and enthalpy are operationally incorporated through Neumann boundary conditions placed on the subgrid fluxes and stresses. Additionally, we incorporate penetrative solar radiation into the Conservative Temperature equation via contributions to its SGS flux vector.
There is no prognostic equation for pressure, with its evolution determined diagnostically. Specifics for its diagnosis depend on the dynamical and/or kinematical approximations. For the unapproximated equation set (A1)-(A5), pressure can be determined by inverting the equation of state after updating , S, and Θ, following the approach used for compressible nonhydrostatic atmospheric models (e.g., Chen et al., 2013). Auclair et al. (2018) also propose a distinct method for ocean modeling. Assuming incompressible Boussinesq kinematics filters acoustic modes to render a three-dimensional Poisson equation whose inversion determines the updated pressure (Marshall et al., 1997). Finally, for hydrostatic flow, pressure is diagnosed as a function of density and depth, with this approach familiar from hydrostatic primitive equation ocean models such as the one pioneered by Bryan (1969). Likewise, for hydrostatic isopycnal layer models, the Montgomery potential takes the role of pressure (e.g., Bleck, 1998).

A7. General Comments on SGSs
Numerical simulations are generally made with discrete grid cells of size L grid ≫ L micro ≈ 10 −3 m. Correspondingly, the equations describing coarse grid simulations involve correlations of nonlinear fluctuations occurring at the finer SGSs. The correlations play a role, sometimes a dominant role, in evolution of the simulated scales of motion. Parameterization of these correlations in terms of simulated motions constitutes an active research problem in theoretical physical oceanography with very important impacts on the integrity of a simulation.
To help in partitioning between resolved and unresolved scales of motion, we here introduce a length scale based on the spatial gradient of the flow field Decomposing velocity fluctuations into Fourier modes leads us to propose that an accurate simulation of velocity fluctuations with length scales L gradient requires a grid length scale That is, we need no less than 2 grid points to accurately represent a simulated fluctuation of length scale L gradient . This is a rather general statement that ignores details of numerical algorithms. However, it accords with the condition proposed by Hallberg (2013) for horizontally resolving baroclinic eddies; by Stewart et al. (2017) for vertically resolving baroclinic modes; and by Chen et al. (2018) for the minimum grid points to accurately capture a wave phase speed at different orders of accuracy and grid staggering. See also Soufflet et al. (2016) for an analysis into the effective resolution associated with specific numerical algorithms.

A8. Comments on Mesoscale Eddy Parameterizations
Mesoscale eddies are commonly parameterized in coarse resolution ocean climate models that do not explicitly resolve transient mesoscale eddies. Following the proposal from Gent and McWilliams (1990) and Gent et al. (1995), these eddies are parameterized by adding an eddy contribution to the velocity that is used to advect tracers. The sum of the resolved velocity and the eddy velocity is termed the residual mean velocity.
In Eulerian geopotential models, the Gent et al. (1995) parameterized eddy-induced velocity is derived from a vector stream function according to where is the locally referenced potential density, GM > 0 is the diffusivity for the eddy-induced velocity, S is the slope of the neutral directions, and the in situ density, , is set to a constant reference value when making the Boussinesq approximation.
Isopycnal layer models implement mesoscale parameterizations by adding a horizontal bolus velocity, u bolus , to the model's resolved horizontal velocity used for seawater mass and tracer transport. The resulting sum is the horizontal residual mean u rm = u + u bolus . (A24) Figure A1. Illustrating how seawater mass and tracer mass evolve within a fluid layer using the vertical Lagrangian-remap method. This layer has one of its sides abutting the ocean bottom as in Figure 2. Mass evolves according to the layer convergence of the horizontal mass flux, −∇ s ·(h u rm ), as well as the convergence of dia-surface mass flux, −Δ s ( w ( . s) ). Tracer evolves similarly along with horizontal and dia-surface flux convergences from diffusion and boundary fluxes, −∇ s · J (h) − Δ s J ( . s) . The two dia-surface fluxes are drawn as vertically oriented fluxes, since when multiplied by the projected horizontal area element they both measure the dia-surface transport crossing the interface. Note that u rm is the horizontal residual mean velocity built from the resolved horizontal velocity plus a parameterized eddy-induced bolus velocity. There is no dia-surface component to the bolus velocity; it is purely horizontal and within layers (see McDougall & McIntosh, 2001, for details). In the Lagrangian portion of the Lagrangian-remap method, the thickness and tracer concentration evolve with w ( . s) = 0, even while maintaining the full subgrid scale fluxes in the tracer equation. The regrid step then rearranges the s-layers according to the prescribed target values. The remap step then estimates the ocean state on the new grid, with remapping equivalent to vertical advective transport using w ( . s) ≠ 0.

10.1029/2019MS001954
The bolus velocity is horizontally divergent, and it contributes to transport within an isopycnal layer. Importantly, it contributes zero diapycnal transport (McDougall & McIntosh, 2001). Hence, it acts to rearrange matter and tracer within isopycnal layers. When following the notions of Gent and McWilliams (1990) and Gent et al. (1995), the bolus velocity is parameterized according to a downgradient diffusion or smoothing of isopycnal heights, which is equivalent to thickness diffusion away from bathymetry for a vertically uniform diffusivity. Vertical Lagrangian-remap models follow the same approach as used by isopycnal models since horizontal advection remains part of the Lagrangian portion of the method. Hence, there is zero dia-surface transport arising from the bolus velocity in vertical Lagrangian-remap models. In Figure A1 we illustrate the terms appearing in the tracer equation as implemented in vertical Lagrangian-remap models, including the bolus velocity.

Appendix B: Forces From Pressure and Subgrid Scales
We here summarize the rudiments of how forces from pressure and subgrid scale stresses are treated in ocean models.

B1. Pressure
Gauss's divergence theorem allows us to write the net pressure force acting on a finite region in two equivalent manners The first expression is the volume integral of the "pressure gradient body force" as numerically realized in Bryan (1969). The second expression is the area integral of the contact pressure force acting in a compressive manner on the grid cell boundary. A contact force formulation leads to a finite volume representation of the pressure force as in Lin (1997) for the atmosphere and Adcroft et al. (2008) for the ocean.
The contact force formulation leads to the notion of form stress, which is the horizontal pressure force per area acting on a sloped surface. (Note that some authors, such as section 3.6 of Vallis, 2017, define form stress as an area average of the form stress defined here.) The form stress acting on the upper side of a surface with depth z = (x, , t) is given by the horizontal vector Hence, on the upper interface of a grid cell the form stress acts in the positive coordinate direction if the interface has a positive slope in that direction. The form stress acting on the lower side is equal in magnitude yet pointed in the opposite direction, as per Newton's law of action/reaction. The net horizontal pressure force acting on a grid cell is thus given by In a hydrostatic fluid, the vertical pressure force acting on the grid cell balances the weight of fluid within the cellẑ with M cell = ∫ cell dV the grid cell mass. The net vertical hydrostatic pressure force acts upward (+ẑ) since the upward pressure acting on the lower interface is greater than the downward pressure acting on the upper interface. Furthermore, this upward hydrostatic pressure force balances the downward (−ẑ) weight of fluid within the cell.

B2. Subgrid Scale Stresses
The dia-surface SGS force acting on the upper and lower surfaces of a finite fluid region can be written The SGS dia-surface stress vector (dimensions of force per area), sgs dia , is commonly parameterized as boundary drag at the ocean bottom (e.g., Deremble et al., 2012) and wind/ice stress at the ocean surface. In the interior it is commonly parameterized as viscous friction due to vertical shear using a viscosity parameterized according to processes such as breaking gravity waves (MacKinnon et al., 2013). In this case the model implementation is generally via the divergence of the stress tensor as written where F sgs dia is the force per volume acting on a grid cell due to the divergence of dia-surface frictional stresses.
The lateral subgrid stress is commonly parameterized via a Hooke's law stress-rate of strain relation (Smagorinsky, 1993) T sgs lateral = C · S, with C the fourth-order viscosity tensor (dimensions L 2 T −1 ). This expression generalizes that for molecular viscosity given by Equation A10 . It is generally designed to enhance grid scale dissipation of the direct (i.e., downscale) cascade of enstrophy, with the cascade induced either through resolved nonlinear flows (via a turbulence cascade) or numerically (via dispersion errors). The most common means for implementing lateral frictional stress is via the divergence of the stress tensor so that we introduce the frictional force vector as Note that the frictional tensor generally has a zero trace when the trace of the full stress tensor is assumed to be captured by pressure. However, it is notable that the mesoscale eddy parameterization from Anstey and Zanna (2017) and Bachman et al. (2018) introduces a subgrid stress tensor that has a nonzero trace, with the trace interpreted as a dynamical modification to pressure. Holloway (1992), Smagorinsky (1993), , Large et al. (2001), Smith and McWilliams (2003), Griffies (2004), Fox-Kemper and Menemenlis (2008), Jochum et al. (2008), Jansen and Held (2014), and Bachman (2017) offer a variety of perspectives on the use of lateral friction operators in numerical ocean models.

Appendix C: Connecting the Integral and Differential Formulations
In this appendix we derive the Leibniz-Reynolds transport theorem and then make use of the theorem to develop the connection between the partial differential equations for a fluid element (strong formulation from Appendix A) and the integral equations for a finite fluid region (weak formulation from section 2).

C1. Leibniz-Reynolds Transport Theorem
We offer a somewhat simplistic derivation of the Leibniz-Reynolds transport theorem that does not make use of the Jacobian methods commonly used in the fluids literature (e.g., Aris, 1962). For this purpose consider the integral of a space-time-dependent function over a one-dimensional spatial domain whose bounds are functions of time Leibniz's rule provides the means to differentiate this integral Generalization to an integral over a three-dimensional domain, , leads to the Leibniz-Reynolds transport theorem

10.1029/2019MS001954
In the second equality we made use of Gauss's divergence theorem to transfer between a boundary integral over  with oriented area elementn d, and a volume integral over the region . We furthermore introduced a shorthand for the velocity of a point on the boundary We emphasize that the boundary term in Equation C3 projects out just the normal component to the boundary velocity; we never make use of the tangential component.
The Leibniz-Reynolds transport theorem is of great practical utility for fluid mechanics since we are often interested in the evolution of fluid properties within an arbitrary moving domain. As shown in the next subsections, it provides the means to connect the weak and strong formulations of the ocean equations as discussed in section 2. Correspondingly, it forms the starting point in our formulation of ocean model dynamical cores making use of finite volume methods as well as arbitrary Eulerian-Lagrangian and Lagrangian methods.

C2. Mass Conservation
Setting = 1 in the Leibniz-Reynolds transport theorem (C3) yields the time change for the region volume d dt whereas = yields the time change for the region mass d dt A region that moves according to the barycentric velocity of fluid elements, (v (b) − v) ·n = 0, maintains constant seawater mass For this result to hold for an arbitrary region of fixed mass requires the mass continuity to be maintained at each space point These arguments exemplify the connection between the weak formulation (integral form) of mass conservation to the strong formulation (differential form). Furthermore, substituting the strong form of mass conservation, Equation C8, into the finite volume budget (C6b) for mass changes over an arbitrary region leads to d dt which is the starting point for the finite volume mass budget in section 2.1. Finally, we note the dual role played by v, which is both the seawater mass flux via the continuity Equation C8 and the linear momentum per volume.
This constraint in turn means that the specific thickness, z s = z∕ s = 1∕s z is single-signed and bounded away from infinity: |z s |<∞. The constraint of a nonvanishing s-stratification ensures that the GVC surfaces are monotonically stacked in the vertical. This property is fundamental to all GVC models. Regions where physical processes can lead to s-surface overturns require either a regridding procedure to ensure nonsingularity or the appending of another submodel to simulate/parameterize processes in the overturning region (e.g., bulk boundary layers appended to isopycnal models as in Hallberg, 2003).
The GVC coordinates are non-orthogonal, thus requiring some care in handling their tensor properties. Although a bit subtle in this regard, they offer great advantages for geophysical modeling. In particular, maintenance of the same horizontal directions as in a geopotential coordinate treatment isolates the gravitational force to the vertical momentum equation. Correspondingly, it means that all GVC models make use of the same horizontal velocity that is aligned perpendicular to gravity This orientation contrasts to that arising in a locally orthogonal coordinate representation, whereby gravity also projects into the lateral directions and we locally rotate velocity components to align with sloped GVC surfaces. As noted by Starr (1945), such a locally orthogonal approach is problematic for geophysical fluid modeling. The reason is that the gravitational acceleration is so dominant in geophysical flows so that it is critical to restrict its appearance to a single component of the momentum equation, with that component reducing to the hydrostatic balance for large-scale flows.

D2. Layer Thickness
Hence, the specific thickness expands or contracts the layer thickness according to stratification of the vertical coordinate surfaces. At the ocean surface, the top interface of the surface grid cell is defined by the free surface, z = (x, , t), whereas the bottom interface of the bottom-most grid cell is defined by the bottom topography, z = −H(x, ) (see Figure 2). When extending our considerations to a finite volume grid cell as in section 3, we consider the cell thickness as defined by Equation 19, which is the average of h over the horizontal extent of the grid cell. Figure D1 shows a geometric expression for the zonal derivative of a function computed using GVCs. This expression allows us to readily compute the transformation between horizontal partial derivatives computed in GVCs to those computed in geopotential coordinates

D3. Partial Derivatives
where x(C) = x(B) follows by construction. The final equality is a finite difference estimate of ( / x) z + ( z/ x) s ( / z). We are thus led to the relation between horizontal derivative operators The time derivative has a similar transformation (D6) Figure D1. Geometric illustration of the horizontal derivative operator on a constant s-surface. As shown here, the zonal derivative is computed by taking the difference of the function along the s-surface and dividing by the horizontal We also show the expression for the slope of a s-coordinate surface, tan , which equals to the change in the depth of the surface moving along the horizontal direction, We understand the need for the second term in each of the constant s derivative operators since we must follow the space-time-dependent s-surfaces when computing the horizontal derivative and the time derivative. Finally, the vertical derivative operators are related through the chain rule We emphasize that the notation, ∇ s , is a mere shorthand for the horizontal derivative operator given by Equation D5. This notation is ubiquitous in the literature on GVCs, including isopycnal and terrain-following models. Even so, as emphasized by Young (2012), there are certain tensor properties that are not carried by ∇ s that have implications for how one treats integrals involving ∇ s acting on a horizontal vector.

D4. Dia-surface Transport
A GVC surface is generally moving. So to diagnose the volume flux of fluid penetrating the surface requires us to subtract the velocity of the surface, v (s) , from the velocity of a fluid element, v. Projecting the component of this relative velocity into the direction of the surface normal measures the amount of fluid penetrating the surface (see Figure D2): volume flux of fluid crossing a moving GVC surface =n · (v − v (s) ).
The dia-surface volume flux is positive if there is motion across the surface in the direction defined byn. There is zero volume flux when the normal component of the surface velocity matches that of fluid elements. We are only concerned with the flux of fluid normal to the surface, with the flux components in the tangent plane requiring dynamical information that is generally unavailable. In the following, we derive a fundamental kinematic relation that connects the flux,n · (v − v (s) ), to the material time derivative of the s-surface.

D4.1. Surface Normal
The surface normal direction can be represented in either of two mannerŝ . The horizontal projection of the surface area element is given by dA = | cos | d, where is the angle between the s-surface and the horizontal.

D4.5. Expressions for the Material Time Derivative
Making use of the expression (D16) for the dia-surface velocity, as well as the partial derivatives (D5)-(D7), renders the equivalent expressions for the material time derivative operator D Dt = We see that the dia-surface velocity, w ( . s) = z s Ds∕Dt, serves a role for GVCs that parallels the vertical velocity component for geopotential coordinates, with the two related by

D5. Dia-surface Subgrid Scale Transport
Analogous to the dia-surface matter transport, we can relate the Cartesian components of a subgrid scale transport, J, to the dia-surface subgrid transport by writing dia-surface SGS transport = J ·n d = J · ∇s | z∕ s| dA ≡ J ( . s) dA.
The SGS flux, J ( . s) , is the "vertical" flux which, when multiplied by the horizontal area element dA = dx d , measures the SGS tracer transport crossing the s-surface in the normal direction.

D6. Vector-Invariant Velocity Equation
Rather than the flux-form horizontal momentum Equation 31, many ocean models implement a discrete version of the vector-invariant velocity equation. They do so for a variety of reasons, one being that the metric terms appearing from expanding the kinetic stress term ∇ s · (h u⊗u) with generalized horizontal coordinates (e.g., spherical coordinates) can be awkward (see section 4.4.1 of Griffies, 2004 for details). Another reason is that the appearance of the vorticity and kinetic energy (see below) facilitate development of energy-enstrophy conserving finite difference methods such as those pioneered by Sadourny (1975) and Arakawa and Lamb (1981). In contrast to the flux-form momentum equation, the vector-invariant velocity equation does not admit a finite volume formulation. Instead, it arises from the continuous velocity equation as derived here.
To derive the vector-invariant velocity equation for a hydrostatic fluid, start from the inviscid horizontal velocity equation Du Dt +ẑ ∧ u = − −1 ∇ z p.
Now expand the material time operator according to Equation D17 so that and express the horizontal self-advection as with K = u · u∕2 and =ẑ · (∇ s ∧ u), the kinetic energy per mass in the horizontal flow and the vertical component to the relative vorticity using GVCs, respectively. These identities then lead to [ u t ] s + w ( . s) u z +ẑ ( + ) ∧ u = −∇ s K − −1 ∇ z p.
As per Equation D5, in GVCs the horizontal pressure gradient is decomposed into two terms with Φ = gz the geopotential (ignoring astronomical tide forcing) and again assuming hydrostatic balance. This decomposition leads to the vector-invariant velocity equation with K + Φ the mechanical energy per mass in the horizontal flow.
Appendix E: Table of Symbols   Table E1 collects the mathematical symbols used in this paper, and Table E2 lists the physical symbols. These tables assist in cross-referencing symbols used in more than one section of this paper.