On Energy Cascades in General Flows: A Lagrangian Application

An important characteristic of geophysically turbulent flows is the transfer of energy between scales. Balanced flows pass energy from smaller to larger scales as part of the well‐known upscale cascade, while submesoscale and smaller scale flows can transfer energy eventually to smaller, dissipative scales. Much effort has been put into quantifying these transfers, but a complicating factor in realistic settings is that the underlying flows are often strongly spatially heterogeneous and anisotropic. Furthermore, the flows may be embedded in irregularly shaped domains that can be multiply connected. As a result, straightforward approaches like computing Fourier spatial spectra of nonlinear terms suffer from a number of conceptual issues. In this paper, we develop a method to compute cross‐scale energy transfers in general settings, allowing for arbitrary flow structure, anisotropy, and inhomogeneity. We employ Green's function approach to the kinetic energy equation to relate kinetic energy at a point to its Lagrangian history. A spatial filtering of the resulting equation naturally decomposes kinetic energy into length‐scale‐dependent contributions and describes how the transfer of energy between those scalestakes place. The method is applied to a doubly periodic simulation of vortex merger, resulting in the demonstration of the expected upscale energy cascade. Somewhat novel results are that the energy transfers are dominated by pressure work, rather than kinetic energy exchange, and dissipation is a noticeable influence on the larger scale energy budgets. We also describe, but do not employ here, a technique for developing filters to use in complex domains.


Introduction
A surprising aspect of nonlinear geophysical flows, namely, that they transfer energy from smaller to larger scales and thus in a sense opposite to that of three-dimensional, isotropic turbulence, has been known for some time (Charney, 1971;Fjortoft, 1953;Kraichnan, 1967). The reason for this behavior has been traced to the quasi-two dimensionality of so-called balanced dynamics, which strongly constrains vortex tube stretching. Scott and Wang (2005) argue for the observation of an upscale cascade in the ocean from analysis of satellite sea surface height measurements, and Scott and Arbic (2007) argue for the same in numerical models. Others have argued that this behavior in wavenumber space is accompanied by a transfer in the frequency domain to lower frequencies (Arbic et al., 2012) and extensions to the combined wavenumber frequency domain have also been found (Arbic et al., 2014).
An upscale cascade dynamically and importantly affects the ocean. For example, the transfer from high frequency mesoscale variability to low frequencies implicates eddies as a mechanism effecting climate variability (Berloff & McWilliams, 1999;Dijkstra & Molemaker, 1999;Holland, 1978). Further, the large-scale circulation receives energy from the large-scale, slowly evolving winds (Ferrari & Wunsch, 2009;Wunsch & Ferrari, 2004), and it is generally thought that this energy is shunted to the mesoscale by geostrophic instabilities. An equilibrated ocean must then find a pathway for the mesoscale to lose its energy, and the so-called upscale cascade presents a hindrance. The return of energy to larger scales where dissipation is vanishingly weak has led to considerable recent interest in mesoscale energy loss.
Given that cross-scale energy exchange in the ocean is central to ocean and climate dynamics, it is important to examine how those transfers can be measured. Most energy exchange diagnoses of both observations and models have proceeded through the use of Fourier transforms (Arbic et al., 2014), which requires that the data be converted into forms consistent with the approach. Speaking primarily about spatial analyses, geographically confined regions of interest are usually defined. Since the data are not spatially periodic, as assumed by Fourier methods, tapering or windowing is usually applied. The subsequent analysis then occurs on the windowed domain, and the results are in some sense representative of an averaged statement of the transfers in the domain. Another complicating factor is that irregularities can appear inside the domain. An example is Bermuda, sitting in the open Atlantic at a location near to the Gulf Stream and the separated North Atlantic jet. Methods, often subjective, are required to handle such areas, frequently involving an interpolation over the problematic regions. Nonetheless, these analyses have usually yielded answers consistent with our theoretical expectations. They do represent spatial averages however, and further understanding of inhomogeneous regions, like the Kuroshio and the Gulf Stream, requires more spatial resolution. This served in part as motivation for Aluie et al. (2018), who examined energy exchanges locally based on a purely spatial filtering approach.
The focus of this paper is to propose a new method for analyzing kinetic energy in models and observations that appears to avoid many of the above-mentioned issues with Fourier methods. It also differs from previous approaches in several ways, primarily in focussing on a Lagrangian-like approach and working directly with the kinetic energy equation. We apply the technique to primitive equation model output describing vortex merger and hence an upscale cascade. The results indicate an upscale cascade, thus demonstrating the viability of the procedure, even though it is not limited to this restrictive setting. The theoretical basis for the procedure is outlined in the next section; the model description and applications appear in the following section. We discuss and summarize the paper in the concluding section, emphasizing new results specific to the energetics of the upscale cascade and outline thoughts for further work. Methods to handle domain irregularities are outlined in Appendix A. An Eulerian version of the present analysis appears in Appendix B. Derivation of comparable formulas for potential energy are discussed in Appendix C, and modifications to the analysis introduced when using isopycnal coordinates are outlined in Appendix D.

