## Abstract

A new approach to parameterizing subgrid-scale processes is proposed: The impact of the unresolved dynamics on the resolved dynamics (i.e., the eddy forcing) is represented by a series expansion in dynamical spatial modes that stem from the energy budget of the resolved dynamics. It is demonstrated that the convergence in these so-called energy modes is faster by orders of magnitude than the convergence in Fourier-type modes. Moreover, a novel way to test parameterizations in models is explored. The resolved dynamics and the corresponding instantaneous eddy forcing are defined via spatial filtering that accounts for the representation error of the equations of motion on the low-resolution model grid. In this way, closures can be tested within the high-resolution model, and the effects of different parameterizations related to different energy pathways can be isolated. In this study, the focus is on parameterizations of the baroclinic energy pathway. The corresponding standard closure in ocean models, the Gent–McWilliams (GM) parameterization, is also tested, and it is found that the GM field acts like a stabilizing direction in phase space. The GM field does not project well on the eddy forcing and hence fails to excite the model’s intrinsic low-frequency variability, but it is able to stabilize the model.

## 1. Introduction

It is crucial that climate models are able to accurately simulate the climate’s internal variability, in addition to the climate’s mean state and externally forced climate changes (IPCC 2013). For example, a correct representation of internal climate variability is needed in climate change detection and attribution studies. Such studies are based on signal-to-noise estimates for which the climate’s intrinsic low-frequency variability (LFV) must be estimated, at least in part, from long control integrations of climate models. Also for climate prediction studies a correct representation of the intrinsic climate variability is crucial such that internally generated sources of predictability can be exploited. Finally, the ability of models to make quantitative projections of changes in climate variability, including the statistics of extreme events under a warming climate, is dependent on an accurate representation of the climate’s internal variability.

The climate’s intrinsic LFV is typically described by large-scale modes of climate variability, which are often either statistical eigenmodes (e.g., EOFs) or dynamical eigenmodes (e.g., linear instability modes; see von Storch and Zwiers 1999; IPCC 2013; Dijkstra 2016). The modes of climate variability are characterized as large scale because they include large spatial structures such as basinwide coupled modes of ocean–atmosphere variability (e.g., El Niño–Southern Oscillation), Rossby wave trains, midlatitude jets and storm tracks, and so on.

In particular with respect to the ocean, a number of LFV modes (i.e., from multiannual to multidecadal time scales) have been described (Deser et al. 2010; Dijkstra 2016). It is clear from observations that multidecadal patterns of sea surface temperature variability exist, such as the Atlantic multidecadal oscillation (Schlesinger and Ramankutty 1994; Kerr 2000) and the Pacific decadal oscillation (Mantua et al. 1997; England et al. 2014). Most of these modes have a particular regional or even global manifestation whose amplitude can be larger than that of human-induced climate change. For example, intrinsic multidecadal variability of the ocean heat content has been held responsible for the relatively low recent trend in the global mean surface temperature anomaly, also referred to as the “global warming hiatus” (Meehl et al. 2011, 2013).

However, care is required when interpreting modes of climate variability since 1) their interpretation depends on how one separates modes of variability from forced changes in the time mean, 2) they may change drastically in space, structure, or probability distribution in response to climate change, and 3) in strongly nonlinear regimes they may be not strictly large scale but the large-scale structures can be entangled with smaller-scale structures such that some modes of climate variability may not be entirely representable in climate models with a too coarse spatial resolution.

Moreover, due to the fact that the real ocean dynamics resides in a highly turbulent regime (with a large Reynolds number leading to a high-dimensional unstable manifold on the attractor) it still remains hard to understand the exact physical mechanisms behind the ocean’s LFV (Berloff and McWilliams 1999a; Hogg et al. 2005; Dijkstra 2016). Plenty of model studies analyzing eddy-resolving ocean models show that LFV in such models is commonplace (Berloff and McWilliams 1999a; Hogg et al. 2005), and it is now known that the collective action of oceanic mesoscale eddies is one of the main drives of the midlatitude LFV (Kwon 2010). But at the same time the strong eddy field can obscure many features of the circulation, making it difficult to agree upon the mechanisms underpinning the variability (Hogg et al. 2005; Dijkstra 2016).

Central questions still need further clarification: Which part of the ocean’s LFV is completely intrinsic to the ocean and which part involves a dynamical coupling to the atmosphere? Which part of the ocean’s intrinsic LFV can be traced back to stationary modes at high viscosity (i.e., low-order bifurcations) and which part represents a genuinely eddy-driven turbulent phenomenon (i.e., physical mechanisms solely active at high Reynolds numbers) (Hogg and Blundell 2006; Berloff et al. 2007; Le Bars et al. 2016; Dijkstra 2016)?

Clarification of these questions is hampered by the fact that computational limitations force most studies on climate variability to employ climate models with ocean components that do not resolve the internal Rossby deformation radius (Hallberg 2013). In these coarse-resolution ocean models (typically operating at a horizontal resolution of 1°) usually deterministic eddy parameterizations are applied that are based on diffusive terms that aim to model the potential and kinetic energy transfer from the mean field to the eddy field. These diffusive eddy parameterizations achieve a reasonable representation of the time-mean effect of the mesoscale eddy field on the time-mean flow (Bryan et al. 2014; Griffies et al. 2015; Viebahn et al. 2016), but they are not able to excite the ocean’s internal LFV observed in eddy-resolving ocean model simulations (Le Bars et al. 2016). Consequently, the estimation of internal variability uncertainty (stemming from the chaotic nature of the system) in climate change detection or projections of climate change is still strongly hampered by model uncertainty (i.e., limitations of a model’s representation of the chaotic nature of the system) in many current climate change studies.

Hence, the search for suitable eddy parameterizations remains a challenging theoretical topic with clear practical dimension. Recently, efforts have been made toward eddy parameterizations that aim to step out of the diffusive parameterization framework and try to represent the eddy effects in terms of stochastic eddy forcing (Berloff 2005c; Grooms and Majda 2013; Porta Mana and Zanna 2014; Verheul et al. 2017). Stochastic climate modeling is based on the concept of scale separation in time (Franzke et al. 2015), namely, that the state vector of the system can be decomposed into fast modes and slow (low frequency) modes such that the time scales of these modes strongly differ. The impact of the fast modes on the slow modes appears as eddy forcing in the equations of motion for the slow modes. The development of stochastic climate models then proceeds by accounting for the effects of the unresolved fast modes in a stochastic fashion.

Moreover, for models formulated in physical space (like most ocean models) the essential difference between a high-resolution model and a low-resolution model is the extent of spatial information. The eddy forcing actually represents the impact of the spatially unresolved (or subgrid scale or small scale) processes on the spatially resolved (or larger scale) processes. Consequently, for models in physical space time-scale separation should imply scale separation in space. That is, the patterns associated with slow variability should exhibit strictly large-scale spatial structures whereas the patterns associated with fast variability should show strictly small-scale spatial structures. Otherwise, the slow modes and the fast modes cannot be disentangled on the low-resolution model grid.

