Sliding or Stumbling on the Staircase: Numerics of Ocean Circulation Along Piecewise‐Constant Coastlines

Coastlines in most ocean general circulation models are piecewise constant. Accurate representation of boundary currents along staircase‐like coastlines is a long‐standing issue in ocean modeling. Pioneering work by Adcroft and Marshall (1998, https://doi.org/10.3402/tellusa.v50i1.14514) revealed that artificial indentation of model coastlines, obtained by rotating the numerical mesh within an idealized square basin, generates a spurious form drag that slows down the circulation. Here, we revisit this problem and show how this spurious drag may be eliminated. First, we find that physical convergence to spatial resolution (i.e., the main characteristics of the flow are insensitive to the increase of the mesh resolution) allows simulations to become independent of the mesh orientation. An advection scheme with a wider stencil also reduces sensitivity to mesh orientation from coarser resolution. Second, we show that indented coastlines behave as straight and slippery shores when a true mirror boundary condition on the flow is imposed. This finding applies to both symmetric and rotational‐divergence formulations of the stress tensor, and to both flux and vector‐invariant forms of the equations. Finally, we demonstrate that the detachment of a vortex flowing past an outgoing corner of the coastline is missed with a free‐slip (zero vorticity) condition at the corner. These results provide guidance for a better numerical treatment of coastlines (and isobaths) in ocean general circulation models.

The choice of lateral boundary condition remains a thorny issue in this regard. It is clear that the "no-slip" condition is the physically founded condition when representing interactions of a fluid with solid boundaries (Richardson, 1973). However it is far from obvious how this fundamental constraint should be translated in OGCMs (Deremble et al., 2011), whose typical spatial resolution (order 10-100 km in the horizontal) does not permit to resolve the viscous boundary layer. In fact, the small aspect ratio of the ocean implies that the viscous boundary layer sees a gently sloping seafloor rather than vertical sidewalls, suggesting that the lateral boundary condition is a numerical requirement lacking physical underpinning. For example, in z-coordinate OGCMs, the sloping seafloor is numerically transformed into a staircase boundary that alternates flat bottom and vertical walls. In this case, combined use of a bottom drag and a no-slip condition at the sidewall may potentially overestimate the frictional boundary effects. In practice, slippery or mixed lateral boundary conditions, such as "free-slip" or "partial slip," are often preferred in coarse OGCMs (Penduff et al., 2007). Deremble et al. (2011) provide some theoretical support for this approach by showing that a free-slip condition yields better agreement with classical turbulence theory than does a no-slip condition. However, they also caution that lateral viscous conditions do not capture all the unresolved physics of the turbulent boundary layer. In general it is not known how to best define boundary conditions along the staircase-like frontiers of OGCMs.
The problem of defining suitable boundary conditions is made even more complex by the (unintended) impacts on the simulated circulation of artificial steps in the lateral boundary. Using an idealized square basin configuration and various orientations of the numerical mesh, Adcroft and Marshall (1998) (hereafter AM98) first pointed out that artificial indentation of a shoreline systematically causes a spurious form drag that slows down the coastal flow. AM98 showed that this drag depends on the numerical formulation of diffusive stresses but always exists, including when a "free slip" boundary condition is implemented. Subsequent investigations within similar configurations showed that the response of the flow to staircase-like coastlines depends on the advection and diffusion schemes (Dupont et al., 2003). These studies also suggested that the spurious drag may persist at higher resolution because coastal steps increase in number as they become smaller (Adcroft & Marshall, 1998;Dupont et al., 2003). Griffiths (2013) argued that the adverse effects of an indented coastline could be handled by implementing the right impermeability condition corresponding to the real coastline. However their findings apply only to traveling Kelvin waves, and it remains unclear how the long-standing issue of spurious form drag should be addressed in OGCMs.
The response to an isolated, large-scale bend in the coastline also deserves attention because real swerves of the shoreline can have pronounced impacts on boundary currents (Magaldi et al., 2008;Warner & MacCready, 2009), and because such swerves may be sharpened in their discrete model representation. Deremble et al. (2016) showed that an outgoing corner in the coastline induces retroflection of a coastal stream. Their findings echo those of Dupont and Straub (2004), who varied the curvature of a wavy wall configuration and found that opposite vorticity filaments were created at the coastal tips and caused detachment of the flow. In both studies, only some numerical formulations enable to capture the expected physical behavior. The reasons for this strong sensitivity to numerical formulations are not fully elucidated.
Here, we address these questions using idealized model configurations. We first reproduce the configuration of AM98 (Section 2) and demonstrate that the circulation is insensitive to the mesh orientation-and related coastline indentation-provided that the model is physically converged to spatial resolution (Section 3). Next, we expose how a true mirror boundary condition on the flow renders the indented coastline as slippery as a straight, free-slip frontier (Section 4). The numerical representation of an isolated, large-scale step in the coastline is Figure 1. Schematic illustrating the coarse-graining of a coastline on a structured mesh. The real coastline (thin brown line) presents an infinite number of details that are too small, relative to the size of grid cells (thin black), to be represented by the numerical model. It must be averaged over at least two grid points to avoid noise generation, resulting in a smooth coastline (red curve). Then the land mask (gray cells) is defined from the projection of the smooth coastline onto the grid, creating artificial abrupt changes in coastline direction. investigated in Section 5 using the configuration of Deremble et al. (2016). We summarize our findings and recommendations in Section 6.

Configuration
We consider the same problem as AM98: a shallow water model with reduced gravity is solved in a square basin of size L = 2,000 km. An anticyclonic wind stress τ = (−τ 0 cos(πy/L), 0) is forcing the active layer. The coordinate system (x, y) has its origin in the lower left corner of the square basin. Equations in vector invariant form read ℎ + div (ℎ ) = 0 (1) where h is the active layer thickness, u = (u, v) represents the horizontal velocity vector, = . (∇ × ) the relative vorticity, k the vertical unit vector, f the Coriolis parameter, g′ the reduced gravity, r the friction coefficient, D ν the diffusion term and ρ 0 the density.
Two formulations of the diffusion term (D ν where ν is the lateral viscosity) are considered, that will lead to different discretisations: the rotational-divergence (hereafter called "rot-div") form (Madec et al., 1991) and the symmetric form (Griffies & Hallberg, 2000). The rot-div form is calculated as ∇( ) − 1 ℎ ∇ × ( ℎ ) where transport divergence χ is defined as = 1 ℎ ( (ℎ ) + (ℎ )) . The symmetric form is expressed instead as 1 ℎ ∇ ⋅ ( ℎ ) (Gent, 1993) A viscous boundary condition (see Section 2.3) is applied at the wall of the basin by enforcing the values of χ and ζ in the "rot-div" form and D S and D T in the symmetric form. In the discrete form, the same viscous boundary condition can thus be written differently depending on the formulation of D ν .
We will also perform simulations using the flux form of the shallow water equations. In this case, Equation 2 rewrites as

Numerical Discretization
We use for the numerical simulations the shallow water option (SWE) introduced in the version 4.2 of the NEMO general circulation model (Madec & NEMO System Team, 2022). The initial state is at rest with h 0 = 500 m the uniform thickness of the active layer. Experiments are integrated over 25 years to achieve a steady state on h. All the simulations presented here use the Leap Frog Robert-Asselin time-stepping scheme (Leclair & Madec, 2009). In contrast, AM98 used the third order Adams-Bashforth (III) scheme. We do not expect the spurious form drag to be sensitive to the order of accuracy of the time-stepping scheme. We performed sensitivity tests with the third order Runge-Kutta method and found no noticeable differences in the equilibrium solutions (not shown).
Following AM98, we employ a Cartesian mesh with uniform resolution. Spatial resolution will be varied from 1/4° to 1/48°. In the reference case, the mesh is aligned with the edges of the square basin, so that the model coastlines are perfectly straight (Figure 2, left). In the rotated cases, the mesh is oriented at some angle (up to 45°) with respect to the physical coastline, so that steps punctuate the model coastlines (Figure 2, right). Hence, Figure 2. Effects of rotating the numerical mesh within an idealized square basin. The mesh is represented by the black lines. Land is shaded in gray, the oceanic domain in white. The red thick line represents the physical coastline.
On the left, the grid is aligned with the basin. On the right, the mesh is rotated so that artificial steps appear in the model shoreline. This figure is adapted from AM98.
the comparison of aligned and misaligned cases allows us to assess potential drag effects of artificial coastal indentation. Both physics (e.g., wind stress, Coriolis parameter) and grid-cell size are kept unchanged when the mesh turns in relation to the physical basin, so that only the shape of the model coastline changes. Physical and numerical parameters are similar to AM98 and listed in Appendix A.
All the simulations solved in vector-invariant form and shown in this study use the potential-enstrophy conserving vorticity scheme (called ENS) (Sadourny, 1975), except in Section 4.4 where we tested the scheme (called EEN) developed by Sadourny which conserves kinetic energy and-provided there is no divergence in the flowpotential enstrophy (Burridge & Haseler, 1977). Tests showed that solutions using the kinetic energy conserving advection scheme (called ENE) (Sadourny, 1975) behave very similarly to the ones shown here with the ENS scheme. Note that AM98 used the "vorticity" scheme given by Bleck and Boudra (1986), which is similar to ENS except for the presence of the vertical scale factor h at the numerator and denominator in ENS (as required to effectively conserve a discrete expression of potential enstrophy, Sadourny, 1975). Additional tests showed that "vorticity" and ENS schemes yield solutions that are slightly different but behave similarly in the presence of staircase-like coastlines. The expression of the various discretization schemes is given in Appendix B.
This configuration assumes an idealistic topography made of vertical walls and a flat bottom. However, the effective topography can be distorted by the traditional estimation of the thickness h at U, V, and F boundary nodes. In the open ocean on a C-grid, h is naturally defined at the center of each cell and calculated as a two-point (or four-point) average at velocity (or vorticity) nodes, for discrete conservation of properties. At the boundary, using the same definition and averaging with masked h can be equivalent to imposing sub-grid-scale topography, and can reinforce topostrophy. Tests under ENS, ENE, and EEN advection schemes showed that simulations are sensitive to the treatment of h at the boundary (not shown), especially with EEN. We chose to calculate boundary h as the masked average of the surrounding masked heights, in order to actually represent vertical walls.

Boundary Condition
The system of Equations 1 and 2 or Equations 1-4 requires two horizontal boundary conditions to be well posed.
The standard conditions on a solid wall consist of the impermeability condition u.n = 0, with n the coast-normal unit vector, and of a slipperiness condition that is a simplified representation of the effects of a viscous boundary layer.
AM98 considered two types of slipperiness condition (hereafter called viscous boundary condition): no-slip and free-slip. No-slip requires the tangential speed to be zero at the boundary which is u.t = 0, with t the coast-tangential unit vector. Combined with the impermeability condition (u.n = 0), no-slip thus entails u = 0 at the border. By contrast, free-slip is the absence of shearing and hence dissipation at the coast. It is defined as the absence of coast-normal shear at the border ∂u/∂n = 0 and can be interpreted as a mirror condition across the border, where virtual flow within land mirrors the oceanic flow.
It is possible to deduce from the viscous boundary condition the value of the vorticity at the boundary. In their Section 2, Verron and Blayo (1996) write the vorticity ζ c at an impermeable boundary regular enough to define in the local frame the vectors (n, t), where n is directed outward the basin: with κ the local curvature of the coastline (note that u is a vector and ∂u/∂n = ∇u ⋅n with ∇u the 2 × 2 Jacobian matrix of u). When the coastline is straight (κ tends to 0), the free-slip boundary condition reduces to ζ c = 0, while the no-slip boundary condition computes ζ c = −(∂u/∂n).t using u.t = 0. When dealing with a numerical model, depending on the numerical grid and discretization schemes, the evaluation of a quantity on a grid-point near the boundary may require the use of another quantity at the boundary. It is for instance the case of ζ c in the non-linear term in Equation 2, and the case of the rot-div formulation of the diffusive term. The boundary conditions therefore provide this information.

Influence of Resolution
The numerical parameters chosen by AM98 are those of an eddy-permitting OGCM: horizontal resolution is 1/4°, corresponding to a grid spacing Δx = Δy = 25 km. At this resolution, neither the internal radius of deformation R nor the boundary layers are properly resolved throughout the basin (Hallberg, 2013). Indeed, in the initial state, the deformation radius is about 35 km at the northern boundary, so that Δx ∼ R. It is therefore expected that increasing spatial resolution while keeping the same values for physical parameters (i.e., viscosity ν and friction r) will give solutions that differ from the reference solution of AM98. Figure 3 shows the steady solution using the vector-invariant form of equations with the ENS advection scheme, the rot-div stress tensor and the free-slip boundary condition. From top to bottom, spatial resolution increases successively from 1/4° to 1/8°, 1/16° and 1/32°. The mesh is either aligned with the coastline (left column) or rotated by 45° (right column). The 45° angle creates an artificially indented coastline, as illustrated in Figure 2. Shading and isolines depict the layer thickness h. In all cases, we obtain an anticyclonic (clockwise) circulation composed of two connected cells: a relatively weak Sverdrup interior intensified at the western boundary in the southern 1,500 km of the domain, and an inertial recirculation sub-gyre confined to the northernmost 500 km. We find that some high resolution simulations vacillate in the eastern part of the recirculation cell; that is, the simulated flow displays an oscillatory behavior in this region over a timescale of about 18 months. Such vacillation is thought to occur in a very restrained parameter range (Holland & Haidvogel, 1981). To ensure consistent comparisons, all shown free-slip solutions are averaged over the final 5 years.
Panels (a) and (b), corresponding to a resolution of 1/4°, reproduce the results of AM98. At this resolution, the inertial recirculation cell extends all the way to the eastern boundary when the mesh is aligned with the physical coastline, whereas it occupies only the western half of the basin when the mesh is rotated by 45°. Hence, the free-slip circulation seems to dwindle as the mesh turns. However, the latter statement is no longer true with a fine spatial resolution. At 1/32°, the two solutions are virtually identical (Figures 3g and 3h). As apparent in the evolution of the shape and maximum of the northern recirculation, both aligned (Figures 3a, 3c, 3e, and 3g) and turned (Figures 3b, 3d, 3f, and 3h) solutions appear to tend toward the same state. Besides, solutions at 1/16° look very similar to the ones obtained at 1/32° meaning that the solutions are physically converged from 1/16° (i.e., the main characteristics of the flow are almost insensitive to the increase of the mesh resolution).
The inertial recirculation sub-gyre is the most sensitive feature to resolution and mesh orientation. Hence, we choose to quantify the model sensitivity to the rotation and resolution of the grid using two diagnostics: the zonal extension of the inertial cell, and its overall intensity (calculated as the maximum of the active layer thickness within this cell). Figure

Preserving 1/4° Staircase Steps
When increasing spatial resolution, the size of the coastal steps decreases as their number increases along the coastline. Is the insensitivity to mesh orientation at high resolution due to the smaller step size? To answer this question, we performed 1/16° simulations with exaggerated coastal indentation (identical to its shape at 1/4° resolution). In this way, we maintain the broken aspect of the shoreline unchanged while reducing the grid spacing. We find that the final free-slip solution ( Figure 5a) is insensitive to the larger steps: it is almost identical to the reference aligned solution ( Figure 3e). Hence, physical convergence allows the simulated circulation to be insensitive to mesh orientation and coastal indentation, irrespective of the size of coastal indents.

Condition for Physical Convergence
The zonal development of the northern recirculation cell results from a non-linear interaction between the sub-gyre and its mirror recirculation induced by the free-slip boundary condition (e.g., Figure 1 of Cessi (1991)).
The northern region is also where the deformation radius is the smallest. This radius is about 35 km at the northern boundary in the initial state and it decreases over time because the westerly wind stress causes upwelling along the northern coast (thus shrinking h). Therefore, the deformation radius near the northern boundary is never properly resolved on a 1/4° mesh. To assess whether resolution of the deformation radius along the northern frontier is key to obtain physical convergence, we performed experiments where the aligned, 1/4° mesh is refined in a narrow northern band. Specifically, we reduce the meridional grid spacing from 1/4° 300 km offshore to 1/16° at the northernmost grid cells. Figure 5b shows that refining the meridional grid spacing near the northern boundary contracts the inertial recirculation sub-gyre (compare with Figure 3a). Local grid refinement thus suffices to bring this recirculation cell closer to the physically converged solution ( Figure 3e) than to the initial 1/4° solution ( Figure 3a). Analogous sensitivity tests with local mesh refinement close to the western, southern or eastern coast showed very little impact on the solution (not shown). Resolving the northern deformation radius, hence the mirror interaction at the north boundary, appears to be the key ingredient for physical convergence. We infer that a minimum of four grid points per deformation radius is necessary.

No-Slip Boundary Condition
All experiments up to here have been conducted using the free-slip boundary condition. The no-slip boundary condition may be expected to generate weaker circulation cells and weaker sensitivity to mesh orientation (AM98). Figure 6 shows no-slip solutions using aligned (left) and 45°-rotated (right) meshes, at 1/4° (top) and 1/16° (bottom) resolution. Under no-slip, there is no large recirculation cell in the northern part of the domain and the zonal transport is two times weaker there. Instead, no-slip solutions converge in time toward an oscillating small gyre nestled in the north-west corner. At 1/4° resolution, a 45° rotation of the mesh causes the small gyre to shift south by about 200 km (Figures 6a and 6b). At 1/16° resolution, this sensitivity vanishes (Figures 6c and 6d). Hence, insensitivity to mesh orientation is again achieved provided that spatial resolution is sufficiently fine. The same conclusion holds when no-slip is applied in the symmetric stress tensor formulation (not shown).  (a) The mesh is oriented at 45° and the spatial resolution 1/16° but coastal steps remain of size corresponding to 1/4° resolution. (b) The mesh is aligned and uniformly at 1/4° except within 300 km of the north coast where the meridional spatial resolution is refined from 1/4° to 1/16° (as illustrated in the right-end panel).

Symmetric Stress Tensor
A major caveat pointed out in AM98 is the extreme sensitivity to mesh orientation of free-slip solutions that resort to the symmetric viscous stress tensor. Figure 7 shows solutions using the free-slip boundary condition and the symmetric stress tensor. When the mesh is aligned with the physical coastline (Figures 7a and 7d), the model behavior is qualitatively similar to what was obtained using the rot-div stress tensor (Figures 3a and 3e). However, when the mesh is rotated by 45° (Figures 7b and 7e), both coarse and fine solutions change starkly and resemble no-slip solutions ( Figure 6) as pointed out in AM98. Hence, the free-slip boundary condition combined with the symmetric stress tensor appears to act as no-slip when the mesh is oriented at 45°, which suggests that its implementation is not suitable. We next examine how to recover a true free-slip condition on the 45°-turned mesh using the symmetric tensor.
When the mesh is not aligned with the physical coastline, the original straight boundaries become broken or indented in the model. To faithfully represent boundary flows, boundary conditions should be written with respect to the original, physical land-ocean frontier. For example, at 45°, the physical shoreline goes through T and F nodes (Figure 8). A free-slip condition is a mirror condition at the coast (e.g., Shchepetkin & O'Brien, 1996). Therefore, we should set virtual flows on "ghost nodes' (AM98), within land grid cells, that are symmetric to the ocean flows with respect to the physical shoreline. These virtual velocities are only used to evaluate the lateral friction term along the border. The example of a western coastline is illustrated in Figure 8. In this case, the mirror condition writes̃+ where the tilde is used to denote inland virtual velocities. It is the same condition as proposed by Griffiths (2013) in the context of Kelvin waves. Both anti-diagonal (D S or ζ) and diagonal (D T or χ) rates of deformation, defined at vorticity and tracer points respectively, are then deduced from this mirror condition. In particular, for a uniform grid spacing Δx (=Δy), we have Hence, D S is not zero at the tips of the coastal steps, but instead doubled compared to the value obtained with zero virtual inland velocities. This is contrary to the traditional implementation of free-slip in the symmetric tensor, which sets D S = 0 at the boundary. D S = 0 is the correct condition when the numerical and physical shorelines perfectly coincide but fails when they are misaligned. Note that zero vorticity condition is ensured at the boundary F nodes with the mirror condition (Equations 6 and 7).
The no-slip condition can be defined following the same rationale, by setting inland virtual velocities as the opposite to their oceanic mirrors: The proposed free-slip and no-slip conditions for a 45°-turned mesh and either stress tensor are summarized in Table 1. Their implementation consists in multiplying the factors defined in Table 1 with the rates of deformation estimated in the trivial case of zero (masked) velocities at the boundary and inland. Note that we could have considered that the physical shoreline goes through (U, V) nodes, as AM98, instead of (T, F) nodes. In this case, the mirror condition requires a slightly more complex interpolation of the virtual velocities.
We implemented the proposed free-slip condition in the symmetric stress tensor and assessed the impact on the equilibrium solution with a 45°-turned mesh (Figures 7c and 7f). In contrast to previous results which relied on  the standard free-slip boundary condition (Figures 7b and 7e), a true freeslip circulation is simulated with a northern recirculation cell that extends roughly to the middle of the basin. The solutions in Figures 7c and 7f are very similar to those previously obtained with the rot-div tensor (Figures 3b  and 3f). Our simulations thus confirm that the traditional way of applying free-slip in the symmetric tensor (D S = 0) is not suitable when the mesh is misaligned with the coastline.
We stress that setting only D S at the tips of steps is insufficient to obtain a true free-slip solution (not shown). Boundary conditions on both D S and D T are necessary to make the indented coastline slippery using the symmetric tensor. In contrast, canceling only ζ in the rot-div tensor proved to be enough to achieve slipperiness; doubling χ only brought minor changes. Since the condition ζ = 0 at the coast was already implemented in the experiments described in Section 3, a correct free-slip behavior was simulated.

Flux-Form Equations
All numerical experiments documented above were solved in vector-invariant form (Equation 2). How slippery are staircase-like coastlines in a model solved in flux form (Equation 4)? Figure 9 shows the steady flux-form solutions on aligned (Figures 9a and 9d) and 45°-turned (Figures 9b, 9c, 9e, and 9f) meshes at 1/4° (top row) and 1/16° (bottom row) resolution, using the rot-div stress tensor. Figures 9a and 9d reveals the same contraction of the northern recirculation cell with increasing resolution as found previously (Figures 3a and 3e) on the aligned mesh. However, when the mesh is turned and the coastline becomes indented, solutions become akin to no-slip: a small gyre is nestled in the northwest corner in an oscillatory steady state (Figures 9b and 9e).
By construction, there is no viscous boundary condition applied in flux-form advection, only the impermeability condition holds. The viscous boundary condition is free-slip and implemented in the rot-div stress tensor. Under such parameters, a coastal step tends to generate filaments of opposite vorticity that cause the coastal current to retroflect (Deremble et al., 2016; see also Section 5). Therefore, solutions solved in flux-form are expected to be sensitive to the presence of steps, as retroflection dynamics hinder along-boundary flow.
To remedy this sensitivity to mesh orientation, we implemented the same mirror condition (Equations 6 and 7) but in the advective trend, by enforcing velocities at the coast to satisfy this condition. The result is shown in Note. The table gives the modifications of quantities used in either tensor to be implemented at the coast, relative to the trivial case of zero virtual velocities. The top row describes quantities of the rot-div tensor while the bottom row describes quantities of the symmetric tensor, as defined in Section 2.1.

Figure 9.
Free-slip solutions solved in flux form using the rot-div stress tensor. Shading and isolines (in black) depict the active layer thickness h (isoline 500 m is thickened). x and y main ticks are 500 km apart. Spatial resolution is 1/4° (top row) and 1/16° (lower row). In the first column (a, d), the mesh is aligned so the coastline is straight. Second and third columns (b, c, e, f) have the mesh rotated at 45°. In the third column, a mirror condition is enforced in the advective trend. Figures 9c and 9f for the 45°-turned mesh. We obtain a sizable inertial recirculation cell in the northern part of the basin, with zonal extensions close to the previous solutions (Figures 3b, 3d, 3f, and 3h and Figures 7c and 7f). Hence, enforcing the free-slip condition in both the diffusive and advective terms makes staircase-like coastlines slippery, including with flux-form equations.

Intermediate Angles
At 1/4° resolution, the fully indented coastline (Figure 3b) leads to a solution closer to the converged solution (Figures 3e-3h) than does the straight coastline (Figure 3a). At intermediate orientations of the mesh (strictly between 0° and 45°), the coastline counts fewer outgoing angles than at 45°, and a solution midway between the aligned and 45°-turned cases might be expected. In reality, intermediate angles generate solutions that depart much more from the converged solution ( Figure 10). Figure 10 uses the exact same numerical choices as Figure 3 except for the mesh orientation, which is either 10° (left column) or 30° (right column). At these orientations, the inertial recirculation cell expands toward the east between 1/4° and 1/16°, then contracts with resolution ( Figure  Why do solutions at intermediate angles differ from their aligned and 45°-turned counterparts? Once the model is converged, dynamics inside the domain must be well captured independently of the orientation of the grid, so that differences are expected to lie at the borders. If the mesh is aligned (Figure 11, left) or turned at 45° (Figure 11, right), the physical coastline coincides with the F nodes of the C-grid and the free-slip condition ζ = 0 exactly matches that of a straight coastline. At intermediate angles however, this condition is inaccurate because the freeslip boundary resulting in the model (red dotted line) is not straight but corrugated (Figure 11, middle). Imposing ζ = 0 on the tips does make artificial steps slippery yet does not achieve the true free-slip condition of a straight coastline. We expect that defining the free-slip condition with respect to a straight boundary would yield the same solutions and convergence rate as with the aligned mesh.
In other words, to accurately represent the free-slip condition, virtual inland velocities should be interpolated as the mirrors of the ocean flows with respect to the straight shoreline. Such a strategy would presumably allow to retrieve a true free-slip behavior also with advection in flux form, and with the symmetric stress tensor, for any orientation of the mesh. We have not endeavored such an interpolation at intermediate angles, for it is tedious and would not generalize to arbitrary (curved) physical shorelines. Instead, we suggest using more general techniques such as immersed boundary methods (Causon et al., 2000;Ketefian & Jacobson, 2011;Kirkpatrick et al., 2003).

Benefit of Wide Stencils
Strictly speaking, the mirror condition should apply in each term of the equations. For example, on the 45° oriented mesh within oceanic cells along the border, this would double the kinetic energy and ∂ t h in Equations 1 and 2 respectively, while the vertical scale factor h and the Coriolis term at the outgoing vorticity points would be mirrored with respect to the cell diagonals. However, our implementation of these conditions did not bring noticeable changes on the 45°-turned solution (not shown), which is already very close to the reference. This result indicates that the sensitivity of circulation within this configuration is controlled primarily by the formulation of advective and diffusive terms, in accord with Dupont et al. (2003).
Motivated by the sensitivity of solutions to the discrete formulation of advection, we investigated the impacts of mesh orientation and coastline indentation using a different advection scheme (EEN). Results are shown in Figure 12 under the free-slip boundary condition. The mesh is progressively turned from left to right (0°, 10°, 30°, and 45°) and spatial resolution increases from 1/4° to 1/8°. First, we find that the aligned and intermediate solutions are physically converged at 1/8° resolution (1/16° not shown), contrary to solutions that resorted to ENS or ENE. Second, all solutions are very similar across mesh orientations, even at 1/4° resolution, with a recirculation cell that extends halfway through the basin. These results are quantified in Figure 4 which shows that EEN-based solutions (in solid lines) are virtually identical as early as 1/8°, and that they converge together to the same state as ENS-based solutions as resolution increases to 1/48°. Vorticity schemes previously used in this study have a 7-point wide stencil, whereas the EEN scheme has a 17-point wide stencil (see Figure B1). We infer that the larger stencil of the EEN scheme effectively smoothens the discontinuity of the coast, making dynamics much less sensitive to the misalignment of the grid with the physical shorelines. These results advocate for the use of schemes having relatively wide stencils, possibly high-order schemes, to minimize spurious effects of staircase-like coastlines.

Dynamics Along a Cornered Coastline
In the previous sections, we explored ways to eliminate spurious effects of artificial steps in model coastlines. However, coastal steps are not always artificial and their effects not necessarily spurious, since real coastlines contain sharp turns that exceed the grid scale and impact boundary currents. For example, a protruding corner in the coastline can cause boundary currents to retroflect (Dengg, 1993), impacting the larger scale circulation (Ansorge & Lutjeharms, 2005;Weeks et al., 2010). Which boundary conditions are most appropriate to model this physical response to a coastal step?
To address this question, we reproduce the configuration of Deremble et al. (2016). The domain is a square basin of 500 km in length, cropped by a 100 km × 250 km land mass at the southwest end, as illustrated by the gray shading in Figure 13. The equations solved are those given in Section 2.1, except that the model is barotropic and excludes wind forcing, bottom friction and the Coriolis effect (parameters are given in Appendix A). The simulation starts with a vortex of negative relative vorticity, standing near the eastern side of the land mass ( Figure 13a). The vortex is then advected up to the corner due to non-linear interaction with the straight coastline. When the Figure 11. Representation on a C-grid of a straight shoreline with the piecewise constant approximation. The black disks are the vorticity (F) nodes and the crosses are the height (T) points. The numerical frontier (red dotted line) connects the vorticity nodes that influence the dynamics. In the aligned case (left) the numerical frontier occupies cell faces and joins F nodes so it coincides with a straight coastline. In a similar way, the numerical frontier in the 45°-turned case (right) goes through cell diagonals, hence intersecting F and T nodes. In the intermediate case (middle), the closest vorticity points to the targeted straight coastline (in blue) are not aligned, resulting in a slithering numerical frontier.

Figure 12.
Free-slip solutions solved in vector-invariant form with the EEN advection scheme and the rot-div stress tensor. Shading and isolines (in black) depict the active layer thickness h (isoline 500 m is thickened). x and y main ticks are 500 km apart. From left to right, the mesh is progressively rotated, at 0°, 10°, 30°, and 45° with respect to the physical coastline. Spatial resolution increases from 1/4° (top) to 1/8° (bottom).
vortex begins to overtake the corner, filaments of positive vorticity are generated at the tip, forcing the retroflection of the flow as explained by Deremble et al. (2016). Figure 13b shows the vortex detaching from the coast along a filament of opposite relative vorticity that stretches northeastward from the tip.

Lateral Boundary Conditions
On a large-scale isolated step, the definitions of the boundary conditions given in Section 2.3 are less straightforward as it is no longer possible to define the local vectors (n, t) at the singularity. A no-slip condition (u = 0) could potentially hold by continuous extension from the adjacent walls to the tip. However, the definition of the free-slip viscous boundary condition becomes unclear as the normal derivative of the flow is unknown at the tip. Furthermore, by bending a regular coastline toward the limiting case of a protruding corner (κ tends to infinity), the vorticity in Equation 5 becomes infinite, suggesting that a singularity in the coastline acts as a source of vorticity in the flow (Deremble et al., 2016). Dengg (1993) argued that the retroflection of a boundary current at a cape is better captured in numerical simulations when a no-slip viscous condition is used instead of free-slip. However, Deremble et al. (2016) described analytically the detachment or retroflection of a vortex under the assumption of no dissipation, with only advection as the driver. They showed that such retroflection requires only impermeability conditions along both side of the walls in the advective trend.
Here we extend this result to several numerical formulations of advection, showing that retroflection is faithfully represented when viscosity ν is set to zero at the F-node of the singularity (Figures 13b, 13e, and 13h). With flux-form equations, there is no required viscous boundary condition in advection, so that only impermeability is imposed and retroflection is simulated (Figure 13b). With vector-invariant equations, retroflection is also captured under either ENS (Figure 13e), ENE (not shown) or EEN (Figure 13h) advection schemes when using the impermeability condition to estimate vorticity at the tip. However, if relative vorticity is enforced to vanish at the tip (ζ c = 0), the retroflection is no longer simulated (Figure 13c). The enforcement of zero vorticity is equivalent to applying a free-slip viscous condition; this result thus matches previous reports that free-slip can suppress retroflection (Dengg, 1993;Deremble et al., 2016;Dupont & Straub, 2004). If instead a no-slip condition is applied at the tip, vorticity filaments change intensity but are still represented (compare Figure 13d with Figures 13b, 13e, and 13h).
When neutralizing the retroflection effect in the advection term (by imposing ζ c = 0 at the tip in the advective trends), it is still possible to recover the filament generation through the viscous dissipation. Indeed, with a non-zero viscosity ν set at the singularity, retroflection is simulated by applying the no-slip viscous boundary condition in the diffusive term (Figure 13f). However the filament of positive relative vorticity is different in shape and intensity (compare Figure 13f with Figures 13b, 13e, and 13h), its characteristics become more sensitive to the chosen value of viscosity ν (not shown), and its physical interpretation as the detachment of the viscous boundary layer (Deremble et al., 2016) is less straightforward.
The rot-div stress tensor is used for diffusion in the solutions discussed above (Figures 13b-13f and 13h). In Figure 13g, we show a solution that uses the symmetric stress tensor combined with the viscosity ν set to zero at the outgoing corner. The solution is quite similar to that obtained using the no-slip viscous condition in the rot-div stress tensor (Figure 13f). This result concurs with those of Section 4 and emphasizes the problematic behavior of the traditional implementation of slipperiness in the symmetric tensor. These findings also explain why Dupont and Straub (2004) and Deremble et al. (2016) obtained unexpected behaviors at coastal tips when using the combination of symmetric tensor and free-slip viscous boundary condition: in effect, this combination produces no-slip on steps.
To summarize, our experiments demonstrate that the detachment of a vortex can be represented in the absence of dissipation by estimating the vorticity at the tip using either impermeable or no-slip conditions. By applying a zero vorticity condition at a coast's singularity, the generation of filaments of opposite vorticity and the retroflection of a current are instead prevented, as explained by Equation 5 and illustrated in Figure 13c. We disentangled the effect of the viscous boundary condition at the singularity under various formulations for advection (vector-invariant and flux-form equations) and diffusion (rot-div and symmetric stress tensors). The results highlight the inadequacy of the usual formulation of free-slip (zero vorticity) to simulate the lateral interaction of a boundary current with a cape, and potentially analogous flow-topography interactions along curved coastlines, as suggested by previous work (Dupont & Straub, 2004).

Implications for Staircase-Like Coastlines
Staircase-like coastlines studied in Sections 3 and 4 using the AM98 configuration can be viewed as a series of small isolated steps. Each step may be expected to have a dynamical effect on the local circulation as described in Section 5.2, and the ensemble of steps may have a cumulative impact on the basin-scale gyres. For example, the behavior of solutions using only the impermeability condition, as it is the case in flux-form equations on the 45°-rotated mesh (Figures 9b and 9e), can be understood by noting that each step works to retroflect the boundary flow and ultimately distort the gyre circulation. Indeed, in these simulations the viscous boundary condition is free-slip and is applied only in the diffusive term. Therefore, the lack of a zonally extended inertial recirculation (Figures 9b and 9e) likely stems from the retroflection-induced by the advection term-of the boundary current on the indented coastline. Interestingly, increasing spatial resolution (while keeping viscosity ν and friction r unchanged) allows the inertial recirculation gyre to grow eastward (compare Figure 9e with Figure 9b). An additional simulation (not shown) at a finer resolution of 1/32° confirms this tendency: the inertial recirculation cell then occupies over half of the basin. We interpret this behavior as the consequence of the reduction of the scale of steps (1/32°) relative to the width of the boundary currents (∼1/4°): the production of vorticity filaments, which feeds upon the discontinuity in the flow field across the tips of steps (Deremble et al., 2016), is damped if steps are too small to generate sizable discontinuities.

Conclusions
We revisited the "staircase problem" highlighted by Adcroft and Marshall (1998), who exposed the existence of a spurious form drag when smooth coastlines are numerically transformed into steps. We reproduced their configuration, which consists of a square closed basin under shallow water dynamics and cyclonic wind forcing, simulated on a Cartesian mesh with varying orientation. We tested various mesh resolutions with many combinations of advection formulations (flux or vector-invariant forms of equations with potential-enstrophy (ENS), energy (ENE) or energy and potential-enstrophy (EEN) conserving schemes developed by Sadourny) and two commonly used viscous stress tensors (rot-div and symmetric formulations).
We first show that the free-slip non-rotated solution is not physically converged at 1/4° resolution, under ENS or ENE with the rot-div stress tensor, but only from 1/16° with the same viscosity and friction parameters as AM98.
By physical convergence we mean the insensitivity of the main characteristics of the flow to further increase of the mesh resolution. Such convergence requires to resolve the inertial dynamics induced by the free-slip (i.e., mirror) boundary condition along the northern coast by having at least four grid points per internal radius of deformation, which is close to 30 km in this region.
In addition, we find that the 45°-rotated free-slip solution is also physically converged from 1/16° resolution and surprisingly looks almost identical to the aligned 1/16° solution, contrary to AM98. The reason is that the free-slip boundary condition (zero vorticity) applied at the tips of coastal indents created by the 45°-rotated mesh exactly stands for a straight shoreline passing through T and F nodes of a C-grid. At intermediate angles of mesh orientation (strictly between 0° and 45°), physically converged solutions are only achieved near 1/32° resolution and depart from the non-or 45°-rotated 1/16° solution. We suggest that the different sensitivity and delayed convergence at intermediate angles stem from inaccurate declaration of the free-slip boundary condition as the resulting numerical frontier is corrugated, not straight.
The above-mentioned results hold in vector-invariant form with the ENS or ENE scheme combined with the rot-div stress tensor. When switching to the symmetric stress tensor or to flux-form advection, rotated solutions no longer converge toward the reference (aligned) solution and instead resemble no-slip solutions. In both cases, applying a true mirror boundary condition with respect to the physical coastline on the 45°-rotated mesh allows to retrieve solutions close to the reference at 1/16° resolution. These results pinpoint the spurious behavior of the traditional implementation of "free-slip" in the symmetric viscous tensor, which works as intended only if the mesh is aligned with the border.
Importantly, using an advection scheme with a larger stencil (EEN in vector-invariant form) makes the freeslip solutions virtually identical from 1/8° for any orientation of the mesh. We infer that larger stencils allow advection schemes and the simulated circulation to become much less sensitive to broken coastlines, providing a practical avenue to mitigate spurious effects of piecewise-constant land-ocean frontiers.
In addition to exposing ways to eliminate spurious effects of artificial coastal steps, we explored the numerical treatment of a single sharp turn in the physical coastline, using the configuration of Deremble et al. (2016). One expected impact of a large protruding corner in the coastline is the generation of vorticity filaments that force the boundary current to retroflect (Deremble et al., 2016). We show that the retroflection is faithfully captured by advection if expressed in flux form, which disallows the use of a viscous condition, or in vector-invariant form, provided an impermeability or no-slip condition is applied at the maxima of curvature. In contrast, using a free-slip (zero vorticity) viscous boundary condition at the tip suppresses the retroflection entirely.
We conclude that staircase-like coastlines can behave like straight and slippery coastlines provided that coastal dynamics are physically resolved and that an accurate mirror boundary condition is used. Naturally, in realistic OGCM configurations, it is not always feasible nor desirable to increase resolution whilst keeping viscosity unchanged to achieve physical convergence. To minimize spurious effects of artificial steps in OGCM boundaries, we thus recommend the use of advection schemes with large stencils (such as the energy-enstrophy conserving (EEN) scheme), combined with free-slip (zero vorticity) boundary conditions and the rot-div viscous stress tensor. If flux-form equations or the symmetric stress tensor are chosen, freeslip along staircase-like coastlines is best implemented with general techniques such as immersed boundary methods. However, when coastal steps represent sharp turns in the real coastline that cause boundary currents to retroflect, perfect slipperiness is no longer desirable as it may suppress retroflection. In this case, application of impermeability or no-slip conditions can better represent physical flow-topography interactions than a free-slip (zero vorticity) boundary condition. Consequently, the best choice of boundary condition may depend on the degree to which the considered boundary steps reflect numerical artifacts versus real topographic features. In realistic configurations, it is thus expected that ideal boundary conditions should vary with location.

Appendix A: Numerical Parameters
In the AM98 configuration, the initial state is at rest with the thickness of the active layer uniformly equal to h 0 = 500 m. Density is ρ 0 = 1,000 kg m −3 and reduced gravity g′ = 0.02 m s −2 . The Coriolis parameter evolves on a beta-plane f(y) = f 0 + βy where f 0 = 0.5 × 10 −4 s −1 and β = 2 × 10 −11 m s −1 so that the internal radius of deformation = √ ′ ℎ0∕ is about 45 km at the mid-basin. Zonal wind stress is τ = −τ 0 cos(πy/L) with τ 0 = 0.2 N m −2 . Uniform friction parameters are considered with r = 10 −7 s −1 the bottom linear friction and ν = 500 m 2 s −1 the lateral viscosity. Uniform mesh resolution (Δx = Δy) and associated time-step for the AM98 configuration are summarized in Table A1. The Asselin filter parameter of the Leap Frog Robert-Asselin time-stepping scheme is ϵ = 10 −1 . In vector-form, gradient of kinetic energy is discretized with a second order centered scheme. In flux form, advection is also discretized with a second order centered scheme while the Coriolis term is discretized with the ENS scheme.
In the single vortex configuration of Deremble et al. (2016), the initial conditions are a flat sea surface and a vortex of negative vorticity placed next to the wall. In the relative frame of reference centered on the vortex, the initial horizontal speeds are given by the azimuthal profile v θ (Lamb-Oseen vortex): with the pseudoradius r 0 = 20 km and the strength of the vortex Γ = −5 × 10 4 m 2 s −1 . The model is barotropic with gravity g = 9.81 m s −2 and the basin is 1 m deep. There is no wind forcing nor bottom friction. Coriolis effects are not considered (f = 0 throughout the basin). Small lateral dissipation is added to ensure numerical stability: viscosity ν = 20 m 2 s −1 . Simulations are ran over 45 days with an uniform spatial resolution of Δx = Δy = 1.25 km and a time-step of 90 s.

Appendix B: Formulation of the Discrete Vorticity Schemes Used in Vector-Invariant Form Simulations
The discrete vorticity schemes considered in this study are defined on the Arakawa C-grid. The grid is staggered so that total vorticity ζ i,j + f i,j , horizontal velocities (u i,j ; v i,j ) and height h i,j variables are arranged as shown in Figure B1a. The total potential vorticity = ( + ) ℎ is needed for vector-invariant form advection so layer thickness at F-nodes h f is deduced from the height h nodes: = ; while vorticity is diagnosed as follows: where ( − , − ) and (δ i , δ j ) are the averaging and differencing operators at the mid point, for example, ℎ = 1 2 (ℎ + ℎ +1 ) and δ i v = v i+1,j − v i,j . The horizontal scale factors e 1t , e 1u , e 1v and e 1f (e 2t , e 2u , e 2v , and e 2f ) are derived analytically at each node from the latitudinal (longitudinal) coordinate; on a uniform regular mesh, they are all equal to Δx = Δy. Figure B1b-B1d represents the discretization at u i,j nodes (blue disks) and illustrates the size of the stencil for each scheme. First, the potential enstrophy conserving scheme ( Figure B1b) (Sadourny, 1975) provides a global conservation of global enstrophy h f q 2 for non-divergent flow. For x and y components of the vorticity term, it writes as: Figure B1. Discrete vorticity schemes on the Arakawa C-grid. (a) Location and indexing of height h nodes (black crosses), horizontal velocity u and v nodes (black circles) and vorticity ζ nodes (black disks) on a single cell. Panels (b-d) illustrate respectively the ENS, ENE, and EEN discretization of vorticity in vector-invariant form. Single arrows (in gray) add themselves; double arrows represent two-point averaging; right-angle cornered arrows symbolize triads j Q ; and rectangles (gray edged) is the factoring of the averaged quantities with the local variables. where (U; V) are the transports across cell faces, for example, V = e 1v h v v. Then, the kinetic energy conserving (ENE) vorticity scheme ( Figure B1c) (Sadourny, 1975) is defined as: Finally, the EEN scheme developed by Sadourny (Burridge & Haseler, 1977) is a member of the family of vorticity schemes derived by Arakawa and Lamb (1981) that conserves kinetic energy and, provided there is no divergence in the flow, potential enstrophy. This scheme relies upon averaging triads of vorticity Q that ultimately widen the stencil up to 17 velocity nodes; instead of 7 in the ENS or ENE schemes as they use a much cheaper two-point averaging. A triad Q is defined as: with (l, m) ∈ I 2 where I = (1/2; −1/2). Each triad combines with the adjacent transport , for example, +1 multiplies with +1 Q −1∕2 1∕2 (in red in Figure B1d). Expressions for EEN scheme summarize as

Data Availability Statement
NEMO code is available at https://forge.nemo-ocean.eu/nemo/nemo. The described version is 4.2. The configurations are available at https://doi.org/10.5281/zenodo.7480139 and the plotting scripts are available at https:// doi.org/10.5281/zenodo.7480159. We also thank David Marshall and two anonymous reviewers for comments that helped improving the manuscript. This work was granted access to the HPC resources of IDRIS under the allocation 2021-A0100107451 made by GENCI. This study was carried out as part of the project 18CP03 (Analyse numérique pour la modélisation de la dynamique océanique à fines échelles) under the auspices of French Ministry of Defense/ DGA, by Shom and l'INRIA.