Overview
Below, we present an analysis of the kinetic energy equation. In brief, our plan is to apply Green's function analysis to the kinetic energy equation. The result relates the final kinetic energy of a parcel to its initial kinetic energy and its history along its (quasi-) Lagrangian trajectory. We then filter this equation in physical space to obtain a statement about a length-scale-based decomposition of the final kinetic energy. Similar filterings of the initial kinetic energy and kinetic energy sources then yield statements about cross-scale energy transfers. Finally, we average the pointwise kinetic energy decomposition to obtain a more representative view of the cross-scale transfers. These steps are spelled out in sequence in the following subsections. The reader wishing to skip mathematical detail can go directly to Equations 36-38, where the results are summarized.
We begin with the hydrostatic equations in Cartesian coordinates, the zonal momentum component of which is where u is zonal velocity and subscripts x, y, z, and t denote partial derivatives with respect to zonal, meridional, and vertical location and time, respectively. The quantity p is pressure, and is viscosity. A partial step to constructing the kinetic energy equation involves multiplying Equation 1 by u, which results in A similar procedure on the meridional momentum equation and adding returns the kinetic energy equation where K = u 2 ∕2 + v 2 ∕2 = |u h | 2 ∕2 and = ∇u h ⋅ ∇u h = |∇u h | 2 is dissipation. We rewrite this as

10.1029/2020MS002090
where R = −u h ⋅ ∇ h p − contains the "sources" and "sinks" for kinetic energy. Given that pressure gradients appear in the momentum equations as a force acting on a fluid particle, the quantity −u h ⋅ ∇ h p represents movement in the direction of a force. We will thus refer to it as "pressure work." The quantity − is simply dissipation of kinetic energy.
Exploiting the hydrostatic balance, pressure work can be rewritten as where b is buoyancy. The last quantity in Equation 5 is the exchange between kinetic and potential energy.
In this paper, we work primarily with pressure work, noting that potential energy exchange is a component of pressure work.

Comparison to Other Approaches
Working directly with kinetic energy, via Equation 4, distinguishes this approach from others. Traditional Fourier analysis proceeds by transforming Equation 1 and multiplying by the momentum transform conjugate û(k) * where k denotes a wavenumber. The result (Ê(k) = (ûû * +vv * )∕2) is usually interpreted as the energy spectrum, although more properly it is related to the momentum spectrum. The relation of Ê(k) to energy comes indirectly from Parseval's equality, i.e., which, while useful, argues this approach is restricted to integral statements. In contrast, the Fourier transform of K isK which, although quite different from Ê(k), still satisfies Parseval's identity while being directly related to kinetic energy.
To probe into the spatial structure of energy transfers, it is natural to remain in the physical domain. In this sense, our approach is like that of Aluie et al. (2018). We differ from them in that they "coarse-grain" Equation 1 and multiply by the coarse-grained velocity to obtain an equation describingK = 0.5(ū 2 +v 2 ), where the bar denotes the coarsening operator. In the process, terms appear in theK equation involving third order quantities likeū(uu x ) −ū(ūū x ), which in the LES community are interpreted as measuring the energy cascade due to wave-mean flow interaction. They describe the effect onK, following a parcel moving atū, of subgrid-scale processes.
Clearly, the coarsening operator is one of the primary distinctions between Aluie et al. (2018) and the present approach. In fact, if we choose the coarsening operator as a Dirac delta function, thenū = u and our approaches merge. The cost of this choice is the wave-mean flow interaction term vanishes, leading to the statement that the unfiltered equations do not cascade. In our view, cascading is broader than wave-mean flow interaction and occurs in unfiltered results, as evidenced in Figure 1. The panels show results from the numerical model discussed in detail later in the paper. For now, we mention they are plots of kinetic energy in a layer and come from Days 20 and 585 of a doubly periodic numerical experiment typified by vortex merger. Both plots are in physical space and neither plot has been coarse-grained nor filtered. Visual comparison suggests that the later plot is characterized by larger length scales than the earlier. Spectra of "normalized" kinetic energy, i.e., of the variable = a(t)K, for both days, taken along the line in Figure 1 (right), appear in Figure 2. The integrated kinetic energies along this line were used to normalize the kinetic energies, i.e., a(t) = , so that the total integral in both plots, ∫ d , was unity. This was done to remove amplitude variations that would interfere with the comparisons of the wavenumber content. The earlier spectrum (red) is characterized by enhanced energies at higher wavenumbers, while the later (blue) shows elevated levels at lower wavenumbers. Clearly a "cascade" has occurred in these unfiltered plots, and it is this behavior we seek to quantify.
Although our procedure differs from that of Aluie et al. (2018), we do not criticize their approach. They adopt a procedure of clear relevance to modelers interested in parameterization. Our approach is meant for the full diagnosis of "eddy-resolving" models.