However, scale separation only holds for regimes in which scales are weakly coupled whereas in turbulent regimes different scales are strongly nonlinearly coupled. The lack of time-scale separation introduces non-Markovian memory effects and complicates the derivation of systematic parameterizations. The lack of scale separation in space implies that the dynamical modes are multiscale patterns both in the horizontal and vertical directions. For example, for the midlatitude ocean gyres it is found that due to the background flow most eigenmodes contain a large variety of scales (Shevchenko et al. 2016). In this case, the LFV is not a single-mode pattern, but rather a coherent pattern phenomenon consisting of a large number of short period phase-related eigenmodes interacting with each other. We note that this can apply to both (high resolution) statistical eigenmodes like EOFs (Gille and Kelly 1996) and linear eigenmodes on a background flow (Shevchenko et al. 2016). Obviously, the small-scale structures of the dynamical modes are not resolvable on a low-resolution model grid.

In this study, we approach the formulation of eddy parameterizations in the following way: First, we define the resolved dynamics and the corresponding instantaneous eddy forcing via spatial filtering (instead of, e.g., temporal averaging) such that we can account for the representation error of the equations of motion on the low-resolution model grid. Second, we represent the impact of the unresolved dynamics on the resolved dynamics (i.e., the eddy forcing) in terms of a series expansion in dynamical spatial modes that stem from the energy budget of the resolved dynamics. These so-called energy modes exhibit strictly large-scale spatial patterns and are equipped with a clear physical interpretation.

In section 2, the resolved dynamics and the related instantaneous eddy forcing are defined. We describe our eddy-resolving ocean model and its LFV in section 2a. Our spatial filtering approach is introduced in section 2b. The corresponding filtered equations of motion and the related eddy forcing terms are presented in section 2c, and in section 2d we analyze the resulting large-scale and small-scale energetics. Subsequently, section 3 deals with developing and testing closures with a focus on the baroclinic energy pathway. We show how we can test parameterizations in a high-resolution model (section 3a) and test the performance of the standard closure of the baroclinic energy pathway in ocean models [i.e., the Gent–McWilliams (GM) parameterization; Gent and McWilliams 1990] in section 3b. Finally, in section 3c we define and test the representation of the eddy forcing in energy modes with a focus on representing the baroclinic energy pathway. We end with a summary (section 4) and discussion (section 5).

## 2. Framework: Eddy forcing of the large-scale flow

Our general starting point is the following (see, e.g., Berloff 2005a): First, an eddy-resolving (ER) model is given (section 2a) in order to obtain a reference solution, say with state vector *ψ*, which contains both the large-scale and eddy components. Second, a non-eddy-resolving (non-ER) model is supposed to have the same general setup as the ER model (e.g., type of governing conservation equations, domain size, and boundary conditions; see section 2c), but the former has a significantly coarser horizontal grid resolution (by a factor of 10 in this study). Consequently, the non-ER model has far fewer degrees of freedom and it can only solve for the large-scale flow evolution. Moreover, the non-ER model may contain additional dynamical terms in the governing conservation equations (e.g., the current deterministic eddy parameterizations) that are supposed to parameterize part of the interactions between large-scale components and (sub)mesoscale eddy components.

Finally, the eddy forcing (EF) is a (not necessarily unique) dynamical term that still needs to be added to the governing conservation equations of the non-ER model at hand such that the non-ER solution, say with state vector , correctly approximates the large-scale structure of the original flow (i.e., of the ER model solution *ψ*). That is, the EF represents interactions between the large-scale flow and eddy fluctuations that are relevant for the large-scale flow evolution. The precise form of the EF depends on 1) the specific definition of the large-scale structure of the original flow (section 2b) and 2) the eddy parameterizations already included in the chosen non-ER model equations (section 3).

### a. Eddy-resolving ocean model exhibiting low-frequency variability

We consider a standard model of idealized ocean dynamics, namely, quasigeostrophic (QG) potential vorticity (PV) equations in a classical double-gyre configuration (see, e.g., Vallis 2006). The fluid-dynamic model describes idealized, wind-driven midlatitude ocean circulation with prescribed density stratification in a flat-bottom square basin with north–south and east–west boundaries. We employ the QG PV conservation equations for two isopycnal layers, representing the simplest description of baroclinically unstable dynamics (Olbers et al. 2012). These are

where the PV of the two isopycnal layers is given by

with the interface displacement , and horizontal velocities given by .

In our numerical model simulations, the flow is driven at the surface by the asymmetric double-gyre zonal wind stress (as, e.g., in Berloff 2005a,c):

where N m^{−2}, and km is the size of the square basin with . The first internal Rossby radius of deformation, , represents a length scale of baroclinic eddies. It is set to *R* = 40 km, a typical value for the midlatitude ocean circulation. We use mean isopycnal layer thicknesses of m and m, such that the mean ocean depth is *H* = 4000 m. We also use typical values for the mean density of seawater kg m^{−3} and the reference Coriolis parameter s^{−1}, such that we have for the meridional variation of the Coriolis parameter m^{−1} s^{−1} and for the reduced gravity m s^{−2}. Finally, we use an eddy-resolving horizontal resolution of 10 km with a correspondingly small lateral viscosity coefficient, m^{2} s^{−1}, as well as no-slip boundary conditions (similar to Berloff 2005a,c). The reference simulation is 500 years long and we analyze daily output.

Figure 1 shows a snapshot (Fig. 1c) and a temporal average (Fig. 1d) of the upper-layer streamfunction [similar to Figs. 1a,c in Berloff (2005a), and Figs. 1a and 2a in Berloff (2005c)]. The upper-ocean time-mean circulation (Fig. 1d) consists of the southern (subtropical) and northern (subpolar) gyres that fill about 2/3 and 1/3 of the basin, respectively, which is consistent with the wind stress pattern. The time-mean flow is characterized by the Sverdrup balance in most parts of the basin. Only in regions related to the pair of the western boundary currents and their eastward jet (EJ) extensions do nonlinear and frictional terms become dominant (Pedlosky 1996). We note that for our specific model setup the boundary currents do not merge with each other but the subpolar gyre enters the subtropical region near the western boundary such that the point of separation from the coast of the subtropical western boundary current is pushed southward relative to the line of zero wind stress curl (similar to Berloff 2005a,c). This is a robust regime that appears at large Reynolds number in the stratified and baroclinically unstable double-gyre flow with no-slip boundary conditions (e.g., Haidvogel et al. 1992; Berloff and McWilliams 1999b; Siegel et al. 2001). In terms of the fluctuations, the basin can be partitioned into the more energetic “western” part, characterized by strong vortices, and the less energetic “eastern” part, dominated by the planetary waves [see Berloff et al. (2002) for details].

The corresponding reservoirs of kinetic energy (KE) and available potential energy (PE) are given by

The two reservoirs are governed by the following conservation equations [obtained by multiplying Eqs. (1) and (2) with followed by global integration]:

where represents the conversion between PE and KE, and the generation of KE [] and the dissipation of KE [] are given by

Figure 2 shows the temporal evolution of the energetics of the reference simulation. The PE (Fig. 2a) exhibits clear cycles of decadal variability. The about 4 times smaller KE also shows cycles of decadal variability, which lags the variability in PE by 1–2 years. The wind energy input (Fig. 2b) is balanced by lateral dissipation , with both showing also significant high-frequency variability on top of low-frequency variability, with higher variance in than in .

