# Cautionary tales from the mesoscale eddy transport tensor

Takaya Uchida, Dhruv Balwada, Quentin Jamet, William K Dewar, Bruno

Deremble, Thierry Penduff, Julien Le Sommer

► To cite this version:

Takaya Uchida, Dhruv Balwada, Quentin Jamet, William K Dewar, Bruno Deremble, et al. Cautionary tales from the mesoscale eddy transport tensor. Ocean Modelling, 2023, 182, pp.102172. 10.1016/j.ocemod.2023.102172. hal-04213098

HAL Id: hal-04213098

https://hal.science/hal-04213098v1

Submitted on 21 Sep 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Cautionary tales from the mesoscale eddy transport tensor

Preprint · September 2022

DOI: 10.31223/X5CW62

CITATIONS

0

READS

121

7 authors, including:

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

Takaya Uchida
Florida State University

39 PUBLICATIONS 269 CITATIONS

SEE PROFILE

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

Dhruv Balwada
New York University

37 PUBLICATIONS 501 CITATIONS

SEE PROFILE

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

Quentin Jamet
Institut national de recherche en informatique et en automatique (Inria)

28 PUBLICATIONS 102 CITATIONS

SEE PROFILE

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

Thierry Penduff
Institute of Environmental Geosciences

169 PUBLICATIONS 4,231 CITATIONS

SEE PROFILE

Some of the authors of this publication are also working on these related projects:

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

"Data Science & Environment", workshop + summer school, 3-7 July 2017, Brest, France View project

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

OCCIPUT View project
This is a revised version of the manuscript submitted for publication in Ocean
Modelling. Please note that, the manuscript is currently under review and has yet to be
formally accepted for publication. Subsequent versions of this manuscript may have
slightly different content. If accepted, the final version of this manuscript will be
available via the ‘Peer-reviewed Publication DOI’ link.

Cautionary tales from the mesoscale eddy
transport tensor

Takaya Uchida, COAPS (tuchida@fsu.edu)

Dhruv Balwada, LDEO (dhruvbalwada@gmail.com)

Quentin Jamet, INRIA (quentin.jamet@inria.fr)

William K. Dewar, FSU (wdewar@fsu.edu)

Bruno Deremble, CNRS (bruno.deremble@univ-grenoble-alpes.fr)

Thierry Penduff, CNRS (thierry.PENDUFF@cnrs.fr)

Julien Le Sommer, CNRS (julien.Lesommer@univ-grenoble-alpes.fr)
Cautionary tales from the mesoscale eddy transport tensor

Takaya Uchida<sup>a,b,*</sup>, Dhruv Balwada<sup>c</sup>, Quentin Jamet<sup>b,d</sup>, William K. Dewar<sup>b,e</sup>, Bruno Deremble<sup>b</sup>, Thierry Penduff<sup>b</sup>, Julien Le Sommer<sup>b</sup>

<sup>a</sup>Center for Ocean-Atmospheric Prediction Studies, Florida State
University, Tallahassee, Florida, USA

<sup>b</sup>Université Grenoble Alpes, CNRS, IRD, Grenoble-INP, Institut des Géosciences de l'Environnement, Grenoble, France

<sup>c</sup>Lamont Doherty Earth Observatory, Columbia University in the City of New York, Palisades, New York, USA

<sup>d</sup>INRIA, ODYSSEY group, Ifremer, Plouzané, France

<sup>e</sup>Department of Earth, Ocean and Atmospheric Science, Florida State
University, Tallahassee, Florida, USA

**Abstract**

The anisotropic mesoscale eddy transport tensor is diagnosed using passive tracers advected online in both an idealized 101-member mesoscale-resolving quasi-geostrophic (QG) double-gyre ensemble, and a realistic 24-member eddying (1/12°) ensemble of the North Atlantic. We assert that the Reynold's decomposition along the ensemble dimension, rather than the spatial or temporal dimension, allows us to capture the intrinsic spatiotemporal variability of the mean flow and eddies. The tensor exhibits good performance in reconstructing the eddy fluxes of passive tracers, here defined as fluctuations about the ensemble thickness-weighted averaged (TWA) mean. However, the inability of the tensor to reconstruct eddy fluxes of QG potential vorticity, which encapsulates the eddy-mean flow interaction, and other active tracers raises the question: To what extent can the diagnosed tensor be applied to inform the parametrization of mesoscale dynamics?

*Keywords:* Ensemble simulation, mesoscale eddy, eddy diffusivity, thickness-weighted average

*Email: tuchida@fsu.edu

Preprint submitted to Ocean Modelling

December 2, 2022
1

# 1. Introduction

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

In the field of computational oceanography, there is never enough computational power to resolve all the spatiotemporal scales of interest, spanning from on the order of $O(1\ \text{m})$ to $O(1000\ \text{km})$ and seconds to centuries respectively. Nevertheless, with the continuous increase in computational power, global fully-coupled ocean simulations with spatial resolution of $O(10\ \text{km})$ have emerged (e.g. Small et al., 2014; Uchida et al., 2017; Chang et al., 2020). This new generation of ocean simulations have improved the representation of oceanic processes, e.g. the oceanic jets in the separated Gulf Stream and Atlantic Circumpolar Current, and the global oceanic heat transport estimates match better with observations (Griffies et al., 2015; Chassignet et al., 2020). Nevertheless, forced ocean simulations with higher model resolution have shown that even at $O(10\ \text{km})$ mesoscale eddies and their feedback onto the jets are not sufficiently resolved (Chassignet and Xu, 2021; Uchida et al., 2019, 2022c; Hewitt et al., 2022).

In order to overcome the insufficient representation of the eddies and their variability, there has been a growing field of research on how to parametrize the effects of eddies in simulations at the grey zone resolution, when the eddy variability is partially resolved. The missing dynamical effects, due to limited resolution, include processes such as energy-backscattering and dissipation of energy and enstrophy (e.g. Bachman et al., 2017; Jansen et al., 2019; Bachman, 2019; Guillaumin and Zanna, 2021; Uchida et al., 2022a, and references therein). In this study, we examine the possibility of capturing the effects of these missing dynamical processes within the framework of eddy transport coefficients, where we represent the sub-grid fluxes of tracers via the gradient flux of the resolved tracer fields. In doing so, we utilized two sets of ensemble simulations: an idealized 101-member quasi-geostrophic (QG) double-gyre ensemble at mesoscale-resolving resolution and a realistic eddying (1/12°) 24-member North Atlantic (NA) ensemble.

In the parametrization literature, it is common to frame the eddies as the sub-grid variability and mean flow as the resolved flow under limited model resolution, where the eddy and mean are defined using a Reynold's decomposition (Bachman et al., 2015). The assumption, and hope, is that the reduced variability in the mean flow, diagnosed by filtering a high-resolution or observed flow field, would mimic the partially resolved variability at coarser model resolution. In a seminal paper, Young (2012) demonstrated that when the equations of oceanic motions are Reynold's decomposed via a thickness-
weighted averaging, the net effect of the eddies onto the mean flow can be represented in terms of eddy Ertel's potential vorticity (PV) fluxes on the thickness-weighted averaged (TWA) mean momentum equations. With this interpretation, the goal of the parametrization problem is to represent the eddy QG or Ertel's PV fluxes in terms of the mean fields (Young, 2012; Marshall et al., 2012; Vallis, 2017). In this study we address the question: Can the eddy fluxes be expressed in terms of a tensor transport model, where the eddy flux is related to the local gradient through a transport tensor? Also, is the transport tensor for active tracers (e.g. PV) similar to a transport tensor for passive tracers? Our approach is to diagnose the transport coefficients using passive tracers outputs and then to examine whether this can be applied to reconstruct eddy fluxes of active tracers. The potential interchangeability between passive (in the sense that they do not affect the dynamics) and active tracers such as PV has its history in that the governing equations for the two take a similar form of an advective-diffusive equation (with the caveat that PV is directly linked to the momentum and buoyancy while passive tracers are not; Killworth, 1997; Wilson and Williams, 2006; Eden and Greatbatch, 2008).

