Representing Eddy Diffusion in the Surface Boundary Layer of Ocean Models With General Vertical Coordinates

The mixing of tracers by mesoscale eddies, parameterized in many ocean general circulation models (OGCMs) as a diffusive‐advective process, contributes significantly to the distribution of tracers in the ocean. In the ocean interior, diffusive contribution occurs mostly along the direction parallel to local neutral density surfaces. However, near the surface of the ocean, small‐scale turbulence and the presence of the boundary itself break this constraint and the mesoscale transport occurs mostly along a plane parallel to the ocean surface (horizontal). Although this process is easily represented in OGCMs with geopotential vertical coordinates, the representation is more challenging in OGCMs that use a general vertical coordinate, where surfaces can be tilted with respect to the horizontal. We propose a method for representing the diffusive horizontal mesoscale fluxes within the surface boundary layer of general vertical coordinate OGCMs. The method relies on regridding/remapping techniques to represent tracers in a geopotential grid. Horizontal fluxes are calculated on this grid and then remapped back to the native grid, where fluxes are applied. The algorithm is implemented in an ocean model and tested in idealized and realistic settings. Horizontal diffusion can account for up to 10% of the total northward heat transport in the Southern Ocean and Western boundary current regions of the Northern Hemisphere. It also reduces the vertical stratification of the upper ocean, which results in an overall deepening of the surface boundary layer depth. Finally, enabling horizontal diffusion leads to meaningful reductions in the near‐surface global bias of potential temperature and salinity.