Figure 1 also shows upper-layer streamfunction anomalies corresponding to a low (Fig. 1e) and a high (Fig. 1f) in PE (see years 52 and 56 in Fig. 2a). These anomaly patterns are similar to those shown in Berloff et al. (2007) (see their Figs. 2 and 4), and demonstrate that the variability is concentrated around the subtropical EJ. More precisely, the decadal transitions are related to coherent meridional shifts and variations of the intensity of the subtropical EJ, likely governed by the nonlinear adjustment of the combined EJ–eddies system [see Berloff et al. (2007) for details].

### b. Flow decomposition into large-scale and eddy components via spatial mode filtering

In this study, the large-scale flow structure is determined by spatial mode filtering. For that the ER model solution *ψ* is expanded in a set of spatial filter modes :

Note that the spatial filter modes are time independent (i.e., nondynamical). The corresponding large-scale (or filtered) component of *ψ* is given by a truncated expansion (equivalent to applying a sharp spectral filter):

The cutoff has to be determined such that the retained spatial filter modes have a consistent representation on the coarse-resolution grid of the non-ER model (see the consistency conditions below). The non-ER model solution (denoting spatial fields on the non-ER model grid by a hat) would then optimally be given by^{1}

More precisely, the specification of the spatial filter modes and the cutoff is guided by the following three consistency conditions:

Scale-content (or image) consistency (SCC): The scale content of the spatial filter modes has to be resolvable by the non-ER model grid resolution in order to avoid aliasing effects. The scale content of a spatial pattern is typically measured by the familiar Fourier modes (i.e., eigenmodes of the Laplacian). The corresponding cutoff is given by the well-known Nyquist criterion, which states that the smallest wavelength included in must not be smaller than twice the grid spacing of the non-ER model. In particular, we note that filtering of the ER model reference solution

*ψ*has to be done on the ER model grid (i.e., by using , and not by using and the injection of*ψ*on the coarse grid) since otherwise aliasing errors occur.Boundary conditions consistency (BCC): The large-scale flow is supposed to be a solution of the non-ER model equations and, hence, has to satisfy its boundary conditions. In the governing model equations, the differential operator with the highest-order derivative (typically related to dissipation) determines the number of boundary conditions that have to be specified. Consequently, the eigenmodes of this differential operator represent a set of spatial modes which are always able to satisfy the boundary conditions (i.e., span the correct function space), and, hence, represent the first choice if the Fourier modes (i.e., eigenfunctions of the Laplacian) cannot satisfy the boundary conditions.

Dynamical consistency (DC): The conservation equations governing the evolution of are obtained by filtering the ER model equations (see section 2c). The non-ER model equations are supposed to represent these equations except that the terms including interactions with eddy components are replaced by eddy parameterizations. For that to hold, the spatial derivatives appearing in the governing equations have to be similar for both and ; that is, the differences in computing dynamical terms on the different grids must be not be significant. Otherwise, the EF would not solely represent the interactions between the large-scale flow and eddy fluctuations that are relevant for the large-scale flow evolution but would also have to compensate for differences simply induced by computing the dynamical budget of the large-scale flow on different grids.

^{2}Consequently, one must generally require .

In our ER model no-slip boundary conditions are applied such that Fourier modes cannot be used as the spatial filter modes (see the BCC condition). Consequently, we use as spatial filter modes the eigenmodes of the bi-Laplacian,

for which no-slip boundary conditions can be prescribed. Note that represents the globally integrated lateral dissipation related to the normalized^{3 } since .

Figure 3 shows selected leading eigenmodes (orthonormalized) of the bi-Laplacian with no-slip boundary conditions computed on the high-resolution (i.e., 10 km) grid. The overall structure (i.e., the scale content) of the is still very similar to the Fourier modes (but note that the quantitative differences are nevertheless global and *not* only localized at the boundary). Computing the eigenspectrum of on both the ER model grid (i.e., 10-km resolution) and the non-ER model grid (i.e., 100-km resolution) enables us to specify a cutoff in accordance with condition DC. Figure 4 shows the corresponding eigenvalues and their relative difference. As a threshold we choose 10% relative difference in globally integrated lateral dissipation, which implies . The corresponding relative difference in globally integrated kinetic energy (also shown in Fig. 4) is about 5%.

Figure 1 also shows the corresponding snapshot (Fig. 1a) and temporal average (Fig. 1b) of the large-scale (i.e., filtered with ) upper-layer streamfunction. The overall structure of the double-gyre circulation is captured by the large-scale flow in both cases. In particular, the separation point of the subtropical western boundary current is exactly recovered. However, local differences are obvious (also in the time-mean patterns); for example, the locations of the local extremes are shifted. Consequently, spatial filtering and temporal filtering are not equivalent.

### c. Conservation equations of the large-scale flow

The conservation equations governing the evolution of the large-scale flow are obtained by applying the filtering operation to the QG PV budget in Eqs. (1) and (2). Filtering and application of the bi-Laplacian obviously commute for the eigenmodes of the bi-Laplacian. However, for no-slip boundary conditions filtering with the eigenmodes of the bi-Laplacian does not commute with both the zonal derivative (i.e., linear beta term) and the Laplacian. Hence, filtering of the governing equations [Eqs. (1) and (2)] leads to the following equations governing the filtered flow^{4}:

where the filtered PV reads

The residual PV fluxes , representing interactions between the large-scale flow and eddy fluctuations that are relevant for the large-scale flow evolution, are given by

with the residual advection of PV and the residual related to the time tendency of relative PV given by

The residual advection of PV can be further decomposed into , with

which are related to residual planetary vorticity advection, residual nonlinear momentum fluxes, and residual buoyancy fluxes (i.e., interface displacements), respectively.

In the following, we focus on a twofold decomposition of the residual PV fluxes into

where represents the part of that is related to the vertical density distribution/layer interaction/interface height/APE, whereas is related to the horizontal eddy PV fluxes. In this study, the main focus will be on (see section 3).

### d. Lorenz energy cycle

The Lorenz energy cycle (LEC) describes the balances of four mechanical energy reservoirs, the large-scale circulation’s kinetic energy ([KE]) and available potential energy ([PE]), the eddy kinetic energy (KE′) and eddy available potential energy (PE′). The four reservoirs are given by

and they are governed by the following conservation equations, which are obtained by multiplying Eqs. (13) and (14) with and global integration:

The respective generation and dissipation terms are given by

and the terms related to energy exchange between the large-scale flow and eddy components read

Note that all LEC terms are instantaneously given due to our spatial (instead of temporal) filtering approach.

Figure 5 shows the different terms of the LEC averaged in time (over 500 years of daily output) and summarizes the time-mean state and variance of the different energy reservoirs and energy pathways (for the reference simulation described in section 2a). The filtered terms and capture 96% and 89% of the full (i.e., unfiltered) and PE values, respectively, implying that and PE are dominated by large-scale structures. In contrast, is very small [1.8% of ], implying that is dominated by small-scale structures. Consequently, almost all of has to be transferred to the eddy field via eddy fluxes. Both conversion terms, and , have the same order of magnitude but dominates (almost twice as large both in temporal average and variance). The two eddy energy reservoirs are of similar magnitude with KE′ being almost 5 times larger than [KE] (capturing 83% of KE).

