Spurious Forces Can Dominate the Vorticity Budget of Ocean Gyres on the C‐Grid

Gyres are prominent surface structures in the global ocean circulation that often interact with the sea floor in a complex manner. Diagnostic methods, such as the depth‐integrated vorticity budget, are needed to assess exactly how such model circulations interact with the bathymetry. Terms in the vorticity budget can be integrated over the area enclosed by streamlines to identify forces that spin gyres up and down. In this article we diagnose the depth‐integrated vorticity budgets of both idealized gyres and the Weddell Gyre in a realistic global model. It is shown that spurious forces play a significant role in the dynamics of all gyres presented and that they are a direct consequence of the Arakawa C‐grid discretization and the z‐coordinate representation of the sea floor. The spurious forces include a numerical beta effect and interactions with the sea floor which originate from the discrete Coriolis force when calculated with the following schemes: the energy conserving scheme; the enstrophy conserving scheme; and the energy and enstrophy conserving scheme. Previous studies have shown that bottom pressure torques provide the main interaction between the depth‐integrated flow and the sea floor. Bottom pressure torques are significant, but spurious interactions with bottom topography are similar in size. Possible methods for reducing the identified spurious topographic forces are discussed. Spurious topographic forces can be alleviated by using either a B‐grid in the horizontal plane or a terrain‐following vertical coordinate.