For practical reasons, a time or spatial mean has often been employed to Reynold's decompose the flow into its eddy and mean in diagnostic studies. However, this comes with its own issues regarding interpretation. For a temporally varying system, such as the real ocean, defining the mean flow via a time mean conflates intrinsic variability of the ocean resulting from the non-linearities in governing equations with the variability in the atmospheric forcing (Aiki and Richards, 2008; Fedele et al., 2021; Uchida et al., 2022b). There are also issues surrounding how to choose the filtering scale, and what model resolution, if any, the scale corresponds to (Bachman et al., 2015). An alternative approach is to define the mean flow via an ensemble mean, where the intrinsic variability is expressed as eddies and the oceanic response to atmospheric forcing as the mean (Sérazin et al., 2018; Leroux et al., 2018; Uchida et al., 2022b). Taking the ensemble outputs, we, therefore, define the eddy-mean flow decomposition along the ensemble dimension and attempt to reconstruct the eddy flux of tracers defined as fluctuations about the ensemble TWA mean. Unlike the spatiotemporal filtering approaches, there is less ambiguity about how the eddy and mean are defined (Jamet et al., 2022), but the connection to parametrization in coarse-resolution models still remains tenuous. Nonetheless, in this study we examine the eddy-mean flow problem under this framework, which is relatively novel in the ocean
76

77

78

79

80

turbulence literature. We also note that because the ensemble dimension is orthogonal to the spatiotemporal dimensions, our Reynold's decomposition is exact. Our approach can be rephrased as us having the lofty long-term goal to parametrize the oceanic intrinsic variability otherwise resolved under sufficient model resolution.

81

82

83

84

85

The paper is organized as follows: We describe the ensemble model configurations and framework of the anisotropic eddy transport tensor in Section 2. The results are given in Section 3 with a discussion on the structural similarity in eddy transport coefficients that emerge between the two ensembles. We conclude and provide some cautionary notes in Section 4.

86

## 2. Methods

87

88

89

We describe the quasi-geostrophic (QG) and realistic North Atlantic (NA) ensembles (sections 2.1 and 2.2), and then provide a brief description of the eddy transport tensor (section 2.3).

90

### 2.1. Quasi-geostrophic ensemble

91

92

93

94

95

96

97

98

99

The model outputs used here are those of Uchida et al. (2021b). For completeness, we provide a brief overview of the configuration here. We use the QG configuration of the Multiple Scale Ocean Model (MSOM; Deremble and Martinez, 2020, hereon referred to as MSQG), based on the Basilisk language (Popinet, 2015), to simulate a three-layer double-gyre flow with a rigid lid and flat bottom. The characteristic length scale of the Rossby radius is prescribed as 50 km and horizontal resolution is ~ 4 km, so we have roughly 12 grid points per radius; our simulation can be considered mesoscale resolving (Hallberg, 2013).

100

101

102

103

104

105

106

107

108

109

110

The model was forced with a stationary wind stress curl without any buoyancy forcing at the surface. A seasonally varying background stratification was prescribed at the first interface but kept stationary at the second interface, which is consistent with the seasonal variability of stratification being confined in the upper few hundred meters in the real ocean (Chelton et al., 1998). The 101 ensemble members are initialized with stream functions slightly perturbed at a single grid point randomly selected in the first layer per member and the surface wind stress and temporally varying background stratification are kept identical amongst the members. We refer the interested reader to Uchida et al. (2021b) for further details on the model configuration and ensemble generation.
In addition, we added four passive tracers to each ensemble member with the governing equation:

$$C_{it} + J(\psi, C_i) = -\mathcal{T}^{-1}(C_i - \dot{C}_i), \quad (1)$$

where $C_i, i = 0, 1, 2, 3$ are the four passive tracers and $\psi$ the stream function. They are relaxed towards a profile orthogonal to each other

$$\dot{C}_0 = \frac{x - L/2}{L/2}, \quad \dot{C}_1 = \frac{y - L/2}{L/2}, \quad \dot{C}_2 = \sin \frac{2\pi x}{L}, \quad \dot{C}_3 = \sin \frac{2\pi y}{L}, \quad (2)$$

with a relaxation time scale of $\mathcal{T} \sim 360$ days for all three vertical layers. $L (= 4000$ km) is the zonal and meridional domain extent, $x = [0, L]$, $y = [0, L]$, resulting in the tracers taking values between $[-1, 1]$. The time scale chosen is similar to previous studies (e.g. Bachman et al., 2020), and meant to be longer than that the typical eddy turn over time scale. In other words, the tracers are passively stirred by the flow realized by the intrinsic variability of each ensemble member but are relaxed towards identical profiles amongst members. The relaxation avoids the tracers from homogenizing (i.e. $\nabla_h C_i \neq 0$ where $\nabla_h$ is the horizontal gradient operator). The tracer profiles were initialized with their respective relaxation profiles at the beginning of the ensemble generation and advected prognostically online.

The QG eddy flux is $\mathbf{J}^{C_i} = J^{C_i 1}\mathbf{i} + J^{C_i 2}\mathbf{j} \stackrel{\text{def}}{=} \overline{u' C'_i}\mathbf{i} + \overline{v' C'_i}\mathbf{j}$ where $(\cdot)$ is the ensemble mean operator and $(\cdot)' \stackrel{\text{def}}{=} (\cdot) - (\cdot)$. $\mathbf{i}, \mathbf{j}$ are the horizontal unit vectors in geopotential coordinates. Snapshots of the eddy kinetic energy (EKE) and tracer variance are shown in the Supplementary Material Fig. S1.

## 2.2. Realistic North Atlantic ensemble

We use the model outputs from the realistic simulations described in Jamet et al. (2019, 2020), and Uchida et al. (2022b), which are 48 air-sea partially coupled ensemble members of the NA ocean at mesoscale-permitting resolution (1/12°) using the hydrostatic configuration of the Massachusetts Institute of Technology general circulation model (MITgcm; Marshall et al., 1997). The modelled domain was configured to wrap around zonally in order to reduce memory allocation in running the simulation. Similar to the QG ensemble, in the latter subset of 24 members, four passive tracers per member were added online using the Relaxation Boundary Conditions (RBCS) package with the relaxation profiles of

$$\dot{C}_0 = \frac{y - L^y/2}{L^y/2}; \quad \dot{C}_1 = \sin \frac{2\pi x}{L^x}; \quad \dot{C}_2 = \frac{z - H/2}{H/2}; \quad \dot{C}_3 = \sin \frac{2\pi y}{L^y}, \quad (3)$$
and the relaxation time scale of $\mathcal{T} = 365$ days. $L^x, L^y$ are the zonal and meridional domain extent respectively and $H$ is the deepest depth in the domain bathymetry. The tracer profiles were again initialized with their respective relaxation profiles at the beginning of the ensemble generation and advected prognostically online.

