<figure><img src="image_1.png" /></figure>

HAL
open science

# Improved Gulf Stream separation through Brinkman penalization

Laurent Debreu, Nicholas Kevlahan, Patrick Marchesiello

<figure><img src="image_2.png" /></figure>

## ► To cite this version:

Laurent Debreu, Nicholas Kevlahan, Patrick Marchesiello. Improved Gulf Stream separation through Brinkman penalization. Ocean Modelling, 2022, 179, pp.102121. 10.1016/j.ocemod.2022.102121. hal-03593222

HAL Id: hal-03593222

https://inria.hal.science/hal-03593222

Submitted on 1 Mar 2022

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Gulf Stream separation through Brinkman penalization

L. Debreu<sup>a,*</sup>, N.K.-R. Kevlahan<sup>b</sup>, P. Marchesiello<sup>c</sup>

<sup>a</sup>Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, LJK, 38000 Grenoble, France

<sup>b</sup>Department of Mathematics and Statistics, McMaster University, Hamilton, Canada

<sup>c</sup>IRD/LEGOS, Toulouse, France

2

3

4

5

## 6 Abstract

The advantage of a smooth representation of bathymetry in terrain-following σ-coordinate ocean models is compromised by the need to avoid numerical errors on steep slopes associated with pressure gradient discretization or spurious diapycnal diffusion. Geopotential z-coordinate models avoid these errors, but greatly underrepresent the interaction of flow with a topographic slope, especially when the bathymetry is underresolved. Hybrid coordinate models are also deficient because it is difficult to find a satisfactory compromise between z and σ coordinates. With volume penalization, we do not seek a compromise, but rather a correction to the usual coordinate systems that realistically recovers continuous and steep bathymetry. The Brinkman volume penalization method studied here is a modified version of the one introduced in Debreu et al. (2020) that simplifies the numerical implementation of the penalization, increases robustness and improves its computational performance for realistic long-term simulations, while preserving accuracy.

We apply this penalization method to the Gulf Stream separation problem that has puzzled modelers for decades. The method improves the representation of the flow-topography interaction and achieves realistic separation of the Gulf Stream at resolutions as coarse as 1/8°. In addition, it provides a tool to separate the effect of eddy activity and topographic slope when changing grid resolution. This has never before been possible because at coarse resolution none of the usual coordinate systems can properly represent a steep continental slope. Our results show that realistic bathymetry is more important than eddy activity in ensuring realistic Gulf Stream separation, even though many recent studies tend to focus on the eddy activity. A steep slope can exert a stabilizing influence that promotes a strong mean slope current with strong inertia that helps it separate from the coast at the topographic curvature of Cape Hatteras.

We anticipate that a successful topographic slope correction will be very valuable to climate models, as their current resolution is far from sufficient to represent western boundary currents (WBCs) using traditional coordinate systems. Our results suggest that a climate model with a 1/4° resolution using volume penalization — and perhaps also some
parameterization of the eddy-mean flow interaction to energize the WBCs — would represent ocean circulation much more realistically.

Keywords: algorithms, penalization, bathymetry, Gulf Stream, ocean modeling

# 1. Introduction

## 1.1. The conundrum of Gulf Stream separation

The Gulf Stream (GS), an example of a Western Boundary Current (WBC), is one of the strongest ocean currents on Earth and a major feature of the global ocean circulation that largely controls the Earth's climate (Lee et al., 2018). There are iconic features of WBCs, for which the mechanisms and modeling capability are still highly debated. One of these is the GS separation from the coast, which was associated early on with the curvature of the coast and the inertia required to overcome the topographic steering (Dengo, 1993; Couvelard et al., 2008; Ezer, 2016; Schoonover et al., 2016). More recently, WBCs have been linked in ocean models to the resolution of mesoscale activity through eddy-mean flow interaction (McWilliams, 2008; Bryan et al., 2007; Chassignet and Marshall, 2008). A mechanistic description of eddies interacting with bathymetry has also emerged (Gula et al., 2015), showing that large isobath curvatures tend to maximize eddy energy transfer, and thus flow-bathymetry interaction. The influence of eddies is somewhat mitigated by dissipation mechanisms, including those provided by surface drag on the atmosphere in air-sea coupled models (Renault et al., 2019). However, flow-bathymetry interaction near isobath curvature remains the primary mechanism.

Failing to properly represent the effect of bottom bathymetry results, in many numerical models, in the GS over-shooting Cape Hatteras, i.e., separating from the coast further north than observed (Chassignet and Marshall, 2008; Schoonover et al., 2016; Ezer, 2016). This is because the bottom pressure torque (not the wind-stress curl) is the main component balancing the planetary vorticity advection (beta term) in the local vorticity budget (e.g., Myers et al. 1996; Couvelard et al. 2008; Schoonover et al. 2016). Therefore, GS separation is related to local, rather than basin-scale, wind-driven dynamics. The implication is that misrepresenting the local interactions between flow and bathymetry

<footnote>*Corresponding author.</footnote>

<footnote>Email addresses: Laurent.Debreu@inria.fr (L. Debreu), kevlahan@mcmaster.ca (N.K.-R. Kevlahan), Patrick.Marchesiello@ird.fr (P. Marchesiello)</footnote>
can have an impact on the entire basin and beyond. It is now accepted that GS separation is more realistic in models using terrain-following σ- or s-coordinates rather than the geopotential z-coordinates (see the model intercomparison study by Schoonover et al. 2016). However, topographic slopes are limited in these models by numerical constraints and a resolution finer than 1/10° is required to avoid excessively smooth bathymetry that results in an unrealistically low bottom pressure torque — high resolution also allows for realistic generation of mesoscale eddies. However, resolutions as fine as 1/50° are recommended for a realistic GS separation (Chassignet and Xu, 2017; Hurlburt and Hogan, 2000). We will demonstrate that penalization permits a similar result at 1/8° resolution.

The advantage of a smooth representation of bathymetry in terrain-following models is impaired by the need to avoid numerical errors over steep bathymetry associated with the pressure gradient discretization (Shchepetkin, 2003) and spurious diapycnal diffusion (Marchesiello et al., 2009; Lemarié et al., 2012). For this reason, some attempts have been made to use hybrid coordinates that transition from terrain following coordinates in shallow regions to z-coordinates at depth (e.g., Ezer and Mellor 2004; Chassignet and Xu 2017). However, these models show the usual deficiencies in GS simulations. Arguably, it is difficult to find a satisfactory compromise that would get the best of each coordinate system over their preferred regions, i.e., the continental shelf, slope, or deep ocean. Here, we are not looking for a compromise, but rather a correction to the classical coordinate systems that realistically recovers a continuous, steep bathymetry (i.e., without the step-like features that produce unrealistic vorticity perturbations). A more accurate and robust vertical coordinate system would be very valuable for climate models, as their resolution is far from sufficient to represent WBCs using z-coordinates. The volume penalization method that performs such a correction is presented below.

## 1.2. Penalization of bathymetry in ocean models

Volume penalization is a simple and elegant way to implicitly impose Neumann (no penetration) and Dirichlet (no slip) boundary conditions in numerical methods for partial differential equations (PDEs) (Arquis and Caltagirone, 1984). In a penalization approach the PDE is solved on a computational domain with simple geometry (e.g., a rectangular domain with a regular grid) and the boundary conditions are imposed approximately by penalizing the solution outside the actual solution domain. Volume penalization has been particularly effective over the past two decades for approximating no-slip boundary conditions in fluid-structure interaction problems (e.g. Angot, 1999;
54

Angot et al., 1999; Kevlahan and Ghidaglia, 2001; Kevlahan and Vasilyev, 2005), but its application to realistic three-
dimensional oceanic problems is new.

55

56