The overall picture is similar to the one found in realistic global ocean models (e.g., von Storch et al. 2012). Common in both the ocean and the atmosphere is that the dominant power pathway is the baroclinic pathway [PE] characterized by a conversion from the large-scale available potential energy to the eddy available potential energy that has about the same magnitude^{5} as the conversion from the eddy potential energy to the eddy kinetic energy. That is, as in the atmosphere, oceanic mesoscale eddies are, to a large extent, generated by baroclinic instability, which is the main mechanism in converting the large-scale available potential energy into the eddy kinetic energy in the ocean. Moreover, and in contrast to the atmosphere, the two conversion terms connected to the large-scale kinetic energy [KE]—that is, and —are directed away from [KE] in the ocean. That is, the two main power pathways in the ocean are [KE] [PE] and [KE] . The oceanic large-scale circulation, being fueled by the winds, converts its kinetic energy into the large-scale available potential energy by Ekman pumping. This conversion substantially facilitates density differences and hence the large-scale available potential energy from which the baroclinic pathway originates. The oceanic large-scale circulation converts also its kinetic energy into the eddy kinetic energy.

Figure 2 shows different terms of the LEC evolving in time. The variability in (Fig. 2a) and (Fig. 2b) is essentially identical to the variability in and , respectively. This means that the low-frequency variability in these fields is indeed large scale, and hence can in principle be adequately captured by a non-ER model. The converse is true for lateral dissipation (Fig. 2b) for which also the variability (next to the time-mean value) of its large-scale component is very small. The KE reservoir (Fig. 2a) represents an intermediate quantity in the sense that its large-scale component [KE] only captures part of the low-frequency variability.

The large-scale wind energy input (Fig. 2b) is balanced by the energy transfer to the eddy components via (Fig. 2c) and (Figs. 2d,e). Most importantly, both and regularly show backscatter, that is, energy transfer from the eddy components to the large-scale components. Moreover, the variances of and [both part of Eq. (27)] are significantly larger than the variances of and . We note that and are highly anticorrelated with a correlation coefficient of −0.89 {they tend to be positively correlated when }, whereas the correlation coefficient of and { and } is 0.42 (0.12).

## 3. Closures for the baroclinic energy pathway

In stratified flows two distinctively different types of energy conversions between large-scale and eddy components exist: the energy conversion involving density perturbations [see Eq. (34)], and the energy conversion solely related to (horizontal) velocity perturbations [see Eq. (33)]. In the temporal average (see Fig. 5), the latter represents a sink of [KE], whereas the former represents a sink of [PE] as part of the baroclinic energy pathway, [PE] PE′ KE′. Instantaneously, both conversion terms can also backscatter, that is, transfer energy from the small-scale components to the large-scale components (Figs. 2c,e). In a non-ER model these two energy transfers have to be adequately modeled. In this study, we focus on closures for , corresponding to in the large-scale PV budget [see Eqs. (34) and (20)], and leave the development of adequate closures for (i.e., ) for future work {note that generally dominates over ; see Fig. 5}.

### a. Testing closures in an eddy-resolving model

To be able to isolate the direct effects of and the performance of corresponding closures we adopt the following approach: For a large-scale flow defined via spatial filtering (section 2b) the corresponding conservation equations (section 2c) can be computed instantaneously from the corresponding eddy-resolving model equations (section 2a). In other words, the non-ER model, Eqs. (13) and (14), can be considered as part of the ER model, Eqs. (1) and (2). To be able to test closures for [Eq. (20)] in an isolated way, that is, without the need to also parameterize , we perform simulations with the ER model Eqs. (1) and (2) in which we employ the following decomposition of the Jacobian *J* at every time step:

where corresponds to the small-scale component^{6} of . Note that can be computed from the large-scale fields and only redistributes large-scale energy but does not contribute to large-scale energy dissipation/generation.

Then a parameterization of , say , can be tested by performing simulations with the ER model Eqs. (1) and (2) and including at every time step the replacement in Eq. (36). That is, the large-scale component of the Jacobian, (which is needed in the non-ER model), is parameterized whereas the small-scale component, , remains explicitly computed. We emphasize that the “true” is always available since we solely perform simulations with the ER model. Hence, quantities like the relative error can be computed at every time step. As demonstrated in the following (sections 3b and 3c), it is by no means a trivial task to parameterize in such a way that the energy level and low-frequency variability of the large-scale flow are captured.

### b. Standard GM parameterization

In general, the GM parameterization is interpreted as the standard downgradient parameterization for the horizontal component of the isopycnal eddy flux (Vallis 2006; Olbers et al. 2012). In a layer model, this corresponds to downgradient diffusion of interface displacement *η*. More precisely, the isopycnal interface PV flux is given by . Assuming a Reynolds decomposition into mean (denoted by an overbar) and eddy (denoted by a prime) components (e.g., via temporal averaging; see, e.g., Pope 2000), the isopycnal interface eddy PV flux is given by , and the GM parameterization reads

where *K*_{GM} is an interfacial diffusivity typically *O*(1000) m^{2} s^{−1}. Finally, the divergence of the interface eddy PV flux (which actually appears in the mean PV budget) becomes

The equivalent to in case the eddy components are defined via spatial filtering is [Eq. (20)]. That is, in our case the GM parameterization reads

such that the unresolved buoyancy fluctuations are represented as local interfacial diffusion. Inserting Eq. (39) into Eq. (34) and assuming a spatially constant *K*_{GM} > 0 we get

Consequently, the GM parameterization represents a sink of [PE] at every instant of time and, hence, excludes any backscatter [Eq. (40) actually corresponds to the kinetic energy of the large-scale baroclinic mode].

#### 1) Constant GM diffusivity

A constant *K*_{GM} can be directly estimated from the energetics of the reference simulation by combining the temporal average (denoted by an overbar) of Eq. (34), shown in Fig. 5, and the temporal average of Eq. (40), such that^{7}

This way the GM parameterization accounts exactly for the time-mean [PE] dissipation, given the reference large-scale flow. For our model results we get a typical value of m^{2} s^{−1}.

Figure 6a shows time series of PE resulting from simulations in which the GM parameterization with a constant *K*_{GM} is employed (blue and green lines). To assure numerical stability *K*_{GM} ≥ 1500 m^{2} s^{−1} is necessary in our model.^{8} The GM parameterization does its job by extracting [PE] from the large-scale flow such that a statistical equilibrium results. However, the low-frequency variability exhibited by the reference simulation (black line) is absent. The dynamics exclusively reside below the PE-level of the low-PE regime of the reference simulation. That is, the low-frequency transitions in phase space to the high-PE regime are suppressed in case the GM parameterization with a constant *K*_{GM} is used. Presumably, backscatter is necessary for the dynamics in order to be able to reach high-PE states.

#### 2) Time-dependent GM diffusivity

For a time-dependent *K*_{GM} the GM parameterization [Eq. (39)] reads

We diagnose *K*_{GM} from the model results via projection on (equivalent to a least squares estimation), that is,