The Depth-Integrated Vorticity Equation
Vorticity diagnostics are an underused tool for interpreting model circulations and offer a description of gyre dynamics that complements textbook theory (Vallis, 2017). A handful of recent papers have used a vorticity budget to diagnose regional and global GCM models (Hughes & de Cuevas, 2001;Le Bras et al., 2019;Le Corre et al., 2020;Schoonover et al., 2016;Yeager, 2015).
To obtain a depth-integrated vorticity budget analytically we start from the vector-invariant form of the momentum equation: Where f is the Coriolis parameter, u h is the ''horizontal'' (parallel to the Earth's surface) velocity vector,  is the vertical divergence of the vertical diffusive momentum fluxes (which relates to the surface momentum fluxes when vertically integrated),  is the horizontal divergence of the horizontal diffusive momentum fluxes, ∇ h is the horizontal gradient operator, and [ ⋅ ] h is the horizontal component of a vector. To derive a depth-integrated vorticity equation, we need to depth-integrate and take the curl of Equation 1. The order of the two operations and any multiplications carried out significantly alters the form and physical meaning of the obtained depth-integrated vorticity equation.
If we choose to depth-integrate the curl of the momentum equation, the pressure gradient vanishes upon taking the curl and bottom vortex stretching represents the interaction of the currents with the sea floor. Both the beta effect and bottom vortex stretching originate from the Coriolis acceleration in Equation 1. In the model, the curl of the single momentum diagnostic associated with the Coriolis acceleration will be responsible for two distinct physical processes. Jagannathan et al. (2021) use this form of vorticity budget to investigate flow interactions with idealized bathymetry. In Section 6.2 we discuss whether the spurious forces identified in this article emerge in this vorticity budget.

of 31
If we choose to take the curl of the depth-averaged momentum equation then sea floor interactions are represented by the JEBAR term (Joint Effect of Baroclinicity and Relief). Cane et al. (1998) and Drijfhout et al. (2013) have questioned the relevance of JEBAR by presenting simple examples in which there is no flow immediately above the bathymetry. In these examples there is trivially no interaction between the flow and the bathymetry, but there is a non-zero JEBAR term.
Throughout this paper we consider the vorticity equation obtained by taking the curl of the depth-integrated momentum equation: (2) Here ζ is the vertical component of the vorticity, τ top is the surface stress due to wind and sea ice, τ bot is the bottom stress due to friction at the sea floor,  is the lateral diffusion of depth-integrated relative vorticity, η is the free surface height, H is the depth of the sea floor, and P b is the pressure at the sea floor. In Equation 2 we omit the contribution from free surface undulations as we assume the time evolution of the free surface is small and we omit atmospheric pressure torques as we are assuming there are no atmospheric pressure gradients above the ocean. The derivation of Equation 2 (including the omitted terms) is presented in Appendix A.
The terms on the right-hand side of Equation 2 are the following: the advection of planetary vorticity; the bottom pressure torque; the surface stress curl; the curl of bottom friction; the lateral diffusion of relative vorticity; and the advection of relative vorticity. The planetary vorticity term in Equation 2 contains contributions from the evolving free surface and surface water fluxes as ∇ h ⋅ U = −∂η/∂t + Q/ρ 0 , where Q is the surface water flux due to evaporation, precipitation, and run-off. In an equilibrated state, the free surface evolution is small and the divergence caused by realistic water fluxes is negligible. Hence, we assume ∇ℎ ⋅ ( ) ≈ where β represents the linear variation of f with latitude and V is the meridional component of the depth-integrated velocity. This formulation is practical as topographic interactions emerge from pressure gradients in the form of the bottom pressure torque and beta effects emerge from the curl of the Coriolis acceleration; the Coriolis acceleration is responsible for one physically meaningful term in the analytic vorticity budget. Equation 2 is also used in Stewart et al. (2021). Jackson et al. (2006) conclude that the leading order balance between bottom pressure torques and the planetary vorticity term in Equation 2 is crucial for steering jets and western boundary currents over bottom topography. They continue to argue that the form of the topographic steering determines if bottom friction is able to modify the geometry of the current.
As a consequence of Stokes' theorem, the area integral of a term from Equation 2 is directly related to the line integral of the depth-integrated forces along the area edge. This is particularly useful when considering area integrals of terms from the vorticity equation and is discussed further in the next sub-section.

Contour Integration Method
All terms in the depth-integrated vorticity equation can be expressed as the curl of a depth-integrated acceleration in the momentum equation: where Ω is a term in the depth-integrated vorticity equation and M is a term in the depth-integrated momentum equation. If we integrate Ω over the area enclosed by a depth-integrated streamline, we can interpret the integral using Stokes' theorem: Where A ψ is the area enclosed by a depth-integrated streamline and Γ ψ is the anticlockwise path along the same streamline. The criteria for selecting the sign in Equation 4 is defined later in this paragraph. The depth-integrated stream function, ψ, only exists if the flow is steady and ∇ h ⋅ U = 0. If a long time-average of a varying flow is taken and the surface water fluxes are sufficiently small, a quasi-streamline can be calculated which approximately follows the circulation. The integral ρ 0 I(ψ) can be interpreted as the work done per unit area by the force associated with M on a fluid column in one circulation of Γ ψ . For a gyre circulating in a clockwise direction, the direction of circulation would be opposite to the conventional anticlockwise direction of Γ ψ . So that the reader does not have to constantly consider the direction of the flow relative to Γ ψ we select the sign in Equation 4 so a positive value of I(ψ) corresponds to a force that is spinning the gyre up.
Analytically, we would expect the planetary vorticity term to vanish upon integration as a consequence of the divergence theorem: where ̂ is the horizontal vector which is normal to the streamline and the depth-integrated velocity. The Coriolis force can still play a role in shaping the streamlines of the circulation but ultimately has no influence on the integrated budget. Although the divergent part of the advection term, ∇ℎ ⋅ ( ∫ ) , has a similar form, we do not expect the same zero integral for this term as the depth-integrated product of ζ and u is not parallel to U in general.
This method has been used in models before. Schoonover et al. (2016) integrated vorticity diagnostics over a limited number of streamlines in the North Atlantic and concluded that the wind stress curl is largely balanced by bottom pressure torques. Stewart et al. (2021) also used this method in an isopycnal model and concluded that wind stress curl is not balanced by bottom pressure torques in general. Stewart et al. (2021) discuss how the integrating area affects the resultant vorticity balances and in their model the wind stress curl is only balanced by bottom pressure torques when integrated over latitude bands. Jackson et al. (2006) note that in their idealized models the vorticity input from bottom friction mostly disappears when integrated over latitude bands but can be significant when integrated over the area enclosed by streamlines.
In Gula et al. (2015), terms in the barotropic vorticity budget are integrated over an entire subdomain covering the South Atlantic Bight to study the interaction between the Gulf Stream and the continental shelf. Flow through the boundaries of the nested model is permitted so we would not expect the planetary vorticity term to vanish upon integration. These integrations show a leading order balance between the planetary vorticity term and bottom pressure torques and suggest that bottom pressure torques are the dominant mechanism for removing the planetary vorticity imported into the subdomain by the Gulf Stream.
It should be noted that Schoonover et al. (2016), Stewart et al. (2021), Jackson et al. (2006), and Gula et al. (2015) all use a terrain-following coordinate in their models but in this article we study the vorticity budget of a z-coordinate model. In Section 6.4 we discuss how the vorticity budget can be affected by the choice of vertical coordinate and how terrain-following coordinates can mitigate spurious Coriolis forces related to the topography.

The Discrete Depth-Integrated Vorticity Equation
In many contemporary ocean GCMs, the discretized model variables are distributed on the C-grid (Mesinger & Arakawa, 1976). The geometry of the C-grid is shown in Figure 1: T points hold scalar information including the divergence of the flow; the U and V points hold the horizontal components of vector quantities including the horizontal velocity, surface stresses, and accelerations in the momentum equation. Values closely related to vorticity are found on F points, this includes the relative vorticity, the streamfunction, and terms in the depth-integrated vorticity equation (Ω). Vertical velocities are located on W points that are directly above and below T points as shown in Figure 1. The Coriolis parameter can be evaluated at any point on the C grid but F point values are used for calculating the Coriolis acceleration in most models that use a vector invariant momentum equation because the relative and planetary vorticity are then evaluated at the same point (see Section 3.2). In this article, f i,j refers to the value of the Coriolis parameter centered on the F point and ( ) refers to the Coriolis parameter centered on the U (V) point as shown in Figure 1.
Every point in the C-grid has an associated cell with a vertical thickness and horizontal width. Throughout this article e 3t is the T cell vertical thickness and e 1t , e 2t are the T cell widths in the i and j direction respectively. The same convention is used for U, V, and F cells also. It should be noted that the values of the F cell thicknesses in this article depend on the scheme used to calculate the Coriolis acceleration (see Section 3.2).
The GCM configurations discussed in this paper use a primitive momentum equation that is a discrete equivalent to the vector invariant momentum equation (Madec et al., 2019). Momentum diagnostics can be combined to represent terms in the analytic momentum equation (Equation 1). The curl of the depth-integrated momentum diagnostics is taken to form a closed discrete vorticity budget that is valid in an unsteady state as the time derivative diagnostic is included. The resultant vorticity diagnostics should closely resemble the terms in the depth-integrated vorticity equation (Equation 2); however, the planetary vorticity diagnostic deviates from the planetary vorticity term in several significant ways.

The Discrete Coriolis Acceleration
The Coriolis acceleration is a product of the Coriolis parameter, f, and the horizontal velocity u h . There are many possible schemes for calculating their product and the choice of scheme affects the quantities that are conserved in the model flow.
Consider the following discrete Coriolis acceleration: where COR x (COR y ) is the x (y) component of the discrete Coriolis acceleration which is centered on a U (V) point; ̃= 1 3 and ̃= 2 3 are volume fluxes; and r n points to one of the four neighboring V or U points. If we depth-integrate the acceleration in Equation 6 and then take the curl, we obtain the following equation for the discrete planetary vorticity term: where PVO is the discrete planetary vorticity term which is centered on an F point. Equation 7 is the discrete calculation of −∇ℎ ⋅ ( U) averaged over the four T cells surrounding the central F point and is therefore closely related to the analytic planetary vorticity term in Equation 2. The Coriolis acceleration given in Equation 6 is not used in C-grid models as it lacks the energy and/or enstrophy conserving properties of other mainstream schemes. However, when studying the discrete depth-integrated vorticity budget it is useful to consider how the used Coriolis accelerations deviate from this reference value as any difference may emerge as a departure from the discrete calculation of −∇ℎ ⋅ ( ) in the vorticity budget.
When using a vector invariant momentum equation, mainstream schemes use multi-point and thickness-weighted averaging of f and u in order to conserve energy and/or enstrophy (Madec et al., 2019). A general form of the discrete Coriolis acceleration under these schemes is: Where a n , b n , and c n are the horizontal locations of three neighboring grid points (not necessarily different) for the nth term of the sum. COR x and COR y are still centered on U and V points respectively. Note that a n is always the location of an F point and c n is always the location of a U or V point. Depending on the scheme, the 3 term can be either a U, V, or F cell thickness so b n is the location of either a U, V, or F point. N is the number of terms in the average which depends on the choice of scheme. Equations 8 and 9 are valid on points near the bathymetry but if b n or c n points to a masked grid point (a point in the bathymetry) then the nth term in the sum is equal to zero.
In this article we consider three popular schemes for calculating the Coriolis acceleration. The energy conserving scheme (ENE) (Sadourny, 1975) conserves total horizontal kinetic energy and uses a four point average (N = 4). The enstrophy conserving scheme (ENS) (Sadourny, 1975) conserves potential enstrophy and has eight terms (N = 8). Finally the energy and enstrophy conserving scheme (EEN) (Arakawa & Lamb, 1981) conserves both horizontal kinetic energy and potential enstrophy and uses a 12 point average (N = 12). Barnier et al. (2006) demonstrates that the choice of scheme can significantly influence the global ocean circulation, especially in areas with strong current-topography interaction.
The explicit forms of the ENE, ENS, and EEN schemes for the Coriolis acceleration are given in Appendix B.
The results in Section 4 and 5 use the EEN scheme; however, in Section 6.1 we argue that all three schemes produce similar spurious forces. This argument is more concise when we use a form of the Coriolis acceleration that is general to the ENE, ENS, and EEN schemes.
We will decompose the general discrete Coriolis acceleration in Equations 8 and 9 by considering variations of f and e 3 around the U and V points: Where fc n ) is the value of the Coriolis parameter centered on the same point as the volume flux. Equation 11 will be applied to COR x (Equation 8) and Equation 12 will be applied to COR y (Equation 9). The α k (b n ) term is of order one and represents the scaling of 3 relative to other local cell thicknesses that only occurs in the EEN scheme. In the EEN scheme, 3 ( ) is an F cell thickness and F cell thicknesses are calculated using: Where masked T cell thicknesses are set to zero. When near bathymetry (masked points), the F cell thickness could be up to four times smaller than the typical unmasked T cell thicknesses surrounding it. The product 3 is the F cell thickness before this scaling is applied and will be more similar to the neighboring T cell thicknesses in Equation 13. This scaling of e 3f is unique to the EEN scheme and therefore α k = 1 in the ENS and ENE cases.
By combining Equations 8-12 we can derive a general decomposition of the Coriolis acceleration: where we have assumed that variations in f and the nonscaled cell thickness, 3 , are small. The x and y components of the Coriolis acceleration have a leading order contribution centered on the U and V point. The leading order term simplifies to the reference Coriolis acceleration in Equation 6 and therefore will resemble − ∇ℎ ⋅ ( ) in the discrete vorticity budget. Equation 14 is valid on points near the bathymetry but if b n or c n points to a masked grid point then the nth term of the entire sum is equal to zero.
The remaining terms may emerge as first order departures from − ∇ℎ ⋅ ( ) in the discrete vorticity budget. The first order contributions are: an f displacement term caused by the difference between the values of f where the volume fluxes are located and the values of f used in the scheme; a topographic effect caused by variations in cell thicknesses; and a coupled f-topographic effect caused by the combined effect of sudden changes in cell thicknesses near masked points and the previously mentioned f displacement term. Note that if α = 1 (true for ENS and ENE) then the f-topographic effect vanishes. The depth-integrated Coriolis acceleration is: Where max and max are the highest unmasked indices in the column and they may vary with horizontal index when z-coordinates are used. The depth-integrated Coriolis acceleration is therefore also sensitive to steps in the bathymetry. This is discussed in the next sub-section.