In Debreu et al. (2020) we introduced a version of volume penalization based on the notion of flow through a porous medium to approximate the complex multi-scale geometry of ocean bathymetry in realistic ocean models. This method is an extension of a previous work by Kevlahan et al. (2015) to approximate the no-slip lateral boundary conditions in a global shallow water model. In this version, the accuracy of the penalization is controlled by two parameters: the porosity α and the permeability ε. The solid regions are defined by a (smoothed) mask. To approach no-slip boundary conditions, both the porosity and permeability are set very small in the solid regions. To approximate non-penetrating conditions, only porosity is used and the friction term associated with permeability is not included.

Kevlahan et al. (2015) showed that the error of the penalization approximation of no-slip boundary conditions is $O(αε^{1/2})$.

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

Volume penalization has several advantages for representing the bottom bathymetry in ocean models. First, because the penalized region does not suffer from the topographic slope limitations imposed by pressure gradient errors, it is possible to use more realistic bathymetry at a given computational grid resolution than with the usual σ-coordinates terrain-following grid. Secondly, because the solid regions are represented implicitly by a mask, it is not necessary to describe the geometry of the solid-fluid interface in detail. Finally, the permeability and porosity parameters can be tuned to model the subgrid-scale structure of bathymetry (Adcroft, 2013), or the composition of the seafloor substrate (e.g., deltas, marshes, flooded urban areas Guinot et al. 2018). Penalization methods are also ideal for dynamically adaptive codes (Hejazialhosseini et al., 2010), which was the original motivation of Kevlahan and Vasilyev (2005) and (Kevlahan et al., 2015).

74

In §3, we describe the penalization method in detail, and propose several improvements to our previous approach.

75

The penalization is then implemented in the Coastal and Regional Ocean COmmunity model (CROCO; Shchepetkin and McWilliams 2005; Debreu et al. 2012) and its performance is validated in realistic simulations of the GS.

76
## 2. Standard model configuration

CROCO is a free-surface, terrain-following coordinate regional ocean model with barotropic-baroclinic mode split-explicit time stepping (Shchepetkin and McWilliams, 2005; Debreu et al., 2012; Soufflet et al., 2016). Here, we use the Boussinesq and hydrostatic approximations, and although a non-Boussinesq and non-hydrostatic solver is also available (Hilt et al., 2020; Marchesiello et al., 2021), it is not required for the mesoscale resolutions used in this paper.

### 2.1. *Model configuration*

The North Atlantic regional domain is a subset of that presented in Renault et al. (2019). The grid encompasses the Florida and GS currents in the western North Atlantic, from 82° W to 45° W and 22° N to 47° N, with three horizontal spatial resolutions: 1/4°, 1/8° and 1/12°. We use 32 generalized σ-coordinate vertical levels (Lemarié et al., 2012) with stretching parameters $h_{cline} = 200$ m, $\theta_b = 2$, and $\theta_s = 7$. The bathymetry is constructed from the Shuttle Radar Topography Mission (SRTM30 plus) dataset based on the 1-min Sandwell and Smith (1997) global dataset and higher-resolution data where available. To avoid aliasing at higher resolution than the computational grid and to ensure the smoothness of topography at the grid scale, a Gaussian smoothing kernel with a width four times the topographic grid spacing is used. In the non-penalized σ-coordinate (or s-coordinate) configuration, pressure gradient errors caused by steep slopes are kept at acceptable levels by applying a local smoothing of the bottom topography so that the r-factor (Beckmann and Haidvogel, 1993) does not exceed 0.25. In addition to the three spatial resolutions, we present a 1/12° simulation with smoother bathymetry interpolated from the 1/4° configuration. This test will allow us to evaluate the effect of resolution that is not associated with a change in topography.

In each case, the model is run for 15 years, using a monthly climatology of atmospheric surface fluxes from COADS (freshwater and heat fluxes) and from QuikSCAT (momentum fluxes). Surface current feedback (CFB) is applied using a stress correction approach for momentum fluxes (Renault et al., 2020). Vertical mixing of tracers and momentum is done with a K-profile parameterization (KPP; Large et al. 1994). At oceanic open boundaries, boundary data of temperature, salinity, surface elevation, and horizontal velocities are taken from the year 2000 monthly-mean Simple Ocean Data Assimilation (SODA) product (Carton and Giese, 2008). The open boundary condition algorithms
consist of an active-passive 2D radiation scheme for the tracers and baroclinic mode and a modified Flather-type scheme for the barotropic mode (Marchesiello et al., 2001). The bottom drag $\tau_b$ is computed assuming that the flow in the bottom boundary layer has a vertical logarithmic profile, defined by a bottom roughness length $z_{0b}$ (m), which takes a constant value of $10^{-2}$ m in our simulations:

$$ (\tau_{bx}, \tau_{by}) = \left( \frac{\kappa}{\log \frac{\zeta}{z_{0b}}} \right)^2 \sqrt{u_b^2 + v_b^2} (u_b, v_b), \quad (1) $$

where $u_b$ and $v_b$ are the bottom currents in the log layer at a height $z$ above the ocean floor (here corresponding to the first vertical grid level) and $\kappa = 0.41$ the von Karman constant. The formulation of the bottom stress given here will be revisited in the penalization section.

## 2.2. Vorticity budget

WBC dynamics can be diagnosed by the barotropic vorticity budget, which is an extension of Sverdrup balance. Following, e.g., Couvelard et al. (2008); Schoonover et al. (2016), the barotropic vorticity budget is computed by taking the rotational of the vertically integrated lateral momentum equations

$$ \frac{\partial \bar{\zeta}}{\partial t} = \frac{J(P_b, h)}{\rho_0} - A - \nabla \cdot (f \mathbf{U}) + \frac{\nabla \times \tau}{\rho_0} - \frac{\nabla \times \tau_b}{\rho_0} + D, \quad (2) $$

where $\partial \bar{\zeta} = (\nabla \times \mathbf{U}) \cdot \hat{\mathbf{z}}$ is the barotropic vorticity, $\mathbf{U}$ the barotropic current, $J$ is the Jacobian operator, $P_b$ is the bottom pressure, $h$ the bathymetry, $A$ the advection torque, $f$ is the Coriolis force, $\tau_b$ is the bottom stress, and $D$ is the viscous torque. In the following, the right-hand-side terms of Eq. 2 are referred to as the bottom pressure torque ($\frac{J(P_b, h)}{\rho_0}$), the advection torque ($-A$), the planetary vorticity advection ($-\nabla \cdot (f \mathbf{U})$), the wind stress curl ($\frac{\nabla \times \tau}{\rho_0}$), the bottom stress curl ($-\frac{\nabla \times \tau_b}{\rho_0}$), and the viscous torque ($D$). The bottom pressure torque is a measure of the topographic steering of the flow (Couvelard et al., 2008). The nonlinear torque represents the advection of vorticity by the mean and eddy flow. The planetary vorticity advection is due to the combined effects of the earth's curvature and rotation ($\beta$ term in the Sverdrup balance). The surface wind stress curl is a top drag curl and can be both a source or dissipation of vorticity. The bottom stress curl is the drag effect of unresolved topography. Finally, the viscous torque represents the vorticity dissipation due to turbulent viscosity. In a regional simulation it has large values only within the numerical sponge layers that implement the open boundaries.
### 3. Description of the penalization method

#### 3.1. Penalization in the primitive equations of the ocean

The volume penalization method investigated here is a modified version of the one introduced in Debreu et al. (2020). The improvements simplify its numerical implementation, but do not affect its accuracy or theoretical basis.

The prognostic equations are written in conservation form using the generalized vertical coordinate system $s(x, z, t)$, with $h = \partial z/\partial s$, as

$$ \frac{\partial \tilde{h}u}{\partial t} + \frac{\partial \tilde{h}u^2}{\partial x} + \frac{\partial \tilde{h}u \Omega}{\partial s} = -\tilde{h} \left( g \frac{\partial \eta}{\partial x} + \frac{1}{\rho_0} \frac{\partial p_h}{\partial x} + \frac{\rho}{\rho_0} g \frac{\partial z}{\partial x} \right) - \frac{1}{\epsilon} \mathbf{1}(x, z) \tilde{h}u, \quad (3) $$