where represents the explicitly computed residual PV flux. That is, *K*_{GM} represents the expansion coefficient of in [see also section 3c, Eq. (44)].

Figure 6a shows the time series of PE resulting from a simulation in which the GM parameterization is employed with *K*_{GM} obtained from projection [i.e., Eq. (43)] at every time step (red line). Now a form of low-frequency variability is indeed excited but the corresponding high-PE regime resides at and below the low-PE regime of the reference simulation (black line), and the PE variability has smaller variance (see also the second row of Table 1). The low-frequency variability actually oscillates around the PE level of the simulation in which the GM parameterization with *K*_{GM} = 1500 m^{2} s^{−1} is employed (blue line). This is consistent with the fact that the time mean of *K*_{GM} obtained from projection is given by m^{2} s^{−1} (see the second row of Table 1 and also the next paragraph). Consequently, also in case of the GM parameterization with a time-dependent *K*_{GM} obtained from Eq. (43) the transitions in phase space to the high-PE regime of the reference simulation are not captured.

Figure 6b shows the estimated pdf of *K*_{GM} computed from Eq. (43) and either employed in the GM parameterization (red) or just diagnosed from the reference simulation (black). Both *K*_{GM} distributions are unimodal and slightly positively skewed. Most striking, however, is that in both cases *K*_{GM} captures a significant amount of negative values. Negative *K*_{GM} values are not consistent with a diffusion model. Hence, the low-frequency variability (red line in Fig. 6a) presumably emerges from the wrong reason, namely, backscatter due to a negative diffusivity. We also note that the temporal average and standard deviation of *K*_{GM} is significantly smaller when only diagnosed from the reference simulation (448 ± 697 m^{2} s^{−1}) than when the GM parameterization is actually applied (1585 ± 1374 m^{2} s^{−1}; see also the second row of Table 1). We discuss this difference in detail in the next section [sections 3c(3) and 3c(4)].

Finally, we emphasize that the relative error of the GM parameterization is about 97% and hence extremely high (see the second row of Table 1). This holds for when the GM parameterization is employed as well as for when the GM parameterization is just diagnosed from the reference simulation [see also Fig. 7a, discussed below in sections 3c(3) and 3c(4)].

### c. Dynamical spatial mode representation of the eddy forcing based on energetics

It is well known that the diffusive closure approach is limited since eddies also act upgradient in geophysical turbulence (Starr 1968; Berloff 2005a), implying energy transfer from the eddy components to the large scale (i.e., backscatter; see Figs. 5c,e). Consequently, instead of aiming for an improved turbulent diffusion closure (e.g., via a spatially/temporally/stochastically varying eddy diffusivity tensor) we seek for additional dynamical large-scale spatial fields (next to the large-scale isopycnal gradient) to represent the eddy forcing more adequately. That is, in order to extend or replace the GM parameterization we think in terms of a dynamical^{9} spatial mode expansion of the eddy forcing,

with time-dependent spatial modes , and evolution coefficients . The GM parameterization in Eq. (42) represents a special case with , , and .

Optimally, the spatial modes can be efficiently obtained from terms of the large-scale flow equations, and the evolution coefficients have clear dynamical or statistical properties such that they may be modeled deterministically or stochastically, . Also a small set of modes should be sufficient in order to assure feasibility. However, dynamical modes are typically constructed via generalized eigenproblems (e.g., linear instability modes; Dijkstra 2005; Berloff 2005b; Shevchenko et al. 2016) or optimization problems (e.g., Lyapunov vectors, CNOPs; Dijkstra 2013; Dijkstra and Viebahn 2015) and, hence, are generally expensive to compute, if at all.

#### 1) Specification of spatial energy modes

In this study, we explore whether spatial fields that stem from the large-scale energetics can suit as dynamical spatial modes [as in Eq. (44)] to parameterize the eddy forcing. More precisely, we focus on large-scale available potential energy budget, Eq. (27), since the eddy forcing related to the baroclinic route, , directly appears therein. The capital letters in Eq. (27) denote globally integrated LEC terms. To express the respective LEC terms as spatially extended PV fields (i.e., as integral kernels of the globally integrated energetics) we use lowercase letters. We then have

with

Note that multiplication of Eq. (45) with and global integration gives Eq. (27). For the two-layer model considered in this study, Eq. (45) corresponds to the PV evolution equation of the large-scale interface displacement (i.e., first baroclinic mode), which is obtained by subtracting Eq. (13) from Eq. (14). Combining Eq. (20) with Eq. (48) gives

such that is solely expressed in terms of large-scale (i.e., filtered) quantities.

This motivates us to consider the following two types of dynamical spatial energy modes

where is related to the temporal change of the APE reservoir at previous time , and is related to the conversion between large-scale APE and KE at previous time . Here *τ* represents the lag relative to the current time *t*. Note that and are not available in a numerical model at time *t* but are given only after the equations of motion are solved. Additionally, the temporal derivatives of the spatial energy modes and can be considered (since e.g., these may improve the convergence behavior of Eq. (44) analogous to a Taylor expansion).

Hence, in terms of a numerical model with a discrete time step the overall set of spatial energy modes reads

The set is obviously infinite. Moreover, the energy modes are generally nonorthogonal. Note that the eddy forcing of the previous time step [i.e., ] is exactly given via the energy fields and [see Eq. (52)]. Consequently, it is essentially the increment of the eddy forcing, , that has to be modeled by Eq. (44) with energy modes.

#### 2) Selection of finite subset of spatial energy modes

To compute a dynamical spatial mode expansion of the eddy forcing as in Eq. (44) in a numerical model one has to select a finite subset of energy modes out of . A detailed analysis of (finding) the optimal subset of energy fields is a topic for future research (see discussion section 5). In this study we investigate the following subsets of :

where represents the lag step size, and *n* determines the cardinality of , given by . Note that for each the contained energy modes vary in time but and *n* are fixed.

In other words, the subspace spanned by is enlarged by increasing *n*, which corresponds to additionally including realizations of the fields further in the past. Enlarging the subspace used to approximate the eddy forcing by field realizations further in the past is a form of delay embedding (Takens 1981). Moreover, it is motivated by the Mori–Zwanzig formalism, which demonstrates that the representation of unresolved physics includes (the estimation of) a memory term that involves the past history of the resolved physics (Wouters and Lucarini 2013; Gottwald et al. 2017). The possible relevance of the flow history for ocean eddy parameterizations has also been pointed out recently by Bachman et al. (2018) in the context of a non-Newtonian fluid mechanics approach to eddy parameterization.

More precisely, in the following sections we investigate the convergence behavior of the following dynamical spatial mode expansion of the eddy forcing:

In this study we consider and h. For completeness we also include the large-scale Jacobian, , in the expansion [see Eq. (52)]. For each choice of *n* and the expansion in Eq. (56) represents a parameterization of the eddy forcing .

The expansion coefficients in Eq. (56) are computed at each model time step by using ordinary least squares with respect to . Analyzing the dynamical and statistical behavior of the expansion coefficients as well as proposing a (possibly stochastic) model for the expansion coefficients in order to build a fully self-consistent closure is a topic for future research (see discussion section 5). Here the aim is to investigate how well the expansion in Eq. (56) approximates (converges to) .