The connection between primitive equations viewed in the thickness-weighted averaged (TWA) framework and quasi geostrophy is that in the latter, the layer thickness $\sigma (\stackrel{\text{def}}{=} \zeta_b)$ only fluctuates on the order of Rossby number ($\zeta$ is the depth of the neutral surface and the subscript ($\widetilde{\cdot}$) denotes derivatives in density coordinates). Hence, it can be argued that QG variables are implicitly TWA and that for a fair comparison, the primitive equations should also be TWA (cf. Marshall et al., 2012). The three-dimensional oceanic motions become quasi two dimensional upon thickness-weighted averaging (Aoki, 2014; Loose et al., 2022), further elucidating the similarity to quasi geostrophy. We proceed in defining our eddy tracer fluxes within the TWA framework, which gives $\mathbf{J}^{C_i} = J^{C_i 1} \bar{\mathbf{e}}_1 + J^{C_i 2} \bar{\mathbf{e}}_2 \stackrel{\text{def}}{=} \widetilde{u''C_i'' \bar{\mathbf{e}}_1} + \widetilde{v''C_i'' \bar{\mathbf{e}}_2}$ where ($\widetilde{\cdot} \stackrel{\text{def}}{=} \bar{\sigma}^{-1}\bar{\sigma}(\widetilde{\cdot})$ and ($\widetilde{\cdot} \stackrel{\text{def}}{=} \widetilde{\cdot} - \widetilde{\cdot}$) (Young, 2012). $\bar{\mathbf{e}}_1$ and $\bar{\mathbf{e}}_2$ are the horizontal unit vectors along the neutral surface, which reduce to $\mathbf{i}$ and $\mathbf{j}$ respectively when the neutral surfaces are flat ($\widetilde{\nabla}_h \zeta = 0$ where $\widetilde{\nabla}_h$ is the horizontal gradient operator in density coordinates). We note that due to the 3rd order direct-space-time (DST) flux-limiter advective scheme used, there are two possible ways to define the eddy flux, i.e. $\widetilde{\mathbf{u}}''C_i'' \approx \bar{\sigma}^{-1}(\bar{\sigma}\mathbf{F}^{C_i}) - \hat{\mathbf{u}}\hat{C_i}$ where $\mathbf{F}^{C_i}$ is the DST advective flux computed by MITgcm. The advective scheme was chosen to be consistent with the schemes used for temperature and salinity. Snapshots of the EKE and tracer variance are given in the Supplementary Material Fig. S2.

In transforming the three-dimensional model outputs in geopotential coordinates to TWA variables, we: i) compute the approximately neutral density surfaces for each member by removing the effect of compressibility, ii) remap the momentum, isopycnal depth and tracers onto the neutral surfaces, iii) thickness weight the variables and take the ensemble mean, iv) apply Reynold's decomposition about the ensemble TWA mean to obtain the eddy terms ($\cdot \cdot''$, and $\cdot$) thickness-weighted average the eddy covariance to get the eddy flux. Further details regarding the coordinate remapping from geopotential to approximately neutral surfaces using the `xgcm` Python package (Abernathey et al., 2021), and the averaging are given in Uchida et al. (2022b).
## 2.3. Eddy transport tensor

As the eddy flux of tracers is generally poorly resolved in global forced and coupled ocean simulations, there has been an effort to parametrize this flux in the quasi-adiabatic interior via a local-gradient based model along neutral surfaces (Redi, 1982; Griffies, 2004; Wilson and Williams, 2006; Holmes et al., 2022)

$$ \mathbf{J}^{C_i} = -\kappa \nabla_h \overline{C_i}, \qquad (4) $$

where $\kappa$ is the scalar eddy diffusivity, hereon referred to as the transport coefficient. The terminology of ‘diffusivity’ has its origin in the works by Redi (1982); Gent and McWilliams (1990) and Gent et al. (1995), but eddies not only diffuse but advect tracers (Smith and Gent, 2004; Griffies, 2004; Haigh et al., 2021b), so we have opted for a more general term.

While it is tempting to directly infer a scalar eddy transport coefficient from equation (4), assuming an isotropic transport coefficient for an anisotropic flow in realistic simulations is a poor approximation (Smith and Gent, 2004; Ferrari and Nikurashin, 2010; Rypina et al., 2012; Fox-Kemper et al., 2013; Kamenkovich et al., 2015, 2020). A more general or appropriate model might be to relate the eddy fluxes to the mean gradients via a transport tensor with four parameters (Plumb and Mahlman, 1987), which has some justification in linear wave theories and mixing-length based models (Bachman and Fox-Kemper, 2013). In both these models, the inherent assumption is that the scalar transport coefficient or transport tensor is only a function of the flow, and is tracer independent.

We, therefore, take the approach of estimating the eddy transport tensor ($\mathbf{K}$) from a least-squares best fit to (Plumb and Mahlman, 1987; Abernathey et al., 2013; Bachman and Fox-Kemper, 2013)

$$ \underbrace{\begin{pmatrix} J^{C01} & J^{C02} \\ J^{C11} & J^{C12} \\ J^{C21} & J^{C22} \\ J^{C31} & J^{C32} \end{pmatrix}}_{J} = \mathbf{G} \cdot \underbrace{\begin{pmatrix} \kappa^{uu} & \kappa^{vu} \\ \kappa^{uv} & \kappa^{vv} \end{pmatrix}}_{\mathbf{K}}, \qquad (5) $$

where for the QG ensemble

$$ \mathbf{G} = - \begin{pmatrix} \overline{C_{0x}} & \overline{C_{0y}} \\ \overline{C_{1x}} & \overline{C_{1y}} \\ \overline{C_{2x}} & \overline{C_{2y}} \\ \overline{C_{3x}} & \overline{C_{3y}} \end{pmatrix}, \qquad (6) $$
and for the realistic NA ensemble

$$ \mathbf{G} = - \begin{pmatrix} \widehat{C}_{0\tilde{x}} & \widehat{C}_{0\tilde{y}} \\ \widehat{C}_{1\tilde{x}} & \widehat{C}_{1\tilde{y}} \\ \widehat{C}_{2\tilde{x}} & \widehat{C}_{2\tilde{y}} \\ \widehat{C}_{3\tilde{x}} & \widehat{C}_{3\tilde{y}} \end{pmatrix}, \qquad (7) $$

respectively. The least-squares fit can be estimated as $K = G^+J$ where $G^+$ is the Moore-Penrose pseudo inverse of $G$ for each data point (Bachman et al., 2015). The rotational components in $J^{C_i}$ are not removed prior to inversion. It is possible to invert equation (5) with just two tracers whose gradients are not aligned with each other, and this will return a unique solution specific to the pair of tracers chosen for the inversion (e.g. Haigh and Berloff, 2021; Sun et al., 2021). Here, we instead focus on estimating a single tensor that works for all possible tracer orientations, even if it does not perfectly reconstruct the eddy flux for any single tracer. We have, thus, kept the system over-determined by using four tracers. The assumption is that by keeping it over determined, $K$ would extract the universal component in the relation between eddy fluxes and gradient flux of the mean tracer fields. From a practical point of view, all tracers should be associated with the same transport coefficients in order to reduce the number of model parameters.

The eddy flux and mean tracer gradient fields for the QG ensemble were coarse grained by 4 × 4 grid-point boxcar filter prior to inverting (5) in order to reduce the computational cost. For the NA ensemble, we first spatially smoothed the eddy flux and mean tracer fields by applying a Gaussian kernel with the standard deviation of 50 km using the gcm-filters Python package (Grooms et al., 2021) as we do not expect the linear model (5) to capture grid-scale features. The eddy flux and mean tracer gradient fields were then further coarse grained by 10 × 10 grid-point boxcar filter in order to reduce the computational cost of inversion. (We demonstrate that other than the significant reduction in computation, the coarse graining has no qualitative impact on our results; Fig. S3. We also note that a boxcar coarse graining commutes with the spatial derivatives and does not result in extra Clark terms. cf. Bachman et al., 2015) Each row in $J$ and $G$ was then normalized by horizontal median of the magnitude of each mean tracer gradient flux (i.e. $\frac{J^{C_i}G^{C_i}}{\text{median}[|G^{C_i}|]}$ where $J^{C_i}$ and $G^{C_i}$ are the smoothed and coarse-grained eddy flux and mean gradient flux of an arbitrary tracer $C_i$) prior to the inversion so that each tracer had roughly equal weighting in inverting equation (5).
## 3. Results

We present results using the first time step of the fifth year, a time when the ensemble spread has converged for both ensembles (Uchida et al., 2021b; Jamet et al., 2019). The QG outputs are instantaneous snapshots while the North Atlantic (NA) outputs are five-day averaged.

### 3.1. Quasi-geostrophic ensemble

The four components in the eddy transport tensor are provided in the top two rows of Fig. 1 with values reaching up to $O(10^4 \text{ m}^2 \text{ s}^{-1})$ in the first and second layer. The bottom layer is quiescent with transport coefficients on the order of $O(10^2 \text{ m}^2 \text{ s}^{-1})$. The diagonal components of the tensor tend to take positive values over the entire domain while the antidiagonal components tend to change signs across the jet centered around $y = 2000 \text{ km}$. This suggests that the eddies are broadly working to dissipate small-scale variance, as one might expect. Focusing on the first layer, $\kappa^{uv}$ tends to be coherently negative south of the jet and positive north of the jet. In the lower two layers, $\kappa^{uv}$ and $\kappa^{vu}$ tend to mirror each other where there is a sign change across the jet but also within each idealized subtropical and subpolar gyre. These cross diagonal terms of the tensor are largely associated with eddy-induced advection, and often act to oppose the mean flow (Marshall, 2011).