$$ \frac{\partial \tilde{h}}{\partial t} + \frac{\partial \tilde{h}u}{\partial x} + \frac{\partial \tilde{h}\Omega}{\partial s} = 0, \quad (4) $$

$$ \frac{\partial p_h}{\partial s} = -h\rho g, \quad (5) $$

$$ \frac{\partial \tilde{h}\rho}{\partial t} + \frac{\partial \tilde{h}u\rho}{\partial x} + \frac{\partial \tilde{h}\Omega\rho}{\partial s} = 0, \quad (6) $$

where the penalized height $\tilde{h} = \phi h$, and $\phi$ is the porosity

$$ \phi(x, z) = \begin{cases} \alpha & \text{in the penalized solid regions,} \\ 1 & \text{in the fluid regions.} \end{cases} \quad (7) $$

Setting $\alpha \ll 1$ approximates a Neumann no-penetration boundary condition at solid horizontal and vertical boundaries.

The last term on the right hand side of equation (3) is the drag due to the permeability of the porous medium, where the mask function $\mathbf{1}$ defining the penalized solid regions is given by

$$ \mathbf{1}(x, z) = \begin{cases} 1 & \text{in the interior of the penalized solid regions,} \\ 0 & \text{in the interior of the fluid regions.} \end{cases} \quad (8) $$

Setting $\epsilon \ll 1$ in addition to $\alpha \ll 1$ approximates a no-slip Dirichlet boundary condition at solid horizontal and vertical boundaries.

To avoid numerical oscillations, a smoothed mask is defined starting from a hyperbolic tangent profile in the vertical direction (see 15).

The volume penalization can therefore approximate either a Neumann no-penetration boundary condition (if the permeability term is not included), or a no-slip Dirichlet boundary condition (using both the porosity and permeability
terms). When a no-penetration boundary condition is desired, it is sufficient to take $\alpha = O(10^{-6})$ since the accuracy of the approximation is $O(\alpha)$. In this case the penalization scheme imposes no stability constraint on the time step. In contrast, when a no-slip boundary condition is implemented, the boundary condition is imposed with an accuracy $O(\alpha\epsilon^{1/2})$. Therefore, both parameters $\alpha$ and $\epsilon$ can be varied to control the accuracy of the penalization. However, the penalization term is stiff and therefore imposes a stability constraint $\Delta t \le \epsilon$ for an explicit Euler method in time. The temporal discretization of the penalization term is discussed in §3.2.2.

## 3.2. Improvements to the original algorithm

In this section we propose some modifications of the discrete implementation of the penalization method to simplify it and improve computational performance compared with the original implementation described in Debreu et al. (2020). The principal changes concern the computation of layer thicknesses and the temporal discretization of the penalization term. In addition, we describe a modification to the traditional barotropic-baroclinic mode splitting that accounts for variable porosity $\phi$.

### 3.2.1. Computation of layer thicknesses

In Debreu et al. (2020), we showed that in order for the method to conserve total energy and, more importantly, to represent correctly the transfer between kinetic and potential energy, the following commutation property must be satisfied,

$$ \frac{\partial}{\partial t} \left( \phi(x, s, t) \frac{\partial z}{\partial s} \right) = \frac{\partial}{\partial s} \left( \phi(x, s, t) \frac{\partial z}{\partial t} \right), \qquad (9) $$

where $s$ is the generalized vertical coordinate. Eq. (9) can be rewritten to give an equation for the change in porous layer thicknesses,

$$ \frac{\partial \tilde{h}}{\partial t} = \frac{\partial}{\partial s} \left( \phi(x, s, t) \frac{\partial z}{\partial t} \right). \qquad (10) $$

In Debreu et al. (2020), Eq. (10) was implemented at the discrete level. However, this imposes a computational overhead by requiring:

(a) Computation of a new porosity function $\phi(x, s, t)$ at each time step (since in 10, the vertical coordinate $s$ depends on time).
(b) Storage of layer thicknesses for diagnostic purposes.

In practice, it is more efficient to compute layer thicknesses using the definition of porosity in the z vertical coordinate system. The thickness of layer k, between two s levels, ($s_{k-1/2}$, $s_{k+1/2}$) is given by

$$ \tilde{h}_k(x, t) = \int_{s_{k-1/2}}^{s_{k+1/2}} \phi(x, s, t)h(x, s, t)ds = \int_{s_{k-1/2}}^{s_{k+1/2}} \phi(x, s, t) \frac{\partial z}{\partial s} ds = \int_{z_{k-1/2}}^{z_{k+1/2}} \phi(x, z, t)dz. $$

As shown in Debreu et al. (2020), property (9) simply expresses the fact that the porosity $\phi$ at a fixed z does not depend on time: $\partial\phi(x, z, t)/\partial t|_z = 0$ or $\phi(x, z, t) = \phi^0(x, z)$. This naturally leads to:

$$ \frac{\partial \tilde{h}_k(x, t)}{\partial t} = \phi^0(x, z_{k+1/2}(x, t)) \frac{\partial z_{k+1/2}(x, t)}{\partial t} - \phi^0(x, z_{k-1/2}(x, t)) \frac{\partial z_{k-1/2}(x, t)}{\partial t} $$

which is the discrete version of (10). As in Debreu et al. (2020), the following two important relations follow.

$$ \tilde{H}(x, t) = \sum_{k=1}^{N} \tilde{h}_k(x) = \int_{-H(x)}^{\eta(x,t)} \phi^0(x, z) dz, $$

and

$$ \frac{\partial \tilde{H}(x, t)}{\partial t} = \phi^0(x, \eta(x, t)) \frac{\partial \eta(x, t)}{\partial t} = \frac{\partial \eta(x, t)}{\partial t}, $$

where at the free surface $z = \eta$, we have $\phi^0(x, \eta(x, t)) = 1$.

In the proposed simplification of the penalization method, at each time step, once the layer interfaces $z_{k+1/2}$ have been computed in the generalized vertical coordinate system, the layer thicknesses are diagnosed using

$$ \tilde{h}_k(x, t) = \int_{z_{k-1/2}(x,t)}^{z_{k+1/2}(x,t)} \phi^0(x, z) dz. \qquad (11) $$

This avoids having to integrate the discrete time evolution equation (10). The new kinematic relation (11) simply requires the specification of the time-independent porosity $\phi^0(x, z)$ defining the solid regions (which is based on the mask $\mathbb{1}(x, z)$ and the porosity parameter $\alpha$).

As in Guinot et al. (2018), the porosity $\phi^0(x, z)$ is first tabulated on a vertical grid $G^0(z^0_{k+1/2}(x))$ with a much vertical higher resolution than the actual computational grid (see Fig. 1). Then, the model porosity is found by integrating the high resolution porosity conservatively between two interfaces of coarser computational grid. In this paper, for simplicity we assume that $\phi^0$ is constant over each layer of the high resolution vertical grid $G_0$: $\phi^0(x, z) = \phi^0_k(x)$, $z^0_{k-1/2}(x) < z < z^0_{k+1/2}(x)$. The details of computation of $\phi^0$ are given in Appendix A. Note that the porosity
on the computational grid includes information about subgrid scale bathymetry encoded as an intermediate value of the porosity $0 < \alpha < 1$ near the fluid-structure interface. This is a kind of homogenization of fine-scale bathymetry on the scale of the computational grid. Note, however, that in this paper, the representation of subgrid bathymetry is not fully exploited. In particular, at the discrete level, the porosity and permeability at the cell interfaces are deduced (by simple averaging) from their cell-averaged values. It would be interesting to consider the ideas of Guinot (2012); Adcroft (2013) in future work.

<figure><img src="image_4.png" /><figcaption>Figure 1: Porosity $\phi^0(x, z)$ defining a higher resolution of the bathymetry. The red lines indicate the high resolution vertical grid $G_0$, the blue lines indicate the actual computational vertical grid. The reference bathymetry is smooth enough to ensure an acceptable pressure gradient error, while penalization is used to provide a more accurate representation of the true bathymetry. Note that we impose the constraint that the penalized area must always lie above the reference bathymetry.</figcaption></figure>