Finally, we contrast the convergence behavior of Eq. (56) in two ways. We first consider the similar dynamical spatial mode expansion:

where the GM term [see Eq. (42)] is additionally included in the expansion. In this way we investigate the impact of the GM field on the convergence behavior.

In addition, we also analyze the convergence behavior of the spatial filter modes as given^{10} in section 2b. That is, we consider the same expansions as in Eqs. (56) and (57) but instead of using the energy modes we use the filter modes (ordered by decreasing eigenvalue/wavenumber). The expansions read

Again, the expansion coefficients are computed at each model time step by using ordinary least squares with respect to .

#### 3) Approximation of the eddy forcing on the reference attractor

In the following we analyze the approximation of the eddy forcing by the different series expansions defined in the previous section [see Eqs. (56)–(59)]. In this section the terms in Eqs. (56)–(59) are diagnosed from the reference simulation (described in section 2a). That is, the replacement (see section 3a) is not applied in the simulation and, hence, the state vector is always on the attractor of the reference ER model.

##### (i) Relative error of eddy forcing

Figure 7a shows the time-mean relative error of the eddy forcing for the filter mode expansion either with GM term [black, Eq. (59)] or without GM term [red, Eq. (58)]. The two curves are almost identical, which demonstrates that the GM term is not able to significantly reduce the relative error of the eddy forcing. In other words, the GM field (as a direction in phase space) is largely orthogonal to the eddy forcing field. For both curves the decrease in relative error (i.e., the slope of the curve) is minimal at the beginning and monotonically increasing with increasing number of filter modes. The value of the relative error for the filter modes is on the order of 10^{−1} and only reaches a very small value [i.e., *O*(10^{−13})] when all filter modes are used. That is, the convergence of Eqs. (58) and (59) is slow.

Figure 7b shows the time-mean relative error of the eddy forcing for the energy mode expansion Eq. (57) with Δ*τ* = 3 h (blue), Δ*τ* = 6 h (black), and Δ*τ* = 12 h (magenta). Note the logarithmic scale on the ordinate. The effect of the GM term on the relative error is again very small such that the curves related to Eq. (56) are indistinguishable from the shown curves. In contrast to the filter modes, the decrease in relative error (i.e., the slope of the curve) is maximal at the beginning and monotonically decreasing with increasing number of energy modes (for comparison the curve of the filter modes is shown by the blue dashed line). With only four energy modes used in Eq. (56) the relative error drops to *O*(10^{−5}) and for Δ*τ* = 3 h a relative error of *O*(10^{−12}) is reached with 30 energy modes. That is, the convergence of Eqs. (56) and (57) is very fast since adding energy modes reduces the order of magnitude of the relative error. Finally, it holds for the reference simulation that the smaller is, the smaller the relative error.

##### (ii) GM diffusivity

Figure 8a shows the time-mean GM diffusivity *K*_{GM} for the filter mode expansion with GM term [blue, Eq. (59)]. The GM diffusivity *K*_{GM} decreases nearly linearly due to the subsequent inclusion of more and more filter modes. However, the value of *K*_{GM} remains on the order of 100 m^{2} s^{−1}. Only when almost all filter modes are included the value of *K*_{GM} becomes small and, hence, the impact of the GM term is insignificant.

Figure 8b shows the time-mean GM diffusivity *K*_{GM} for the energy mode expansion with GM term [Eq. (57)] with Δ*τ* = 3 h (blue), Δ*τ* = 6 h (black), and Δ*τ* = 12 h (magenta). The behavior of *K*_{GM} resembles the behavior of the relative error (Fig. 7b). Note again the logarithmic scale on the ordinate. Including energy modes drastically reduces the value *K*_{GM}, that is, by orders of magnitude. With only four energy modes used in Eq. (57) the value of *K*_{GM} drops to m^{2} s^{−1}, indicating that the GM term is essentially without impact. Finally, it holds for the reference simulation that smaller correlates with smaller *K*_{GM}.

#### 4) Approximation of the eddy forcing in the presence of error perturbations

In this section we analyze the approximation of the eddy forcing by the different series expansions defined by Eqs. (56)–(59). At each time step the corresponding replacement is performed (see section 3a), resulting in a different simulation for each parameterization (e.g., series expansion). Consequently, error perturbations due to the approximate representation of the eddy forcing are introduced in each simulation, and hence the state vector can be pushed away from the attractor of the reference simulation. If the parameterization of the eddy forcing is accurate enough it can compensate for the error perturbations and can keep the system within or near the attractor of the reference simulation. On the other hand, if the parameterization of the eddy forcing is not accurate enough then the respective model will exhibit a different attractor.

##### (i) Relative error of eddy forcing

Figure 7a shows the time-mean relative error of the eddy forcing for the filter mode expansion either with the GM term [blue, Eq. (59); see also Table 1] or without the GM term [magenta, Eq. (58)]. The relative error for the simulations with error perturbations is slightly smaller than for the reference simulation (black and red curves). Nevertheless, the overall behavior is very similar to the reference simulation: The blue and magenta curves are nearly identical, which indicates that the GM term is not able to significantly reduce the relative error of the eddy forcing. For both curves the decrease in relative error (i.e., the slope of the curve) is minimal at the beginning and monotonically increasing with increasing number of filter modes. The convergence of Eqs. (58) and (59) is slow since the value of the relative error is on the order of 10^{−1} and only reaches a very small value [i.e., *O*(10^{−13})] when all filter modes are used.

A crucial point in Fig. 7a is that the filter mode expansion without the GM term [magenta, Eq. (58)] leads to a model blow-up if fewer than 10 filter modes are used (the magenta curve only starts when the number of modes = 10). On the other hand, the filter mode expansion with GM term [blue; Eq. (59)] leads to stable model simulations for any number of filter modes (see also Table 1). Hence, the effect of the GM term becomes clearer: the GM term cannot not significantly reduce the relative error of the eddy forcing but it can stabilize the model. In dynamical systems terms the GM term acts as a stabilizing direction in phase space. That is, the GM term cannot direct the system’s state along the attractor (it cannot excite the intrinsic low-frequency variability transitions in phase space as done by unstable directions) but it mainly keeps the system from diverging.

Figure 7b shows the time-mean relative error of the eddy forcing for the energy mode expansion Eq. (57) with Δ*τ* = 3 h (red; see also Table 1), Δ*τ* = 6 h (green), and Δ*τ* = 12 h (cyan). Note the logarithmic scale on the ordinate. The effect of the GM term on the relative error is again very small such that the curves related to Eq. (56) are indistinguishable from the shown curves. On the other hand, the stabilizing effect of the GM term also appears for the energy modes: for the application of Eq. (56) (i.e., energy mode expansion without GM term) with only two energy modes the model blows up whereas for the application of Eq. (57) (i.e., energy mode expansion with GM term) the model is stable.

