Cautionary tales from the mesoscale eddy transport tensor

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 ed-dying (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 ﬂow and eddies. The tensor exhibits good performance in re-constructing the eddy ﬂuxes of passive tracers, here deﬁned as ﬂuctuations about the ensemble thickness-weighted averaged (TWA) mean. However, the inability of the tensor to reconstruct eddy ﬂuxes of QG potential vorticity, which encapsulates the eddy-mean ﬂow interaction, and other active tracers raises the question: To what extent can the diagnosed tensor be applied to inform the parametrization of mesoscale dynamics?


Introduction
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 m) to O(1000 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 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 (Gri es et al., 2015;Chassignet et al., 2020).Nevertheless, forced ocean simulations with higher model resolution have shown that even at O(10 km) mesoscale eddies and their feedback onto the jets are not su ciently resolved (Chassignet and Xu, 2021;Uchida et al., 2019Uchida et al., , 2022c;;Hewitt et al., 2022).
In order to overcome the insu cient representation of the eddies and their variability, there has been a growing field of research on how to parametrize the e↵ects of eddies in simulations at the grey zone resolution, when the eddy variability is partially resolved.The missing dynamical e↵ects, 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 e↵ects of these missing dynamical processes within the framework of eddy transport coe cients, 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 e↵ect 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 coe cients 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 a↵ect 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-di↵usive 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 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 su cient model resolution.
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 coe cients that emerge between the two ensembles.
We conclude and provide some cautionary notes in Section 4.

Methods
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).

Quasi-geostrophic ensemble
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).
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: where C i , i = 0, 1, 2, 3 are the four passive tracers and the stream function.
They are relaxed towards a profile orthogonal to each other with a relaxation time scale of T ⇠ 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.r h C i 6 = 0 where r 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

Realistic North Atlantic ensemble
We use the model outputs from the realistic simulations described in Jamet et al. (2019Jamet et al. ( , 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 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 (Young, 2012).e 1 and e 2 are the horizontal unit vectors along the neutral surface, which reduce to i and j respectively when the neutral surfaces are flat ( e r h ⇣ = 0 where e r h is the horizontal gradient operator in density coordinates).We note that due to the 3 rd order direct-space-time (DST) flux-limiter advective scheme used, there are two possible ways to define the eddy flux, i.e. [ where 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

Eddy transport tensor
As the eddy flux of tracers is generally poorly resolved in global forced and coupled ocean simulations, there has been an e↵ort to parametrize this flux in the quasi-adiabatic interior via a local-gradient based model along neutral surfaces (Redi, 1982;Gri es, 2004;Wilson and Williams, 2006;Holmes et al., 2022) where  is the scalar eddy di↵usivity, hereon referred to as the transport coe cient.The terminology of 'di↵usivity' has its origin in the works by Redi (1982); Gent and McWilliams (1990) and Gent et al. (1995), but eddies not only di↵use but advect tracers (Smith and Gent, 2004;Gri es, 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 coecient from equation ( 4), assuming an isotropic transport coe cient 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., 2015Kamenkovich et al., , 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 coe cient 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 (K K K) from a least-squares best fit to (Plumb and Mahlman, 1987;Abernathey et al., 2013;Bachman and Fox-Kemper, 2013) where for the QG ensemble and for the realistic NA ensemble respectively.The least-squares fit can be estimated as the Moore-Penrose pseudo inverse of G G 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 Berlo↵, 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 overdetermined by using four tracers.The assumption is that by keeping it over where J J J C i and G G 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).

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.

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 m 2 s 1 ) in the first and second layer.The bottom layer is quiescent with transport coe cients on the order of O(10 2 m 2 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 km.This suggests that the eddies are broadly working to dissipate small-scale variance, as one might expect.Focusing on the first layer,  uv tends to be coherently negative south of the jet and positive north of the jet.In the lower two layers,  uv and  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 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 where h•i 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 |" C i j | computed as 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 coe cients are not shown for the middle and bottom layer as it is evident from Fig. 1 that the flux and reconstruction are completely decorrelated.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., 2015Bachman et al., , 2020)).However, it is important to note that such a good fit, although not perfect (Figs.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 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.
bottom two rows of Fig. 1, we see that while the reconstruction r h q • K K 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).

Realistic North Atlantic ensemble
Following the QG analyses, we now present the four components of the tensor K K K from the NA ensemble.The magnitude of the tensor components are similar to the QG ensemble, reaching up to O(10 4 m 2 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 m 2 s 1 ) in the abyssal ocean (bottom four panels in Fig. 2).
We remind the reader that we define the eddy flux as following the convention of the TWA framework, the reconstructed eddy flux 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 coe cients 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.
Tracer Entire domain 10 N-30 N r C 0 1 0.918 0.956 In examining the reconstruction, we limit it to regions where the neutral surface is deeper than 150 m and to grids upon the 10 ⇥ 10 coarse graining included no land cells in order to minimize the e↵ect 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 C i and J C i reconstructed 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 normed relative errors for passive tracers generally decrease when they are taken over the subdomain of 10 N-30 N.
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.
Tracer Entire domain 10 N-30 N 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 ⇥ reconstructed and J S reconstructed capture the sign structure south of 30 N in J ⇥ 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.

Conclusion and discussion
In this study, we have diagnosed the anisotropic mesoscale eddy transport tensor K K 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 Berlo↵, 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 K 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 S S = (K K K + K K K T )/2 where K K K T is the transpose).The physical significance of the eigenvalues is that the eddies tend to di↵use the tracers in the direction of the eigenvectors with an associated transport coe cient of each respective eigenvalue (Abernathey et al., 2013;Haigh et al., 2021a).We find a robust 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 coe cients 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 ( uu ,  vv ) tend to be positive while ( uv ,  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 2 nd -order centered scheme.We also diagnosed the tensor with the eddy fluxes defined as C i but the performance of reconstruction deteriorated compared to when J C i def = [ u 00 C 00 i (not shown).We speculate that the linear gradient flux model ( 5) is unable to capture the non-linearities in the fluxlimiter 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 coe cients in the top layer diagnosed from a subset of 12, 25 and 50 members in  (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 K K K as our interest has been in whether this tensor can be used to approximate the total eddy PV flux itself (F # in 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 K 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 a↵ect 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.
where (•) is the ensemble mean operator and (•) 0 def = (•) (•).i, 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.
and the relaxation time scale of 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 thicknessweighted averaged (TWA) framework and quasi geostrophy is that in the latter, the layer thickness ( def = ⇣ b) only fluctuates on the order of Rossby number (⇣ is the depth of the neutral surface and the subscript f (•) denotes 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 e↵ect 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 (•) 00 , and v) 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 inUchida et al. (2022b).
determined, K K 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 coe cients 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 J J and G G G was then normalized by horizontal median of the magnitude of each mean tracer gradient flux (i.e. (J J J C i j would return positive values if the reconstruction had the correct sign as the actual flux even if it had an o↵set.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. S3and S4), suggests that the transport tensor model with four free parameters is a good model to represent how the eddies flux di↵erent 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.
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 coe cients 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).
Fig. S11.While noisier, the transport coecients 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 di↵erent eddy transport coe cients.In other words, passive and active tracers have di↵erent 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 (r 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 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.

Figure 1 :
Figure1: The eddy transport tensor K K 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 J C2 and its reconstruction r h C 2 • K K 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.