### 3.2.2. Temporal discretization of the penalization term

In Debreu et al. (2020), a simple explicit formulation of the penalization term was used. For simplicity, barotropic-baroclinic mode splitting was avoided (see subsection 3.2.4) which required using a small time step to satisfy the stability condition based on the external wave propagation speed $c = \sqrt{gH}$. This small time step allowed the use of
a similarly small (and highly accurate) penalization coefficient in comparison with the time scale of other physical processes, since $\epsilon \le \Delta t$.

Considering only the vertical direction, the one-dimensional penalized equation is

$$ \frac{\partial u}{\partial t} = F(u) - \frac{1}{\epsilon} \mathbb{1}(z)u. \qquad (12) $$

A splitting method can be used to integrate (12) in two steps,

$$ 1) \frac{\partial u}{\partial t} = F(u), \quad 2) \frac{\partial u}{\partial t} = -\frac{1}{\epsilon} \mathbb{1}(z)u. $$

A provisional value of $u$ at time $(n + 1)$ is first obtained:

$$ u^{n+1, \star} = u^n + \Delta t F(u), $$

and then the penalization term is included, either explicitly,

$$ u^{n+1} = \left(1 - \frac{\Delta t}{\epsilon} \mathbb{1}(z)\right) u^{n+1, \star}, \qquad (13) $$

or implicitly,

$$ u^{n+1} = \frac{1}{1 + \frac{\Delta t}{\epsilon} \mathbb{1}(z)} u^{n+1, \star}, $$

or with exact time integration,

$$ u^{n+1} = e^{-\frac{\Delta t}{\epsilon} \mathbb{1}(z)} u^{n+1, \star}. $$

The interpolating penalization (Rasmussen et al., 2011) consists in starting from the explicit version (13) with the specific choice of $\epsilon$ such that $\mu = \Delta t / \epsilon = 1$ (i.e. the minimum stable $\epsilon$) which simply leads to,

$$ u^{n+1} = (1 - \mathbb{1}(z)) u^{n+1, \star}. \qquad (14) $$

This implementation choice will be made in the following.

In practice, the smoothed version of the mask function is defined as a function of $r$, the signed distance to the solid-fluid interface scaled by the local grid size. $r$ is defined by $r = -(z - z_{bottom}) / \Delta z$ and is positive in the solid region and negative in the fluid region. As an example, in Debreu et al. (2020), the smoothed mask function $\mathbb{1}_S(r)$ was given by an hyperbolic tangent:

$$ \mathbb{1}_S(r) = \frac{1}{2} (1 + \tanh \lambda(r - r_0)), \qquad (15) $$
where the values of $\lambda$ and $r_0$ were such that $1_S(-1/2) = 1/25$ and $1_S(1/2) = 2/5$. Following Rasmussen et al. (2011), Figure (2) shows $u^{n+1}/u^{n+1,*}$ for different values of $\mu$ and for different temporal discretizations (explicit, implicit, interpolating). The main advantage of an explicit scheme with $\mu = 1$ (i.e. the above interpolating penalization) is full control of the profile of the penalized velocity near the fluid-solid boundary, which is directly linked to the shape of the smoothed mask function. On the contrary, an implicit scheme requires adjusting both the value of $\mu$ and the shape of the original mask function in order to obtain the desired velocity profile. In Figure (2), large values of $\mu$ are required for the solid region to be fully penalized (i.e. zero velocity). This advantage, however, comes at the price of an extension of the penalized domain to the fluid region (i.e. the effective position of the boundary is shifted slightly inwards). This could be corrected by increasing the value of $\lambda$ in the definition of the smooth mask function (15).

<figure><img src="image_5.png" /><figcaption>Figure 2: Different time discretizations of the penalization term for various values of the parameter $\mu = \Delta t/e$. $r = 1$ corresponds to only one vertical grid spacing $\Delta z$. The exact solution $u^{n+1}$ would be a Heaviside (step) function $H(-r)$.</figcaption></figure>

### 3.2.3. Penalization and bottom stress parameterization

With the interpolating penalization given by (14), choosing the shape of smoothed mask function is identical to the choice of the penalization itself. Therefore, in the following, we do not differentiate between the two. The shape of the smoothed mask function (15) and its coefficients $\lambda, r_0$ are somewhat arbitrary. This is particularly true for the
induced penalization value in the first cell above the ocean bottom. For realistic applications, a reasonable choice could be for this value to coincide with the one given by the vertical logarithmic profile and associated bottom stress (see subsection 2.1). Doing this, we ensure that in the limiting case where the penalization domain is empty (the bottom of the discrete domain matches the true ocean bottom), we recover the solution of the unpenalized equations.

At the discrete level, assuming a basic Euler time stepping, the bottom friction is implemented in the numerical model as

$$ (\tilde{h}u)^{n+1} = (\tilde{h}u)^n - \Delta t \tau_{bx}, \quad (16) $$

while the penalization term is implemented as described above by

$$ (\tilde{h}u)^{n+1} = (\tilde{h}u)^n(1 - \mathbb{1}(z)). \quad (17) $$

The procedure is then as follows: the bottom stress $\tau_{bx}$ is computed using the velocity at the last ocean level and at an altitude $z = z_b$ (or $r = r_b$) corresponding to this last ocean level. Then, in order for the two formulations (16) and (17) to match, we impose

$$ \mathbb{1}_S(r_b) = \min\left(\frac{\Delta t}{(\tilde{h}u)^n} \tau_{bx}, 1\right), \quad (18) $$

where $\tau$ is given by (1). Note that this includes a clipping of the penalization coefficient to ensure that its value is less than one. In the original unpenalized equations, the same clipping is applied to the bottom stress $\tau_{bx}$. The general expression of the penalization is then given by:

$$ \begin{aligned} \mathbb{1}_S(r) &= \mathbb{1}_S(r_b) + (1 - \mathbb{1}_S(r_b)) \tanh(\lambda(r - r_b)) && \text{for } r \ge r_b \\ &= 0 && \text{for } r < r_b \end{aligned} \quad (19) $$

where, by construction, $r_b = -(z_b - z_{bottom})/\Delta z \in ]-1, 0]$. The penalization coefficient thus changes from $\mathbb{1}(z_b)$ at the last ocean level to 1 inside the solid region. The maximum slope of the smooth mask function is given by $\lambda(1 - \mathbb{1}(r_b))$, and naturally decreases when the bottom stress (and thus the penalization term) is already large (close to 1) at the last ocean level. In the realistic experiments of Section 4 $\lambda$ will be taken as in Debreu et al. (2020), $\lambda \approx 1.38629$.

### 3.2.4. Mode splitting and barotropic-baroclinic correction with porosity

In contrast to the initial idealized experiments of Debreu et al. (2020), the present realistic numerical experiments take advantage of barotropic-baroclinic mode splitting. We first recall the natural definition of the penalized barotropic
238 velocity component

$$
\bar{u}(x,t) = \frac{1}{\tilde{H}(x,t)} \int_{z=-H(x)}^{z=\eta(x,t)} \phi(x,z) u(x,z,t) \, dz, \quad \tilde{H}(x,t) = \int_{z=-H(x)}^{z=\eta(x,t)} \phi(x,z) \, dz,
$$

239 or in s coordinates by

$$
\bar{u}(x,t) = \frac{1}{\tilde{H}(x,t)} \int_{s=-1}^{s=0} \tilde{h}(x,s,t) u(x,s,t) \, ds, \quad \tilde{H}(x,t) = \int_{s=-1}^{s=0} \tilde{h}(x,s,t) \, ds.
$$

240 In discrete form, the penalized barotropic velocity is

$$
\bar{u}(x,t) = \frac{1}{\tilde{H}(x,t)} \sum_{k=1}^{N} \tilde{h}_k(x,t) u_k(x,t), \quad \tilde{H}(x,t) = \sum_{k=1}^{N} \tilde{h}_k(x,t).
$$