Taking the diagnosed tensor, we examine the reconstruction of the eddy flux of $C_2$, i.e. $J_{\text{reconstructed}}^{C_2} \stackrel{\text{def}}{=} -\nabla_h \overline{C_2} \cdot \mathbf{K}$. The middle two rows of Fig. 1 show that the performance of reconstruction is good across all three layers. The spatial correlations between the actual eddy flux and its reconstruction are higher than 0.9 for all four tracers (Table 1). The spatial correlation was computed as

$$r_{ij}^{C} \stackrel{\text{def}}{=} \frac{\sum \left[ (J_{ij}^{C} - \langle J_{ij}^{C} \rangle) \left( J_{\text{reconstructed}}^{C} - \langle J_{\text{reconstructed}}^{C} \rangle \right) \right]}{\sqrt{\sum (J_{ij}^{C} - \langle J_{ij}^{C} \rangle)^2} \sqrt{\sum \left( J_{\text{reconstructed}}^{C} - \langle J_{\text{reconstructed}}^{C} \rangle \right)^2}}, \quad (8)$$

where $\langle \cdot \rangle$ is the horizontal spatial mean, $j = 1, 2$ correspond to the zonal and meridional orientation, and the summation is taken over the entire domain of interest.

We also exhibit the spatial median of normed relative errors $|\varepsilon_{ij}^C|$ computed as

$$\varepsilon_{ij}^{C} \stackrel{\text{def}}{=} \frac{J_{ij}^{C} - J_{\text{reconstructed}}^{C}}{J_{ij}^{C}}, \quad (9)$$
<table><caption>Table 1: Spatial correlation for the four passive tracers and QGPV between the eddy flux and its reconstruction for the zonal and meridional orientation and all three layers. The PV coefficients are not shown for the middle and bottom layer as it is evident from Fig. 1 that the flux and reconstruction are completely decorrelated.</caption><thead><tr><th>Tracer</th><th>Top</th><th>Middle</th><th>Bottom</th></tr></thead><tbody><tr><td>r<sup>C01</sup></td><td>0.952</td><td>0.989</td><td>0.999</td></tr><tr><td>r<sup>C02</sup></td><td>0.957</td><td>0.990</td><td>0.996</td></tr><tr><td>r<sup>C11</sup></td><td>0.980</td><td>0.966</td><td>0.981</td></tr><tr><td>r<sup>C12</sup></td><td>0.999</td><td>0.999</td><td>0.997</td></tr><tr><td>r<sup>C21</sup></td><td>0.939</td><td>0.978</td><td>0.993</td></tr><tr><td>r<sup>C22</sup></td><td>0.949</td><td>0.966</td><td>0.973</td></tr><tr><td>r<sup>C31</sup></td><td>0.992</td><td>0.969</td><td>0.980</td></tr><tr><td>r<sup>C32</sup></td><td>0.999</td><td>0.999</td><td>0.997</td></tr><tr><td>r<sup>q1</sup></td><td>0.844</td><td>—</td><td>—</td></tr><tr><td>r<sup>q2</sup></td><td>0.906</td><td>—</td><td>—</td></tr></tbody></table>

in Table 2. This allows us to quantify the normalized error in addition to the sign; $\epsilon^{Cij}$ would return positive values if the reconstruction had the correct sign as the actual flux even if it had an offset. When the reconstruction has the wrong sign, the errors would become negative. The spatial median of the normed errors are all positive with values smaller than unity indicating that the reconstruction overall predicts the correct sign of the eddy flux. The spatial fields given in the Supplementary Material Fig. S4 also support the good performance.

As the inversion was done as a least squares fit to (5), the good reconstruction of $J^{C_i}$ may not come as a surprise and is consistent with previous studies (e.g. Abernathey et al., 2013; Bachman et al., 2015, 2020). However, it is important to note that such a good fit, although not perfect (Figs. S3 and S4), suggests that the transport tensor model with four free parameters is a good model to represent how the eddies flux different tracers along neutral surfaces. This model is able to separate the contribution of the flow from any dependence on the orientation of tracer gradients well, and this model of reduced complexity can be a target for parametrizations that are able to represent the eddy fluxes of any arbitrary passive tracer.

One may now ask the question: Can the tensor diagnosed from passive tracers be applied to active tracers, here chosen as QG potential vorticity (PV), which is the sole active tracer in quasi geostrophy? Looking at the
<table><caption>Table 2: Spatial median of the normed relative error for the four passive tracers and QGPV between the eddy flux and its reconstruction for the zonal and meridional orientation and all three layers. The PV errors are not shown for the middle and bottom layer for the same reasons as correlation.</caption><thead><tr><th>Tracer</th><th>Top</th><th>Middle</th><th>Bottom</th></tr></thead><tbody><tr><td>|ε<sup>C01</sup>|</td><td>0.0440</td><td>0.0190</td><td>0.0122</td></tr><tr><td>|ε<sup>C02</sup>|</td><td>0.0622</td><td>0.0517</td><td>0.0169</td></tr><tr><td>|ε<sup>C11</sup>|</td><td>0.141</td><td>0.0746</td><td>0.0252</td></tr><tr><td>|ε<sup>C12</sup>|</td><td>0.0596</td><td>0.0153</td><td>0.00721</td></tr><tr><td>|ε<sup>C21</sup>|</td><td>0.0770</td><td>0.0326</td><td>0.0189</td></tr><tr><td>|ε<sup>C22</sup>|</td><td>0.112</td><td>0.0838</td><td>0.0280</td></tr><tr><td>|ε<sup>C31</sup>|</td><td>0.135</td><td>0.0994</td><td>0.0373</td></tr><tr><td>|ε<sup>C32</sup>|</td><td>0.0596</td><td>0.0210</td><td>0.0112</td></tr><tr><td>|ε<sup>q1</sup>|</td><td>1.161</td><td>—</td><td>—</td></tr><tr><td>|ε<sup>q2</sup>|</td><td>0.945</td><td>—</td><td>—</td></tr></tbody></table>

bottom two rows of Fig. 1, we see that while the reconstruction $-\nabla_h \bar{q} \cdot K$ captures some features in the top layer (cf. Table 1) and in the regions away from the jet in the middle layer, the reconstruction performs poorly for the bottom layers and in the jet in the middle layer, with the spatial correlations being smaller than 0.1. The normed relative errors are also quite large being on the order of $O(1)$ (Table 4; Fig. S5).

## 3.2. Realistic North Atlantic ensemble

Following the QG analyses, we now present the four components of the tensor $K$ from the NA ensemble. The magnitude of the tensor components are similar to the QG ensemble, reaching up to $O(10^4 \text{ m}^2 \text{ s}^{-1})$ (top left four panels in Fig. 2). Also similar to the QG ensemble, the diagonal components tend to take positive values while the antidiagonal components tend to change signs across the separated Gulf Stream (GS). The separated GS can be identified in the top right panel in Fig. 2 where the ensemble-mean depth associated with the neutral surface shoals around 35°N. We also show the vertical transects along 300°E. The patterns in sign persist over depth with the order of magnitude decreasing towards $O(10^2 \text{ m}^2 \text{ s}^{-1})$ in the abyssal ocean (bottom four panels in Fig. 2).

We remind the reader that we define the eddy flux as $J^{C_i} \stackrel{\text{def}}{=} \tilde{u}''C''_i$, and following the convention of the TWA framework, the reconstructed eddy flux
<table><caption>Table 3: Spatial correlation for the four passive tracers and temperature and salinity between the eddy flux and its reconstruction for the zonal and meridional orientation along the neutral surface shown in Fig. 2. The correlation coefficients are shown for when they are diagnosed over the entire horizontal domain and only between 10°N-30°N both excluding the hatched regions in Figs. 3 and 4.</caption><thead><tr><th>Tracer</th><th>Entire domain</th><th>10°N-30°N</th></tr></thead><tbody><tr><td>r<sup>C01</sup></td><td>0.918</td><td>0.956</td></tr><tr><td>r<sup>C02</sup></td><td>0.917</td><td>0.971</td></tr><tr><td>r<sup>C11</sup></td><td>0.935</td><td>0.959</td></tr><tr><td>r<sup>C12</sup></td><td>0.917</td><td>0.936</td></tr><tr><td>r<sup>C21</sup></td><td>0.975</td><td>0.948</td></tr><tr><td>r<sup>C22</sup></td><td>0.983</td><td>0.969</td></tr><tr><td>r<sup>C31</sup></td><td>0.924</td><td>0.967</td></tr><tr><td>r<sup>C32</sup></td><td>0.931</td><td>0.977</td></tr><tr><td>r<sup>Θ1</sup></td><td>0.344</td><td>0.303</td></tr><tr><td>r<sup>Θ2</sup></td><td>0.494</td><td>0.353</td></tr><tr><td>r<sup>S1</sup></td><td>0.181</td><td>0.390</td></tr><tr><td>r<sup>S2</sup></td><td>0.308</td><td>0.128</td></tr></tbody></table>