of 20
The most prevalent mesoscale eddy parameterizations generally consist of two parts: (a) eddy-diffusive transport, where tracers are diffused along isopycnal surfaces (or surfaces of constant neutral density) using a down-gradient approach (Redi, 1982;Solomon, 1971) and (b) eddy-advective transport, where an additional advection of tracers by the eddy-induced velocity acts to flatten isopycnals, thereby reducing potential energy (Gent & McWilliams, 1990;Gent et al., 1995;Griffies, 1998). However, in the surface and bottom boundary layers vigorous microscale turbulence is induced by several processes, including breaking surface and internal waves, destabilizing buoyancy fluxes (e.g., Taylor & Ferrari, 2010), and boundary stresses (e.g., Thomas, 2005;Thomas & Ferrari, 2008). The resultant mixing of the boundary layer stratification makes the notion of neutral surfaces tenuous. Large-scale turbulent motions are thus able to break the constraint of following any "minimal-disturbance surface" (e.g., Fox-Kemper et al., 2013), and their direction instead becomes governed by the boundary itself. As argued by Tréguier et al. (1997) and Ferrari et al. (2008Ferrari et al. ( , 2010, it is for these reasons that as a boundary is approached, the eddy fluxes rotate to become parallel to the boundary rather than acting along neutral surfaces. In the context of eddy parameterizations, the parameterized fluxes thus must be oriented in a purely horizontal direction near the surface and along bathymetric contours near the bottom. Ideas for applying horizontal mesoscale eddy fluxes in the surface boundary layer were first proposed by Tréguier et al. (1997). Following these ideas, horizontal eddy diffusion in the upper boundary layer was implemented in the Modular Ocean Model versions 3 and 4 (MOM3 and MOM4, Griffies et al., 2005;Pacanowski & Griffies, 2000). Subsequently, Ferrari et al. (2008) derived a new eddy parameterization where the diabatic nature of the eddy fluxes could be retained near ocean boundaries. In this parameterization, eddy-induced velocity and diffusion are set parallel to the boundary within the turbulent boundary layer, while in the ocean interior eddy fluxes still occur along neutral planes. The two regimes are matched in the so-called transition layer, where the fluxes are progressively tilted from being aligned with the neutral slope to being purely horizontal as the boundary layer is approached. This method was implemented and tested in a climate model, leading to improvements in the solution when compared to the results using the traditional tapering approach (Danabasoglu et al., 2008).
Imposition of such near-surface eddy fluxes across neutral planes is now a common practice in ocean general circulation models (OGCMs) with geopotential vertical coordinates, where horizontal becomes synonymous with along-layer, such as the Parallel Ocean Program version 2 (POP2, Danabasoglu et al., 2012), the Modular Ocean Model versions 4 and 5 (MOM4 and MOM5, Griffies, 2012;Griffies et al., 2005), the Finite element/ volume Sea ice-Ocean Model version 1.4 (FESOM, Wang et al., 2014), and the Nucleus for European Modeling of the Ocean version 3.6 (NEMO, Madec et al., 2017). In these models the transition from neutral to dianeutral physics is handled by a near-boundary eddy flux parameterization (Danabasoglu et al., 2008;Ferrari et al., 2008). This approach has become a fundamental element of coupled Earth system models that use these OGCMs. However, the recent trend toward OGCMs with general (or Lagrangian) vertical coordinates, such as the MOM version 6 (MOM6, Adcroft et al., 2019) and the Model for Prediction Across Scales-Ocean (MPAS-O, Petersen et al., 2018), makes a similar implementation of diabatic near-surface eddy fluxes much more complex because the coordinate surfaces can be tilted with respect to the horizontal. Furthermore, this tilt is determined dynamically, and it evolves in time. Although applying neutral diffusion in such models is also non-trivial, recently Shao et al. (2020) have developed a neutral diffusion operator that is appropriate for this class of OGCMs. Therefore, only the horizontal diffusive part remains missing.
To fill this gap, we present a method for applying horizontal diffusive tracer fluxes due to mesoscale eddies within the surface boundary layer of general vertical coordinate ocean models. To make the method useable for climate simulations, we ensure that it satisfies the following requirements: (a) the total tracer content is conserved; (b) tracer monotonicity is preserved (i.e., it does not create new tracer extremes); and (c) it does not lead to a significant increase in the computational cost.
This manuscript is organized as follows. A description of the method is presented in Section 2. In Section 3 we use an idealized test case to provide a proof of concept for how this method operates in conjunction with the neutral diffusion scheme developed by Shao et al. (2020). In Section 4 we explore the effects of including surface horizontal diffusion in global forced ocean-sea-ice simulations configured using two different vertical coordinate systems within the Community Earth System Model version 2 (CESM2) framework. A summary and discussion are given in Section 5.
MARQUES ET AL.

Description of the Method
In this section, we present an algorithm for applying horizontal eddy diffusion in the surface boundary layer of models using Lagrangian vertical coordinates. This algorithm accounts for the effects of mesoscale diffusive fluxes whose direction is constrained by the geometry of the ocean surface. That is, unlike eddy fluxes in the ocean interior that are applied along neutral-density surfaces, the following method imposes fluxes that are purely horizontal.
We will denote the parameterized diffusive flux of an arbitrary, averaged tracer ϕ as The nature of the averaging is not crucial for describing the method and is deliberately left unspecific here; we only require that the average is taken over sufficiently many samples of ϕ so that there is statistical convergence (i.e., adding another sample will not meaningfully change the average). The parameterized flux is directed down the mean gradient of ϕ using a symmetric diffusion tensor, K H , and the fluxes are set to be purely horizontal by choosing where x, y, and z are the zonal, meridional, and vertical (positive upward) directions, respectively, t is time, and κ H (x, y, z, t) is a user-specified scalar diffusion coefficient that varies in time and space. Note that K H is equivalent to the Redi (1982) tensor in the limit where the neutral slopes are zero and there is no dianeutral flux. For simplicity of presentation, we will assume an isotropic diffusion that is directed along the model coordinates; a more general, anisotropic prescription would replace κ H and the identity matrix in the upper left minor of K H with a symmetric matrix consisting of unequal diffusivities (Bachman et al., 2020;Smith, 1999).
Given that layer thicknesses can vary between two horizontally adjacent cells in a Lagrangian vertical coordinate model, applying commonly used rotation operators to ensure that the fluxes transition from neutral to horizontal is not possible. To overcome this issue, we propose a method that uses regridding/remapping techniques before applying diffusive fluxes.

Step 1: Regridding/Remapping
The first step is to define a new vertical grid, which we refer to as the Horizontal Boundary Diffusion grid (HBD hereafter), using layer interfaces from the native vertical grid, the boundary layer depth (BLD), and the maximum depth (H) from two adjacent water columns ( Figure 1a). In this manuscript, BLD is based on a critical Richardson number value, following the definition given in the K-profile parameterization (KPP) for vertical mixing (Large et al., 1994). However, in practice, a MLD can be used instead of BLD to define the depth of the surface diabatic layer in models that do not carry/use a BLD (see Danabasoglu et al. (2008)). Using the above-mentioned information, a rescaled height (z*) coordinate (Adcroft & Campin, 2004;Stacey et al., 1995) is constructed by combining layer thicknesses (h) and BLDs from both columns, starting at z* = 0 ( Figure 1b). Duplicate interface values are removed and the maximum depth (z max ) in the HBD grid is set to where BLD max = max(BLD L , BLD R ) is the deeper BLD of the two columns, H min = min(H L , H R ) is the shallower total depth, and the subscripts L and R refer to the left and right columns, respectively. The HBD grid has enough resolution to correctly represent tracer concentrations in both columns as well as the tracer fluxes (at velocity points) between the two columns. It is not possible to define a unique HBD grid that works globally because An example of the HBD grid, which is constructed by combining the distance between interfaces (i.e., layer thicknesses) and the boundary layer depths from both profiles shown in panel (a). The first interface in the HBD grid is set to z = 0 (green circle) and the last interface (z max ) is calculated using BLD L , BLD R , and the depths of the left and right columns (H L and H R , respectively) as given by Equation 3.

of 20
interfaces can move with time. Therefore, we define a new HBD grid for each pair of water columns at every tracer time step. The HBD grid is stored in two 3-dimensional arrays defined at the zonal and meridional velocity points of the Arakawa C-grid (Arakawa & Lamb, 1977). The number of vertical layers in these arrays (NZ) is set to NZ = (2 × NK) + 2, where NK is the number of vertical layers in the native grid. This gives enough points to represent the HBD grid even when all the thicknesses from two adjacent columns are distinct and the BLD in one of the columns reaches the maximum depth of the ocean.
Once the HBD grid is defined, tracer profiles from each neighboring cell are remapped onto this grid using a remapping scheme. The remapping must be both conservative and monotonic. Conservative remapping is necessary to preserve the integrated value of scalar concentration to within machine precision. This is particularly important for long-duration climate simulations running for centuries or millennia, where the accumulation of spurious tracer content can significantly influence the solution. Monotonic remapping assures that no overshoots or new extrema are created. This is crucial for ocean variables that must be bounded (e.g., seawater salinity) because the lack of monotonicity can lead to undesirable effects, such as triggering nonphysical convective adjustments. The reader is referred to White and Adcroft (2008) for an overall description of high-order remapping schemes, including the piecewise parabolic method used throughout this manuscript. Note that, once a tracer is remapped to the HBD grid, "horizontal" becomes synonymous with "layer-wise," and correctly orienting the diffusive fluxes becomes much more straightforward. Figure 2a shows two adjacent water columns where the tracer field has been remapped to the HBD grid shown in Figure 1b. This figure is used to describe the steps listed below.

Finding Vertical Indices Containing the Boundary Layer Depths
A key priority is to ensure that the horizontal diffusion in the boundary layer tapers smoothly to neutral diffusion in the ocean interior, with the transition point occurring at the shallowest BLD. To this end, the first step in this part of the algorithm is to find the vertical indices of the layers containing BLD in both columns. For the scenario shown in Figure 2a, these are k = 9 and 4 for the left and right columns, respectively. These vertical indices are then compared and the minimum (k min ) and maximum (k max ) indices are identified, that is, k min = 4 and k max = 9. Note that, by definition, k max is always the vertical index of the deepest layer in the HBD grid. Between k min and k max , we impose a transition zone where the strength of the horizontal fluxes is gradually tapered to zero as one moves downward (toward larger k). Figure 2a shows the presence of a transition layer, but this layer can be absent in certain situations (e.g., when the BLD is the same in both columns, i.e., when k min = k max ).

Computing Fluxes at Each Layer
For each vertical layer index k on the HBD grid, the diffusive flux (here presented for the zonal direction, with a meridional flux obtained analogously) is calculated at the interface between the two columns as where h(k) is the layer thickness, ϕ R (k) and ϕ L (k) are tracer concentrations on the right and left cells, respectively, and Δx and Δy are the width of the tracer cell in the zonal and meridional directions, respectively. Fluxes are calculated for each layer beginning at the surface layer (k = 1) and ending at the layer bounded by the deepest BLD (k = k max ). This is illustrated by the black arrows in Figure 2a. Note that, assuming κ H > 0, the diffusive fluxes are always layerwise down-gradient.

Tapering the Fluxes in the Transition Layer
If there is a transition layer, that is, k max > (k min + 1), the diffusive horizontal flux decays linearly between k min and k max . In this zone, Equation 4 becomes where FT x (k) is the diffusive flux in the transition layer. The nondimensional horizontal tapering coefficient α H is where H T is the distance between the bottom interfaces of the layers with indices k min and k max , z is the mean depth of the layer (the depth of the tracer point), and BLD min = min(BLD L , BLD R ). The linear decay imposed in Equation 5 ensures that the fluxes are at full strength at z = −BLD min and become zero at z = −(BLD min + H T ). Similarly, a nondimensional neutral tapering coefficient α N is defined as Unlike α H which is defined at the cell centers of the HBD grid, α N is defined at interfaces of the native grid and then used to scale the neutral diffusion fluxes. To guarantee conservation, this is done by performing a simple 4-point average within each layer bounded by neutral surfaces between two water columns. Because of this, it is not possible to assure that neutral diffusive fluxes are zero above BLD min , but the non-zero values should be limited to the layer where BLD min is located (see Section 3). Figure 2b shows the horizontal and neutral tapering coefficients for the scenario depicted in Figure 2a. There is no physically based justification for applying a linear profile for these coefficients, and this choice was made because this is the simplest way to implement tapering. The same linear approach was used in Danabasoglu et al. (2008) and we do not expect any significant changes in the results if a different decay (e.g., quadratic) is applied.

Remaping Fluxes Onto the Native Grid
At this point in the algorithm, the thickness-weighted horizontal fluxes calculated in Section 2.2 must be transferred onto the model's native vertical coordinate. These fluxes are vertically extensive quantities (i.e., a quantity i.e., already volume weighted with respect to the vertical axis) and so the regridding/remapping approach used to transform the tracer concentrations in Section 2.1 is no longer appropriate. Instead, we apply a one-dimensional conservative remapping that ensures that vertical integrals with the same depth extent are conserved between the two vertical coordinates. From now on we follow the convention that (˙) represents the function specific to the model's native vertical grid and variables without the dot represent the function in the HBD grid.
A thickness-weighted flux is constant between the interfaces of the discretized HBD vertical grid and so this constraint is satisfied by a simple binning approacḣ where ̇( ) is the thickness-weighted diffusive flux in layer n on the model's native grid, k is a layer index on the HBD grid, a(k) is the fraction of that layer which falls within the depth range spanned by layer n, and F(k) 6 of 20 is the diffusive flux on the HBD grid. This method only has first-order accuracy but seems to be sufficient in prac tice (see the experiments in Section 3). Higher-order methods can be constructed by adding additional integral constraints (e.g., smoothness criteria based on vertical gradients).
In this procedure, the target layers in the destination grid are the layer thicknesses at cell interfaces, which is where fluxes are computed. To avoid non-zero fluxes in the presence of vanished layers, the layer thicknesses at cell interfaces are computed using the harmonic mean of two neighboring thicknesses at tracer points. Once the fluxes are remapped onto the target grid, a flux limiter must be applied to avoid up-gradient fluxes and maintain monotonicity. The maximum diffusive flux between two adjacent cells (̇m where ̇( ) and ̇( ) are the volume of the left and right cells at vertical level n, respectively. The non-dimensional constant c is set to 0.2 in the simulations shown in Sections 3 and 4, which represent the maximum fraction of the tracer concentration that can be isotropically diffused to each neighboring cell. That is if a cell has an initial concentration of one and this cell is connected to four neighboring cells with zero initial tracer concentration, the final (i.e., after steady state) tracer concentration in all cells will be 0.2. The diffusive flux is then limited as follows: • If ̇×̇m ax < 0 , ̇= 0 ; • If ̇×̇m ax > 0 and ̇>̇m ax , ̇=̇m ax ; • If ̇×̇m ax > 0 and ̇<̇m ax , ̇ is not modified. The last step is to add the diffusive fluxes to the tracer tendency and advance in time.

Summary of the Algorithm and Additional Remarks
The main steps for the algorithm are depicted in Figure 3. It starts with the tracer concentration in a pair of adjacent water columns on the native model grid (Figure 3a). The first step is to construct the HBD grid using layer thicknesses and the boundary layer depths from both columns, and then conservatively and monotonically remap the initial tracer concentration onto this grid (Figure 3b). The horizontal tracer fluxes are computed on this grid by simply calculating the diffusive fluxes within each layer according to Equation 4. These fluxes are then remapped to the velocity points on the native grid. After applying a flux limiter following Equation 9, the fluxes are added to the tracer tendencies at each point and the tracer field is advanced in time. This procedure is repeated for every pair of water columns in both zonal and meridional directions and at every tracer time step.

Proof of Concept Using Idealized Simulations
The algorithm presented in Section 2 is implemented in MOM6 (Adcroft et al., 2019), which uses a vertical Lagrangian-remap algorithm that enables general vertical coordinates (see Section 4.1 for additional details about this model). To show how the algorithm described in the previous section works in practice, we consider a . The neutral diffusion algorithm is described in Shao et al. (2020), and it is modified here to only act below z min ; neutral diffusion fluxes linearly increase from zero at z min to full strength at z max (see Section 2.2.3). Both, horizontal and neutral diffusion coefficients are set to 1,000 m 2 s −1 .
The horizontal domain is 200 × 100 km in the x and y directions, respectively, with a constant grid spacing Δx = Δy = 100 km. The east and west boundaries are closed and the north and south boundaries are periodic. The ocean bottom is flat and the maximum depth (H max ) is 500 m. The initial potential temperature field, θ(z, x), is defined as where Δθ = 15°C. A zonal gradient is imposed by setting θ surf (x) = 20°C at x = 50 km and θ surf (x) = 19°C at x = 150 km. We also define a salinity field, S(z), which acts as a passive tracer and only varies vertically: where ΔS = 1 psu and S surf = 35 psu. The initial conditions for both θ and S are shown in Figures 4a and 4b, respectively.
A linear equation of state is applied so that where ρ is the in situ density, ρ 0 = 1,035 kg m −3 is the reference density, and ∂ θ ρ = −0.255 kg m −3 °C −1 . This leads to a constant isopycnal slope of 3.3 × 10 −4 everywhere (Figure 4a).  We now evaluate layer-integrated tracer tendency (concentration times thickness) profiles due to horizontal and neutral diffusion after one tracer time step. These profiles are taken at x = 50 and 150 km. Layer-integrated tendencies are used because MOM6 is formulated to conserve the vertically extensive tracer content. The algorithm is not going to give identical tendencies for HBD-Z and HBD-H because the tracer values, and therefore the diffusive fluxes, are dictated by the layer thicknesses. Instead, the goal of the algorithm is to ensure that diffusion is purely horizontal near the surface and shifts from horizontal to neutral within a transition layer, regardless of the vertical coordinate. In HBD-Z, the left and right tendencies are symmetric because layer thicknesses do not vary horizontally in the z* grid. However, in HBD-H the tendencies are symmetric only near the surface where layer thicknesses do not vary horizontally, and non-symmetric where layer thicknesses along the same vertical index vary between the left and right columns (−300 < z < −45 m). This non-symmetric aspect of the HBD scheme is a consequence of the conservative remapping applied when transferring the fluxes from the HBD grid 9 of 20 onto a native vertical grid that has horizontally varying thicknesses (Section 2.3). We should note that the magnitude of the column-integrated heat (or any other tracer) content is always identical in the left and right columns.
In HBD-Z and HBD-H, the tendency in heat content from horizontal diffusion decays linearly within the transition zone and then vanishes below the deepest BLD (z = −300 m). By design, the tendency in heat content from neutral diffusion is zero throughout the entire water column in both experiments (Figure 5b). Similarly, the tendency in salt content due to horizontal diffusion is also zero everywhere (Figure 5c). This is despite the fact that in HBD-H, S can vary horizontally for a fixed vertical index because thicknesses are not the same in the left and right columns. However, diffusive fluxes are computed after the tracers are remapped to the HBD grid, where ∂ x S = 0.
The tendency in salt content due to neutral diffusion is shown in Figure 5d. Again, the left and right profiles differ because of the difference in layer thicknesses in each experiment. However, a linear pattern within the transition zone occurs in both experiments. Notice that the right column in HBD-H (purple line) has a non-zero value within the surface zone (z > −100 m). As explained in Section 2.2.3, this is because α N is defined at the interfaces of the native grid, and BLD min = 100 m falls between two interfaces. Nevertheless, only the layer that is intercepted by BLD min has a non-zero value. To satisfy the condition for iso-neutrality, the diffusive flux due to neutral diffusion can also be "non-local" (non-symmetric), as in Figure 5d.

Model Description
Global simulations are performed using the CESM2 framework  with active ocean and sea-ice components.
The ocean model is MOM6 (Adcroft et al., 2019), which is the same model used in Section 3. MOM6 uses an Arbitrary-Lagrangian-Eulerian (ALE) algorithm in the vertical, which allows the application of any vertical coordinate system (e.g., geo-potential, isopycnal, terrain-following, or any combination of them). Errors due to remapping are minimized via high-order accurate reconstructions (White & Adcroft, 2008;White et al., 2009). Unless otherwise stated, MOM6's dynamical core has been configured following Adcroft et al. (2019). Below is a brief description of the sub-grid-scale parameterizations employed in the present study. We attempt to keep these parameterizations and their settings as close as possible to what has been used in recent simulations with POP2 within CESM (e.g., Danabasoglu et al., 2020;Tsujino et al., 2020). We emphasize that the parameters and choice of physics in MOM6 for CESM are continuously evolving and, therefore, the description below reflects the configuration employed when the experiments presented here were conducted.
The KPP parameterization for vertical mixing (Large et al., 1994) is incorporated via the Community ocean Vertical Mixing (CVMix, Griffies et al., 2015) framework. In addition to accounting for the mixing in the surface boundary layer, KPP is also used to represent the vertical mixing in the ocean interior due to convection, double-diffusion, and vertical shear of the horizontal velocity. The latitude-dependent diffusivity due to internal wave mixing defined in Danabasoglu et al. (2012) is included. Energy dissipation from tidally induced breaking internal gravity waves is represented using the scheme developed by Simmons et al. (2004).
The restratifying effects of baroclinic eddies in the mixed layer are represented using the parameterization developed by Fox-Kemper et al. (2008) as implemented by Fox-Kemper et al. (2011). The MOM6 implementation of this scheme has been modified as described in Adcroft et al. (2019). We set the frontal length scale to 1 km and the efficiency coefficient to 0.0625. The MLD used in this scheme is calculated based on a potential density change of 0.03 kg m −3 from the surface. To ensure that restratification of the deepest mixed layer is persistent, a running-mean filter with a time scale of 5 days is applied to MLD.
In addition to the near-surface horizontal eddy diffusion scheme that is the focus of this manuscript, mesoscale eddies are represented by activating two additional schemes in the tracer equation. The first scheme follows the ideas of Gent and McWilliams (1990), where available potential energy is removed from the large scale by flattening isopycnals (hereafter GM). This scheme is implemented using the stream function formulation of Ferrari et al. (2010). By following what is commonly done in layer models (e.g., Bleck, 2002), the associated MARQUES ET AL.

10.1029/2023MS003751
10 of 20 eddy-induced transport is applied as a bolus velocity. To avoid the problems associated with layer thickness diffusion described by Holloway (1997), the scheme is implemented as an interface height diffusion. The second scheme applies the diffusive mixing of tracers along neutral directions following the idea of Solomon (1971) and Redi (1982). The implementation of the neutral diffusion algorithm in MOM6 is described in Shao et al. (2020) and we chose the continuous reconstruction option for the present study. As mentioned in Section 3, we have modified this scheme to taper up the fluxes between the minimum and maximum boundary layer depths (see Section 2.2.3).
The mesoscale tracer diffusivity (κ) is computed as where α is a non-dimensional constant set to 0.07, H is the total thickness of the water column, and N 2 = ∂ z b and M 2 = |∇ h b| are the vertical and horizontal buoyancy (b) gradients, respectively. The depth-averaged mesoscale eddy kinetic energy (e) is computed using a prognostic equation (Jansen et al., 2015). Equation 13 is based on the geometric formalism of Marshall et al. (2012), except that we employ the eddy kinetic energy instead of the full (kinetic plus potential) eddy energy. The eddy kinetic energy field is initialized by assuming an instantaneous balance between the bottom friction and the baroclinic source terms (Equations 2 and 3 in Jansen et al. (2015)), which yields a simple algebraic expression for the eddy kinetic energy that is based on the stratification parameters and bottom drag coefficient. The eddy kinetic energy prognostic equation is then iterated forward in time to predict the evolution of e and hence the diffusivity (κ). The same diffusivity field given by Equation 13 is used by the neutral and horizontal diffusion schemes. The diffusivity used in the GM parameterization is κ GM = κ × EBT(z), where a vertical structure is added via the equivalent barotropic mode, EBT(z) (Adcroft et al., 2019). At the surface both κ and κ GM have the same values.
Viscous terms are added to the horizontal momentum equation using both Laplacian and biharmonic operators with coefficients set via the eddy kinetic energy scheme. Momentum is extracted via a quadratic drag law with a constant bottom friction coefficient C d = 3 × 10 −3 . The non-linear equation of state for seawater defined by Wright (1997) is applied.
The sea-ice component is CICE Version 5.1.2 (CICE5, Bailey et al., 2018), with the improvements listed in Danabasoglu et al. (2020). With the exception of the horizontal grid, CICE5 has been configured following the description for the CESM-POP model provided in Tsujino et al. (2020).

Experimental Design
The four global forced simulations conducted here are summarized in Table 1. They differ in terms of the vertical coordinate system (hybrid or z*, the same coordinates employed in Section 3), number of vertical layers (NK), whether neutral diffusion is applied throughout the entire water column or just below BLD min , and whether the horizontal diffusion scheme is enabled. The z* vertical coordinate has 65 layers with Δz = 2.5 m down to 10-m depth. The vertical resolution follows a hyperbolic tangent function where Δz increases to 250 m at ∼3,000-m depth, remaining constant until the bottom is reached. The hybrid vertical coordinate has 75 layers and is a combination of z* near the surface and potential density (referenced to 2,000 dbar) elsewhere. The depth of transition between z* to isopycnal is shallowest in the tropics (∼50 m) and deepens toward high latitudes (∼1,200 and 2,000 m in the Southern and Northern Hemispheres, respectively).
The simulations start from rest with the initial potential temperature and salinity fields from the January-mean climatology of the World Ocean Atlas 2018 (WOA18, Locarnini et al., 2018;Zweng et al., 2019). The sea surface salinity is restored to the monthly climatology of the upper 10-m averaged salinity from WOA18 using a piston velocity of 0.1667 m/day, corresponding to ∼300 days over 50 m.
Both the sea-ice and ocean models use the same tripolar horizontal grid with a nominal resolution of 2/3° and equatorial refinement of 1/4°. Bottom topography and coastlines are derived from the ETOPO1 data set. The minimum and maximum depths are set to 10 and 6,000 m, respectively. The simulations are forced using the JRA55-do v1.3 data set (Tsujino et al., 2018) and the total integration time is one forcing cycle (1958-2018; a total of 61 years). Unless otherwise stated, the results presented in the next section have been averaged over the last 31 years of the integrations.

Results
In this section, we compare the simulations listed in Table 1 focusing on the impacts of the near-surface horizontal diffusion scheme outlined in Section 2 on several climate-relevant oceanic metrics.

Winter-Mean Surface Boundary Layer Depth
We start by evaluating how BLD is modified when near-surface horizontal diffusion is included. We will focus on the winter-mean values because this is when BLDs are the deepest and the effects of adding horizontal diffusion are more pronounced; recall that the horizontal diffusion scheme is only acting within BLD.
The winter-mean BLDs are shown in Figure 6. The choice of a vertical coordinate has a strong effect in BLD, with experiments using the z* coordinate (Figures 6a and 6b) showing overall shallower values than experiments using the hybrid coordinate (Figures 6d and 6e) that are reflected in the respective global-mean BLDs. However, for either choice of vertical coordinate system, adding horizontal diffusion deepens BLDs almost everywhere (Figures 6c and 6f). The effect is more pronounced in experiments using the hybrid coordinate (the maximum difference is ∼130 m, Figure 6f) versus in experiments using z* (the maximum difference is ∼53 m, Figure 6c).
Differences in the time-averaged diffusion coefficient between experiments with and without horizontal diffusion are relatively small and cannot account for these differences in BLD (Appendix A). The deepening of BLDs in experiments with horizontal diffusion is more pronounced in regions that tend to have relatively deep winter-mean values, such as in the Labrador, Greenland, and Norwegian Seas as well as in the Southern Ocean (Figures 6c and 6f).

Potential Temperature and Salinity
We next show winter-mean vertical profiles of potential temperature, salinity, and the square of buoyancy frequency (computed online as 2 = − −1 0 , where g is the gravitational acceleration), focusing on two regions where BLD differences are large in comparison to the respective control simulations: region #1 is located in the Labrador Sea (55°-60°N; 50°-55°W, the black box in Figure 6f) and region #2 is located in the Southern Ocean (55°-60°S; 100°-120°W, the red box in Figure 6f).
In the Labrador Sea, including near-surface horizontal diffusion leads to a slightly better representation of θ (Figures 7a and 7b) and S (Figures 7c and 7d) profiles when compared to the WOA18 winter-mean climatology. Indeed, in general, the bias reduction is relatively small for θ and S profiles from the z* case (Figures 7a and 7c), as well as for the θ profile from the hybrid experiment ( Figure 7b). The improvement is more evident when comparing S profiles of the hybrid experiments below ∼75 m (Figure 7d). Another consequence of applying horizontal diffusion is the overall slight reduction in the vertical stratification of the upper ocean within this region, which occurs in both hybrid and z* cases (Figures 7e and 7f).
Near-surface horizontal diffusion does not reduce biases everywhere and for all tracers. In the Southern Ocean, adding horizontal diffusion increases the biases in the winter-mean S profiles (Figures 8c and 8d) when compared to the WOA18 climatology, regardless of the vertical coordinate system. As for the Labrador Sea, the bias in the winter-mean θ profiles (Figures 8e and 8f) is also slightly reduced because of horizontal diffusion. Most importantly, horizontal diffusion still leads to an overall reduction in the vertical stratification of the upper ocean (Figures 8e and 8f). We note that BLD is a measure of the depth over which turbulent boundary layer eddies can penetrate before becoming stable relative to the local velocity and buoyancy. Therefore, the above results suggest that the deepening of BLDs in experiments with horizontal diffusion enabled (HBD-Z and HBD-H) is a consequence of the reduction in the vertical stratification.
So far we have seen that in the two regions considered above adding near-surface horizontal diffusion can either increase or decrease biases in θ and S. To further assess the broader effects of horizontal diffusion on these tracers, we show the time evolution of global bias profiles of θ ( Figure 9) and S (Figure 10), focusing on the 12 of 20 upper 1,500 m of the water column. CTRL-Z and CTRL-H show an overall warming bias in θ between 100 and 800 m (Figures 9a and 9b). Toward the end of the simulation, the maximum bias in CTRL-Z (∼1°C, Figure 9a) is larger than the maximum bias in CTRL-H (∼0.7°C, Figure 9b). Enabling near-surface horizontal diffusion leads to an overall reduction in the θ bias with respect to the control experiments (Figures 9c and 9d). The bias reduction is relatively small in the z* cases (up to ∼2%, Figure 9c) when compared to the hybrid cases (up to ∼10%, Figure 9d).

of 20
In terms of global salinity biases, CTRL-Z and CTRL-H show a fresh bias between 0 and 300 m (Figures 10a  and 10b). The overall bias pattern is again stronger in CTRL-Z (∼−0.12 psu) when compared to CTRL-H (∼−0.07 psu). When near-surface horizontal diffusion is included, this upper-ocean salinity bias is reduced by up to ∼20% in both HBD-Z and HBD-H (Figures 10c and 10d).

Global Northward Heat Transport
Finally, we explore the impacts of the near-surface horizontal diffusion scheme on ocean heat transport. Figure 11a shows the total global northward heat transport, which represents the sum of the advection, neutral, and horizontal diffusion components. Overall, experiments using the hybrid coordinate (CTRL-H and HBD-H) display a stronger poleward heat transport in both hemispheres. The effects of the horizontal diffusion scheme on the total northward heat transport are not noticeable in Figure 11a because advection alone accounts for the majority of the total transport almost everywhere (Figure 11b). Contributions from both neutral ( Figure 11c) and horizontal (Figure 11d) diffusion become notable in the Southern Ocean and in the Western boundary current regions of the Northern Hemisphere. These are the regions that tend to display relatively large eddy Figure 7. Winter-mean (January, February, and March) vertical profiles of potential temperature (θ, panels (a, b)), salinity (S, panels (c, d)), and the square of buoyancy frequency (N 2 , panels (e, f)) averaged within a region in the Labrador Sea (55°-60°N; 50°-55°W, see the black box in Figure 6f for location) and averaged over years 31-61. Corresponding wintermean profiles taken from the WOA18 climatology are shown in black for comparison. The dashed and solid gray lines in panels (a and b) represent the time-averaged boundary layer depth for cases with and without horizontal diffusion, respectively. Only the upper 300 m of the water column is shown.

10.1029/2023MS003751
14 of 20 diffusivities ( Figure A1) and deep BLDs ( Figure 6). Despite the fact that horizontal diffusion is only active within BLD, its contribution is comparable to that of neutral diffusion regardless of the coordinate system (compare Figures 11c and 11d). The transport due to horizontal diffusion is not sensitive to the choice of vertical coordinate (Figure 11d). However, because the advective transport is overall stronger in experiments with a hybrid coordinate, the contribution from horizontal diffusion to the total transport is stronger (up to 10%) in HBD-Z when compared to HBD-H (up to 5%).

Summary and Discussion
We have developed a new algorithm for applying horizontal eddy diffusion within the surface boundary layer of general vertical coordinate system ocean models. The algorithm uses regridding/remapping techniques to represent tracer profiles in a geopotential vertical coordinate, where horizontal fluxes are easily calculated and then remapped back to the native grid. Combined with a neutral diffusion operator appropriate for this class of models (e.g., Shao et al., 2020), the algorithm can be used to represent the transition from neutral to dianeutrally oriented diffusive fluxes in ocean models. The effect is equivalent to what is achieved via the near-boundary eddy flux parameterization (Danabasoglu et al., 2008;Ferrari et al., 2008) that has been implemented in OGCMs with geopotential vertical coordinates (e.g., Danabasoglu et al., 2012;Griffies, 2012;Griffies et al., 2005;Wang et al., 2014).
The algorithm was implemented in a general vertical coordinate OGCM (MOM6) and we have run a set of forced global experiments using the CESM2 framework to assess the effects of including horizontal diffusion in two different vertical coordinate systems (z* and hybrid). Horizontal diffusion leads to a reduction of the vertical stratification within the surface boundary layer of certain regions, which results in an overall deepening of BLD. This feedback is more pronounced in the hybrid experiments, where the winter-mean BLD can be more than 100 m deeper over a 31-year period. The maximum winter-mean BLD deepening in the z* experiment with horizontal diffusion is approximately half (up to 50 m) of what is achieved with the hybrid experiment. Furthermore, including horizontal diffusion results in an overall reduction of the near-surface global biases in potential temperature and salinity, regardless of the coordinate system. Horizontal diffusion is also necessary to ensure a physically consistent suite of mesoscale eddy parameterizations, particularly in the near-surface region where eddy fluxes are known to be large (e.g., Robbins et al., 2000).
The northward heat transport due to horizontal diffusion can account for 5%-10% of the total heat transport in the Southern Ocean and in the Western boundary current regions of the Northern Hemisphere. This is approximately half of the heat transport due to neutral diffusion in these regions. In the experiments conducted here, the spatially varying mesoscale eddy diffusivity (κ), as defined in Equation 13, does not vary with depth. However, observations and theory suggest that κ is surface intensified and strongly influenced by vertical stratification (Groeskamp et al., 2020). Therefore, the contribution from horizontal diffusion to the total heat transport is likely to change when employing more sophisticated depth-dependent κ schemes.
We have implemented horizontal diffusion within the surface BLD and we defined a transition zone based on BLDs from two neighboring cells. In this zone, the horizontal diffusion fluxes are tapered down while  Figure 2b). We chose linear tapering as the simplest approach and following previous work (Danabasoglu et al., 2008). Previous studies have defined the transition zone using the concept of a layer being intermittently exposed to strong turbulent mixing (Danabasoglu et al., 2008;Ferrari et al., 2008). However, the effect of including such a transition zone in a climate model was negligible (Danabasoglu et al., 2008). We have adopted a different definition than that of Danabasoglu et al. (2008) for practical reasons. Our approach assumes that the total diffusive fluxes (the combined contribution from horizontal and neutral fluxes) vary continuously within the transition zone and this assures that diffusion is applied at all layers.
The algorithm presented here is both conservative and monotonic and, therefore, it is suitable for use in climate studies. The averaged computational cost of the algorithm is ∼4% of the total integration time when employing 3 tracers. This is less than one-quarter of the cost of the tracer advection scheme and also of the fastest neutral diffusion scheme (method 3) proposed by Shao et al. (2020). The construction of the HBD grid accounts for ∼2.5% of the total integration time and is the most expensive part of the algorithm. The HBD grid is constructed at every time step using the layer thicknesses of each pair of adjacent columns (see Section 2.1). This step does not depend on the number of tracers used and the HBD grid is stored in memory to make the algorithm computationally efficient when multiple tracers are employed (e.g., in biogeochemical applications).
This study focuses on applying horizontal eddy diffusion within the surface boundary layer of general vertical coordinate system ocean models using regridding/remapping techniques and tapering. Although the method was developed and applied to represent diffusive fluxes only within the surface boundary layer, a similar approach can be used for other ocean boundaries.