The overall behavior of the relative error for the simulations with error perturbations (red, green, cyan) is similar to the results of the reference simulation (blue, black, magenta). That is, the relative error decreases much faster (adding energy modes reduces the order of magnitude of the relative error) than for the filter modes (shown for comparison by the blue dashed line). However, because of the induced error perturbations the decrease in relative error is weaker than for the reference simulation. For example, for 30 energy modes and Δ*τ* = 3 h the relative error is *O*(10^{−5}) instead of *O*(10^{−12}) for the reference simulation. Moreover, the impact of is more complicated than for the reference simulation. Roughly speaking, if fewer than 20 energy modes are used in Eq. (56) or (57) then the relative error is slightly smaller for larger whereas if more than 20 energy modes are used then the situation of the reference situation is reencountered (i.e., the smaller the smaller the relative error).

##### (ii) GM diffusivity

Figure 8a shows the time-mean GM diffusivity *K*_{GM} for the filter mode expansion with GM term [red, Eq. (59); see also Table 1]. The behavior is largely similar to the results of the reference simulation (blue), namely, the GM diffusivity *K*_{GM} decreases nearly linearly due to the subsequent inclusion of more and more filter modes. However, the value of *K*_{GM} is significantly larger (about one order of magnitude) than when diagnosed from the reference simulation. This in accordance with the interpretation of the GM term as a stabilizing direction in phase space because in the presence of error perturbations (driving the system away from the attractor) the eddy forcing will project more on stable directions (driving the system back to the attractor). In other words, in the presence of error perturbations the GM term has work to do.

Figure 8b shows the time-mean GM diffusivity *K*_{GM} for the energy mode expansion with GM term [Eq. (57)] with Δ*τ* = 3 h (red; see also Table 1), Δ*τ* = 6 h (green), and Δ*τ* = 12 h (cyan). The behavior is largely similar to the results of the reference simulation (blue, black, magenta): including energy modes drastically reduces the value *K*_{GM}, that is, by orders of magnitude. It also largely holds that the smaller , the smaller *K*_{GM}. On the other hand, the value of *K*_{GM} is larger than when diagnosed from the reference simulation (i.e., the stabilizing direction projects on the error perturbations). Nevertheless, the value of *K*_{GM} is still significantly smaller [ for only four energy modes] compared to the values typically used in ocean models [].

##### (iii) Time series of potential energy

Figure 9 shows time series of PE related to simulations employing the filter mode expansion with GM term [red, Eq. (59); see also Table 1] and the energy mode expansion with GM term and Δ*τ* = 3 h [blue, Eq. (57); see also Table 1]. For comparison the time series of the PE of the reference simulation is shown in black.

The energy mode expansion exhibits monotonic and fast convergence behavior in terms of PE (i.e., low-frequency variability). If only two energy modes are used (Fig. 9d) the PE variability is still significantly different from the reference PE. Intense low-frequency variability is present but it is situated between the low-PE regime of the reference simulation and another very-low-PE regime. Already with four energy modes in the expansion the high-PE regime of the reference simulation is regularly reached (not shown). But the low-PE regime is still bit lower than for the reference case. For six or more energy modes (Figs. 9f,h,j) the PE variability of the reference simulation appears to be essentially recovered.

As expected, the situation is different for the filter mode expansions. The convergence behavior is nonmonotonic. Even for 20 filter modes (Fig. 9i) the PE variability is significantly different from the reference simulation. When using filter modes it appears to be difficult to reach the high-PE regime of the reference simulation. Either the PE variance is significantly smaller than for the reference simulation (Figs. 9c,g) or the low-PE regime is lower than for the reference simulation (Figs. 9e,i). This is also visible in Table 1.

## 4. Summary

The three key points of this study can be summarized as follows: First, we propose a new approach to parameterizing subgrid-scale processes. In this approach the impact of the unresolved dynamics on the resolved dynamics (i.e., the eddy forcing) is represented by a series expansion in dynamical spatial modes stemming from the energy budget of the resolved dynamics. More precisely, the so-called energy modes are directly obtained from the equations of motion of the resolved flow by identifying the integral kernels that lead to the different reservoir, generation, dissipation, and conversion terms in the large-scale energy budget. Hence, the energy modes exhibit strictly large-scale patterns and they are equipped with a clear physical interpretation in terms of energetics. Convergence toward the eddy forcing is accomplished via delay embedding by including additional realizations of these fields further in the past. We also note the relation to the Mori–Zwanzig formalism, which indicates that the representation of unresolved physics needs to include a memory term that involves the past history of the resolved physics. For the two-layer QG ocean model considered in this study, we demonstrate that the convergence of a series expansion in the energy modes is orders of magnitude faster than the convergence of a series expansion in Fourier-type modes. That is, the eddy forcing can be accurately approximated with a very limited number of energy modes, which enables a feasible parameterization.

Second, we explore a novel way to test parameterizations in models. The resolved dynamics and the corresponding instantaneous eddy forcing are defined via spatial filtering, which accounts for the representation error of the equations of motion on the low-resolution model grid. In this way closures can be tested within the high-resolution model. Whereas in low-resolution models all energy pathways between large-scale and eddy components must be parameterized simultaneously, testing parameterizations in the high-resolution model offers the possibility to isolate the effects of a single parameterization (related to a single energy pathway) while the other large-scale eddy energy conversions are correctly computed. For the two-layer QG ocean model considered in this study, we focus on parameterizations of the baroclinic energy pathway while the barotropic energy pathway is correctly computed by the high-resolution model.

Third, we test the standard closure of the baroclinic energy pathway in the ocean components of state-of-the-art climate models [i.e., the Gent–McWilliams (GM) parameterization with a scalar diffusivity] in the high-resolution QG ocean model considered in this study. It turns out that the GM field steers trajectories along a stabilizing direction in phase space. That is, the GM field does not project well on the eddy forcing (it exhibits a very high relative error) and fails to excite the model’s intrinsic low-frequency variability (i.e., it is not able to propagate the model’s state along the correct attractor e.g., along an unstable direction). The GM field mainly stabilizes the model. That is, if the representation of the eddy forcing is very inaccurate (e.g., small number of modes used in expansion) the GM term performs the necessary dissipation of available potential energy such that the model does not diverge.

## 5. Discussion

Finally, we elaborate on open issues of this study and related future research directions.

### a. Self-consistent closure of the baroclinic energy pathway

A closure of the baroclinic energy pathway is self-consistent if it does not involve the actual (“true”) baroclinic eddy forcing. However, in this study we still use for the computation of expansion coefficients (i.e., the coefficients that appear in a spatial mode expansion) via ordinary least squares. Determining a self-consistent closure of the baroclinic energy pathway is related to three intricate and intimately related issues: 1) determining the optimal subset of energy fields [see Eq. (54)], 2) diagnosing the corresponding expansion coefficients, and 3) proposing a (possibly stochastic) self-consistent model for the expansion coefficients. The choices made with respect to these three issues can have an effect on the accuracy of the approximation (as indicated in this study by the different choices for ), the computational cost and complexity of the model, the regularity of the expansion coefficients, and the uniqueness and hence physical interpretation of the series expansion in energy modes.

For example, a problem related to these issues and well-known in statistics and machine learning is the issue of overfitting versus underfitting or the bias/variance trade-off. Low-bias approaches can usually give accurate representation of the data but produce large variances. In contrast, models with higher bias produce lower variances but less accurate representations. Regularization methods introduce bias into the regression solution that can reduce variance considerably. In this way the behavior of the expansion coefficients becomes simpler and easier to model but the approximation becomes less accurate.