The Influence of Model Level Steps on the Coriolis Acceleration
In this section, we present a toy configuration that highlights how model levels can influence the discrete Coriolis acceleration. The configuration is shown in Figure 2. The configuration has two model levels, three U-grid points in the i direction, two in the j direction, and a rigid lid. The points in the upper level are surrounded by unmasked points, we assume the grid is regular, and cell widths are the same in the i and j direction. We also assume an f-plane so f does not vary.
The configuration has a step bathymetry and a current running alongside it. The current has no y component so v = 0 everywhere and therefore COR x = 0 at all points. The lower limb of the current decelerates by an amount u 1 and as a consequence of incompressibility a vertical velocity is induced which accelerates the upper current by u 1 .
Under these assumptions, the discrete Coriolis acceleration does not vary between the ENE, ENS, and EEN schemes and is: Figure 2. A toy model demonstrating how model levels influence the discrete Coriolis acceleration. A horizontal plan is shown for the upper and lower level as well as a view of the depth-integrated fields divided by the cell thickness Δz. Single arrows represent prescribed velocities; double arrows represent calculated Coriolis accelerations; and shaded cells represent bottom topography. Accelerations on the lower level are masked to prevent the velocity field from evolving into a flow that would violate the no penetration boundary condition. The central F point is marked by a cross and is where the depthintegrated vorticity is generated.
Which is effectively f multiplied by the four point average of u.
In the upper layer, the Coriolis accelerations, located on the V points marked by red triangles in Figure 2, are: In the lower layer, the Coriolis accelerations are set to zero as they lie on masked V points. The V points are masked to prevent accelerations into the topography that would violate the no penetration boundary condition. The depth-integrated Coriolis accelerations are: Where Δz is the constant cell thickness. It should be noted that u 1 vanishes when calculating the depth-integrated velocities but remains in the depth-integrated acceleration. The depth-integrated Coriolis acceleration depends on more than the depth-integrated velocities.
When we take the curl of the depth-integrated accelerations, we can see how a depth-integrated vorticity is generated: Where Δx is the constant cell width. Note that this value is located on the central F point shown in Figure 2.
The pressure gradient, lateral diffusion term (unless no-slip boundary conditions are used), and the horizontal advection term are ambiguous on the masked velocity points at the edge of the bathymetry (e.g., the V points in the upper right diagram of Figure 2). An explicit momentum balance cannot be resolved and the Coriolis acceleration is the only non-zero and unambiguous acceleration into the sea floor. There should be no net acceleration into the bathymetry or else the no penetration boundary condition would be violated, so all accelerations that are incident on bathymetry are masked and set to zero. The masking of all accelerations can be interpreted as the addition of a spurious term to the discrete Coriolis acceleration. This spurious force is of unclear physical origin and is not realistic as it is localized to grid points that lie near model level steps. We can think of the result in Equation 21 as either the curl of this spurious force or as a form of spurious vortex stretching that takes place on F points near model level steps (cf. Bell, 1999).

Decomposing the Planetary Vorticity Term
In Section 3.2 we concluded that the discrete Coriolis acceleration used in mainstream schemes contained spurious contributions caused by f displacement, variations in cell thicknesses, and a coupled f-topographic effect. In Section 3.3 we demonstrated how spurious contributions from model level steps exist in the depth-integrated discrete Coriolis acceleration. The four found spurious contributions have the potential to emerge in the planetary vorticity diagnostic which is calculated by taking the curl of the depth-integrated Coriolis acceleration: We can therefore express the planetary vorticity diagnostic as the sum of five components: is also needed. The calculations are listed in Table 1 along with the components from Equation 23 they include. By linearly combining the fields from each calculation we can isolate each component in Equation 23. The f-topographic component is calculated by finding the difference between the complete planetary vorticity diagnostic and the sum of the four other components; therefore the five components add up to the complete planetary vorticity diagnostic by construction.