241 In a standard primitive equation implementation of barotropic-baroclinic mode splitting, the two-dimensional
242 barotropic velocity field, which is derived from the time integration of the depth averaged equations, is used to correct
243 the three-dimensional baroclinic velocity field. Because the vertical dependence of the error between the barotropic
244 velocity and the vertically integrated baroclinic velocities is not known, this is usually done by adding a vertically
245 constant correction to the three-dimensional velocity. The corrected velocity $u_k^{n+1}$ is then derived from the provisional
246 velocity $u_k^{n+1,*}$ as

248

$$
u_k^{n+1} = u_k^{n+1,\star} + \left[ \bar{u}^{n+1} - \frac{1}{\bar{H}} \sum_k \tilde{h}_k u_k^{n+1,\star} \right].
$$

247 When a volume penalization method is used, this is of course not efficient since the correction is applied even
248 in the fully penalized solid regions (where the mask is equal to one), leading to spurious non-zero velocities. We
249 therefore propose to take into account the mask of the penalized domain by modifying the correction according to the
250 following relation

252

259

250

$$
u_k^{n+1} = u_k^{n+1,\star} + \beta_k(1 - \mathbb{1}_k) \left[ \bar{u}^{n+1} - \frac{1}{\bar{H}} \sum_{k=1}^{N} \tilde{h}_k u_k^{n+1,\star} \right].
$$

251 The coefficient $\beta_k$ is chosen to ensure that the corrected depth-integrated three-dimensional velocity is, as required,
252 equal to the barotropic velocity $\bar{u}^{n+1}$. This leads to the following constraint

252

$$
\frac{1}{\tilde{H}} \sum_{k=1}^{N} \beta_k (1 - \mathbb{1}_k) \tilde{h}_k = 1. \qquad (20)
$$

253 We then simply compute a set of weights $\beta_k$ which minimizes $\sum_k (\beta_k - 1)^2$ under the constraint (20). This leads to the
254 following computation of $\beta_k$:

254

$$
\beta_k = 1 + \frac{\sum_{k=1}^{N} \mathbb{1}_k \tilde{h}_k}{\sum_{k=1}^{N} [(1 - \mathbb{1}_k) \tilde{h}_k]^2} (1 - \mathbb{1}_k) \tilde{h}_k. \quad (21)
$$
In practice, this method slightly enhances the barotropic-baroclinic correction in the fluid region $(1-1_k)\beta_k = 1+\epsilon$, $\epsilon \ll 1$, but does not alter the 3D velocities in the solid region $(1-1_k)\beta_k \approx 0$.

# 4. Results: comparing $\sigma$-coordinates and penalized simulations

In this section we assess the qualitative and quantitative accuracy of the penalization algorithm by comparing simulations of the GS using penalized (denoted $\sigma_p$ in the following) and $\sigma$ vertical coordinates at resolutions $1/4^\circ$, $1/8^\circ$ and $1/12^\circ$. We also present results computed at resolution $1/12^\circ$ using smoother bathymetry interpolated from the $1/4^\circ$ standard configuration. This provides a basis for discussing the sensitivity of the penalization method to the choice of reference bathymetry, independent of grid resolution. When available, observational AVISO satellite data products are also included.

## 4.1. Penalized configurations

To construct the penalization we need to define the base (or reference) bathymetry ($h_{\text{base}}$) and “true” bathymetry ($h_{\text{true}}$). At a given resolution, the base bathymetry is the smooth $\sigma$-coordinate bathymetry ($r_{\text{max}} = 0.25$), while the “true” bathymetry is the base bathymetry, but without smoothing. (To avoid aliasing and ensure smoothness at the grid scale, Gaussian smoothing is applied to the true bathymetry to obtain the base bathymetry used in the $\sigma$-coordinate simulations.) The penalized area is the difference between the true and base bathymetries. We impose the constraint that $h_{\text{base}} \ge h_{\text{true}}$. This is ensured by increasing $h_{\text{base}}$ where needed.

## 4.2. Sensitivity to resolution

Figures 3 and 4 presents the 5-year mean sea surface height (SSH) and sea surface temperature (SST) for simulations at different resolutions in the standard $\sigma$-coordinate and penalized cases. In the standard $\sigma$-coordinate case, an overshooting of the GS is clearly evident at $1/4^\circ$ resolution, but improves with resolution $1/8^\circ$ and is accurate only at $1/12^\circ$. The $1/12^\circ$ solution is similar to observations, and correctly shows a clear separation at Cape Hatteras (Renault et al., 2019). We note also a loop over the Charleston Bump in the upstream region between the Florida straight and Cape Hatteras. This loop is described in the literature (Gula et al., 2015), but not represented in the AVISO product because it is too coarse. The GS path and its recirculation are also realistic at this resolution.
The above results are typical for terrain-following coordinate models applied to the GS system (Schoonover et al., 2016). An important question raised by the present study is whether penalization can improve the solution at a resolution coarser than 1/12°. Penalization does indeed appear to improve realism at lower resolutions since the 1/4° penalized run is significantly closer to the 1/12° σ-coordinate run than the 1/4° σ-coordinate run. The 1/8° penalized run is very close to the 1/12° σ-coordinate run, with a realistic GS separation. On the other hand, the 1/12° penalized run is similar to the 1/12° σ-coordinate run, which is expected if resolution convergence has been achieved. Nevertheless, the penalized run produces a more realistic GS: finer and stronger than the σ-coordinate run.

The results for the mean SST are consistent with those of the mean SSH. In the 1/4° penalized run, the GS appears to have two temperature signatures, one of which overshoots, while the other separates at Cape Hatteras. The separation is not present at all in the 1/8° σ-coordinate run. The 1/12° σ-coordinate run appears to have converged, although the penalized run again produces a better defined SST signature in the mean.
<figure><img src="image_6.png" /><figcaption>Figure 3: Mean sea surface height in CROCO simulations at different resolutions for the standard case with terrain-following ($\sigma$) coordinates (top) and penalization ($\sigma_p$) (below). The third row shows the observational AVISO product for comparison.</figcaption></figure>
<figure><img src="image_7.png" /></figure>