### b. Self-consistent closure of the barotropic energy pathway

To make the equations for the large-scale flow [Eqs. (13) and (14)] completely self-consistent one also has to specify a self-consistent closure of the barotropic energy pathway (i.e., ). The standard closure of the barotropic energy pathway is lateral viscous dissipation with an enhanced “eddy” viscosity coefficient. Similar to the GM parameterization the lateral viscosity parameterization suffers from the lack of backscatter (see Fig. 2). But an adequate closure of the energy exchange between large-scale and eddy components is necessary in order to be able to perform low-resolution model simulations exhibiting eddy-driven low-frequency variability. One option is to proceed in a way similar to this study: exploring whether spatial fields that stem from the large-scale kinetic energy budget can suit as dynamical modes to parameterize the eddy forcing .

### c. Dynamical systems analysis of the large-scale flow in the turbulent regime

As soon as adequate closures for both the baroclinic and the barotropic energy pathways are available it is in principle possible (i.e., feasible due to low model resolution) to analyze the dynamics of the large-scale flow in the turbulent regime in a systematic way. In the case of deterministic closures this is related to the existence of multiple equilibria, stability properties, bifurcations, and chaotic attractors (Dijkstra 2005). In the case of stochastic closures the investigation will be from the perspective of random dynamical systems, which is related to stochastic bifurcations (i.e., changes in the probability density function), pullback attractors, and invariant measures (Dijkstra 2013). We note that in case of low model resolutions a whole set of numerical techniques to investigate transitions in stochastic dynamical systems becomes feasible (Dijkstra et al. 2016). For example, it becomes possible to numerically solve the stochastic partial differential equations (SPDEs) via dynamical mode expansions (Sapsis and Lermusiaux 2009) and to investigate the interaction of external noise forcing with internal nonlinear variability in the turbulent regime (Sapsis and Dijkstra 2013).

### d. Comparison with other approaches to eddy parameterization

In this study, we compared turbulence closures based on energy modes with the GM eddy parameterization approach. We focused on a positive and spatially constant *K*_{GM} because it is straightforward to diagnose (e.g., not entering issues around rotational eddy fluxes), and, more importantly, because a spatially homogenous *K*_{GM} is still regularly applied in state-of-the-art realistic ocean models. On the other hand, the estimation and performance of a spatially inhomogeneous (and possibly tensor-valued) *K*_{GM} remains a crucial topic (Eden et al. 2007, 2009; Viebahn and Eden 2010). The relation between energy modes and the GM parameterization, as well as other approaches to eddy parameterization (Porta Mana and Zanna 2014; Jansen and Held 2014; Bachman et al. 2018), will hopefully be further elucidated in future studies.

### e. More realistic ocean model configurations

The ocean model considered in this study is situated at the more idealized end in the hierarchy of ocean models. Several features and processes must be included in order to make the details more realistic. These include higher vertical resolution, diabatic terms like buoyancy forcing and buoyancy sinks, and realistic topography and coastlines. We are currently extending our results to a three-layer model including realistic topographic interactions. Eventually, one also has to consider the primitive equations in order to be able to investigate global realistic ocean models. The corresponding energy budgets are more complicated but detailed analyses are becoming available nowadays (von Storch et al. 2012; Wu et al. 2017; Jüling et al. 2018).

### f. Climate model simulations subject to intrinsic (eddy driven) low-frequency variability

Finally, when adequate closures for the energy pathways in realistic ocean models are available then long-period low-resolution climate model simulations exhibiting eddy-driven low-frequency variability become possible. This is crucial since then issues related to anthropogenic climate change (forced variability) versus intrinsic low-frequency variability (internal variability) can be addressed in a statistically significant manner.

## Acknowledgments

The authors thank Inti Pelupessy, Alexis Tantet, and Fred Wubs for helpful discussions. This research is funded by the Netherlands Organization for Scientific Research (NWO) through the Vidi project “Stochastic models for unresolved scales in geophysical flows.” HD acknowledges support by the Netherlands Earth System Science Centre (NESSC), financially supported by the Ministry of Education, Culture and Science (OCW), Grant 024.002.001.

## REFERENCES

*Nonlinear Physical Oceanography*. Springer, 532 pp.

*Nonlinear Climate Dynamics*. Cambridge University Press, 367 pp.

*Nonlinear and Stochastic Climate Dynamics*, C. Franzke and T. J. O’Kane, Eds., 209–240.

*Climate Change 2013: The Physical Science Basis.*T. F. Stocker et al., Eds., Cambridge University Press, 1535 pp.

*Ocean Dynamics*. Springer, 704 pp.

*Ocean Circulation Theory*. Springer, 455 pp.

*Turbulent Flows*. Cambridge University Press, 771 pp.

*Physics of Negative Viscosity Phenomena*. McGraw-Hill, 256 pp.

*Dynamical Systems and Turbulence*, D. Rand and L.-S. Young, Eds., Lecture Notes in Mathematics Series, Vol. 898, 366–381.

*Atmospheric and Oceanic Fluid Dynamics*. Cambridge University Press, 745 pp.

*Statistical Analysis in Climate Research*. Cambridge University Press, 499 pp.

## Footnotes

© 2019 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).

^{1}

Identity with respect to time evolution is meant in a statistical/dynamical sense.

^{2}

Note that this criterion is expressed here with respect to the equations in physical space. With respect to the equations in modal/wavenumber space it says that the constant interaction coefficients (obtained by computing the amplitude equations of the individual modes) related to the resolved large-scale modes should not significantly change whether computed from the high-resolution or low-resolution representation of the modes.

^{3}

Orthonormalized in the streamfunction norm (equivalent to PE norm), .

^{4}

Note that there is a subtlety here: We assume that the terms , , , and do not project on modes that lie outside the subspace defined by the filter cutoff . Of course, every model discretized and stepped forward in physical space (and not directly in modal/wavenumber space) suffers from the fact that energy can be transferred to small-scale modes that cannot be adequately represented on the spatial grid (leading, e.g., to aliasing). However, since we use the eigenmodes of the frictional term (which typically represents the most small-scale patterns) we expect to essentially remain within the subspace spanned by the large-scale modes (defined via the filter cutoff ).

^{5}

In our model setup the two conversion terms are identical in the time mean (see Fig. 5) due to the absence of buoyancy sources/sinks.

^{6}

Note that since we use a sharp filter.

^{7}

Note that this estimation is not affected by rotational eddy fluxes since it is not computed on the level of fluxes [like Eq. (37)] but on the level of dynamical terms appearing in the PV budget.

^{8}

Note that the model blows up if is simply set to zero, consistent with the fact that acts as a sink of time-mean [PE] (see Fig. 5).

^{9}

Dynamical modes are time-dependent and budget-based in the sense that their computation explicitly involves the governing conservation equations (see, e.g., Dijkstra 2016). In contrast, for example, statistical modes (e.g., EOFs) are data based and not budget based.

^{10}

Loosely speaking, these are Fourier-type modes. More precisely, the spatial filter modes are eigenmodes of the bi-Laplacian in this study.