Contour Integration on a C-Grid
Calculating the curl on a C-grid is consistent with Stokes' law applied to an F cell, and integrating ∇ × M ⋅ k over several adjacent F cells is equivalent to a line integral of M around them (see Figure 3). As the streamfunction, ψ, is defined on F points we can approximate that the area enclosed by a streamline is a collection of interior F cells and that the area integral of vorticity diagnostics is the line integral of model accelerations around them. This is an approximation as we are assuming that the streamline follows the rectangular edges of the interior F cells but the resultant error is minimized if we first interpolate the points onto a sufficiently fine grid. The asymptotic value the contour integral tends toward as the interpolation resolution is increased should be free of area error. This method is applied to all contour integrals presented in Sections 4 and 5. Any non-topographic contributions to the contour integral that remains after the interpolation will be described as a numerical beta effect.
A numerical beta effect can emerge from −∇ h ⋅ (f U)| i,j even after being interpolated onto a fine grid as the divergence is calculated over the four T cells that surround the central F point (see Equation 7). When the internal F points are summed within the contour, the local domains for calculating the grid point divergences will overlap meaning the resultant area integral will not satisfy the divergence theorem in general. Overlapping local domains are a requirement of the C-grid's horizontal geometry. In Section 6.3 we highlight how the divergence calculation on a B-grid only requires information from a single tracer cell. The local domains for calculating the divergence do not overlap when integrating on the B-grid and the associated numerical beta effect will not emerge.