of $C_1$ becomes $J_{\text{reconstructed}}^{C_1} \stackrel{\text{def}}{=} -\tilde{\nabla}_h \widehat{C_1} \cdot \mathbf{K}$. In examining the reconstruction, we limit it to regions where the neutral surface is deeper than 150 m and to grids upon the $10 \times 10$ coarse graining included no land cells in order to minimize the effect of diabatic mixing (here parametrized by the K-Profile Parametrization; Large et al., 1994). As we see from Fig. 3, the reconstruction of eddy passive tracer fluxes is good generally across the entire three dimensional domain with the spatial correlation again higher than 0.9 (Table 3). The largest disagreements between $J_{\text{reconstructed}}^{C_1}$ and $J_{\text{reconstructed}}^{C_2}$ can be seen in the separated GS region where eddy activity and vertical fluctuations of the neutral surface are vigorous (Uchida et al., 2022b). The horizontal spatial correlation improves for the quiescent gyre interior slightly for all four passive tracers.

Similar to the QG ensemble, we also show the normalized relative error from the NA ensemble in Table 4. The magnitude of errors are smaller than unity over most of the domain (Fig. S6 in Supplementary Material). As expected from the improvement seen in the spatial correlations, the normalized relative errors for passive tracers generally decrease when they are taken over the subdomain of 10°N-30°N.
<table><caption>Table 4: Spatial median of the normed relative errors for the four passive tracers and temperature and salinity between the eddy flux and its reconstruction for the zonal and meridional orientation along the neutral surface shown in Fig. 2. The errors are shown for when they are diagnosed over the entire horizontal domain and only between 10°N-30°N both excluding the hatched regions in Figs. 3 and 4.</caption><thead><tr><th>Tracer</th><th>Entire domain</th><th>10°N-30°N</th></tr></thead><tbody><tr><td>|ε<sup>C01</sup>|</td><td>0.0257</td><td>0.0141</td></tr><tr><td>|ε<sup>C02</sup>|</td><td>0.0208</td><td>0.00593</td></tr><tr><td>|ε<sup>C11</sup>|</td><td>0.0118</td><td>0.0257</td></tr><tr><td>|ε<sup>C12</sup>|</td><td>0.0108</td><td>0.00847</td></tr><tr><td>|ε<sup>C21</sup>|</td><td>0.0372</td><td>0.0324</td></tr><tr><td>|ε<sup>C22</sup>|</td><td>0.0726</td><td>0.0582</td></tr><tr><td>|ε<sup>C31</sup>|</td><td>0.0334</td><td>0.00343</td></tr><tr><td>|ε<sup>C32</sup>|</td><td>0.0261</td><td>0.00153</td></tr><tr><td>|ε<sup>Θ1</sup>|</td><td>0.343</td><td>0.286</td></tr><tr><td>|ε<sup>Θ2</sup>|</td><td>0.428</td><td>0.331</td></tr><tr><td>|ε<sup>S1</sup>|</td><td>0.234</td><td>0.183</td></tr><tr><td>|ε<sup>S2</sup>|</td><td>0.314</td><td>0.222</td></tr></tbody></table>

We may now also examine if we can use the tensor to reconstruct the eddy temperature and salinity fluxes, variables which were not included in the inversion of (5) and are active tracers. Overall the level of reconstruction is poorer than that of passive tracers particularly in the separated GS region; the discrepancy north of 30°N extends vertically down to ~1000 m (bottom two rows in Fig. 4). It is interesting, however, that there seems to be some utility of the tensor in the quiescent gyre interior (top three rows of Fig. 4; Table 3); $J^\Theta$ reconstructed and $J^S$ reconstructed capture the sign structure south of 30°N in $J^\Theta$ and $J^S$ respectively. The normed errors are also on the order of $O(0.1)$ indicating some level of reconstructing the eddy flux particularly in the quiescent eastern part of the gyre (cf. Fig. S7 in Supplementary Material). Given the QGPV results from the previous section, we did not attempt the reconstruction for Ertel’s PV.

## 4. Conclusion and discussion

In this study, we have diagnosed the anisotropic mesoscale eddy transport tensor **K** using passive tracer outputs from two sets of ensemble simulations: an idealized 101-member quasi-geostrophic (QG) double-gyre ensemble and
a realistic 24-member ensemble of the North Atlantic (NA). In decomposing the eddies and mean flow, we have chosen to take the averaging operator over the ensemble dimension rather than the often employed spatial or temporal dimension (e.g. Balwada et al., 2020; Bachman et al., 2020; Kamenkovich et al., 2020; Haigh and Berloff, 2021) with the aim of capturing the intrinsic variability of eddy transport. We have also investigated the tensor in the thickness-weighted averaged (TWA) context, which to our knowledge, is the first study to do so. While we have only utilized one time slice of output, we do not expect the performance of inverting for $K$ to qualitatively vary over time in the quasi-adiabatic interior of the ocean; in Fig. S8 in the Supplementary Material, we exhibit the time series of the spatial median of the major and minor eigenvalues of the symmetric part of the QG tensor ($S = (K + K^T)/2$ where $K^T$ is the transpose). The physical significance of the eigenvalues is that the eddies tend to diffuse the tracers in the direction of the eigenvectors with an associated transport coefficient of each respective eigenvalue (Abernathey et al., 2013; Haigh et al., 2021a). We find a robust seasonal cycle repeating itself (larger values during late summer and smaller values during late winter) with a lag to the seasonality in the prescribed background stratification. When examining the climatological summer average of transport coefficients diagnosed at each time step and averaged over five years, they exhibit less spatial structure but with the same tendencies as Fig. 1, i.e. the diagonal components are largely positive and antidiagonal components change sign across the separated jet (Fig. S9).

The diagnosed tensor shows good performance in reconstructing the eddy fluxes of passive tracers from the gradient flux of mean passive tracer fields, which were weakly restored with a one-year relaxation timescale. The agreement between the spatial patterns in the transport coefficients emerging from the QG and NA ensemble in respect to the position of the jet is also comforting; they have similar orders of magnitude and the diagonal components ($\kappa_{uu}, \kappa_{vv}$) tend to be positive while ($\kappa_{uv}, \kappa_{vu}$) tend to have opposite signs across the jet and change signs within each gyre. This partially justifies our assumption that there is a universal eddy transport tensor, which is able to represent the eddy flux across passive tracers (at least for passive tracers with similar source/sink terms; cf. Fig. S10). As noted in Section 2.2, there are two possible ways to define the eddy tracer fluxes as soon as the advective scheme becomes more complicated than a 2nd-order centered scheme. We also diagnosed the tensor with the eddy fluxes defined as
$$ \mathbf{J}^{C_i} \stackrel{\text{def}}{=} \bar{\sigma}^{-1} (\overline{\sigma \mathbf{F}^{C_i}}) - \hat{\mathbf{u}} \hat{C}_i $$ but the performance of reconstruction deteriorated compared to when $\mathbf{J}^{C_i} \stackrel{\text{def}}{=} \mathbf{u}'' \overline{C_i''}$ (not shown). We speculate that the linear gradient flux model (5) is unable to capture the non-linearities in the flux-limiter advective scheme used in our simulations for all passive and active tracers.

For studies using ensemble simulations, whether the ensemble size is large enough to achieve statistical convergence is always in question. We exhibit the QG transport coefficients in the top layer diagnosed from a subset of 12, 25 and 50 members in Fig. S11. While noisier, the transport coefficients diagnosed from 25 members already demonstrate similar values and spatial patterns to the top left panels in Fig. 1. The histograms also align more closely to the distribution from 101 members and the tails of histogram shorten as the size of subset increases (bottom four panels in Fig. S11).