Figure 4: Mean sea surface temperature for the standard case with standard σ-coordinates (top) and penalization ($\sigma_p$) (below). The bottom figure shows the Satellite-based 5-km climatology (1985-2013) from NOAA Coral Reef Watch (https://coralreefwatch.noaa.gov).
Figure 5 shows the resolution sensitivity of the model for the mean eddy kinetic energy (EKE) with respect to a gridded EKE Drifter AVISO product (Laurindo et al., 2017). Since mesoscale eddies arise from the mixed barotropic and baroclinic instability of the mean flow (Gula et al., 2015), it is not surprising to find maximum EKE values along the mean GS path. At 1/4° the σ-coordinate shows a weak overshooting of the GS. It provides very little EKE and the EKE maximum is located too far north. The solution is improved at 1/8° and, in particular, at 1/12° where a relatively narrow region of EKE extends from Cape Hatteras. The penalized simulations already demonstrate clear improvement at 1/4°. The offshore EKE patch is still weak, but stronger and much better positioned than in the 1/4° σ-coordinates run. The 1/8° penalized run is close to the 1/12° σ-coordinate run, while the 1/12° penalized run shows an interesting narrowing of the EKE patch off Cape Hatteras. This is consistent with the well-known high stability of the GS trajectory at this location (Renault et al., 2019), and thus appears to be an improvement over the 1/12° σ-coordinate run.

Of particular interest is that, in the 1/8° σ-coordinate run, the GS overshoot is associated with an excess of EKE along the US east coast prior to separation. The excess energy is reduced in the 1/12° run and in the 1/8° and 1/12° penalized runs. Our interpretation is that the smoothing of the continental slope in the σ-coordinate runs releases the steering effect of bathymetry on the slope current, which can then more easily develop instabilities. Therefore, the steep slope has the effect of trapping a strong, but more stable, current which retains more inertia when it reaches the Cape Hatteras curvature. Gula et al. (2015) showed that submesoscale eddies can increase the sensitivity of the GS to topographic curvature, but this effect may not be strong enough to compensate for the loss of topographic strain. Therefore, smooth bathymetry produces a weaker mean current that can still produce excess EKE, but will overshoot due to lack of inertia. Penalization can effectively rectify the process unless the resolution is too coarse for an effective interaction between the eddies and the mean flow to produce a strong GS.
<figure><img src="image_8.png" /><figcaption>Figure 5: Mean eddy kinetic energy (EKE) in CROCO simulations at different resolutions for the standard case with terrain-following ($\sigma$) coordinates (top) and penalization ($\sigma_p$) (below). The bottom row shows the observational Drifter EKE product for comparison.</figcaption></figure>

Figure 6 shows the model sensitivity to resolution for surface currents past the separation point, compared with a 20-year climatology (1993-2012) along the Oleander transect (Rossby et al., 2014). This comparison with observational climate data underlines the alternating currents between the north-eastward GS and the south-westward slope current. In the $\sigma$-coordinate model, only the 1/12° run can reproduce the transect structure. In contrast, the position of the GS is already realistic in the penalized run at 1/4°. At this resolution in the $\sigma$-coordinate run (and to a large extent also at 1/8°), the GS is stuck to the coastal margin.
<figure><img src="image_9.png" /><figcaption>Figure 6: Mean surface velocity across the Oleander transect for the standard case with terrain-following (σ) coordinates (top) and penalization (σp) (middle). The bottom figure shows the mean near-surface velocity as measured by the Oleander (1993-2012) with shipboard ADCP (from Rossby et al. 2014).</figcaption></figure>

## 4.3. Sensitivity to bathymetry smoothing

In this section we attempt to separate the effect of resolution on eddy-mean flow interaction from its effect on flow-bathymetry interaction. We compare results for the σ-coordinate and penalized runs using standard 1/12° bathymetry and bathymetry interpolated from the 1/4° σ-coordinate configuration. For the penalized run, the reference bathymetry hbase is interpolated from the 1/4° σ-coordinate configuration. In this way we separate the effect of simulation resolution from the effect of bathymetry resolution. At the same time, this experiment illustrates the sensitivity of the
penalization procedure to the reference bathymetry.

Figure 7 (top rows) compares SSH from a standard 1/12° σ-coordinate run with a similar 1/12° run using bathymetry interpolated from the 1/4° configuration. The difference between the two 1/12° simulations is due to bathymetry smoothing, while the difference between the 1/12° run using interpolated 1/4° bathymetry and the 1/4° resolution run (Figure 3) is due to the effect of better mesoscale resolution. Comparing the 1/12° run using interpolated 1/4° bathymetry and the 1/4° run, we see that better mesoscale resolution considerably improves the GS path, but is not sufficient to prevent overshooting. In other words, the flow-bathymetry interaction is improved by the feedback of eddies on the mean flow, but using smooth 1/4° bathymetry in a 1/12° run still leads to an underestimation of the topographic strain.

We now consider similar results for penalized runs with smoother reference bathymetry interpolated from the 1/4° configuration. The results compare well with the penalized run using steeper 1/12° bathymetry. The GS is a bit stronger past the separation point, but no other difference is visible for the mean SSH. This suggests that the penalized results are relatively insensitive to the size of the penalized area. This result increases our confidence in the accuracy and performance of the penalization.

A similar analysis of EKE (not shown) confirms the previous results. The 1/12° runs using interpolated 1/4° bathymetry are closer to the standard 1/8° runs. Therefore, the improvement generally observed in the models when increasing the resolution is not due only to stronger eddy activity, but at least equally to a more realistic representation of the continental slope. This is true for both the σ-coordinate and penalized runs, but for the penalized run a resolution of 1/8° is sufficient for a realistic GS separation.
<figure><img src="image_10.png" /></figure>

Figure 7: Mean SSH for the for terrain-following ($\sigma$) coordinates (top) and penalization ($\sigma_p$) (bottom). The left column shows results using the usual 1/12° bathymetry for each case and the right column shows results using smooth bathymetry interpolated from the respective 1/4° configurations.

### 4.4. Vorticity diagnostics

The terms in the right hand side of the barotropic vorticity budget (Eq. 2) are the bottom pressure torque (BPT), nonlinear torque (ADV), planetary vorticity advection (BETA), wind stress curl (WSC), bottom drag curl (BSC), and viscous torque (D). The latter is zero outside of sponge layers in CROCO.

In the ocean interior, away from steep slopes, the Sverdrup vorticity balance between the wind stress curl and BETA term would prevail. However, with steep bathymetry near coastal margins or elsewhere, topographic stretching expressed by BPT becomes dominant in balancing BETA (Couvelard et al., 2008) — The Munk viscous torque is an avatar of BPT in simple models for closing the basin-scale vorticity budget. This balance between BPT and BETA is particularly appropriate in the GS region and the vorticity budget presented in Table 1 confirms others (Schoonover et al., 2016; Ezer, 2016).

The budget terms are integrated over the GS region, defined as in (Schoonover et al., 2016) between the 1- and 30-Sv contours, upstream of the separation. It shows that the wind stress curl and the bottom stress curl play only a
<table><caption>Table 1: The GS barotropic vorticity budget ($10^{12} m^3/s^2$) averaged over the GS region, defined between the 1- and 30-Sv contours, upstream of the separation.</caption><thead><tr><th></th><th>BPT</th><th>BETA</th><th>WSC</th><th>BSC</th><th>ADV</th></tr></thead><tbody><tr><td>σ, 1/8°</td><td>0.93</td><td>-1.78</td><td>0.05</td><td>0.24</td><td>0.54</td></tr><tr><td>σ, 1/12°</td><td>1.83</td><td>-2.48</td><td>0.09</td><td>0.27</td><td>0.29</td></tr><tr><td>σ, 1/12° smooth</td><td>1.13</td><td>-2.44</td><td>0.05</td><td>0.30</td><td>0.96</td></tr><tr><td>penalized, 1/8°</td><td>1.77</td><td>-2.65</td><td>0.08</td><td>0.48</td><td>0.26</td></tr><tr><td>penalized, 1/12°</td><td>2.24</td><td>-2.75</td><td>0.08</td><td>0.45</td><td>0.04</td></tr><tr><td>penalized, 1/12° smooth</td><td>2.16</td><td>-2.90</td><td>0.07</td><td>0.34</td><td>0.09</td></tr></tbody></table>

minor role, while the nonlinear torque is a secondary, but significant, player. The sensitivity to resolution is particularly interesting (we do not present the 1/4° results as the GS region is more difficult to define in this case). The topographic effect grows with increasing numerical resolution from 1/8° to 1/12°. There is no sign of convergence at 1/12°, but the topographic slope is steep enough at this resolution to drive a realistic GS separation.

When penalization is used, the vorticity budget is generally similar, but BPT is boosted. At 1/8°, BPT is already very close to the 1/12° standard model budget. Also, the rate of increase in BPT between 1/8° and 1/12° is lower than for the σ-coordinates run. This suggests the topographic effect converges faster with decreasing resolution when using penalization. These numbers confirm quantitatively the interpretation drawn from the SSH maps: penalization produces a physically correct GS separation at lower resolutions than σ-coordinates. With penalization, only a slight overshoot is observed at a 1/4° resolution, and this is entirely corrected at 1/8°. Ocean circulation is less sensitive to resolution when using penalization for vertical coordinates than when using σ-coordinates. The remaining sensitivity is due to more active mesoscale activity at higher resolutions.

The planetary vorticity advection (BETA) term in the penalized runs is almost as strong at 1/8° as it is at 1/12°. This is consistent with the fact that the mean slope current has more inertia at 1/8° than in the standard case, which favors its separation from the coast. The increased inertia of the slope current probably results from the steering and stabilizing effect of steeper bathymetry. This stabilizing effect, already suggested by the EKE maps (Section 4.2), is confirmed by the ADV term (mostly representing mean eddy advection over the slope region), which tends to decrease
with resolution and is similar in the penalized 1/8° run and the σ-coordinate 1/12° run.

When using interpolated 1/4° bathymetry in the 1/12° σ-coordinate run, the result is largely degraded, with BPT losing about 40% of its value. In addition, the GS overshoots even more than in the 1/4° penalized run, despite stronger mesoscale activity. On the contrary, using a smoother base bathymetry $h_{\text{base}}$ in the penalized 1/12° simulation does not significantly change the value of BPT, the structure of GS, or its separation.

Note that the bottom friction appears higher in the penalized run than in the standard run. This may be an artefact, as the value of the bottom friction term is particularly sensitive to the definition of the penalized solid region using smooth mask functions (see Sec. 3.2.2). The exact diagnosis of bottom friction is not essential to the purpose of this paper, but it should be addressed in the future.

# 5. Conclusion

The advantage of a smooth representation of bathymetry in terrain-following models is compromised by the need to avoid numerical errors on steep slopes associated with pressure gradient discretization or spurious diapycnal diffusion. Geopotential coordinate models avoid these types of numerical errors, but greatly under-represent the interaction of flow with a topographic slope, especially when the resolution of the slope is insufficient. Hybrid σ-z coordinate models are also deficient because it is difficult to achieve a satisfactory compromise between these two errors.

With volume penalization, we do not seek a compromise between the usual vertical coordinate systems. Instead, we seek a correction to these systems that realistically recovers a continuous and steep bathymetry and avoids step-like features that lead to unrealistic small-scale vorticity perturbations. Moreover, since the penalization approximates flow through a porous medium, the associated permeability and porosity parameters can be tuned to model the intensity of bottom drag, the structure of the bathymetry at the subgrid scale, or the composition of the seafloor substrate. Thus, the penalization approach to modeling bathymetry increases the realism of many aspects of ocean modeling.

The volume penalization method studied here is a modified version of the one introduced in Debreu et al. (2020). The main changes concern the calculation of the sublayer thicknesses in the penalized region and the temporal discretization of the penalization term in the momentum equations. In addition, the traditional splitting of barotropic and baroclinic modes is modified to account for the variable porosity. These modifications generally simplify the
numerical implementation of the penalization, bring more robustness and improve its computational performance for a realistic long-term simulation, while preserving its accuracy.

In this study we applied the volume penalization method to the Gulf Stream (GS) separation problem that has perplexed ocean modelers for decades. Penalization improves the representation of the flow-bathymetry interaction, allowing realistic separation of the GS even at relatively coarse resolutions. In addition, it provides a tool to separate the effects of eddy activity and topographic slope when changing grid resolution. This has never before been possible, because at coarse resolution none of the usual coordinate systems can properly represent a steep continental slope (e.g., neither z-level systems nor the terrain following σ-coordinates). Our results show that realistic bathymetry is more important than eddy activity in ensuring a realistic GS separation. This is in contrast to many recent studies, which tend to focus on eddy activity as the most important factor. A steep slope can exert a stabilizing influence that promotes a strong mean slope current with strong inertia that helps its separation from the coast at the topographic curvature of Cape Hatteras. Therefore, with penalization, the resolution required for a correct representation of the GS trajectory can be coarser than the usual recommendation for resolutions as fine as 1/50° (Chassignet and Xu, 2017). We present penalized results that show a realistic GS structure at 1/8° and a separation that already occurs at 1/4°.

We anticipate that a successful topographic slope correction will be especially valuable to climate models, as their current resolution is far from sufficient to represent western boundary currents (WBCs) using z-level coordinate systems. The error of climate models in the Gulf Stream separation and its path plays an important role in affecting the atmospheric circulation in the North Atlantic (Keeley et al., 2012), and also in remote regions via a planetary Rossby wave response (Lee et al., 2018). Our results suggest that a climate model with a 1/4° resolution using volume penalization — and perhaps also some parameterization of the eddy-mean flow interaction to energize the WBCs — would produce much more realistic simulations of ocean circulation. The associated transport of heat and salt through the global ocean would benefit greatly from a better representation of these large currents. Another area of application in the coastal ocean is more realistic representation of the seafloor substrate, for example, the vegetation in deltas and marshes, or the structures of roads and buildings in flooded urban areas. This could be achieved by carefully optimizing the porosity and permeability parameters of the penalization model.
6. Acknowledgements
---

We appreciate the support of the GdR CROCO and thank Lionel Renault for useful discussions. NKRK was
funded by an NSERC Discovery Grant.

References
---

Adcroft, A., 2013. Representation of topography by porous barriers and objective interpolation of topographic data. Ocean Modelling 67, 13–27.

Angot, P., 1999. Analysis of singular perturbations on the Brinkman problem for fictitious domain models of viscous flows. Mathematical Methods in the Applied Science 22, 1395–1412.

Angot, P., Bruneau, C.H., Fabrie, P., 1999. A penalization method to take into account obstacles in viscous flows. Numerische Mathematik 81, 497–520.

Arquis, E., Caltagirone, J.P., 1984. Sur les conditions hydrodynamiques au voisinage d'une interface milieu fluide - milieu poreux : application à la convection naturelle. C. R. Acad. Sci. Paris II 299, 1-4.

Beckmann, A., Haidvogel, D.B., 1993. Numerical simulation of flow around a tall isolated seamount. part i: Problem formulation and model accuracy. Journal of Physical Oceanography 23, 1736–1753.

Bryan, F.O., Hecht, M.W., Smith, R.D., 2007. Resolution convergence and sensitivity studies with north atlantic circulation models. part i: The western boundary current system. Ocean Modelling 16, 141–159.

Carton, J.A., Giese, B.S., 2008. A reanalysis of ocean climate using simple ocean data assimilation (SODA). Monthly Weather Review 136, 2999–3017.

Chassignet, E.P., Marshall, D.P., 2008. Gulf stream separation in numerical ocean models. Ocean Modeling in an Eddying Regime, 39–61.

Chassignet, E.P., Xu, X., 2017. Impact of horizontal resolution (1/12° to 1/50°) on gulf stream separation, penetration, and variability. Journal of Physical Oceanography 47.

Couvelard, X., Marchesiello, P., Gourdeau, L., Lefèvre, J., 2008. Barotropic zonal jets induced by islands in the southwest pacific. Journal of Physical Oceanography 38, 2185–2204.

Debreu, L., Kevlahan, N.R., Marchesiello, P., 2020. Brinkman volume penalization for bathymetry in three-dimensional ocean models. Ocean Modelling 145, 101530. URL: http://www.sciencedirect.com/science/article/pii/S146350031930174X, doi: https://doi.org/10.1016/j.ocemod.2019.101530.

Debreu, L., Marchesiello, P., Penven, P., Cambon, G., 2012. Two-way nesting in split-explicit ocean models: Algorithms, implementation and validation. Ocean Modelling 49, 1–21.

Dengo, J., 1993. The problem of gulf stream separation: A barotropic approach. Journal of Physical oceanography 23, 2182–2200.

Ezer, T., 2016. Revisiting the problem of the gulf stream separation: on the representation of topography in ocean models with different types of vertical grids. Ocean Modelling 104, 15–27.
Ezer, T., Mellor, G.L., 2004. A generalized coordinate ocean model and a comparison of the bottom boundary layer dynamics in terrain-following and in z-level grids. Ocean Modelling 6, 379-403.

Guinot, V., 2012. Multiple porosity shallow water models for macroscopic modelling of urban floods. Advances in Water Resources 37, 40-72.

Guinot, V., Delenne, C., Rousseau, A., Boutron, O., 2018. Flux closures and source term models for shallow water models with depth-dependent integral porosity. Advances in Water Resources 122, 1-26. URL: http://www.sciencedirect.com/science/article/pii/S0309170818300484, doi:https://doi.org/10.1016/j.advwatres.2018.09.014.

Gula, J., Molemaker, M.J., McWilliams, J.C., 2015. Gulf stream dynamics along the southeastern us seaboard. Journal of Physical Oceanography 45, 690-715.

Hejazialhosseini, B., Rossinelli, D., Bergdorf, M., Koumoutsakos, P., 2010. High order finite volume methods on wavelet-adapted grids with local time-stepping on multicore architectures for the simulation of shock-bubble interactions. J. Comput. Phys. 229, 8364-8383. doi:10.1016/j.jcp.2010.07.021.

Hilt, M., Auclair, F., Benshila, R., Bordois, L., Capet, X., Debreu, L., Dumas, F., Jullien, S., Lemarié, F., Marchesiello, P., Nguyen, C., Roblou, L., 2020. Numerical modelling of hydraulic control, solitary waves and primary instabilities in the strait of gibraltar. Ocean Modelling 151, 101642.

Hurlburt, H.E., Hogan, P.J., 2000. Impact of 1/8° to 1/64° resolution on gulf stream model-data comparisons in basin-scale subtropical atlantic ocean models. Dynamics of Atmospheres and Oceans 32, 283-329.

Keeley, S., Sutton, R., Shaffrey, L., 2012. The impact of north atlantic sea surface temperature errors on the simulation of north atlantic european region climate. Quarterly Journal of the Royal Meteorological Society 138, 1774-1783.

Kevlahan, N., Ghidaglia, J.M., 2001. Computation of turbulent flow past an array of cylinders using a spectral method with brinkman penalization. Eur. J. Mech./B 20, 333-350.

Kevlahan, N., Vasilyev, O., 2005. An adaptive wavelet collocation method for fluid-structure interaction at high Reynolds numbers. SIAM J. Sci. Comput. 26, 1894-1915.

Kevlahan, N.K.R., Dubos, T., Aechtner, M., 2015. Adaptive wavelet simulation of global ocean dynamics using a new Brinkman volume penalization. Geoscientific Model Development 8, 3891-3909. doi:http://10.5194/gmd-8-3891-2015.

Large, W.G., McWilliams, J.C., Doney, S.C., 1994. Oceanic vertical mixing: A review and a model with a nonlocal boundary layer parameterization. Reviews of Geophysics 32, 363-404.

Laurindo, L.C., Mariano, A.J., Lumpkin, R., 2017. An improved near-surface velocity climatology for the global ocean from drifter observations. Deep Sea Research Part I: Oceanographic Research Papers 124, 73-92.

Lee, R., Woollings, T., Hoskins, B., Williams, K., O'Reilly, C., Masato, G., 2018. Impact of gulf stream sst biases on the global atmospheric circulation. Climate Dynamics 51, 3369-3387.

Lemarié, F., Kurian, J., Shchepetkin, A.F., Molemaker, M.J., Colas, F., McWilliams, J.C., 2012. Are there inescapable issues prohibiting the use of terrain-following coordinates in climate models? Ocean Modelling 42, 57-79.

Marchesiello, P., Auclair, F., Debreu, L., McWilliams, J., Almar, R., Benshila, R., Dumas, F., 2021. Tridimensional nonhydrostatic transient rip
483 currents in a wave-resolving model. Ocean Modelling 163, 101816.

484 Marchesiello, P., Debreu, L., Couvelard, X., 2009. Spurious diapycnal mixing in terrain-following coordinate models: the problem and a solution. Ocean Modelling 46, 156–169.

485 Marchesiello, P., McWilliams, J.C., Shchepetkin, A., 2001. Open boundary conditions for long-term integration of regional oceanic models. Ocean Modelling 3, 1–20.

488 McWilliams, J., 2008. The nature and consequences of oceanic eddies. Washington DC American Geophysical Union Geophysical Monograph Series 177, 5–15.

490 Myers, P.G., Fanning, A.F., Weaver, A.J., 1996. Jebar, bottom pressure torque, and gulf stream separation. Journal of Physical Oceanography 26, 671–683.

492 Rasmussen, J.T., Cottet, G.H., Walther, J.H., 2011. A multiresolution remeshed vortex-in-cell algorithm using patches. Journal of Computational Physics 230, 6742–6755. URL: http://www.sciencedirect.com/science/article/pii/S0021999111003032, doi:https://doi.org/10.1016/j.jcp.2011.05.006.

495 Renault, L., Marchesiello, P., Masson, S., McWilliams, J.C., 2019. Remarkable control of western boundary currents by eddy killing, a mechanical air-sea coupling process. Geophysical Research Letters 46, 2743–2751.

497 Renault, L., Masson, S., Arsouze, T., Madec, G., McWilliams, J.C., 2020. Recipes for how to force oceanic model dynamics. Journal of Advances in Modeling Earth Systems 12, e2019MS001715.

499 Rossby, T., Flagg, C.N., Donohue, K., Sanchez-Franks, A., Lillibridge, J., 2014. On the long-term stability of gulf stream transport based on 20 years of direct measurements. Geophysical Research Letters 41, 114–120.

501 Sandwell, D.T., Smith, W.H., 1997. Marine gravity anomaly from Geosat and ERS 1 satellite altimetry. Journal of Geophysical Research: Solid Earth (1978–2012) 102, 10039–10054.

503 Schoonover, J., Dewar, W., Wienders, N., Gula, J., McWilliams, J.C., Molemaker, M.J., Bates, S.C., Danabasoglu, G., Yeager, S., 2016. North atlantic barotropic vorticity balances in numerical models. Journal of Physical Oceanography 46, 289–303.

505 Shchepetkin, A., 2003. A method for computing horizontal pressure gradient force in an ocean model with a non-aligned vertical coordinate. Journal of Geophysical Research 108.

507 Shchepetkin, A.F., McWilliams, J.C., 2005. The Regional Oceanic Modeling System (ROMS): A split-explicit, free-surface, topography-following-coordinate oceanic model. Ocean Modelling 9, 347–404.

509 Soufflet, Y., Marchesiello, P., Lemarié, F., Jouanno, J., Capet, X., Debreu, L., Benshila, R., 2016. On effective resolution in ocean models. Ocean Modelling 98, 36–50.
## 511 Appendix A. Computation of high resolution discrete porosity $\phi_{k_p}^0$

We first define a high resolution vertical grid with $N^{\text{poro}} \gg N$ layers, between $z = -h$ and $z = 0$.

$$z_{k_p}^{\text{poro}} = -h + h \frac{p}{N^{\text{poro}}}. $$

Then for each high resolution level $k_p \in [1, \dots, N^{\text{poro}}]$, we find the layer $k$ of the computational grid which includes

<figure><img src="image_11.png" /><figcaption>Figure A.8: The computational level $k$ containing the high resolution interface $z_{k_p}^{\text{poro}}$. The mask function associated with high resolution level $k_p$ is smoothed on the scale of computational level $k$ where it is located.</figcaption></figure>

513

514

515

$z_{k_p}^{\text{poro}}$ (see Figure A.8), and compute its thickness $\Delta z_{k_p}$. We can then compute the smoothed mask function associated with the high resolution layer $k_p$ as

$$\mathbb{1}_{k_p} = 1 - \frac{1}{2} \left( 1 + \tanh \left( \frac{1}{2} \lambda (r_{k_p} - r_0) \right) \right),$$

516

where

$$r_{k_p} = - \frac{\frac{1}{2}(z_{k_p-1} + z_{k_p}) + h_{\text{true}}}{\Delta z_{k_p}}.$$

517

Note that the smoothing scale is the thickness of the computational grid $\Delta z_{k_p}$. Finally, the porosity of layer $k_p$ is

$$\phi_{k_p}^0 = 1 + \mathbb{1}_{k_p}(\alpha - 1).$$

518

519

The porosity of the computational layer $k$ is found by integrating the high resolution porosities $\phi_{k_p}^0$ included in the computational layer of thickness $\Delta z_{k_p}$.