Details of the Configuration
The first experiment in this article is an idealized double gyre configuration based on the GYRE PISCES reference configuration in NEMO. The GYRE PISCES reference configuration has been used for a wide range of experiments (Lévy et al., 2010(Lévy et al., , 2015Perezhogin, 2019;Ruggiero et al., 2015). The domain is a closed rectangular basin which is 3,180 km long, 2,120 km wide, and is rotated at an angle of 45° relative to the zonal direction. The basin exists on a beta plane where f varies linearly around its value at ∼30°N.
The model has a regular 122 × 82 grid that is aligned with the rotated basin. The horizontal resolution is equivalent to a 1/4° grid at the equator and the configuration has 31 model levels. Two forms of bathymetry are used in this section. The FLAT configuration has a fixed depth of 4.5 km and no partial cells are used. The SLOPED configuration has a linear slope that extends from the North West side of the basin and spans half the basin (see Figure 4a). The maximum depth of the SLOPED configuration is 4.5 km, the minimum depth is 2 km, and partial cells are used to represent the slope.
The circulation is forced by sinusoidal analytic profiles of surface wind stress and buoyancy forcing. The wind stress is zonal and only varies with latitude so that the curl changes sign at 22°N and 36°N (see Figure 4b). The wind stress profile is designed to spin up a subpolar gyre in the north, a subtropical gyre in the south, and a small recirculation also emerges in the bottom corner. The net surface heat flux takes the form of a restoring to a prescribed apparent temperature. Further details about the buoyancy forcing can be found in Lévy et al. (2010). The wind stress and buoyancy forcing varies seasonally in a sinusoidal manner.
The model uses a free slip condition on all boundaries except at the bottom where a linear friction drag is applied. A simplified linear equation of state is used with a thermal expansion coefficient of a 0 = 2 × 10 −4 kg m −3 K −1 , and a haline coefficient of b 0 = 7.7 × 10 −4 kg m −3 psu −1 . Horizontal and biharmonic diffusion of momentum is implemented with a diffusivity of 5 × 10 10 m 4 s −1 . Biharmonic diffusion of tracers along isopycnals is implemented with a diffusivity of 10 9 m 4 s −1 (Lemarié, Debreu, et al., 2012;Madec et al., 2019).
The model is spun up for 60 years and the experiment was run for an additional 10 years with monthly mean outputs. A steady state is not required for the diagnostics to be valid as the time derivative term is present in the vorticity budget. A time step of 10 min is used for the model integration.
The EEN vorticity scheme is used for consistency with all analysis discussed in Section 3 and the results from the Weddell Gyre in Section 5. The EEN method calculates F cell thicknesses using the method described by Equation 13 and we therefore expect sudden changes in the F cell thickness near the domain edge for both the FLAT and SLOPED configurations.

Methods
Momentum diagnostics are calculated for every time step and the discrete vorticity diagnostics are calculated by depth-integrating the momentum diagnostics and taking the curl. The resultant diagnostics are time-averaged over the 10 year experimental period. The extensive time-averaging will influence the advection vorticity diagnostic as there is an added contribution from the eddy vorticity flux.
For contour integration, the vorticity diagnostics and depth-integrated stream function are then linearly interpolated onto a regular 1/12° grid. This is to minimize errors caused by the difference between the true enclosed streamline area and the total area of the enclosed F cells. Interpolation beyond 1/12° resolution makes little difference to the results, suggesting that the area error has been significantly suppressed.
For 1,001 values of ψ, closed streamline contours are identified using a marching squares algorithm from the scikit-image package (Van Der Walt et al., 2014). Streamlines that are near the recirculation gyre (south of 20°N) are ignored in this experiment and for some values of ψ no closed streamlines could be found. For each closed streamline found, the vorticity diagnostics are integrated over the area enclosed; this is equivalent to calculating I(ψ) in Equation 4 over many values of ψ. The freshwater fluxes mean that ∇ h ⋅ U ≠ 0 even in a steady state and an exact stream function cannot be calculated. To test how closely the calculated streamlines follow the circulation we integrate the positive quantity | 0 (∇ℎ ⋅ )| over the same enclosed areas to estimate the magnitude of the error caused by the divergent flow. The maximum value of f is used as f 0 and the largest contour integral of | 0 (∇ℎ ⋅ )| is 1.84 m 3 s −2 which is substantially smaller than the leading contour integrals presented in the next sub-section. In addition to this test we used an elliptical solver to calculate the Helmholtz decomposition of the depth-integrated velocity field (e.g., Marshall & Pillar, 2011); using the streamlines from the incompressible component does not change the results presented in the next sub-section.
Multiple closed contours can be found for the same value of ψ so an additional contour constraint is needed to ensure I(ψ) is single-valued. In this experiment we always choose the contour that spans the largest area which minimizes the influence of small pocket circulations that are not a part of the gyre. Closed streamlines that run along the edge of the domain can be hard to identify so a discontinuity in I(ψ) near ψ = 0 is expected as the largest detected contours will suddenly become pocket circulations as ψ approaches zero.

Results
The depth-integrated streamfunction from the FLAT and SLOPED configurations is shown in Figure 5. The vorticity of the depth-integrated velocity field is shown in Figure 6. In both configurations a subtropical and  subpolar gyre can clearly be identified and a small recirculation gyre can be found in the Southernmost corner. The subtropical gyre circulation is clockwise and the subpolar gyre circulation is anticlockwise.
In the FLAT configuration the subtropical gyre has a transport of 65 Sv and the subpolar gyre has a transport of 18 Sv. In the SLOPED configuration the subtropical gyre has a transport of 38 Sv and the subpolar gyre has a transport of 14 Sv. We note that the sloped bathymetry significantly alters the form of the subtropical gyre streamlines.
The depth-integrated vorticity diagnostics of the FLAT and SLOPED configuration are shown in Figures 7 and 8 respectively alongside the decomposition of the planetary vorticity diagnostic introduced in Section 3.4. In the FLAT configuration we note that the non-linear advection of vorticity and the planetary vorticity diagnostic have the largest grid point values (∼10 −9 ms −2 ) near the western boundary currents of both gyres. The wind stress curl is one order of magnitude smaller (∼10 −10 ms −2 ) but changes sign less frequently within the gyre regions. We see that the planetary vorticity diagnostic is almost entirely a result of the beta effect (Figures 7g and 7h). We note that the contribution from varying cell thicknesses in the FLAT configuration is non-zero and localized to the edge (Figure 7j) where the EEN Coriolis scheme artificially shrinks F cell thicknesses near masked points.
In the SLOPED configuration ( Figure 8) the advection and planetary vorticity diagnostics are still large but have an elongated structure similar to the SLOPED streamlines in Figure 5b. The bottom pressure torque is significant and is localized to the sloped region ( Figure 8b). The planetary vorticity diagnostic has a more complex decomposition as the influence of varying cell thicknesses extends beyond the edge of the domain and model level steps also contribute (Figure 8k).
The integrals of the vorticity diagnostics over areas enclosed by streamlines are shown in Figure 9 and Figure 10 for the FLAT and SLOPED configurations respectively as well as the integrals of the planetary vorticity diagnostic components. Example streamline contours are also shown. In these figures ψ > 0 describes the subtropical gyre and ψ < 0 describes the subpolar gyre. The subtropical and subpolar gyres circulate in the opposite direction but the sign of the integration results are adjusted so that positive integrals correspond to forces that spin the gyres up.
In the FLAT configuration we see that the subtropical and subpolar gyre are entirely driven by wind stress curl. For area integrations that envelop most of the subtropical gyre (small and positive values of ψ), the wind stress curl is largely balanced by the advection of relative vorticity. This balance implies a net import of positive vorticity into the gyre. The imported vorticity cannot originate from the subpolar gyre as the advection of relative vorticity plays no role in spinning the subpolar gyre down. The streamlines at the exterior of the gyre envelop both cells (maxima in ψ) of the subtropical gyre so the advection of vorticity between the cells is not a contribution to the signal. The imported vorticity must originate from the recirculation gyre in the southernmost corner. In the subtropical gyre interior (large positive values of ψ) the wind stress curl is largely balanced by the curl of bottom friction, matching the balance proposed by Niiler (1966).
The planetary vorticity diagnostic is significant in both of the FLAT gyres and is the dominant drag for the subpolar gyre. For area integrals that include the exterior of either gyre (small values of ψ), the integrated planetary vorticity diagnostic is a combination of a numerical beta effect originating from the discrete calculation of − ∇ℎ ⋅ ( ) and the influence of partial F cells that are artificially created by the EEN scheme. At the interior of both gyres (large values of ψ) the numerical beta effect is the only component.
In the SLOPED configuration we see that both the subtropical and subpolar gyre are almost entirely driven by wind stress curl. There is no dominant force spinning the gyres down. Advection, bottom pressure torques, lateral diffusion, bottom friction, and planetary vorticity all make a similar contribution to spinning the gyres down. The planetary vorticity diagnostic is similarly mixed as both the numerical beta effect and partial cells make up the signal. The gyres in the SLOPED configuration appear to be an intermediate case between a topographically steered gyre and an advective regime.
Spurious forces that emerge from the discrete Coriolis acceleration are significant in idealized models with and without variable bathymetry and appear to have a large influence on gyre circulations. In the next sub-section we see if these forces are also significant in a realistic global model.

Details of the Configuration
We now consider a more realistic configuration based on the NEMO global model with realistic forcing and bathymetry. In this experiment, we use an ocean-ice global configuration that is similar to that described in Storkey et al. (2018) but based on NEMO version 4. The global grid is based on the 'ORCA' family of grids within the NEMO framework (Madec et al., 2019). In this article we only consider the configuration using the ORCA025 grid (1/4° horizontal resolution at the equator). Most of the model bathymetry for ORCA025 is derived from the ETOPO1 data set (Amante & Eakins, 2009). Bathymetry on the Antarctic shelf is based on IBSCO (Arndt et al., 2013) and has been smoothed by three applications of a first order Shapiro filter. The bathymetry is represented in z-coordinates by partial cells (Barnier et al., 2006). Surface forcing is taken from the CORE2 surface forcing data set (Large & Yeager, 2009) and includes contributions from sea ice. The bathymetry is shown in Figure 11a.
The model uses a free slip lateral boundary condition with a quadratic drag along the bottom boundaries and the TEOS-10 equation of state (McDougall & Barker, 2011). Biharmonic diffusion of momentum is implemented and acts along model level surfaces with a diffusivity that varies with local horizontal grid spacing (Willebrand et al., 2001). Laplacian diffusion of tracers is implemented and acts along isopycnal surfaces with a diffusivity that also varies with local horizontal grid spacing. The EEN vorticity scheme is used again for consistency with the analysis in Section 3 and the results in Section 4.

Methods
The methods used for calculating the depth-integrated streamfunction, vorticity diagnostics, and contour integrals are identical to those described in Section 4.2. We study the area including and surrounding the Weddell Gyre in the model (see Figure 11) and consider the time-averaged fields over a typical year. The stream function is interpolated onto a regular 1/12° grid and closed contours are identified for 201 values of ψ. Interpolating beyond 1/12° resolution makes little difference to the results, suggesting that any area errors have been  significantly suppressed. We test how closely the calculated streamlines follow the circulation by integrating the positive quantity | 0 (∇ℎ ⋅ )| over the same enclosed areas to estimate the magnitude of the error caused by the divergent flow. The maximum value of |f| is used as f 0 and the largest contour integral of | 0 (∇ℎ ⋅ )| is 19.52 m 3 s −2 which is substantially smaller than the leading contour integrals presented in the next sub-section. In addition to this test we used an elliptical solver to calculate the Helmholtz decomposition of the depth-integrated velocity field; using the streamlines from the incompressible component does not change the results presented in the next sub-section.
As we are studying a one gyre system we choose to only identify contours where ψ > 0. This effectively filters out the vorticity budget of closed circulations in the Antarctic Circumpolar Current. The sign of the integration results are adjusted so that positive integrals correspond to forces that spin the Weddell Gyre up.

Results
The depth-integrated streamfunction of the Weddell Gyre is shown in Figure 11b and it can be seen that the Weddell Gyre has a transport of 60 Sv. The streamlines follow the isobaths closely suggesting the circulation is largely constrained by the bathymetry. The vorticity of the depth-integrated velocity field is shown in Figure 12. The depth-integrated vorticity diagnostics are shown in Figure 13. The fields shown in Figure 13 have been smoothed using 25 point nearest neighbor averaging over a local 5 × 5 grid. The contribution from model level steps (Figure 13k) has not been smoothed to show that it is localized to specific lines where the number of model levels change. The combined effect of the wind stress and stress due to sea ice are shown in Figure 13e. With realistic topography and forcing, the grid point values of depth-integrated vorticity diagnostics are very noisy (even when smoothed) with the exception of the surface stress curl. This highlights how important it is to integrate the vorticity diagnostics when interpreting them. For individual grid points we see that the planetary vorticity diagnostic is made up of contributions from the beta effect, partial cells, and a significant contribution from model level steps. The beta effect is the most coherent of the contributions and is mostly negative in the western limb of the gyre where V > 0 and positive in the eastern limb where V < 0. As expected, the contribution from model levels steps is localized to areas where the number of model levels change.
Unlike in the double gyre model, bottom friction appears to be small and incoherent in the Weddell Gyre region and is unlikely to have any significant influence on the vorticity budget. The total time tendency (Figure 13d) is non-zero in this vorticity budget suggesting that the model is not in a completely steady state; however, the grid point values are only significant in the Drake Passage and are noisy.
The integrals of the depth-integrated vorticity diagnostics over areas enclosed by streamlines are shown in Figure 14 alongside integrations of the planetary vorticity components. We see that the Weddell Gyre is almost entirely spun up by the wind stress curl. The stress due to sea ice (marked by hatching in Figure 14a) and the advection of relative vorticity also help to spin the Weddell Gyre up. The advective contribution is caused by vorticity exchange at the interface between the Weddell Gyre and the ACC.
Bottom pressure torques and lateral diffusion play a notable role in spinning the Weddell Gyre down but the planetary vorticity diagnostic is the most significant contribution. Looking at the decomposition of the planetary vorticity diagnostic we see that the signal is mostly determined by changes in model level and the remainder is determined by variations in cell thickness. This suggests that the Weddell Gyre is almost entirely spun down by topography due to the combined effect of bottom pressure torques and the planetary vorticity diagnostic, but the majority of the gyre's interaction with the sea floor is spurious. This conclusion is true in both the interior and exterior of the gyre.
The results in Figure 14 are concerning as they suggest that the vorticity input from the realistic surface stresses is largely balanced by spurious topographic accelerations. In an area of the ocean with such strong bathymetric  feratures, it is not surprising that topographic forcing is important but we would expect the topographic accelerations from a realistic bathymetry to also be realistic. Instead, the dominant component of the topographic forcing is a spurious acceleration that is localized to discrete lines where the number of model levels change (see Figure 13k) and arises from the masking of the non-topographic Coriolis acceleration. This suggests that the partial cell representation of the sea floor is not providing realistic topographic forcing in the Weddell Gyre region.

Discussion
We have shown that the vorticity dynamics of both highly idealized and realistic gyre configurations are greatly influenced by spurious forces that emerge from the discrete Coriolis force and the step-like representation of bathymetry. In the idealized double gyre configuration (Section 4) the spurious force is a combination of numerical beta and topographic effects that are present in both the FLAT and SLOPED configuration. In the realistic Weddell Gyre (Section 5) the spurious force is the dominant drag and is entirely determined by model level steps and partial cells. In this section we discuss possible methods to mitigate these spurious forces.

Alternative Vorticity Schemes
The results presented in Sections 4 and 5 both use the EEN vorticity scheme and it is tempting to dismiss the spurious forces as an artifact of the selected scheme. The analysis in Section 3.2 is general for three popular schemes: EEN, ENE, and ENS. The methods and decomposition used in this article are applicable under any scheme where the Coriolis acceleration can be expressed in the form of Equations 8 and 9. Results from the SLOPED double gyre configuration using the different schemes are presented in Appendix C and the vorticity budgets are qualitatively similar. Spurious topographic forces and the numerical beta effect are still significant. It therefore seems that switching between the available vorticity schemes will not alleviate the spurious signal. It is possible that a new scheme could be formulated which is designed to significantly reduce the spurious forces, but that will most likely require abandoning the conserved quantities that characterize the existing schemes.

Alternative Depth-Integrated Vorticity Equations
In Section 2.1 we derived a depth-integrated vorticity equation by taking the curl of the depth-integrated momentum equation and we calculated the model vorticity diagnostics using the equivalent discrete method. As discussed in Section 2.1, there are alternative formulations of the depth-integrated vorticity equations with different physical meanings. An accurate model should be able to represent all forms of the depth-integrated vorticity budget so switching between formulations does not alleviate any spurious forces, but it is interesting to see if any of the spurious contributions in this article can spill over into other vorticity budgets.
If we derive a continuous depth-integrated vorticity equation by depth-integrating the curl of the momentum equations then the Coriolis acceleration emerges in the vorticity budget as: where u t and u b are the horizontal velocities at the free surface and sea floor respectively. When compared with Equation 2 we can see that the planetary vorticity term has an additional topographic and free surface term. The second term on the right hand side of Equation 24 describes a vortex stretching acting on the vertical velocity induced by the bottom topography. In configurations with no variable bathymetry and small variations in the free surface, the order of taking the curl and depth-integrating no longer affects the vorticity budget so the non-topographic spurious forces identified in this article will remain in either formulation.
To calculate the discrete curl of a horizontal vector field near the bathymetry we need to make an assumption about how the along-slope component varies as it approaches the edge of the domain. We can assume either a free slip or no slip boundary condition by using a ghost point that mirrors the location of the closest grid point into the bathymetry. For a free slip boundary condition the ghost point value matches the closest grid point value, F ∥ ; for a no slip boundary condition the ghost point value will be the negative of the closest grid point value, −F ∥ . A partial slip boundary condition also exists where the value of the ghost point will be between −F ∥ and F ∥ .
Let us return to the simple flow introduced in Section 3.3 and illustrated in Figure 2 but this time when we calculate the planetary vorticity diagnostic we will calculate the curl of the Coriolis acceleration on each model level and then depth-integrate. For the lower level, the horizontal flow is entirely in the x direction so there is a zero along-slope component of the Coriolis acceleration near the bathymetry (F ∥ = 0). This means that if a free slip, no-slip, or partial slip boundary condition are used the ghost point value will be zero and the curl of the Coriolis force (centered on the purple cross in Figure 2) will be zero in all three cases. As all vorticity generation takes place in the upper level, the planetary vorticity diagnostic is the same if we take the curl before or after depth-integrating (Equation 21) and the effect of model level steps can exist in either vorticity budget.
The result of Equation 21 can be interpreted as a vortex stretching acting on the vertical velocity that is induced by the change in horizontal velocity u 1 (see Figure 2). The vertical velocity seems unlikely to originate from topographic upwelling as there is no flow in the y direction. This fact combined with the ambiguity of ∇H at model level steps means we would advise caution before comparing the discrete vortex stretching that originates from model level steps to the analytic vortex stretching in Equation 24.

The B-Grid
Altering the grid geometry can significantly change the behavior of model forces. To highlight this we consider how the Coriolis force behaves on the B-grid. The B-grid excels at representing geostrophic flows as u, and v are located on the same vector point. The streamfunction and relative vorticity are located on the tracer point as shown in Figure 15.
On the B-grid the Coriolis acceleration is simply: The Coriolis acceleration does not rely on multi-point averaging or thickness weighting of f so numerical contributions do not emerge in the grid point acceleration.
On the B-grid u and v lie on the same point so they share the same mask. This means that non-zero Coriolis accelerations are never masked near model level steps and the depth-integrated Coriolis acceleration is a function of the depth-integrated velocities only: We therefore conclude that the spurious force caused by model level steps on the C-grid (see Section 3.3) is not present on the B-grid. The corresponding planetary vorticity diagnostic is equal to −∇ h ⋅ (fU)| i,j calculated over a single tracer cell.
Calculating the curl on a B-grid is consistent with Stokes' law applied to a tracer cell but the vector information is found on the corners of the cell. As the stream function is defined on the tracer point we can approximate that the area enclosed by a streamline is a collection of interior tracer cells. Similarly to the C-grid case in Section 3.5 this is an approximation as we are assuming that the streamline follows the rectangular edges of the interior tracer cells so interpolation may be required to remove any significant area error. Unlike the C-grid case, the planetary vorticity diagnostic is equal to −∇ h ⋅ (f U)| i,j calculated over a single tracer cell. Therefore, the area integral of the planetary vorticity diagnostic will satisfy the divergence theorem applied to the internal tracer cells. It seems likely that this discrete integral may vanish on a sufficiently fine grid but further investigation with idealized and realistic streamlines is needed.
Using the B-grid would remove all of the spurious topographic forces identified in this article. This highlights how a model circulation's interaction with the sea floor is significantly affected by the grid geometry.

Terrain-Following Coordinates
The spurious topographic effects found in this article are a consequence of how bottom topography is represented in z-coordinates. In the Weddell Gyre especially we see how model level steps can create large spurious contributions to the depth-integrated vorticity budget.
Terrain-following coordinates (or σ-coordinates) are an alternative form of vertical coordinate where the vertical resolution adjusts with the bottom topography so that the same number of model levels are present in all fluid columns (Song & Haidvogel, 1994). σ-coordinates are used in Stewart et al. (2021), Schoonover et al. (2016), and Jackson et al. (2006) and have the advantage of removing spurious terms that emerge from model level steps. The forms of the EEN, ENE, and ENS vorticity schemes are unchanged when using terrain-following coordinates so the horizontal variations in cell thicknesses could still cause a spurious signal.
Terrain-following coordinates are not used widely in climate models because of the difficulty in calculating accurate horizontal pressure gradients (near the equator), advection, and isoneutral tracer diffusion. A full discussion of the current advantages and limitations of terrain following coordinates can be found in Lemarié, Kurian, et al. (2012).

Isopycnal Coordinates and the Vertical Lagrangian-Remap Method
In isopycnal C-grid models, where density is used as a vertical coordinate, cell thicknesses still vary and in models with many density layers the model levels are free to incrop to the sea floor. The forms of the EEN, ENS, and ENE schemes are unchanged when using density coordinates so the spurious signals in the planetary vorticity diagnostic seem to be possible. In configurations where density layers infrequently incrop to the sea floor, the effect of model level steps will be significantly suppressed as the grid is approaching the limit of a terrain-following coordinate system.
In C-grid models that use the vertical Lagrangian-remap method (Adcroft et al., 2019;Bleck, 2002) the vertical coordinate evolves with the flow and is then conservatively remapped onto a target grid (see Griffies et al. (2020) for a review). The forms of the EEN, ENS, and ENE schemes are unchanged when using this method. If the target coordinate grid still has horizontal variations in cell thicknesses and incrops with the sea floor, we would expect there to be spurious topographic interactions with the sea floor. It is possible that in areas of topographic upwelling the effect of model level steps could be reduced as Coriolis accelerations near the bathymetry are elevated by the vertical motion and are partially projected onto unmasked points when remapped onto the target grid.

Summary
The depth-integrated vorticity budget is a valuable tool for identifying important model forces in gyre circulations. Vorticity diagnostics can be integrated over the area enclosed by streamlines to identify forces responsible for spinning the gyre up and down. By considering how the vorticity budget is represented on a C-grid with step-like bathymetry we identified spurious forces that emerge from the representation of bottom topography and the discrete Coriolis acceleration. Model level steps and partial cells produce two distinct spurious topographic forces. In the absence of bottom topography, it is shown that the discrete planetary vorticity term does not generally vanish when integrated over the discrete area enclosed by a streamline. This suggests that a spurious non-topographic force, described as a numerical beta effect, is also present.
We first studied the vorticity budget of an idealized double gyre configuration with analytic geometry, forcing, and two bathymetry options. The FLAT variant has a constant depth and the SLOPED variant has a linear slope that extends over half the domain. The subtropical gyre of the FLAT configuration is non-linear at the exterior (wind stress curl balanced by advection) and is in a Stommel (1948) regime in the interior (wind stress curl balanced by friction). The FLAT subpolar gyre is spun up by wind stress curl and mostly spun down by spurious forces found in the planetary vorticity diagnostic. Spurious forces are significant in both FLAT gyres and are a consequence of the numerical beta effect and partial F cells that are artificially introduced by the EEN vorticity scheme. Artificial partial F cells would not be present in the ENS or ENE vorticity schemes.
The vorticity budget of the SLOPED gyres features bottom pressure torques and an increased influence of partial cells on the planetary vorticity diagnostic. The SLOPED subtropical gyre is an intermediate case between a topographically steered gyre and a non-linear circulation. The SLOPED subpolar gyre is driven by wind stress curl but spun down by the combined effect of bottom pressure torques and spurious interactions with the topography via partial cells. This first case study highlighted how spurious terms can dominate a vorticity budget in idealized configurations with and without variable bathymetry.
The second case study was the Weddell Gyre in a global model where the forcing and geometry are more realistic. By studying the vorticity budget of the Weddell Gyre we conclude that the model circulation is mostly spun up by wind stress curl and spun down by the combined effect of bottom pressure torques and spurious interactions with the topography. The largest of the topographic forces spinning the Weddell Gyre down is the spurious and unrealistic force caused by model level steps.
Switching to alternative vorticity schemes is not effective at reducing spurious contributions to the vorticity budget. By presenting a general form of the discrete Coriolis acceleration we are able to quickly conclude that the topographic and non-topographic spurious forces will remain under all three vorticity schemes and any other scheme that uses this general form. The influence of model level steps is a direct consequence of the C-grid geometry when using vertical coordinates that intersect the bathymetry and is relatively insensitive to the choice of vorticity scheme.
Altering the geometry of the discretization is an effective method for reducing spurious topographic forces. The B-grid is better at representing the Coriolis force and it is not possible for model level steps or partial cells to influence the Coriolis acceleration. Model level steps and their influence on the Coriolis acceleration can be avoided altogether by using terrain-following coordinates.
The B-grid and terrain-following coordinates have their own unique limitations and it is unclear how much the identified spurious forces corrupt circulation variables such as the gyre transport. It is possible that the spurious forces are inadvertently performing the role of one or more real ocean processes that are required for accurate simulations. If a combination of non-spurious forces can fully account for the spurious forces found in this article then the identified problem is purely diagnostic in nature. Otherwise, any part of the spurious forcing that cannot be accounted for by non-spurious forces should be considered as a numerical error. This numerical error could be small but may also accumulate under specific conditions and corrupt model circulations. The spurious cooling (Hecht, 2010) that occurs when a dispersive advection scheme is used with the Gent and McWilliams (1990) eddy parametrization highlights the dangers of ignoring numerical errors.
It is also possible that other model forces contain spurious contributions that have not been uncovered in this article. These contributions could be significant and may have the potential to cancel the spurious effects found in this article. When looking at the integrated diagnostics in Figures 9, 10 and 14 we see that usually the only model force with an opposite contribution to the Coriolis force that is large enough to cancel the found spurious effects is the surface stress. It seems unlikely that the surface stress contains spurious contributions that are closely tied to bathymetry and the Coriolis parameter.
It is important for the ocean modeling community to continue developing new ways of representing bathymetry and we hope that vorticity budgets and the diagnostic method presented in this article will provide a valuable tool for assessing and quantifying representations of the sea floor in current and future ocean models.

Appendix A: Deriving the Depth-Integrated Vorticity Equation
Here we derive the depth-integrated vorticity equation (Equation 2) including the omitted contributions from surface undulations and atmospheric pressure torques. We start from the vector invariant form of the momentum equation, Which has already been introduced in Section 2.1. To derive the depth-integrated vorticity equation, we must first depth-integrate the equation and then calculate the vertical component of the curl. In this appendix, we consider how each term in Equation A1 is transformed by this operation. When depth-integrating the time derivative term in Equation A1, we must respect the time dependency of the free surface, η. We therefore use the Leibniz integration rule, Where the second term on the right hand side of Equation A2 is the contribution from free surface undulations.
The non-linear term in Equation A1 can be rewritten as, The non-linear term emerges as the advection term in the depth-integrated vorticity equation and we note that, Similarly the curl of the depth-integrated Coriolis acceleration is the planetary vorticity term, When depth-integrating the pressure gradient in Equation A1, we must respect the x and y dependency of the sea floor and the free surface. We therefore use the Leibniz integration rule, Where P a is the atmospheric pressure at the free surface. The second term on the right hand side of Equation A6 is the atmospheric pressure torque.
The surface forcing term in Equation A1 emerges as the difference between the curl of the top and bottom stresses, And the diffusion term emerges as  , we note that each term in the ENE and ENS forms can be written in the general form of Equations 8 and 9 as 1 =̃∕ 3 and 2 =̃∕ 3 . In the ENE and ENS cases 3 ( ) = 3 ( ) in Equations 8 and 9.
In the EEN vorticity scheme, the x and y components of the Coriolis acceleration are: where F NE , F NW , F SE , and F SW are thickness-weighted triads of the Coriolis parameter: Where ̃= ∕ 3 using the EEN definition of e 3f shown in Equation 13.
To calculate the planetary vorticity diagnostic we take the curl of the depth-integrated Coriolis acceleration using Equations 15 and 22. In general the resulting equation of the vorticity diagnostic is very difficult to interpret. We only present the form of the planetary vorticity diagnostic for the EEN scheme on a grid with no partial cells or model level steps as it is used to derive the numerical beta effect in Section 3.5:

Appendix C: Alternative Vorticity Schemes in the Double Gyre Model
In this section we present various integrations of the SLOPED double gyre configuration using different vorticity schemes: EEN, ENS, and ENE. All other aspects of the experiment are as described in Section 4.1. The results are shown in Figure C1. The vorticity budget is qualitatively similar between the three cases as well as the decompo sition of the planetary vorticity diagnostic. It should be noted that the circulations do differ as the transports vary and the separation points of the western boundary currents change. Figure C1. Stacked area plots showing the integrals of depth-integrated vorticity diagnostics for the SLOPED configuration (time-averaged) using the EEN, ENE, and ENS vorticity schemes. Positive values correspond to a force that spins the subtropical (ψ > 0) or subpolar (ψ < 0) gyre up. A decomposition of the planetary vorticity diagnostic integrals are given on the right (b,d,f). Figure D1. Stacked area plots showing the integrals of depth-integrated vorticity diagnostics (time-averaged) for the FLAT configuration without using interpolated fields. Positive values correspond to a force that spins the subtropical (ψ > 0) or subpolar (ψ < 0) gyre up. (b) Shows the area integrals of the planetary vorticity diagnostic and its components. The vorticity budget and decomposition are qualitatively similar to that shown in Figure 9.