Kinetic Energy at a Single Point
We now apply Green's function analysis, which leads eventually to the construction of kinetic energy at a point and begins by multiplying Equation 4 by some space-time-dependent function G.
Straightforward calculus leads to Again, appealing to straightforward Green's function methods, we suppose G satisfies where is the usual Dirac delta function centered on the so-called source points x o , t o and̂is a viscosity-like quantity. It is not necessarily equal to , although in this paper, we will always discuss cases wherê= . Boundary conditions are required to solve Equation 10. As is typical for Green's functions, the boundary conditions will always be homogeneous; however, how they will be distributed between specified values of G and its normal derivatives on the boundary is yet to be determined.
Equation 10 is almost a simple advection-diffusion equation, differing only because viscosity appears multiplied by −1, i.e., −̂< 0, as opposed tô> 0 occurs on the right-hand side. In fact, such a result is standardin Green's function analyses, where the equation solved by G is the adjoint of the equation governing K, and is handled by solving Equation 10 backwards in time. Instead of an initial condition, a "posterior" condition on G is applied, and integration proceeds to temporally earlier values. For this case, we will always choose a posterior condition of G(t ) = 0. Notationally, G = G(x, t; x o , t o ), expressing the dependence of G on both the "source" points, x o , t o , and the "observation" points, x, t.
Now, we perform an integral of Equation 9 over a volume V of the observation space, x, and over the time where S is the surface bounding the volume V. At this point, we exploit the available flexibility in boundary conditions. Equation 11 describes how various kinetic energy sources contribute to the final kinetic energy at a point. From a physical perspective, we expect that kinetic energy transport across the boundary S can contribute to the final kinetic energy inside S. This can occur by kinetic energy advection and viscous exchange across S. Both of these appear in Equation 11 multiplied by G. We therefore choose to apply ∇G ⋅ ⃗ n = 0, i.e., homogenous Neumann boundary conditions, to G, thus eliminating the underlined term from Equation 11 and preserving the aforementioned physical effects.
Suppose first V is a closed domain. With this, and a rearrangement of terms, we obtain as the boundary fluxes all vanish.
Equation 12 illustrates how the kinetic energy at location x o and time t o has been constructed. The first term on the right-hand side represents the "initial" kinetic energy, K(x, 0) that ends up at (x o , t o ). The volume defining the impact of this K(x, 0) contribution depends on the development of G from t o . The last term involves the internal fluid sources of kinetic energy and how they have contributed to K along the trajectory that eventually ends up at x o . Equation 12 provides a full deconstruction of the kinetic energy in the basin, regardless of spatial inhomogeneities in the flow.
If V is an open domain or a subdomain of a full, closed basin, Equation 12 is augmented by advective and diffusive fluxes of kinetic energy across the subdomain boundary, to the extent they are involved as measured by G. These vanish identically in the presence of doubly periodic conditions, which will describe our application in this manuscript. But even for the most general case, boundary contributions can be made small by choosing a domain large enough that advection and diffusion do not have time to reach the boundaries. We will proceed here ignoring boundary contributions.
The clearest interpretation of this formula occurs if G, defined in Equation 10, measures purely Lagrangian trajectories, which occurs in the limit̂= 0. For the remainder of this paper, we will continue witĥ= both for simplicity and because our numerical application will employ a very small value for . We have also checked that our results are insensitive tô< . It is in the nearly Lagrangian nature of Equation 12 that the present approach somewhat resembles that of Nagai et al. (2015) and Shakespeare and Hogg (2017). K. Srinivasan has also alerted us to early contributions by Kraichnan (Kraichnan, 1964(Kraichnan, , 1967 in which Green's function-based analysis of turbulence was employed. There are similarities in our approaches, but the objectives our studies differ. Kelley et al. (2013) also explored the coarse-graining decomposition in Aluie et al. (2018) in a Lagrangian setting with a focus on Lagrangian Coherent Structures.

Cross-Scale Contributions to Kinetic Energy
We now compute how energy is transferred between various length scales. These exchanges become more complicated when K(x o , t o ) is strongly inhomogeneous and anisotropic, as is typical in the ocean. In this and the next subsections, we will build to a formal statement of cross-scale transfers in several steps, proceeding from a statement valid at a single point to one averaged over a specified volume.
As a first step, consider a "filter" applied to K(x o , t o ), i.e., where f L is a filter associated with a length scale, L, and fixed at location x 1 . What the filter might be is a question deserving discussion, because the filtering will be carried out in physical space. For realistic problems, physical space typically contains a lot of irregularity: examples are complex lateral boundaries and major topography that can extend through domains at given depths and reach up to the surface. Any useful filter must contend with these features. Guidance about the filters can be obtained from the simplest case of an unbounded domain without interior obstacles. In that case, as we argue below, the obvious choices are trigonometric functions. At a more fundamental level, these work because they are a collection of functions forming a complete set due to their being the eigenfunction, eigenvalue solutions to the Helmholtz equation in an unbounded domain. We argue in Appendix A basis functions constructed from the Helmholtz equation with homogeneous conditions on S, the boundary of the domain V, will for complex domains provide a collection of functions with the desired properties. We stress here, however, that Appendix A is not exploited in the present paper but serves as a procedure in need of further examination.
But whatever the filter is, an important point is that the filter works on the "source points" (x o ) of Green's function. Operating on Equation 12, we obtain where We can similarly filter Equation 10 to show We want the filtering to decompose K(x 1 , t o ) completely into scale-dependent information, i.e.,

Employing Equation 13 in Equation 17
, it is simple to show where is the Dirac delta function (This property of the basis set is essential to their role in filtering. We will revisit it in Appendix A when developing more general basis sets.). A classic result from the theory of generalized functions is where the index i denotes dimensions and Π implies multiplication. Working now with one of the dimensions, suppressing the subscript i, breaking the integral into real and imaginary parts, and using symmetry implies 10.1029/2020MS002090 Last, we employ the relation between wavenumber and wavelength, k = 1∕L, to obtain which effectively defines our original filter, f L , for a single spatial dimension To filter over a group of length scales L x , L y , and L z in three dimensions, one employs the product of Equation 22 appropriate to each dimension. The quantity K L (x 1 , t o ) represents the energy resident in K(x 1 , t o ) between the length scales L and L + L.

Band Pass Filtering of Final Kinetic Energy
The second step in formulating cross-scale transfers is to represent the energy at x 1 resident between scales of L j and L k by suitable integration over length scales. For example, the energy at scales L j and greater becomes (again working in one dimension) is the well-known "sinc" function, often found in signal detection applications. The remainder of the energy at shorter lengths is Analogous to Equation 15, we define and find that it is governed by This more general formulation anticipates that the length-scale-dependent basis functions needed to probe the scale-dependent structure of K in general domains will not typically be trigonometric functions. See Appendix A for more discussion. We will frequently refer to L k L in the following as a "filter." The construction of the kinetic energy at x 1 between L j and L k can be written Note that Equation 29 is a Lagrangian statement of kinetic energy exchange, in that the modified Green's function Γ to use this formalism to render an Eulerian statement about kinetic energy exchange. This is discussed in Appendix B. We continue here with the Lagrangian formalism, such having been less discussed in the literature (but see Kelley et al., 2013).

Band Pass Filtering of Kinetic Energy Sources
The third step addresses the cross-scale exchanges involved in the creation of L k L . Consider the "initial" kinetic energy, appearing as the first term on the right-hand side of Equation 29. Exploiting Equation 18, this quantity can be rewritten where the interval L = 0 to L = ∞ has been broken into N + 1 separate intervals and, equating L j to L i and Equation 31 expresses how the scale-dependent decomposition of the initial kinetic energy contributes to the kinetic energy between L j and L k found in the final state. To the extent that the quantity does not vanish for L i , L i + 1 outside the range L j to L k , there has been cross-scale transfer of kinetic energy. This calculation can be repeated for each band L i to L i + 1 to obtain a complete breakdown of the cross-scale transfers from the initial data to the final. Further, this computation can be performed for any time between the initial and final times, thus providing a time series of cross-scale kinetic energy exchanges. Similar decompositions into scale-dependent contributions can be formed from the boundary terms and the source terms in R.
These lead to a complete reconstruction of the kinetic energy at any given scale from large-and small-scale structures occurring along the quasi-Lagrangian fluid path which, in turn, are measures of energy transfers across scales. Examining the time series of these terms identifies critical points in time when transfers are at their greatest. This information, along with knowledge of the spatial Green's functions structures at those times, identifies the particular events associated with the transfers. A relatively complete overall picture of the energetics of the fluid can be developed.

Area Averaging
At this point, cross-scale transfers have been formally identified, with perhaps the strongest restriction being that they apply to the development of kinetic energy at a given point. However, it might prove useful to have such a statement averaged over a given area. For example, for the case of merging vortices, one might be interested in an averaged statement of the construction of the energy averaged over the area of the merged vortex. We will pursue this goal in the next section. For now, we note that this once again reduces to an operation on the source points of Green's functions and can be obtained following the procedures described above.
Suppose that we want a simple "top-hat" averaging of the results over a given domain, such as would be obtained according to where Π Q is the top-hat function defined by

10.1029/2020MS002090
with A a square centered on x = 0 of area Q 2 . Following the steps outlined previously leads to where L n+1 L n is the analog to Equation 29 for the energy source terms and we have assumed the boundary We will mostly work with the "modified" Green's function Φ L k L in this paper and will call L k Q,L a sampling function.
Although not pursued here, a similar approach can be applied to the potential energy equation. The details appear in Appendix C. Recall kinetic energy connects to potential energy through pressure work.
To summarize, Equation 36 equates the scale-dependent content of the kinetic energy between lengths L j and L k , averaged over the A (see Equation 34), to contributions from the initial kinetic energy, pressure work, and dissipation.

Application to Merging Vortices
We now apply the previous methodology to a classical GFD problem, namely, the merger of like-signed vortices. We choose this as a first test case because this scenario is usually identified with the upscale cascade of rotating flows. In this sense, we know what the answer should be, broadly speaking, and so can judge the utility of the procedure. In addition, we can examine what new information is provided by this approach.
For this purpose, we used the recently developed MOM6 model (https://github.com/NOAA-GFDL/MOM6). This model includes a number of modern features: most importantly, it can be deployed in a purely isopycnal configuration. There are relatively minor modifications to the kinetic energy equation in this setting, all of which are detailed in Appendix D. Perhaps the most significant details are that Green's function equation is modified by the appearance of layer thickness in the viscous term, and the velocity field must be treated as compressible.
Our model consisted of a four layer f -plane fluid, with three subsurface interfaces, in a flat bottomed basin of 4,000 m of total depth. The free surface and third layer interface were initially taken as flat, and the second layer, bounded by the first and second interfaces, was seeded with an organized distribution of 200-m thick anomalies. All initial velocities were assumed to vanish. As expected, after a very short time adjustment involving gravity wave radiation, the fluid followed a classical path of like-signed vortex merger and opposite signed vortex propagation toward a final state of a single cyclonic-anticyclonic pair. Along the way, smaller vortices are swallowed by larger ones and a recognizable upscale energy cascade takes place. Examples of the fluid state at Days 50, 135, and 220 in Year 2 of the simulation appear in Figure 3 in which this sequence is evident.
Numerical details of the simulation appear in Table 1; here, we emphasize that the resulting internal deformation radii, 34, 19, and 14 km, are representative of the open ocean. The experiment was simulation for a total of 4 years, from which we focus on the analysis of Days 50-220 from the second year of the simulation. The sequence in Figure 3   Closing budgets as we are attempting to do here requires high accuracy, especially in the solution of the adjoint equation. MOM6 does not have an adjoint option, so the following procedure was adopted. Output from the merger experiment of horizontal velocity components and thickness were cataloged at intervals of 1 day, from a simulation that used time steps of 300 s. The computational domain was doubly periodic and f -plane, so we developed a spectral model for the adjoint equation in MATLAB. We adopted a comparable 300-s time step for the adjoint equation which, of course, required knowledge of the model state at comparable temporal resolution. This was obtained for any given one 1-day interval by using 6 days of MOM6 data centered on that day to interpolate via cubic splines the velocities and thickness to the required adjoint model times. The interpolated velocity and thickness data were used to compute kinetic energy and dissipation at the mass points of the model grid, and the pressure work at that point was estimated as the residual in the kinetic energy equation. Our estimates were checked against direct estimates of the pressure work which could be computed using the daily archived fields, and the two were found to agree to a high degree of accuracy. In this way, we developed time series of the various inputs to the kinetic energy equation that balanced to machine precision, which in turn allowed us to balance our adjoint-based kinetic energy budgets.
As is evident in Equation 37, the evolving modified Green's functions are related to Lagrangian measures of fluid displacement. This point is emphasized in Figure 4 which compares Green's functions corresponding to the beginning and the end of our experiment, i.e., to Days 50 (left) and 220 (right). Comparatively speaking, the left-hand side looks considerably more disorganized than the right-hand side. This reflects the strongly interactive nature of the evolving fields, which moves particles where L j = 0 km, L k = 1,536 km, and A o is the domain area. Intermediate lengths were set to (L 2 → L 13 ) = 1,200, 1,000, 800, 700, 600, 500, 450, 400, 350, 300, 250, and 150 km. The adjoint problem appropriate to each of these sampling functions was solved to obtain fourteen different sets of Φ L i+1 L i functions, and the kinetic energy, pressure work, and dissipation fields were filtered using each of the associated filtering functions. This resulted in a 14 × 14 matrix for each of the kinetic energy and energy input fields, detailing how the various scale-dependent quantities contributed to the final kinetic energy at a given time. We will refer to the results in terms of the bands associated with the initial Φ L i+1 L i functions, as spelled out in Table 2.
For the purposes of this paper, we consider the transition from Days 50 to 220, i.e., over the entire interval of the vortex merger. We compare initial and final kinetic energy decompositions as functions of length scale in Figure 5. Note that some of the values are not positive. This is a consequence of Equation 17, which insures that when summed all filtered contributions converge on the value of the energy in physical space. What is guaranteed to be positive is the final sum. Contributions from individual length-scale bands can contribute negative values to achieve this goal.
The red line is L i o of the initial kinetic energy within the square seen in Figure 3c. Although relatively flat, it does exhibit a peak in the range L = 400-600 km which reflects the two vortices present at Day 50. While indicative of the early energy distribution, a more telling decomposition comes from that found by subsampling the field according to the Day 50 Green's function (see Figure 4). This appears in the green line, which shows a weak dip in values at small length scales, followed by a steady growth. The final, blue, line is the decomposition associated with the Day 220 state inside the square in Figure 3c. It is clear that, compared to both other curves, the kinetic energy has both grown considerably in amplitude, particularly at large L, which is consistent with an upscale cascade.

Journal of Advances in Modeling Earth Systems
A more detailed view of the cross-scale exchanges is shown in Figure 6. Here, we plot the decomposition of the large-scale (Band 2), intermediate-scale (Band 8), and small-scale (Band 12) kinetic energy in the upper three panels. Each plot consists of 43 separate values. The first value is the band-passed kinetic energy at Day 220. This is followed in the next 14 slots by the initial kinetic energy contributions from the various bands. Those values are distinguished by a magenta axis. The next 14 slots are the various band-passed contributions from the pressure work (cyan axis), and the final 14 slots (blue axis) are the contributions from dissipation. Each group of 14 contributions is also separated by a vertical black line. The solid black horizontal lines in each slot are the values of the contributions, and the red line starting from the second slot shows the cumulative sum of the values. As shown in Equation 29, the sum of all the inputs should match the final band-passed kinetic energy; equivalently, the final red value at Slot 43 should match the black line in the first slot in amplitude. Note that this holds in all three plots. Finally, the three dashed green vertical lines mark the band corresponding to the band of the final kinetic energy for each of the initial kinetic energy, pressure work, and Figure 6. Three kinetic energy decompositions corresponding to Bands 2 (upper), 8 (next upper), and 12 (next lower). Vertical black lines limit the final kinetic energy in the band, the contributions to that band from the 14 bands of the initial kinetic energy, pressure work, and dissipation, respectively. Horizontal black lines denote values of the quantity, and the red line is the cumulative sum, moving to the right, of the various contributions. The abscissa denotes band number, and the ordinate is energy in m 2 /s 2 . The bottom plot connects the band number on the y-axis to the length scales within that band, appearing on the x-axis in km. dissipation. Bands left of the green line correspond to larger length scales and those to the right to shorter length scales.
The bottom panel connects the band number on the y-axis with the length-scale range of that band on the x-axis.

Large Scales
For example, consider Figure 6a. The first slot is occupied by which is Band 2 kinetic energy in the final state, with a value of roughly 5.5 × 10 −3 m 2 /s 2 . A plot of the physical space structure of this quantity appears in Figure 7d. The kinetic energy field at Day 220 from which this filtered content is extracted appears above it in Figure 7b. The next 14 slots, from Slots 2 to 15, are occupied by is the initial kinetic energy content contained in band n (see Table 2). The dashed green line marks the third slot, which is the second band of the initial kinetic energy. This is the same band appearing in the first slot, and a physical space picture of this transfer term appears in Figure 7c. Day 50 kinetic energy field from which this term is extracted by Equation 42 appears above it in Figure 7a. It is clear from comparing the black line in the green dashed slot to the black line in the first slot from Figure 6a that this band has grown in kinetic energy over the experiment, from an initial value of 1.3 × 10 −3 m 2 /s 2 to roughly 5.5 × 10 −3 m 2 /s 2 . The values in the remaining slots in the kinetic energy zone represent the contributions from shorter length scales in Figure 8. The distribution of dissipation at Day 135 appears on the left. Note that it is concentrated into small areas. Dissipation as it affects Band 2 appears on the right. Although this is a large-scale mode, dissipation is still a major input to the kinetic energy development, as seen in Figure 6.
the initial kinetic energy to Band 2 final kinetic energy structure. The next four values in the initial kinetic energy zone are positive, indicating a transfer of smaller scales upward. The remaining slots are filled with very small values. The cumulative red line indicates the total transfer to Band 2 final kinetic energy from initial kinetic energy is roughly 3.3 × 10 −3 m 2 /s 2 , with the biggest single contributor being local Band 2. In summary, this plot suggests an upscale kinetic energy to kinetic energy cascade.
The slots from 16 to 29 are occupied by is the pressure work performed in band n at time t of the calculation. The green dashed line falls in Slot 17 because that band corresponds to the band appearing in the first slot. The values in the remaining slots represent contributions to Band 2 final kinetic energy from other bands in the pressure work occurring along the fluid trajectory. As can be seen from the red lines, pressure work is considerably more active across the bands than kinetic energy. The green dashed slot contains a negative value, indicating a net loss of kinetic energy by the pressure work in this band. However the next several slots are positive, indicative of an upscale cascade of energy due to pressure work, and in total, these transfers are considerably stronger than the in-band loss. In total, pressure work elevates kinetic energy by 3.8 × 10 −3 m 2 /s 2 and is roughly 3-4 times larger than the kinetic to kinetic transfer, with most of the transfers coming from smaller scales. This result is somewhat surprising, as most previous studies have focussed on kinetic to kinetic transfers. Contributions from other sources, such as pressure work, have not received as much attention. Here, they are dominant.
Last, we come to the dissipative contributions, appearing in Slots 30-43. Here, the surprise is that the net dissipative effects acting on Band 2 are relatively large. Over all scales, the net effect is slightly less than −1.3 × 10 −3 m 2 /s 2 , although the total is dominated by the in-band, Band 2, dissipative impact of slightly more than −1.3 × 10 −3 m 2 /s 2 . The reason this is unexpected is that Band 2, from 1,536 to 1,200 km, is a large-scale band, and dissipation is dominantly a small-scale process. This is certainly true as well in the present simulation, as is seen in Figure 8a of the modeled dissipation field at Day 135. However, dissipation is broadly distributed and occurs at locations governed by the large scale, so that when filtered, dissipation imprints importantly on the evolution. This appears in Figure 8b which is the dissipation field at Day 135 filtered through Band 2.

Intermediate and Small Scales
Band 8, in Figure 6b, loses kinetic energy over the simulation, as can be seen by comparing the black line in the first slot with that in the first green dashed slot. However, the scale of this loss is close to two orders of magnitude smaller than those described for Band 2. Essentially, this is an inert band in regard to kinetic to kinetic exchanges. Across the bands, initial kinetic inputs to final kinetic energy are quite small, at a Figure 9. Time series of the pressure work input to kinetic energy. The pressure effect on Band 2 of Band 2 appears on the left and on Band 2 due to Band 3 on the right. Note the difference in scale on the y-axis. The net pressure work between these two cases is opposite, with the Band 3 to Band 2 case leveling off late in the event, while the Band 2 to Band 2 case continues to drain kinetic energy. net value of 3 × 10 −6 m 2 /s 2 . Note however that pressure work is quite active across all bands. Larger scale bands, to the left of the green dashed line in the pressure work zone, tend to be positive. This is indicative of a down-scale transfer of energy from larger scales to Band 8. However, smaller scales tend to be negative in amplitude and overall comparable in size to the contributions from the larger scales. The net effect of pressure work across all bands is negative, at a value of −1.7 × 10 −5 m 2 /s 2 , indicative of a net transfer to small scales. Band 8, with respect to pressure work, acts like a pass-through, allowing energy from larger scales to move through to smaller scales, while adding a small amount to the net. In this intermediate-scale band, the principal signal is a down-scale contribution to kinetic energy.
Finally, Band 12, in Figure 6c, is also a band of relatively weak energy transfers, if still stronger than Band 8. Overall, there is a stronger kinetic energy loss at these small scales, roughly −5.5 × 10 −5 m 2 /s 2 , most of which, −4.0 × 10 −5 m 2 /s 2 is accounted for by losses to larger length scales. This is indicative of a kinetic to kinetic upscale transfer. Once again, it is seen that the pressure work to kinetic energy transfers is far more active. While the net effect is relatively small, −1 × 10 −5 m 2 /s 2 , it is the result of large transfers of both signs across the length-scale spectrum. A clear pattern does not emerge. Losses and gains to other length-scale bands of both larger and small character occur. Perhaps the safest interpretation is that the weak values involved in this band indicate that it responds erratically to the conditions imposed on it by the energetically dominant larger length scales.
Time histories of two of the pressure work energy transfers is shown in Figure 9, i.e., those for Band 2 to Band 2 and for Band 3 to Band 2. These are among the most active pressure work bands seen in Figure 6and, despite their similarity in length scales, show opposite trends. The Band 2 to Band 2 transfer appears on the left. Early in the experiment, pressure work adds to kinetic energy, but the trends reverses around Day 120, and pressure work depletes kinetic energy strongly to Day 220. The last 100 days are thus associated with a large-scale tendency for flow to be directed toward high pressures and away from low pressures. Figure 9b, the pressure work contribution of Bands 3 and 2, shows the opposite trend. Pressure work persistently builds kinetic energy through most of the experiment, except for the period from Days 100 to 120. Toward the end, the net levels to a constant value, equivalent to the cessation of pressure work contributions between bands. The pressure effect on Band 2 of Band 2 appears on the left and from Band 3 to Band 2 on the right. The former case is dominated by negative values near regions of cyclones and reflects the tendencies at these largest scales for cyclones to merge. The distribution on the right-hand side is much more evenly distributed between positive and negative values and indicates a weaker pressure work associated with a more symmetrized, smaller scale structure.
These tendencies are explained in part by examining the anatomy of the pressure work contributions. For this, we show in Figure 10 the pressure work field as it appears to Band 2 at Day 200, for Bands 2 and 3. On the left is the pressure work as seen by Band 2, having been filtered through Band 2. Note that the structure is dominated by relatively large-scale negative values. These are zones dominated by cyclones and hence low pressures. The negative pressure work values indicate flows away from the low-pressure centers, which is a flow field designed to bring the lows together. As a result, the negative pressure work at largest scales is a signal of vortex merger. In contrast, the pressure work filtered through Band 3 as it appears to Band 2 appears on the right. Here, the distribution consists much more evenly of both highs and lows of comparable magnitude. Day 200 is relatively late in the merger event and is a point where on these scales of Band 3; the fields are becoming relatively symmetrized. As a result, the net pressure work from this band is quite small, consistent with the late parts of the time series in Figure 9 (right), where the time history of integrated pressure work has leveled off.

Summary
A new method for computing cross-scale energy transfers in fluid flows has been described and applied to a well-known example. The technique has both Lagrangian and Eulerian forms. In this paper, we mostly examine the former method, which exploits the Lagrangian nature of kinetic energy evolution to relate kinetic energy at any one location and time to its flow history. The Lagrangian history is computed by solving an adjoint equation obtained by a classical Green's function approach. The scale content of the final kinetic energy structure can be obtained by solving a filtered version of the adjoint equation and projecting the kinetic energy field onto the result. In addition, the same filter used on Green's function equation can be used to diagnose the scale-dependent inputs to the kinetic energy evolution, such as pressure work and dissipation, in a relatively straightforward way. This leads rather naturally to a quantification of how energy from one band of length scales transfers to another during the nonlinear evolution of the flow. The procedure also allows for a diagnosis of the nature of the energy transfers, e.g., kinetic to kinetic or kinetic to pressure work, as well as examination of how the transfers are structured in time and space. While we have here focussed on kinetic energy, the same approach can be applied to potential energy evolution.
We have applied the procedure to the classic problem of vortex merger in a stratified fluid. The well-known MOM6 model was deployed in a four layer isopycnal configuration and with an initial condition consisting of a sequence of undulations in second layer thickness. The system evolution proceeds with like-signed vortices merging and opposite-signed vortices pairing up and propagating. We analyzed a particular merger event in detail.
The subsequent analysis showed the expected result that kinetic energy exhibits an upscale cascade; however, the procedure also illustrates several somewhat novel aspects about merger. Perhaps most significantly, the primary energy transfers involve pressure work to kinetic energy exchanges, rather than kinetic to kinetic energy exchanges. This result complements several published studies of upscale energy transfers which have focussed on Fourier transforms of the nonlinear advection terms in the momentum equations by emphasizing the role played by pressure work. Dissipation also emerged as a surprisingly large effect on larger scale energy exchanges. This was due in part to the tendency for the larger scale flow to organize regions of small-scale dissipation so that the large-scale projection of dissipation was substantial. Also, several of the intermediate ranges in an integral sense were relatively inert, in that initial and final kinetic energies were often quite similar; however, those same ranges were found to be very dynamically active and acted as pass-throughs, receiving potential energy from smaller scales and passing them along to larger scales. Finally, the diagnoses identify the critical points in time when exchanges are at their peak and allows for an examination of the events associated with the exchanges. Spatial information about the nature of the exchanges is also provided and allows for a diagnosis of the involved mechanisms.
The example problem considered here is highly idealized, with an isotropic setting, doubly periodic boundaries and an unobstructed geometry. The promise of the technique, however, is that it potentially generalizes in a straightforward way to realistic ocean simulations, in which the flows are highly anisotropic and inhomogeneous, and the domains are often geometrically complex. The examination of such regions, such as the Gulf Stream, the Kuroshio, and other major ocean currents, is of importance to the understanding of the global ocean energy cycle. It will be of great interest to see what new light this method might cast on such regions.  insures they are also suitably normalized. Further classical analysis shows that arbitrary functions in V can be represented by an appropriate superposition of the eigenfunctions (in a least-squares sense) where the coefficients of the representation are computed by exploiting orthonormality.
In regard to the length scale, L, we would associate with any particular eigenfunction; the form of the Helmholtz equation suggests We will see in multiple dimensions that this simple relationship is problematic.
While for irregular domains with islands and topography, the solutions cannot generally be written down, they can be developed in relatively straightforward numerical ways. As an example, consider the domain shown in Figure A1, consisting of a single circular island inside of an otherwise square domain. The Laplacian operator is discretized using a 5-point stencil. Periodic boundary conditions are placed on the outer edge of the domain, and the eigenfunctions are required to vanish on the island boundary. These are conditions consistent with the orthogonality of the resulting Helmholtz solutions and are entirely reflected in the details of the stencil coefficients. In principle, the coefficients can be formed into a nx × ny by nx × ny square matrix, where nx and ny are the lengths of the domain in the east-west, north-south directions, respectively. For the present calculation, nx = n = 70 and grid spacing dx, d = 5m is isotropic. In fact, this does not account for the island, inside of which are several of the regularly spaced grid points but which represent locations inaccessible to the fluid. It is necessary to remove those points, resulting in slightly smaller matrix, the regular, pentadiagonal structure of which is then interrupted. This matrix can be fed into a standard eigenvalue extraction routine; we used the function "eig" in MATLAB. The gravest and 10th-most gravest eigenmodes appear in Figure A2; these are associated with eigenvalues k 2 1 ≈ 10 −4 m −2 on the left and k 2 10 ≈ 1.4 × 10 −3 m −2 , or equivalently L k ≈ 200 m (left) and L 10 ≈ 52 m (right). It is easy to associate a length scale of 200 m with the former, but the latter also measures comparable length scales while being spatially considerably more complex. It is thus not straightforward to simply connect the eigenvalue size and the associated length scale. This is due to the multiple spatial dimensions involved in the Figure A2. The first (left) and tenth (right) eigenmodes of the Helmholtz equation for the domain in Figure A1. The island location appears in black. Clearly, the gravest mode captures the broadest possible scales of the domain and can be associated with the length scale L = 200 m. The higher mode, however, while more complex spatially, is also measuring comparable length scales of roughly 200 m, in spite of being associated with a much higher eigenvalue.  Figure A1. The island location appears in black (in all panels), and the delta function appears just to the south of mid island. The γ 100 sampling function (top right). Note that it does not penetrate the island, instead sampling only the active fluid. The dashed black lines denote the locations of the transects appearing in the lower panels. The north-south transect through γ 100 appears on the lower left. The black lines denote the island location at this longitude. The east-west transect through γ 100 appears lower right. Note the qualitative resemblance of the structure to a sinc function. The distance between neighboring high lobes is roughly 75 m and consistent with the rough length-scale estimate of 90 m associated with the eigenvalue. analysis and arises in the unbounded domain using simple trigonometric functions as well. For example, high wavenumbers in the meridional direction coupled to low wavenumbers in the zonal direction can be associated with eigenvalues very similar to comparable wavenumbers of intermediate values distributed symmetrically in the two directions. We recommend a different way to associate length scales to the analysis below.
We now demonstrate that these functions satisfy constraint (Equation 18). Exploiting orthonormality, we anticipate ∑ n a n n (x) = ( where a n = n (x o ) The result of evaluating this formula using the present eigenfunctions appears in Figure A3 (top left), where the source point x o = 50 and 150 m, i.e., just south of the island. Clearly, we have reproduced to within domain discretization a Dirac delta function. The analysis in the main part of the paper, however, emphasized the use of integrals of filters between two length scales. We present in Figure A3 (top right) the result of summing all the eigenfunctions below a given eigenvalue. This corresponds to Equation 28, in particular to L o for the 100th eigenvalue (k 2 100 ≈ .012). Cross sections of the sampling function on the dashed lines in Figure A3b appear in the two lower panels. Note that the sampling function 100 avoids the island, as seen in Figure A3 (bottom left), while maintaining a qualitative resemblance to the sinc function ( Figure A3, bottom right). The length scale L associated with the eigenvalue is roughly 90 m, which corresponds well to the length scale of the sampling function, particularly as suggested by the distance separating the neighboring highest positive lobes seen in Figure A3. As suggested by this example, we recommend associating length scales L with sums of eigenfunctions corresponding to ranges of the eigenvalues.
Clearly, this technique, while promising, is in need of considerable study.

Appendix B: Eulerian Kinetic Energy Exchanges
This paper deals primarily with the Lagrangian version of the kinetic energy equation, somewhat in a manner reminiscent of Nagai et al. (2015) and Shakespeare and Hogg (2017). However, an Eulerian statement can be made as well, and the formulation is outlined here. Starting with Equation 29, the final time t o is chosen to be very small, such that so Equation B1 can be rewritten Using the chain rule converts Equation B3 to Finally, substituting with Equation 26 and integrating by parts lead to The left-hand side is recognized as the filtered version of the time rate kinetic energy change, and the right-hand side breaks it up into advective, diffusive, pressure work and dissipative contributions. The filter L n+1 L n is independent of time, so the left-hand side of Equation B5 can be rewritten Thus, the Eulerian version of the present analysis concerns the time rate of change of the band-passed kinetic energy. The kinetic energy appearing in the advective and diffusive terms can be similarly decomposed, which in turn yield formulae expressing the exchange of kinetic energy between bands involving spatial operations on the filter. For example, the advective contribution becomes After some algebra and modifying the order of integration In this Eulerian view, advective contributions appear explicitly. This procedure is somewhat like that described in Aluie et al. (2018), although we work directly in kinetic energy rather than the momentum equations. Last, this result can also be obtained by filtering the Eulerian kinetic energy equation directly.
where b is buoyancy and the velocity is fully three dimensional. For a linear equation of state, the quantity wb can be related to h = −bz via h t + ∇ ⋅ uh = −wb + ∇ ⋅ ∇h + 2 b z (C2) where is diffusivity.
For a nonlinear equation of state, the choice for h of "dynamic enthalpy" (Young, 2010) where o is a reference density, g gravity, P = − o gz is static pressure, and P o is a reference surface pressure, leads to a generalized form of Equation C2. This equation involves a lot of terms and is not written down here.
Note that Equation C2 can be rewritten as where = −wb + 2 b z represents the "sources" of potential energy. The form of Equation C4 is identical to that of Equation 4, so the previous analysis and results all apply to potential energy. If the fluid is modeled assuming identical diffusivity and viscosity (i.e., a Prandtl number of 1), the analysis in this paper can be applied to the sum of kinetic and potential energy. In any case, potential energy can always be analyzed separately from kinetic and the two connected through their common quantity, wb.