Given the diagnosed transport tensor from passive tracers, we have further examined whether it can be carried over to inform the parametrization of active tracer fluxes. However, when applying the tensor to reconstruct the eddy fluxes of active tracers, here QG potential vorticity (PV), our results suggest that passive and active tracers have significantly different eddy transport coefficients. In other words, passive and active tracers have different relations between the eddy fluxes and mean fields and this likely stems from the fact that for QGPV, there is no restoring force, similar to what we prescribed for the passive tracers, particularly below the first layer; the gradient of mean QGPV had very little spatial structure in the separated jet region as opposed to the passive tracers (Figs. S12 and S13). Without any restoring for passive tracers, their concentrations would completely homogenize over time ($\nabla_h C_i = 0$) rendering the gradient flux model (5) useless. Our emphasis on the eddy PV flux is because it encapsulates the energy backscattering onto the mean flow (Young, 2012; Marshall et al., 2012; Uchida et al., 2021a). Unlike previous studies (e.g. Mak et al., 2016; Haigh et al., 2021a), we have not split the eddy flux of passive tracers into its divergent and rotational component in diagnosing $\mathbf{K}$ as our interest has been in whether this tensor can be used to approximate the total eddy PV flux itself ($\mathbf{F}^\# \in \text{Young, 2012}$), and not its divergence, that forces the mean momentum. This being said, our least-square approach does fairly well at reconstructing also the eddy flux divergence of passive tracers in the quiescent gyre interior where the relative errors are small (Figs. S13 and S14). Machine learning methods or further generalizations to the framework may provide a pathway forward in finding
the relation for active tracers (e.g. Zanna and Bolton, 2020; Guillaumin and Zanna, 2021; Frezat et al., 2021; Lu et al., 2022).

Nonetheless, the development of a prognostic and physically consistent **K** will likely benefit biogeochemical modelling since biogeochemical tracers are passive (cf. Jones and Abernathey, 2019; Uchida et al., 2020). We end on the note that for the realistic NA ensemble, the tensor shows some skill in reconstructing the eddy temperature and salinity fluxes in the gyre interior. While temperature and salinity are not purely passive as they affect the dynamics via the hydrostatic pressure, the partial reconstruction may suggest that they have less of a direct role on the dynamics compared to PV.

### Data availability statement

The Jupyter notebooks used for the analyses of the QG and NA ensembles are available via Github (https://github.com/roxyboy/qg_eddy_diffusivity and https://github.com/roxyboy/TWA-eddy-diffusivity respectively; DOIs will be added upon acceptance of the manuscript).

### Declaration of competing interest

The authors declare no conflict of interest.

### Acknowledgements

We thank the editor Dr. Fangli Qiao along with three anonymous reviewers for their comments. This study is a contribution to the 'Assessing the Role of forced and internal Variability for the Ocean and climate Response in a changing climate' (ARVOR) project supported by the French 'Les Enveloppes Fluides et l'Environnement' (LEFE) program. This work was supported via National Science Foundation (NSF) grants OCE-1829856, OCE-2023585, OCE-2123632 and the French 'Make Our Planet Great Again' (MOPGA) program managed by the Agence Nationale de la Recherche under the Programme d'Investissement d'Avenir, reference ANR-18-MPGA-0002. The latter two grants served as the primary support for Uchida and partially for Jamet. Balwada acknowledges support from the NSF grant OCE-1756882. We acknowledge high-performance computing resources for generating the QG ensemble on Occigen, which was maintained by CINES with the reference A0090112020. The open-source software for MSOM can
be found at github.com/bderembl/msom. It was developed as a module of Basilisk (available at basilisk.fr). High-performance computing resources on Cheyenne (doi:10.5065/D6RX99HX) used for running the NA ensembles were provided by NCAR's Computational and Information Systems Laboratory, sponsored by NSF, under the university large allocation UFSU0011. Uchida personally thanks the Oceanography with an Alpine Flair workshop held in Grenoble 2022, which has led to fruitful exchanges. The model diagnostics were executed on the Florida State University and Université Grenoble Alpes clusters.

References

Abernathey, R.P., Busecke, J., Smith, T., Banihirwe, A., Fernandes, F., Bourbeau, J., Cherian, D., Dussin, R., Swanson-Hysell, N., Constantinou, N., Ponte, A., et al., 2021. xgcm: General circulation model postprocessing with xarray. URL: https://xgcm.readthedocs.io/en/latest/, doi:10.5281/zenodo.3634752.

Abernathey, R.P., Ferreira, D., Klocker, A., 2013. Diagnostics of isopycnal mixing in a circumpolar channel. Ocean Modelling 72, 1–16. doi:10.1016/j.ocemod.2013.07.004.

Aiki, H., Richards, K.J., 2008. Energetics of the global ocean: The role of layer-thickness form drag. Journal of Physical Oceanography 38, 1845–1869. doi:10.1175/2008JP03820.1.

Aoki, K., 2014. A constraint on the thickness-weighted average equation of motion deduced from energetics. Journal of Marine Research 72, 355–382. doi:10.1357/002224014815469886.

Bachman, S.D., 2019. The GM+E closure: A framework for coupling backscatter with the Gent and McWilliams parameterization. Ocean Modelling 136, 85–106. doi:10.1016/j.ocemod.2019.02.006.

Bachman, S.D., Fox-Kemper, B., 2013. Eddy parameterization challenge suite I: Eady spindown. Ocean Modelling 64, 12–28. doi:10.1016/j.ocemod.2012.12.003.

Bachman, S.D., Fox-Kemper, B., Bryan, F.O., 2015. A tracer-based inversion method for diagnosing eddy-induced diffusivity and advection. Ocean Modelling 86, 1–14. doi:10.1016/j.ocemod.2014.11.006.
490

Bachman, S.D., Fox-Kemper, B., Bryan, F.O., 2020. A diagnosis of anisotropic eddy diffusion from a high-resolution global ocean model. Journal of Advances in Modeling Earth Systems 12, e2019MS001904. doi:10.1029/2019MS001904.

491

492

493

494

Bachman, S.D., Fox-Kemper, B., Pearson, B., 2017. A scale-aware subgrid model for quasi-geostrophic turbulence. Journal of Geophysical Research: Oceans 122, 1529-1554. doi:10.1002/2016JC012265.

495

496

497

Balwada, D., Smith, K.S., Abernathey, R.P., 2020. Measuring eddy driven transport in a zonally inhomogeneous flow, in: Ocean Sciences Meeting 2020, AGU.

498

499

500

Chang, P., Zhang, S., Danabasoglu, G., Yeager, S.G., Fu, H., Wang, H., Castruccio, F.S., Chen, Y., Edwards, J., Fu, D., et al., 2020. An unprecedented set of high-resolution earth system simulations for understanding multi-scale interactions in climate variability and change. Journal of Advances in Modeling Earth Systems 12, e2020MS002298. doi:10.1029/2020MS002298.

501

502

503

504

505

Chassignet, E.P., Xu, X., 2021. On the importance of high-resolution in large-scale ocean models. Advances in Atmospheric Sciences, 1-14doi:10.1007/s00376-021-0385-7.

506

507

508

Chassignet, E.P., Yeager, S.G., Fox-Kemper, B., Bozec, A., Castruccio, F., Danabasoglu, G., Horvat, C., Kim, W.M., Koldunov, N., Li, Y., et al., 2020. Impact of horizontal resolution on global ocean-sea ice model simulations based on the experimental protocols of the Ocean Model Intercomparison Project phase 2 (OMIP-2). Geoscientific Model Development 13, 4595-4637. doi:10.5194/gmd-13-4595-2020.

509

510

511

512

513

514

Chelton, D.B., DeSroeke, R.A., Schlax, M.G., El Naggar, K., Siwertz, N., 1998. Geographical variability of the first baroclinic Rossby radius of deformation. Journal of Physical Oceanography 28, 433-460. doi:10.1175/1520-0485(1998)028<0433:GVOTFB>2.0.CO;2.

515

516

517

518

Deremble, B., Martinez, E.M., 2020. MSOM: Multiple Scale Ocean Model. URL: https://github.com/bderembl/msom, doi:10.5281/zenodo.4669909.

519

520

521

Eden, C., Greatbatch, R.J., 2008. Towards a mesoscale eddy closure. Ocean Modelling 20, 223-239. doi:10.1016/j.ocemod.2007.09.002.

522
523

Fedele, G., Penduff, T., Pierini, S., Alvarez-Castro, M.C., Bellucci, A., Masina, S., 2021. Interannual to decadal variability of the Kuroshio extension: analyzing an ensemble of global hindcasts from a dynamical system viewpoint. Climate Dynamics, 1–18doi:10.1007/s00382-021-05751-7.

524

525

526

527

Ferrari, R., Nikurashin, M., 2010. Suppression of eddy diffusivity across jets in the Southern Ocean. Journal of Physical Oceanography 40, 1501–1519. doi:10.1175/2010JP04278.1.

528

529

530

Fox-Kemper, B., Lumpkin, R., Bryan, F.O., 2013. Lateral transport in the ocean interior, in: International Geophysics. Elsevier. volume 103, pp. 185–209. doi:10.1016/B978-0-12-391851-2.00008-8.

531

532

533

Frezat, H., Balarac, G., Le Sommer, J., Fablet, R., Lguensat, R., 2021. Physical invariance in neural networks for subgrid-scale scalar flux modeling. Physical Review Fluids 6, 024607. doi:10.1103/PhysRevFluids.6.024607.

534

535

536

537

Gent, P.R., McWilliams, J.C., 1990. Isopycnal mixing in ocean circulation models. Journal of Physical Oceanography 20, 150–155. doi:10.1175/1520-0485(1990)020<0150:IMIOCM>2.0.CO;2.

538

539

540

Gent, P.R., Willebrand, J., McDougall, T.J., McWilliams, J.C., 1995. Parameterizing eddy-induced tracer transports in ocean circulation models. Journal of Physical Oceanography 25, 463–474. doi:10.1175/1520-0485(1995)025<0463:PEITTI>2.0.CO;2.

541

542

543

544

Griffies, S.M., 2004. Fundamentals of ocean climate models. Princeton university press.

545

546

Griffies, S.M., Winton, M., Anderson, W.G., Benson, R., Delworth, T.L., Dufour, C.O., Dunne, J.P., Goddard, P., Morrison, A.K., Rosati, A., et al., 2015. Impacts on ocean heat from transient mesoscale eddies in a hierarchy of climate models. Journal of Climate 28, 952–977. doi:10.1175/JCLI-D-14-00353.1.

547

548

549

550

551

Grooms, I., Loose, N., Abernathey, R., Steinberg, J., Bachman, S.D., Marques, G., Guillaumin, A.P., Yankovsky, E., 2021. Diffusion-based smoothers for spatial filtering of gridded geophysical data. Journal of Advances in Modeling Earth Systems, e2021MS002552URL:

552

553

554
https://gcm-filters.readthedocs.io/en/latest/, doi:10.1029/2021MS002552.

Guillaumin, A.P., Zanna, L., 2021. Stochastic-deep learning parameterization of ocean momentum forcing. Journal of Advances in Modeling Earth Systems 13, e2021MS002534. doi:10.1029/2021MS002534.

Haigh, M., Berloff, P., 2021. On co-existing diffusive and anti-diffusive tracer transport by oceanic mesoscale eddies. Ocean Modelling 168, 101909. doi:10.1016/j.ocemod.2021.101909.

Haigh, M., Sun, L., McWilliams, J.C., Berloff, P., 2021a. On eddy transport in the ocean. Part I: The diffusion tensor. Ocean Modelling 164, 101831. doi:10.1016/j.ocemod.2021.101831.

Haigh, M., Sun, L., McWilliams, J.C., Berloff, P., 2021b. On eddy transport in the ocean. Part II: The advection tensor. Ocean Modelling 165, 101845. doi:10.1016/j.ocemod.2021.101845.

Hallberg, R., 2013. Using a resolution function to regulate parameterizations of oceanic mesoscale eddy effects. Ocean Modelling 72, 92-103. doi:10.1016/j.ocemod.2013.08.007.

Hewitt, H., Fox-Kemper, B., Pearson, B., Roberts, M., Klocke, D., 2022. The small scales of the ocean may hold the key to surprises. Nature Climate Change doi:10.1038/s41558-022-01386-6.

Holmes, R.M., Groeskamp, S., Stewart, K.D., McDougall, T.J., 2022. Sensitivity of a coarse-resolution global ocean model to a spatially variable neutral diffusivity. Journal of Advances in Modeling Earth Systems 14, e2021MS002914. doi:10.1029/2021MS002914.

Jamet, Q., Dewar, W.K., Wienders, N., Deremble, B., 2019. Spatiotemporal Patterns of Chaos in the Atlantic Overturning Circulation. Geophysical Research Letters 46, 7509-7517. doi:10.1029/2019GL082552.

Jamet, Q., Dewar, W.K., Wienders, N., Deremble, B., Close, S., Penduff, T., 2020. Locally and remotely forced subtropical AMOC variability: A matter of time scales. Journal of Climate 33, 5155-5172. doi:10.1175/JCLI-D-19-0844.1.
586

Jamet, Q., Leroux, S., Dewar, W.K., Penduff, T., Le Sommer, J., Molines, J.M., Gula, J., 2022. Non-local eddy-mean kinetic energy transfers in submesoscale-permitting ensemble simulations. Journal of Advances in Modeling Earth Systems 14, e2022MS003057. doi:10.1029/2022MS003057.

587

588

589

590

Jansen, M.F., Adcroft, A., Khani, S., Kong, H., 2019. Toward an energetically consistent, resolution aware parameterization of ocean mesoscale eddies. Journal of Advances in Modeling Earth Systems 11, 2844–2860. doi:10.1029/2019MS001750.

591

592

593

594

Jones, S.C., Abernathey, R.P., 2019. Isopycnal mixing controls deep ocean ventilation. Geophysical Research Letters 46, 13144–13151. doi:10.1029/2019GL085208.

595

596

597

Kamenkovich, I., Berloff, P., Haigh, M., Sun, L., Lu, Y., 2020. Complexity of mesoscale eddy diffusivity in the ocean. Geophysical Research Letters, e2020GL091719doi:10.1029/2020GL091719.

598

599

600

Kamenkovich, I., Rypina, I.I., Berloff, P., 2015. Properties and origins of the anisotropic eddy-induced transport in the North Atlantic. Journal of Physical Oceanography 45, 778–791. doi:10.1175/JPO-D-14-0164.1.

601

602

603

Killworth, P.D., 1997. On the parameterization of eddy transfer Part I. Theory. Journal of Marine Research 55, 1171–1197. doi:10.1357/0022240973224102.

604

605

606

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–403. doi:10.1029/94rg01872.

607

608

609

Leroux, S., Penduff, T., Bessières, L., Molines, J.M., Brankart, J.M., Sérazin, G., Barnier, B., Terray, L., 2018. Intrinsic and atmospherically forced variability of the AMOC: insights from a large-ensemble ocean hindcast. Journal of Climate 31, 1183–1203. doi:10.1175/JCLI-D-17-0168.1.

610

611

612

613

Loose, N., Bachman, S., Grooms, I., Jansen, M., 2022. Diagnosing scale-dependent energy cycles in a high-resolution isopycnal ocean model. Journal of Physical Oceanography doi:10.1175/JPO-D-22-0083.1.

614

615

616

Lu, Y., Kamenkovich, I., Berloff, P., 2022. Properties of the lateral mesoscale eddy-induced transport in a high-resolution ocean model: Beyond the

617
flux-gradient relation. Journal of Physical Oceanography doi:10.1175/JPO-D-22-0108.1.

Mak, J., Maddison, J.R., Marshall, D.P., 2016. A new gauge-invariant method for diagnosing eddy diffusivities. Ocean Modelling 104, 252–268. doi:10.1016/j.ocemod.2016.06.006.

Marshall, D.P., Maddison, J.R., Berloff, P.S., 2012. A framework for parameterizing eddy potential vorticity fluxes. Journal of Physical Oceanography 42, 539–557. doi:10.1175/JPO-D-11-048.1.

Marshall, J., 2011. A pedagogical example of Eulerian and Lagrangian flows induced by Rossby wave rectification in an ocean basin. Journal of Marine Research 69, 647–657. doi:10.1357/002224011799849318.

Marshall, J., Hill, C., Perelman, L., Adcroft, A., 1997. Hydrostatic, quasi-hydrostatic, and nonhydrostatic ocean modeling. Journal of Geophysical Research: Oceans 102, 5733–5752. URL: https://mitgcm.readthedocs.io/en/latest/, doi:10.1029/96JC02776.

Plumb, R., Mahlman, J., 1987. The zonally averaged transport characteristics of the GFDL general circulation/transport model. Journal of the atmospheric sciences 44, 298–327. doi:10.1175/1520-0469(1987)044<0298: TZATCO>2.0.CO;2.

Popinet, S., 2015. A quadtree-adaptive multigrid solver for the Serre-Green-Naghdi equations. Journal of Computational Physics 302, 336–358. doi:10.1016/j.jcp.2015.09.009.

Redi, M.H., 1982. Oceanic isopycnal mixing by coordinate rotation. Journal of Physical Oceanography 12, 1154–1158. doi:10.1175/1520-0485(1982)012<1154:0IMBCR>2.0.CO;2.

Rypina, I.I., Kamenkovich, I., Berloff, P., Pratt, L.J., 2012. Eddy-induced particle dispersion in the near-surface North Atlantic. Journal of Physical Oceanography 42, 2206–2228. doi:10.1175/JPO-D-11-0191.1.

Sérazin, G., Penduff, T., Barnier, B., Molines, J.M., Arbic, B.K., Müller, M., Terray, L., 2018. Inverse cascades of kinetic energy as a source of intrinsic variability: A global OGCM study. Journal of Physical Oceanography 48, 1385–1408. doi:10.1175/JPO-D-17-0136.1.
650

Small, R.J., Bacmeister, J., Bailey, D., Baker, A., Bishop, S., Bryan, F., Caron, J., Dennis, J., Gent, P., Hsu, H.m., et al., 2014. A new synoptic scale resolving global climate simulation using the Community Earth System Model. Journal of Advances in Modeling Earth Systems 6, 1065–1094. doi:10.1002/2014MS000363.

651

652

653

654

655

Smith, R.D., Gent, P.R., 2004. Anisotropic Gent–McWilliams parameterization for ocean models. Journal of Physical Oceanography 34, 2541–2564. doi:10.1175/JP02613.1.

656

657

658

Sun, L., Haigh, M., Shevchenko, I., Berloff, P., Kamenkovich, I., 2021. On non-uniqueness of the mesoscale eddy diffusivity. Journal of Fluid Mechanics 920. doi:10.1017/jfm.2021.472.

659

660

661

Uchida, T., Abernathey, R.P., Smith, K.S., 2017. Seasonality of eddy kinetic energy in an eddy permitting global climate model. Ocean Modelling 118, 41–58. doi:10.1016/j.ocemod.2017.08.006.

662

663

664

Uchida, T., Balwada, D., Abernathey, R.P., McKinley, G.A., Smith, K.S., Lévy, M., 2019. The contribution of submesoscale over mesoscale eddy iron transport in the open Southern Ocean. Journal of Advances in Modeling Earth Systems 11, 3934–3958. doi:10.1029/2019MS001805.

665

666

667

668

Uchida, T., Balwada, D., Abernathey, R.P., McKinley, G.A., Smith, K.S., Lévy, M., 2020. Vertical eddy iron fluxes support primary production in the open Southern Ocean. Nature communications 11, 1–8. doi:10.1038/s41467-020-14955-0.

669

670

671

672

Uchida, T., Deremble, B., Dewar, W.K., Penduff, T., 2021a. Diagnosing the Eliassen-Palm flux from a quasi-geostrophic double gyre ensemble, in: EarthCube Annual Meeting, NSF. URL: https://earthcube2021.github.io/ec21_book/notebooks/ec21_uchida_etal/notebooks/TU_05_ Diagnosing-the-Eliassen-Palm-flux-from-a-quasi-geostrophic-double-gyre-ensemble.html, doi:10.5281/zenodo.5496375.

673

674

675

676

677

678

679

Uchida, T., Deremble, B., Penduff, T., 2021b. The seasonal variability of the ocean energy cycle from a quasi-geostrophic double gyre ensemble. Fluids 6, 206. doi:10.3390/fluids6060206.

680

681
682

Uchida, T., Deremble, B., Popinet, S., 2022a. Deterministic model of the eddy dynamics for a midlatitude ocean model. Journal of Physical Oceanography 52, 1133–1154. doi:10.1175/JPO-D-21-0217.1.

683

684

685

Uchida, T., Jamet, Q., Dewar, W.K., Le Sommer, J., Penduff, T., Balwada, D., 2022b. Diagnosing the thickness-weighted averaged eddy-mean flow interaction from an eddying North Atlantic ensemble: The Eliassen-Palm flux. Journal of Advances in Modeling Earth Systems 14, e2021MS002866. doi:10.1029/2021MS002866.

686

687

688

689

690

Uchida, T., Le Sommer, J., Stern, C., Abernathey, R.P., Holdgraf, C., Alber, A., Brodeau, L., Chassignet, E.P., Xu, X., Gula, J., et al., 2022c. Cloud-based framework for inter-comparing submesoscale-permitting realistic ocean models. Geoscientific Model Development 15, 5829–5856. doi:10.5194/gmd-15-5829-2022.

691

692

693

694

695

Vallis, G.K., 2017. Atmospheric and oceanic fluid dynamics. 2 ed., Cambridge University Press.

696

697

Wilson, C., Williams, R.G., 2006. When are eddy tracer fluxes directed downgradient? Journal of Physical Oceanography 36, 189–201. doi:10.1175/JPO2841.1.

698

699

700

Young, W.R., 2012. An exact thickness-weighted average formulation of the boussinesq equations. Journal of Physical Oceanography 42, 692–707. doi:10.1175/JPO-D-11-0102.1.

701

702

703

Zanna, L., Bolton, T., 2020. Data-driven equation discovery of ocean mesoscale closures. Geophysical Research Letters, e2020GL088376doi:10.1029/2020GL088376.

704

705
<figure><img src="image_8.png" /><figcaption>Figure 1: The eddy transport tensor **K** diagnosed by inverting (5), eddy flux of $C_2$ and PV and their reconstruction from the QG ensemble for all three layers. The four components of the tensor are exhibited in the top two rows, the eddy flux $\mathbf{J}^{C_2}$ and its reconstruction $-\nabla_h \bar{\mathbf{C}}_2 \cdot \mathbf{K}$ in the middle two rows, and the tensor applied to PV ($q$) in the bottom two rows. The eddy PV fluxes in the bottom two layers are two orders of magnitude smaller than the top layer.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 2: The eddy transport tensor $K$ diagnosed by inverting (5) from the realistic North Atlantic ensemble on January 3, 1967. The top left four panels show the four components along the neutral surface associated with the ensemble-mean depth $\bar{\zeta}$ shown in the right top panel. The bottom four panels exhibit the vertical transect of the four components along 300°E.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 3: The eddy flux of tracers $C_1$ and $C_2$ and their reconstruction from the realistic North Atlantic ensemble on January 3, 1967. The top three rows present the eddy flux $\mathbf{J}^{C_i}$, reconstruction $-\nabla_h \tilde{C}_i \cdot \mathbf{K}$, and the difference between the two respectively. The hatches indicate regions where the neutral surface is shallower than 150 m. The bottom two rows exhibit the vertical transect of the flux and reconstruction along 300°E.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 4: The eddy flux of potential temperature and practical salinity, and their reconstruction from the realistic North Atlantic ensemble on January 3, 1967. The top three rows present the eddy flux $\mathbf{J}^\Theta, \mathbf{J}^S$, reconstruction $-\nabla_h(\hat{\Theta}, \hat{S}) \cdot \mathbf{K}$, and the difference between the two respectively. The hatches indicate regions where the neutral surface is shallower than 150 m. The bottom two rows exhibit the vertical transect of the flux and reconstruction along 300°E.</figcaption></figure>