††thanks: Co-first author††thanks: Co-first author

# Spatio-temporal spectral transfers in fluid dynamics

Avik Mondal<EMAIL_ADDRESS>Department of Physics, University of Michigan,
Ann Arbor, MI 48109, USA  Andrew J. Morten  Department of Physics, University
of Michigan, Ann Arbor, MI 48109, USA  Brian K. Arbic  Department of Earth and
Environmental Sciences, University of Michigan, Ann Arbor, MI 48109, USA
Glenn R. Flierl  Department of Earth, Atmospheric, and Planetary Sciences,
Massachusetts Institute of Technology, Cambridge, MA 02139, USA  Robert B.
Scott  Département de Physique et Laboratoire des Mathémathiques Bretagne
Atlantique, Université de Bretagne Occidental, CNRS, Brest 29238, France
Joseph Skitka  Department of Earth and Environmental Sciences, University of
Michigan, Ann Arbor, MI 48109, USA

(August 28, 2024)

###### Abstract

Motivated by previous work on kinetic energy cascades in the ocean and
atmosphere, we develop a spatio-temporal spectral transfer tool that can be
used to study scales of variability in generalized dynamical systems. In
particular, we use generalized time-frequency methods from signal analysis to
broaden the applicability of frequency transfers from theoretical to practical
applications such as the study of ocean or atmosphere data or simulation
output. We also show that triad interactions in wavenumber used to study
kinetic energy and enstrophy cascades can be generalized to study triad
interactions in frequency or wavenumber-frequency. We study the effects of
sweeping on the locality of frequency transfers and frequency triad
interactions to better understand the locality of spatio-temporal frequency
transfers. As an illustrative example, we use the spatio-temporal spectral
transfer tool to study the results of a simulation of two-dimensional
homogeneous isotropic turbulence. This simulated fluid is forced at a well-
defined wavenumber and frequency with dissipation occurring at both large and
small scales, making this one of the first studies of “modulated turbulence”
in two dimensions. Our results show that the spatio-temporal transfers we
develop in this paper are robust to potential practical problems such as low
sampling rates or nonstationarity in time series of interest. We anticipate
that this method will be a useful tool in studying scales of spatio-temporal
variability in a wide range of fluids applications as higher resolution data
and simulations become more widely available.

##  I Introduction

The cascade model plays a fundamental role in turbulence theory. According to
the model, conserved quantities such as energy or enstrophy in two-dimensions
are transferred locally between nearby spatial scales starting from some
initial spatial scale and ending at some faraway dissipation scale. The range
of scales over which the cascade takes place is called the inertial range,
within which forcing and dissipation have negligible direct effect. While
energy and enstrophy transfers are known to not be perfectly local [1, 2, 3,
4, 5], the cascade framework has proven useful for predictions of wavenumber
spectrum $\mathcal{E}(k)$ of energy in the inertial range. Wavenumber spectra
predictions were derived by Kolmogorov [6, 7], Obukhov [8, 9], Onsager [10]
for three-dimensional turbulence and by Kraichnan [1], Leith [11], Batchelor
[12] for two-dimensional turbulence, and verified experimentally (less so for
the two-dimensional case) and numerically many times over.

An important component of the cascade picture is the spectral transfer, which
we shall refer to as the _spatial_ spectral transfer to avoid ambiguity. The
spatial spectral transfer quantifies the time rate of change of some spatial
spectral quantity, such as the wavenumber spectra of kinetic energy
$\mathcal{E}(k,t)$ or enstrophy $\mathcal{Z}(k,t)$, due to terms in the
equation of motion. For example, an energy spectral budget could take the form

$\partial_{t}\mathcal{E}(k,t)=\sum_{n}\mathcal{A}_{n}(k,t),$ (1)

where $t$ is time and $k$ is isotropic wavenumber. In this paper we will refer
to each term $\mathcal{A}_{n}(k,t)$ in Eq. 1 as a spatial spectral transfer.

Each term $\mathcal{A}_{n}(k,t)$ can be thought of as transferring energy into
or out of $\mathcal{E}(k,t)$ at wavenumber $k$. For example, an external
forcing term might transfer energy into the system, and a dissipation term
might transfer energy out of the system in the form of heat. In the cascade
picture, the spatial spectral transfer due to the nonlinear advection term
plays a central role, as it is the only term that can move energy between
scales. Indeed, sometimes the terminology “spectral transfer” is reserved
exclusively for the nonlinear advection term. However, we find it useful to
think of each term $\mathcal{A}_{n}(k,t)$ as a spectral transfer, whether the
energy transfer is into the system, out of the system, or between spatial
scales.

In contrast with spatial spectral transfers, temporal and spatio-temporal
spectral transfers, which would quantify the transfer of energy in frequency
or wavenumber-frequency space respectively, have had little success in terms
of theoretical predictions. For example, an analogous derivation of spectral
slopes in frequency space would be problematic, because nonlinear transfers
between time scales are expected to be less local than between spatial scales
if small-scale structures are swept by large-scale structures [13]. However,
while temporal and spatio-temporal spectral transfers may have limited
theoretical applications, they can still be used as practical diagnostics.
Temporal and spatio-temporal spectral transfers have already proven to be
useful diagnostics in a variety of contexts. Scott and Arbic [14] used spatial
spectral transfers to identify and examine an inverse cascade of kinetic
energy present in a simple two-layer, quasigeostrophic model. Arbic  _et al._
[15] diagnosed temporal spectral transfers in the ocean using altimetric
measurements of sea surface height, output from a realistic general ocean
circulation model, and the output from a simple two-layer quasi-geostrophic
model. Arbic  _et al._ [16] diagnosed spatio-temporal spectral transfers using
similar data sets. Both studies were motivated by the possibility of linking
the well-known nonlinear inverse cascade of energy towards larger spatial
scales with low-frequency variability in the ocean.

Müller  _et al._ [17] quantified spatio-temporal spectral transfers due to
internal gravity waves in two global ocean models at different resolutions.
They found that the internal gravity wave frequency spectrum in the higher
resolution simulation agreed better with spectra calculated from observations.
Pan  _et al._ [18] and Skitka  _et al._ [19] used spatio-temporal spectral
transfers to study internal gravity wave energy transfers in high-resolution
regional ocean models. Sérazin  _et al._ [20] applied spatio-temporal spectral
transfers to the long-duration output from a global ocean model, finding that
advective transfers move energy from high-frequency frontal Rossby waves
towards lower-frequency westward-propagating mesoscale eddies. Morten  _et
al._ [21] used spatio-temporal spectral transfers to study kinetic energy
dynamics in a simulation of single-layer, shallow-water beta-plane quasi-
geostrophic turbulence. More recently, spatio-temporal and temporal transfers
have been used to study interactions between the ocean and the atmosphere.
O’Rourke  _et al._ [22] compared spatio-temporal transfers due to advection
with transfers due to wind stress in the oceanic component of a coupled ocean-
atmosphere model. Martin  _et al._ [23] employed temporal transfers of kinetic
energy in the atmosphere and the ocean to understand the sources of
variability at different timescales of each fluid in a quasi-geostrophic,
fully coupled ocean-atmosphere simulation. Martin  _et al._ [24] extended this
framework by computing a temperature variance budget to diagnose the sources
of variability in an idealized ocean-atmosphere system as a function of
frequency.

The purpose of this paper is to provide a theoretical framework to aid in the
interpretation of temporal and spatial-temporal spectral transfers calculated
from theory, simulation, or data. Previous studies have derived spatio-
temporal spectral transfers in theoretical contexts. In related work, Chiu
[25] derived kinetic energy spectral equations in the frequency domain for
large-scale atmospheric motions. Sheng and Hayashi [26, 27] independently
derived similar equations and applied them to global atmospheric simulations.
Elipot and Gille [28] also derived spectral equations to study the frequency
components that dominate the wind energy input into the Ekman layer in the
Southern Ocean. All of these studies interpret the frequency spectra
diagnostic as a balance among terms. For example, Chiu [25] notes that while
the spectral budget in wavenumber space deals with the time rate of change of
the spectrum, the analogous equation in frequency space contains no time rate
of change and is therefore simply a balance among various co-spectra. However,
we show that a time rate of change can indeed be incorporated into the
spectral budget in the frequency domain, resulting in an improved
interpretation.

To further aid in the interpretation of temporal and spatio-temporal spectral
transfers, we investigate the effect of a mean flow. When there is a mean
flow, we may invoke Taylor’s hypothesis [29, 30, 31, 32, 33, 34, 35], which
states that small-scale structures are advected without significant distortion
by the large-scale mean flow. In other words, at sufficiently small scales the
turbulent structures appear to translate uniformly at the mean velocity
$\bm{U}$. Under this assumption, there is a mapping between spatial and
temporal Fourier modes given by $\omega=\bm{k}\cdot\bm{U}$, where $\omega$ is
frequency. We calculate this mapping in the framework of spatio-temporal
transfers by approximating a mean flow with a Galilean transformation of a
fluid’s velocity field and comparing this boosted velocity field to a
reference velocity field.

Analysis of Taylor’s hypothesis brings to light questions about the locality
of spatio-temporal transfers as a result of a mean flow. To further study this
locality, we introduce spatio-temporal triad interactions, which are a
generalization of the wavenumber triad interactions introduced by Kraichnan
[1]. We also discuss other uses and properties of these spatio-temporal triad
interactions.

To illustrate the use of spatio-temporal spectral transfers, we apply the
diagnostic to a simple fluid system: a modified version of the incompressible
two-dimensional Navier–Stokes equation [for reviews of two-dimensional
turbulence see 36, 37, 38, 39]. Motivated by the utility of diagnosing
_spatial_ spectral transfers in systems where energy and enstrophy are
injected within a narrow range of wavenumber [recent examples include 40, 41,
42, 43, 44, 45, 46], we choose to study the _spatio-temporal_ spectral
transfers in a system where energy and enstrophy are injected within a narrow
range of wavenumber _and frequency_. As such, we apply a forcing that is
approximately sinusoidal in space and time.

The effect of a sinusoidal modulated forcing has been investigated for three-
dimensional turbulence [47, 48, 49, 50, 51, 52]. In the three-dimensional case
resonant frequencies were predicted using the static structure functions of
Effinger and Grossmann [53]. To our knowledge no analogous result has been
calculated for two-dimensional turbulence, which is complicated by the dual
cascade of energy and enstrophy. The motivation for studying two-dimensional
turbulence is to connect our results with prior studies of spatio-temporal
spectral transfers calculated using oceanic data and output of realistic ocean
models [15, 16]. In those studies, for example, one region experienced an
energy transfer to smaller time scales, an unexpected result that may be
partly explained by the idealized numerical investigations used here.

A second motivation for the numerical investigation is to test the accuracy of
the temporal spectral transfers as a diagnostic for datasets with limited
temporal resolution or duration. Researchers studying large oceanic datasets
often have to contend with both problems. For example, satellite altimeter
measurements of sea surface height used to calculate spectral transfers suffer
from both limited temporal resolution and duration, such that the data cannot
resolve some relevant dynamical time scales. Because our numerical simulations
produce data that resolve all dynamical time scales, we can comprehensively
study the effects of limited record duration, limited temporal resolution, and
temporal detrending in a way that is not necessarily possible with realistic
data. Our investigation shows that temporal spectral transfers can be accurate
over a range of frequencies under certain conditions even when the data set is
severely limited.

In Sec. II, we derive temporal and spatio-temporal spectral transfers for a
general equation of motion. Because data sets are typically detrended in
practical applications, we also show how a detrending operation may be
incorporated into the spectral transfers. If the correct detrending operation
is applied, the spectral budget remains exact.

In Sec. III, we assume that the equation of motion is two- or three-
dimensional forced-dissipated Navier-Stokes, which has a nonlinear advection
term. We derive temporal and spatio-temporal triad interactions and show how a
mean flow affects the spectral transfers and triad interactions.

In Sec. IV, we calculate spatial, temporal, and spatio-temporal transfers
using the output of a simulation of the incompressible two-dimensional Navier-
Stokes equation. This simulation is a simplified version of the one used in
Morten  _et al._ [21]. We begin Sec. IV by introducing the details of this
simplified simulation. We then diagnose the spatial, temporal, and spatio-
temporal spectral transfers and fluxes for specific simulation outputs. We
show how various limitations of the data or simulation output may affect the
diagnosis of the transfers. Such limitations include poor temporal resolution,
inadequate duration of the dataset, and the existence of a trend. We show that
for our two-dimensional simulations, the temporal transfers and fluxes are
fairly robust when applied to imperfect data. We also show how the temporal
spectral transfers change in time starting with a fluid at rest, diagnose the
effects of detrending, and look for evidence of sweeping in the spatio-
temporal spectral transfers.

###  I.1 Notation

Throughout this paper, we employ a useful notational convention: a single
function name may be used more than once to refer to several different
functions, with each unique function distinguished by the dimensions and units
of the input parameters. For example, $\mathcal{A}_{n}(\bm{k},t)$,
$\mathcal{A}_{n}(k,t)$ and $\mathcal{A}_{n}(k)$ are three different spatial
spectral transfers. Similarly, we will later derive three different spatio-
temporal and temporal spectral transfers:
$\mathcal{A}_{n}(\bm{k},\omega,\tau)$, $\mathcal{A}_{n}(k,\omega,\tau)$, and
$\mathcal{A}_{n}(\omega,\tau)$. The distinction implied by the notation is
that $\mathcal{A}_{n}(\omega)$ is $\mathcal{A}_{n}(\bm{k},\omega)$ with the
$\bm{k}$ dependence integrated out. The angular dependence is integrated out
of $\mathcal{A}_{n}(\bm{k},\omega,\tau)$ to get
$\mathcal{A}_{n}(k,\omega,\tau)$. Spectral energy densities
$\mathcal{E}(\cdot)$ and spectral fluxes $\Pi_{\mathcal{A}_{n}}^{>}(\cdot)$
will also follow the same notation. This convention greatly reduces the number
of symbols that the reader would otherwise need to remember. However,
sometimes we want to leave out a parameter to avoid repeatedly writing a
ubiquitous parameter. In this case, we will introduce a function with
parameters after a semicolon but then neglect to write these parameters in
future uses of this function. This does not imply that the function’s
dependence on the neglected parameter is integrated out. We simply leave them
out for terser notation. For example, $\mathcal{E}(\bm{k},\omega,\tau;T)$ is
the same function as $\mathcal{E}(\bm{k},\omega,\tau)$, but with the
dependence on the parameter $T$ made explicit. Additionally, we will sometimes
consider functions and their corresponding transforms (i.e various types of
Fourier transforms or filtered functions). For such functions, we will use
diacritics such as tildes, overbars, carons, and circumflexes depending on
which transform we are applying to a function. We avoid such diacritics with
spectral transfers $\mathcal{A}$, spectral fluxes $\Pi_{\mathcal{A}_{n}}^{>}$,
and spectral energy densities $\mathcal{E}$ because we only consider these in
spectral space. When multiple transforms are made, we will use multiple
diacritics.

##  II Spectral transfers for general system

We derive temporal and spatio-temporal spectral transfers in a general way so
that they may be used in a wide variety of fluid dynamical applications. To be
widely applicable, the initial equation of motion should be as general as
possible. Therefore, we consider a general equation of the form

$\partial_{t}g(\mathbf{x},t)=f(\mathbf{x},t),$ (2)

where $g$ is a general scalar field and $f$ is the flux of $g$. $g$ can be
used to refer to the components of vector fields. The function $f$ should be
considered a finite sum of terms, such as forcing, dissipation, and nonlinear
advection. It will be useful to write this explicitly as

$f(\bm{x},t)=\sum\limits_{n}A_{n}(\bm{x},t),$ (3)

where each $A_{n}$ contributes to the time derivative of $g$. We derive
temporal and spatio-temporal spectral transfers with Eq. 2 as the starting
point. Later on, we will specifically consider two- and three-dimensional
fluid dynamical equations that contain a nonlinear advection term.

It may help to motivate the use of Eq. 2 by considering some examples. We
could start with the incompressible Navier–Stokes equations with general
dissipation and forcing terms:

$\displaystyle\partial_{t}\bm{u}(\bm{x},t)=-(\bm{u}\bm{\cdot}\bm{\nabla})\bm{u}-$
$\displaystyle\bm{\nabla}\Pi+\bm{d}[\bm{u}]+\bm{f}(\bm{x},t),$ (4)
$\displaystyle\bm{\nabla}\bm{\cdot}\bm{u}$ $\displaystyle=0,$ (5)

where $\bm{u}(\bm{x},t)$ is the Eulerian velocity field, $\Pi$ is pressure
normalized by the density of water, $\bm{d}[\cdot]$ is a fairly general linear
dissipation operator

$\bm{d}[\bm{u}]:=-\sum_{n}\nu_{n}(-\nabla^{2})^{n}\bm{u},$ (6)

with dissipation coefficients $\nu_{n}$, and $\bm{f}(\bm{x},t)$ is an external
force. We could start with the spatial Fourier transform of Eq. 5:

$\displaystyle\partial_{t}\widetilde{\bm{u}}(\bm{k},t)=-\widetilde{(\bm{u}\bm{\cdot}\bm{\nabla})\bm{u}}-$
$\displaystyle\widetilde{\bm{\nabla}\Pi}+\widetilde{\bm{d}[\bm{u}]}+\widetilde{\bm{f}}(\bm{k},t),$
(7) $\displaystyle\bm{k}\bm{\cdot}\widetilde{\bm{u}}(\bm{k})$
$\displaystyle=0,$ (8)

where “tilde” indicates a two- or three-dimensional spatial Fourier transform
over a periodic domain. Saltzman [54] uses a version of Eq. 7 to study large
scale atmospheric turbulence. In this paper, Eq. 7 will be the starting point
for Sec. III. If instead, we were interested in looking at individual terms in
the multi-scale gradient expansion of Eyink [55, 56], we could begin with a
low-pass filtered equation [57, 58]:

$\displaystyle\partial_{t}\overline{\bm{u}}(\bm{x},t,\ell)=-(\overline{\bm{u}}\bm{\cdot}\bm{\nabla})\overline{\bm{u}}-\bm{\nabla}\overline{\Pi}$
$\displaystyle-\bm{\nabla}\bm{\cdot}\bm{\tau}+\bm{d}[\overline{\bm{u}}]+\overline{\bm{f}}(\bm{x},t,\ell),$
(9) $\displaystyle\bm{\nabla}\cdot\overline{\bm{u}}$ $\displaystyle=0$ (10)

where

$\overline{A}(\bm{x},t,\ell):=\int
d^{n}\bm{r}G_{\ell}(\bm{r})A(\bm{x}+\bm{r}),$ (11)

$\bm{\tau}:=\overline{\bm{u}\bm{u}}-\overline{\bm{u}}\;\overline{\bm{u}}$, and
$G_{\ell}$ is a low-pass filter kernel that keeps only length scales larger
than $\ell$. We could also apply a spatial wavelet transform to Eq. 5 or start
with some other primitive fluid equation.

We have great flexibility in both the choice of the primitive equation and the
choice of spatial transform.

###  II.1 Spatial spectral transfers

Spatial spectral transfers are relatively simple to derive. Combining Eq. 2
and Eq. 3 and taking the spatial Fourier transform of the resulting equation,
we can write

$\partial_{t}\widetilde{g}(\bm{k},t)=\sum_{n}\widetilde{A}_{n}(\bm{k},t).$
(12)

We can then multiply Eq. 12 by $\widetilde{g}^{*}(\bm{k},t)$ and take its real
part to obtain the spectral budget,

$\partial_{t}{\textstyle\frac{1}{2}}|\widetilde{g}(\bm{k},t)|^{2}=\sum_{n}\mathcal{A}_{n}(\bm{k},t),$
(13)

where

$\mathcal{A}_{n}(\bm{k},t):=\mbox{Re}\left[\widetilde{g}^{*}(\bm{k},t)\widetilde{A}_{n}(\bm{k},t)\right],$
(14)

and the superscript $*$ denotes complex conjugate. For each term
$\widetilde{A}_{n}(\bm{k},t)$ in the equation of motion,
$\mathcal{A}_{n}(\bm{k},t)$ is the corresponding spatial spectral transfer. If
$g$ were a vector rather than a scalar, then we would use a dot product with
$\widetilde{g}^{*}(\bm{k},t)$. The imaginary part of the product of Eq. 12 and
$\widetilde{g}^{*}{(\bm{k},t)}$ would tell us about the evolution of the phase
of the $\widetilde{g}(\bm{k},t)$. We discuss the meaning of this imaginary
part in the Appendix A.

As a concrete example, if the equation of motion were Eq. 7 (as will be the
case in Sec. III), then we would have
$\widetilde{g}=\widetilde{\bm{u}}(\bm{k},t)$, and the spatial spectral budget
would involve the time rate of change of

$\mathcal{E}(\bm{k},t):={\textstyle\frac{1}{2}}\left|\widetilde{\bm{u}}(\bm{k},t)\right|^{2},$
(15)

where $\mathcal{E}(\bm{k},t)$ is the energy in wavevector mode $\bm{k}$ at
time $t$. Motivated by this example, we refer to
${\textstyle\frac{1}{2}}\left|\widetilde{g}(\bm{k},t)\right|^{2}$ as the
“energy.” The spatial spectral budget Eq. 13 tells us the time rate of change
of energy in spatial mode $\bm{k}$ due to each term in the equation of motion.
In other words, $\mathcal{A}_{n}(\bm{k},t)$ is the spatial spectral transfer
corresponding to the term $\widetilde{A}_{n}(\bm{k},t)$ in the equation of
motion.

###  II.2 Spatio-temporal spectral transfers

One might expect to be able to derive temporal and spatio-temporal spectral
transfers analogous to the derivation of spatial spectral transfers in Sec.
II.1. However, temporal spectral transfers are fundamentally different because
of the temporal derivative in Eq. 2. While the temporal derivative commutes
with spatial detrending, spatial Fourier transforms, and other spatial
transforms, the temporal derivative does not in general commute with temporal
transforms. In this section, we discuss the difficulties with deriving
spectral transfers in frequency space, many of which have been observed in
previous works. Furthermore, we will demonstrate how to correctly calculate
such spectral transfers.

Applying a temporal Fourier transform to Eq. 2 results in the replacement of
the time-derivative by $i\omega$, giving

$i\omega\check{\widetilde{g}}(\bm{k},\omega)=\check{\widetilde{f}}(\bm{k},\omega)=\sum_{n}\check{\widetilde{A}}_{n}(\bm{k},\omega).$
(16)

In this paper, the “caron” operator, $\check{\hskip 7.22743pt}$, denotes a
temporal Fourier transform.

Multiplying by the complex conjugate
$\check{\widetilde{g}}^{*}(\bm{k},\omega)$ and taking the real part gives

$0=\sum_{n}\mbox{Re}\left[\check{\widetilde{g}}^{*}(\bm{k},\omega)\check{\widetilde{A_{n}}}(\bm{k},\omega)\right]$
(17)

because multiplying any complex function by its complex conjugate and
$i\omega$ results in an entirely imaginary function. Eq. 17 looks similar to
the spatial spectral budget given by Eq. 13, except that there is no longer a
temporal derivative that leads to the interpretation of each term as a
spectral transfer (i.e. as a rate of change). All that remains is a balance of
terms. This is indeed what Sheng and Hayashi [26] found when calculating
atmospheric kinetic energy budgets.

If we instead take the imaginary part we obtain

$\omega\left|\check{\widetilde{g}}(\bm{k},\omega)\right|^{2}=\sum_{n}\mbox{Im}\left[\check{\widetilde{g}}^{*}(\bm{k},\omega)\check{\widetilde{A_{n}}}(\bm{k},\omega)\right].$
(18)

The resulting Eq. 18 resembles an energy budget and has been used in several
other papers [59, 28]. However, the right-hand-side of Eq. 18 does not say
anything about the transfer of energy in spectral space. We will show in
Appendix A that it instead tells us about the phase dynamics of the time
series of $g$. Our above usage of the Fourier transform, which resembles the
derivation of a spatial-spectral transfers, cannot in general be used for
practical applications for several reasons.

One of the issues is that the Fourier transform is only defined for functions
that decay to zero sufficiently rapidly as $t\rightarrow\pm\infty$. Time
series from data or simulation do not necessarily explicitly decay to zero
during the timescales of observation. This issue also applies to calculating
the Fourier transforms of the autocorrelation function of $g(\bm{x},t)$ to
exploit the Wiener-Khinchin theorem using the method of Chiu [25]. A Fourier
transform of data that is not set to zero at the beginning and the end of the
time window of observation will have ringing that will obscure the true
spectral behavior of a function.

Furthermore, we are often interested in time series that are not statistically
stationary in time. This usually manifests as a time signal having a time
average over a finite window that changes significantly in time. Such time
signals will have a periodic component and a trend. The Fourier transform of a
time series that includes a trend can be corrupted due to potential
aperiodicities.

The above considerations motivate the need for a more careful derivation of
spatio-temporal transfers. We turn to generalizations of the Fourier transform
that work for statistically non-stationary, aperiodic time series, the so-
called time-frequency representations [60]. Time-frequency representations
$\mbox{TFR}_{f}(\omega,\tau)$ are functions derived from some time signal
$f(t)$ of two variables: a frequency or scale $\omega$; and a “central time”
$\tau$, which is typically the center of some time window of interest. Such
time windows are usually defined by the range over which data of interest is
available. At a fixed $\tau$, a time-frequency representation gives a
transformed spectrum of a time series as a function of frequency or scale with
respect to a chosen set of orthonormal basis functions. These transformations
can solve the problems of aperiodicity and nonstationarity mentioned
previously in this section. See Hlawatsch and Auger [60] for details on
general bilinear time-frequency representations. In this paper, we primarily
focus on the spectrogram, the most commonly used time-frequency
representation. We briefly discuss the scalogram in Appendix C.

The spectrogram is the modulus squared of the short-time Fourier transform
(STFT). For any function $f(t)$, the STFT is

$\displaystyle\widehat{f}(\omega;\tau,T)$
$\displaystyle:=\int_{\tau-T/2}^{\tau+T/2}\sigma(t-\tau;T)f(t){e^{-i\omega(t-\alpha\tau)}dt},$
(19)

where $\sigma(t-\tau;T)$ is a taper function of width $T$ centered at
$t=\tau$. We use the circumflex to denote an STFT. We assume
$\sigma(t-\tau;T)$ to be identically zero outside of the range
$\tau-T/2<t<\tau+T/2$, so one can also define the integration to be over the
range $(-\infty,\infty)$ without changing the result of the integral.

Figure 1: Candidate time series $f(t)$ for STFT with taper function
$\sigma(t)$ within a specified time window, denoted here in white. The grey-
shaded region is not considered in an STFT calculation of $f(t)$ with chosen
central time $\tau$ and window width $T$. The black solid curve corresponds to
an artificial time series consisting of white noise with standard deviation
$\frac{V}{2}$ and periodic functions of magnitude $V$. The red dotted curve
corresponds to the taper function $\sigma(t)$ that guarantees that the time
series decays to $0$ at the boundaries of the specified time window. We
emphasize that $T$ and $\tau$ are time variables, $T$ sets a timescale while
$\tau$ sets the central time location of the specified window.

The STFT is effectively a Fourier transform applied to a tapered signal with
the central position $\tau$ of the taper function made explicit. Whenever the
choice of taper function $\sigma$ must be made explicit, we use the notation
$\widehat{f}[\sigma](\omega,\tau)$. We show an example of a potential time
series with a Hann function as the taper function in Fig. 1. While many
candidate functions for $\sigma(t)$ are referred to as windows (Hann window,
Tukey window, etc.), we refer to them as taper functions in this paper to
avoid confusion with the time windows centered at $\tau$. Note that because we
apply STFT’s with a finite $T$, the case of “no taper” is effectively the same
as setting $\sigma(t)$ as a rectangular function:

$\sigma_{R}(t-\tau;T)=\Theta\Big{(}t-\tau+\frac{T}{2}\Big{)}-\Theta\Big{(}t-\tau-\frac{T}{2}\Big{)}.$
(20)

In the above equation, $\Theta$ is the Heaviside step function.

The parameter $\alpha$ in Eq. 19 takes one of two values: $\alpha=0$ or
$\alpha=1$. The standard definition of STFT has $\alpha=0$ with basis
functions $e^{-i\omega t}$ centered at $t=0$, which corresponds to complex
exponential functions that do not move with the time window. It will be useful
to consider the $\alpha=1$ case, which has Fourier basis functions
$e^{-i\omega(t-\tau)}$ centered at $t=\tau$, which corresponds to complex
exponential functions that do move with the time window. The two cases differ
only by a multiplicative factor of $e^{i\omega\tau}$ that commutes with the
integral in Eq. 19.

Applying the STFT to Eq. 12 gives

$\int_{-\infty}^{\infty}\sigma(t-\tau;T)\left(\partial_{t}\widetilde{g}(\bm{k},t)\right)e^{-i\omega(t-\alpha\tau)}dt=\sum_{n}\widehat{\widetilde{A_{n}}}(\bm{k},\omega,\tau).$
(21)

Integrating by parts and transforming $t$-derivatives into $\tau$-derivatives
gives

$\partial_{\tau}\widehat{\widetilde{g}}(\bm{k},\omega,\tau)+(1-\alpha)i\omega\widehat{\widetilde{g}}=\widehat{\widetilde{f}}(\bm{k},\omega,\tau).$
(22)

Multiplying Eq. 22 by $\frac{1}{2}\widehat{\widetilde{g}}^{*}$ and then adding
the product of $\frac{1}{2}\widehat{\widetilde{g}}$ and the complex conjugate
of Eq. 22 gives

${\textstyle\frac{1}{2}}\partial_{\tau}|\widehat{\widetilde{g}}(\bm{k},\omega,\tau)|^{2}=\mbox{Re}[\widehat{\widetilde{g}}^{*}\widehat{\widetilde{f}}].$
(23)

Expanding $f$ as a sum of terms gives the temporal spectral budget:

${\textstyle\frac{1}{2}}\partial_{\tau}|\widehat{\widetilde{g}}(\bm{k},\omega,\tau)|^{2}=\sum_{n}\mathcal{A}_{n}(\bm{k},\omega,\tau).$
(24)

where

$\mathcal{A}_{n}(\bm{k},\omega,\tau):=\mbox{Re}[\widehat{\widetilde{g}}^{*}(\bm{k},\omega,\tau)\widehat{A_{n}}(\bm{k},\omega,\tau)].$
(25)

Unlike the balance obtained in Eq. 17, the spectral budget Eq. 24 contains a
$\tau$-derivative of the energy. The $\tau$-derivative of the energy
${\textstyle\frac{1}{2}}|\widehat{\widetilde{g}}(\bm{k},\omega,\tau)|^{2}$
tells us how the energy in mode $(\bm{k},\omega)$ changes as the central time
changes (i.e. as the window of the STFT advances). The interpretation of the
temporal spectral budget Eq. 24 is therefore analogous to that of the spatial-
spectral budget Eq. 13 with the $t$-derivative replaced by a
$\tau$-derivative. For a statistically homogeneous time series with a large
enough time window, the $\tau$-derivative should act as a residual. If all the
appropriate terms are included in Eq. 12, the $\tau$-derivative should equate
to zero. If there are stochastic flux terms in the equation of motion, then
one would expect the average of the $\tau$-derivative over a sufficiently
large ensemble would be zero. Additionally, one may find that the
$\tau$-derivative is non-zero if the chosen time window is too small to
capture the long period contributions to the input signal; we see examples of
this in Figs. 10 and 11. The spectral budget Eq. 24 can be derived using the
autocorrelation function of a time series in the manner of Chiu [25] instead
of the above method; we show this derivation in detail in Appendix B. In the
STFT defined in Eq. 19, we did not incorporate any means of dealing with
nonstationary functions; we dedicate the next subsection to this topic.

###  II.3 Effects of detrending

In practice, the standard method of dealing with nonstationarity is detrending
[61]. Detrending subtracts a chosen function from a time series, derived from
a best-fit of a function with the time series. In previously published papers
by some of the authors, detrending consists of computing the linear trend of
the time series and subtracting it from the time series [24, 23, 21, 15, 16].
Removing this best-fit line removes a square wave or sawtooth waveform from
the time series. This removes ringing artifacts from the Fourier transform
that corrupt the entire frequency range of a time series. We emphasize that
simply subtracting a linear trend is distinct from applying any type of low or
high-pass filters which can be used to separate time series into “slow” and
“fast” parts. Separating time scales does not entirely remove ringing effects
because ringing corrupts the spectra at all frequencies.

In this section we show how to derive an exact spectral budget, giving the
$\tau$ rate of change of the energy in the detrended signal:

$\frac{1}{2}\partial_{\tau}|\widehat{\widetilde{g_{detrend}}}(\bm{k},\omega;\tau)|^{2}=\mbox{Re}[\widehat{\widetilde{g_{detrend}}}^{*}\widehat{\widetilde{f_{detrend}}}].$
(26)

Here, the subscript “detrend” indicates the application of an operator that
removes the trend from a function. Eq. 26 is not obvious and in fact depends
on the details of how the detrending is formulated and incorporated with the
STFT. In this section, we allow for detrending operations beyond subtracting a
linear trend. We derive the requirements necessary for a generalized
detrending operation to give a spectral budget of the form Eq. 26.

Given a function $h(t)$, we calculate its trend over a finite window of width
$T$ centered at time $\tau$. Thus, the trend will depend on the position of
the window, so it will have the form $h_{trend}(t,\tau)$. For example, a
linear trend takes the form

$h_{trend}(t,\tau)=c_{0}(\tau)+c_{1}(\tau)(t-\tau),$ (27)

where $c_{0}(\tau)$ and $c_{1}(\tau)$ are the fit coefficients, which depend
on the position of the window. For more general trends, because we only
calculate spectrograms for finite time windows, we can define a trend function
as:

$h_{trend}(t,\tau):=\sum_{n}c_{n}(\tau)\varphi_{n}(t-\tau)$ (28)

where the $\varphi_{n}(t-\tau)$ refer to members of an orthogonal basis that
move with the position of the window. We use angle brackets to denote the
inner product for basis functions $\varphi$:

$\innerproduct{\varphi_{n}(t^{\prime})}{\varphi_{m}^{*}(t^{\prime})}_{t^{\prime}}:=\int_{-T/2}^{T/2}dt^{\prime}\varphi_{n}(t^{\prime})\varphi_{m}^{*}(t^{\prime}).$
(29)

For simplicity, we consider trends constructed from bases with unweighted
inner products; for example, for polynomial trends of order $b$ one can use
the first $b$ Legendre polynomials.

We use Eq. 28 to define the detrending operation:

$h_{detrend}:=h(t)-h_{trend}(t,\tau).$ (30)

Because it is commonly used, we consider in detail the linear least squares
method for detrending. The coefficients $c_{n}$ are determined by minimizing
the squared residual as shown in Appendix D. The resulting best-fit
coefficients are

$c_{n}(\tau)=\sum_{m}(M^{-1})_{nm}\innerproduct{h(t^{\prime}+\tau)}{h(t^{\prime}+\tau)}{\varphi_{m}^{*}(t^{\prime})}_{t^{\prime}},$
(31)

In Eq. 31, $M$ is a square invertible matrix independent of $t$ and $\tau$.
Plugging Eq. 31 into Eq. 28 gives the trend $h_{trend}(t,\tau)$ in terms of
$f(t)$:

$h_{trend}(t,\tau)=\sum_{n}\varphi_{n}(t-\tau)\sum_{m}(M^{-1})_{nm}\innerproduct{h(t^{\prime}+\tau)}{\varphi_{m}^{*}(t^{\prime})}_{t^{\prime}}.$
(32)

With Eq. 32, we can construct trends of any form that can be derived with
linear least squares and an orthogonal basis with an unweighted inner product.
To derive a detrended spectral budget, we can combine Eq. 2 and Eq. 26 to get:

$(\partial_{t}g)_{detrend}(\bm{x},t)=f_{detrend}(\bm{x},t).$ (33)

We separately calculate $(\partial_{t}g)_{trend}$, $\partial_{t}g_{trend}$ and
$\partial_{\tau}g_{trend}$ in Appendix D and combine the results to obtain

$(\partial_{t}g)_{trend}=(\partial_{\tau}+\partial_{t})g_{trend}(t,\tau),$
(34)

from which it easily follows that

$(\partial_{t}g)_{detrend}=(\partial_{\tau}+\partial_{t})g_{detrend}(t,\tau).$
(35)

In other words, applying the time-derivative _before_ removing the trend gives
the same result as applying the operator $(\partial_{\tau}+\partial_{t})$
_after_ removing the trend.

In addition to satisfying Eq. 35, linear least squares satisfies linearity,

$(f_{1}+f_{2})_{detrend}=(f_{1})_{detrend}+(f_{2})_{detrend},$ (36)

which is useful because the spectral budget can remain a sum of terms before
and after a linear least squares detrending operation.

To continue the derivation of the spectral budget, we would like to apply the
STFT after the detrending operation. However, some care is required in
applying the STFT after detrending. Whereas the standard STFT takes as input a
function dependent only on the time variable $t$, after detrending we actually
have a function dependent on $t$ and window time position $\tau$. It is
natural to equate the window position $\tau$ used by the detrending operation
with the window position $\tau$ used by the STFT, which is why we have already
used $\tau$ in both situations. Given a function of the form $f(t,\tau)$ as
input, we therefore define the generalized STFT as

$\displaystyle\mbox{STFT}_{f}(\omega,\tau):=\int_{\tau-T/2}^{\tau+T/2}\sigma_{T}(t-\tau)f(t,\tau){e^{-i\omega(t-\alpha\tau)}dt}.$
(37)

The change is subtle, but it is important to understand that the input to the
STFT depends on the central time when we incorporate detrending. We use the
same notation as for the standard STFT, since the only change is in the form
of the input function.

Combining Eqns. 19, 22, and 35, we find that the generalized STFT of the
detrended derivative of a function $g$ is

$\displaystyle\widehat{\left(\partial_{t}g\right)_{detrend}}(\omega,\tau)=\partial_{\tau}\widehat{g_{detrend}}(\omega,\tau)+(1-\alpha)i\omega\widehat{g_{detrend}}(\omega,\tau).$
(38)

Applying Eq. 38 to the general equation of motion Eq. 12 gives

$\partial_{\tau}\widehat{g_{detrend}}(\omega,\tau)+(1-\alpha)i\omega\widehat{g_{detrend}}=\widehat{f_{detrend}}.$
(39)

Multiplying by $\widehat{g_{detrend}}^{*}$ and taking the real part gives

${\textstyle\frac{1}{2}}\partial_{\tau}\left|\widehat{g_{detrend}}(\omega,\tau)\right|^{2}=\mbox{Re}\left[\widehat{g_{detrend}}^{*}\widehat{f_{detrend}}\right].$
(40)

Expanding $f(t)$ as a sum and using the linearity condition Eq. 36 gives the
exact spectral budget:

${\textstyle\frac{1}{2}}\partial_{\tau}\left|\widehat{g_{detrend}}(\omega,\tau)\right|^{2}=\sum_{n}\mathcal{A}_{n}(\omega,\tau),$
(41)

where

$\mathcal{A}_{n}(\omega;\tau):=\mbox{Re}\left[\widehat{g_{detrend}}^{*}\widehat{(A_{n})_{detrend}}\right].$
(42)

Thus, we have shown that the temporal spectral budget remains balanced even
when a detrending operation is incorporated. While linear least squares was
used in the previous derivation, any detrending operation that satisfies Eq.
35 and Eq. 36 will also result in an exact equality in the spectral budget.
For example, any linear filter acting upon windowed data works as well. That
is, $h_{trend}$ can take the form of a convolution of the windowed signal with
an impulse response function. This can be expressed as
$h_{trend}(t,\tau):=h(t)*_{t}(\sigma(t-\tau)f(t))$, where $*_{t}$ indicates
convolution over the variable $t$, and where $h(t)$ is the impulse response
function (i.e. inverse Laplace transform of the transfer function).
Importantly, the taper function $\sigma(t-\tau)$ must again depend on $t$ and
$\tau$ specifically through the combination $t-\tau$.

##  III Spectral kinetic energy transfers in fluid dynamics

We now turn our attention to temporal and spatio-temporal spectral kinetic
energy transfers in fluid dynamics applications. We derive for spatio-temporal
and temporal transfers the nonlinear triad interactions already identified for
spatial transfers in Kraichnan [1]. Additionally, we examine how flow features
such as a mean flow affect the relationship between spatial, temporal, and
spatio-temporal spectral transfers. We also derive the effects of isotropic
sweeping on the temporal triad interactions, particularly with regard to the
locality of the triad interactions. Many of these results will be applicable
to dynamics in terms other than kinetic energy; we briefly discuss this in
Appendix F.

Our starting point is the spatially Fourier transformed two- or three-
dimensional Navier–Stokes equations defined by Eq. 7, although for the later
discussion of triad interactions we need only assume the existence of a
nonlinear advection term. The spatio-temporal spectral equation of motion
takes the form

$\displaystyle\partial_{\tau}\mathcal{E}(\bm{k}$
$\displaystyle,\omega,\tau;T)=\sum_{n}\mathcal{A}_{n}(\bm{k},\omega,\tau;T)$
$\displaystyle=\mathcal{N}(\bm{k},\omega,\tau;T)+\mathcal{D}(\bm{k},\omega,\tau;T)+\mathcal{F}(\bm{k},\omega,\tau;T),$
(43)

where

$\mathcal{E}(\bm{k},\omega,\tau;T)={\textstyle\frac{1}{2}}|\widehat{\widetilde{\bm{u}}}(\bm{k},\omega,\tau;T)|^{2}$
(44)

is the kinetic energy density;

$\mathcal{N}(\bm{k},\omega,\tau;T):=-\mbox{Re}\left[\widehat{\widetilde{\mathbf{u}}}^{*}\cdot\left(\widehat{\widetilde{\bm{\nabla}\Pi}}+\widehat{\widetilde{\bm{u}\cdot\bm{\nabla}\bm{u}}}\right)\right]$
(45)

is the transfer due to nonlinear advection;

$\mathcal{D}(\bm{k},\omega,\tau;T):=-\mbox{Re}\left[\sum_{m}\nu_{m}(-k^{2})^{m}\mathcal{E}(\bm{k},\omega,\tau;T)\right]$
(46)

is the transfer due to dissipation; and

$\mathcal{F}(\bm{k},\omega,\tau;T):=\mbox{Re}\left[\widehat{\widetilde{\bm{u}}}^{*}\cdot\widehat{\widetilde{\bm{f}}}\right]$
(47)

is the transfer due to forcing.

###  III.1 Doppler Shifting by a Mean Flow

We show how the spatio-temporal and temporal transfers change under a Galilean
boost, which approximates the imposition of a mean flow. Such a transformation
is important because many flows of practical interest either exhibit a strong
mean flow or are thought to exhibit isotropic sweeping, in which the smaller-
scale structures are swept without distortion by larger-scale structures. The
following calculation will not precisely model a mean flow when there is a
forcing term that sets a preferred frame of reference. In that case, there
could still be some regions of $(\bm{k},\omega)$-space where some of the
spectral transfers are independent of the spatio-temporal structure of the
forcing.

Suppose we have two frames of reference in which the Eulerian velocity fields
are $\bm{u}^{(1)}$ and $\bm{u}^{(2)}$, and that $\bm{V}$ is the (constant)
velocity of the second frame with respect to the first. The velocity fields
and external forcing in the two frames are then related by

$\displaystyle\bm{u}^{(2)}(\bm{x},t)=\bm{u}^{(1)}(\bm{x}-\bm{V}t,t)+\bm{V},$
(48) $\displaystyle\bm{f}^{(2)}(\bm{x},t)=\bm{f}^{(1)}(\bm{x}-\bm{V}t,t).$
(49)

As a means of setting notation, the energy budgets in frames $(1)$ and $(2)$
take the form

$\displaystyle\frac{\partial}{\partial\tau}\mathcal{E}^{(m)}(\bm{k},\omega,\tau)=\mathcal{N}^{(m)}(\bm{k},\omega,\tau)$
$\displaystyle+\mathcal{D}^{(m)}(\bm{k},\omega,\tau)$
$\displaystyle+\mathcal{F}^{(m)}(\bm{k},\omega,\tau),$ (50)

where $m=1,2$ labels the reference frame.

Using the STFT with or without detrending, it can be shown that each of the
terms in Eq. 50 transform according to

$\displaystyle\mathcal{A}^{(2)}(\bm{k},\omega,\tau)$
$\displaystyle=\mathcal{A}^{(1)}(\bm{k},\omega+\bm{V}\bm{\cdot}\bm{k},\tau),$
(51)

where $\mathcal{A}=\mathcal{N},\mathcal{D},\mathcal{F}$, or
$\partial_{\tau}\mathcal{E}$. The transformation of each term is simply a
$\bm{k}$-dependent shift in the frequency variable. The transfer in mode
$(\bm{k},\omega)$ in frame (2) corresponds with the transfer in mode
$(\bm{k},\omega+\bm{V}\cdot\bm{k})$ in frame (1).

After integrating over all spectral angles in $\bm{k}$-space in order to
obtain $\mathcal{A}(k,\omega,\tau)$, or after summing over all $\bm{k}$ to
obtain $\mathcal{A}(\omega,\tau)$, we obtain the transformation rules

$\displaystyle\mathcal{A}^{(2)}(k,\omega,\tau)$ $\displaystyle:=\int
d\Omega_{d}\,k^{d-1}\mathcal{A}^{(2)}(\bm{k},\omega,\tau)$ (52)
$\displaystyle=\int
d\Omega_{d}\,k^{d-1}\mathcal{A}^{(1)}(\bm{k},\omega+\bm{V}\bm{\cdot}\bm{k},\tau;T),$
(53) $\displaystyle\mathcal{A}^{(2)}(\omega,\tau)$
$\displaystyle:=\sum_{\bm{k}}\mathcal{A}^{(2)}(\bm{k},\omega,\tau)$ (54)
$\displaystyle=\sum_{\bm{k}}\mathcal{A}^{(1)}(\bm{k},\omega+\bm{V}\bm{\cdot}\bm{k},\tau)$
(55)

where the integrals in Eq. 52 and Eq. 53 are over all angles in
$d$-dimensional $\bm{k}$-space.

According to Eq. 53 and Eq. 55, after summing over multiple wavevectors
$\bm{k}$ the spectral content in mode $\omega$ in frame (2) corresponds to a
sum of the spectral content over a range of wavevector-frequency modes in
frame (1). In general, the transfers $\mathcal{A}^{(1)}(\bm{k},\omega,\tau)$
contribute substantially to the integral or sum only for wavenumbers
$k<k_{max}$ for some wavenumber $k_{max}>0$. With $\omega$ fixed, the
frequencies $\omega^{\prime}:=\omega+\bm{V}\cdot\bm{k}$ that contribute
substantially to the sum will lie in the range $\omega-
k_{max}V<\omega^{\prime}<\omega+k_{max}V$. Thus, the effect of a mean flow is
to “spread” out each frequency mode into nearby frequency–wavenumber modes. We
assert that this result, represented specifically by Eq. 53 and Eq. 55, is
identical to the spectral broadening discussed by Tennekes [13] and a
generalization of the broadening shown by Wilczek and Narita [62], who
implemented a Gaussian distributed sweeping velocity in addition to a mean
velocity to show how $\bm{V}\cdot\mathbf{k}$ acts as a Doppler shift for
frequency spectra.

###  III.2 Triad interactions

Spatial triad interactions, $\mathcal{T}(\bm{k},\bm{p},\bm{q})$, give the rate
of energy injection into wavevector mode $\bm{k}$ due to nonlinear interaction
with modes $\bm{p}$ and $\bm{q}$ [63, 64, 65, 1]. The sum of all the triad
interactions is the total transfer due to nonlinear advection. In this
section, we define spatio-temporal and temporal triads interactions using
asymmetric and symmetric diagrams defined with the STFT without detrending.
Asymmetric diagrams were introduced by Skitka  _et al._ [19] to describe the
individual energy exchanges between pairs of modes within a triad taken in
isolation. They have an associated scalar value of energy transfer. Asymmetric
diagrams describe a portion of the full interaction between three spectral
modes; in particular, they define a direction of energy transfer between two
distinct modes. Instead of asymmetric diagrams, one can also decompose triad
interactions into symmetric diagrams. Symmetric diagrams also have an
associated scalar value of energy transfer. They describe the total energy
injection into one mode from the two other modes in the triad interaction. We
will use the notation $\mathcal{T}_{u}$ to represent asymmetric diagrams and
$\mathcal{T}_{s}$ to represent symmetric diagrams. More explicitly, asymmetric
spatio-temporal diagrams,
$\mathcal{T}_{u}(\bm{k},\bm{p},\bm{q},\omega,\omega_{\bm{p}},\omega_{\bm{q}})$
give the rate of energy flow from wavevector-frequency mode
$(\mathbf{q},\omega_{q})$ into wavevector-frequency mode
$(\mathbf{k},\omega)$, where $\mathbf{p}=\mathbf{k}-\mathbf{q}$ and
$\omega_{p}=\omega-\omega_{q}$. Temporal diagrams
$\mathcal{T}_{u}(\omega,\omega_{\bm{p}},\omega_{\bm{q}})$ give the rate of
energy flow from frequency mode $\omega_{\bm{q}}$ to $\omega$ with
$\omega_{p}=\omega-\omega_{q}$. Symmetric spatio-temporal diagrams,
$\mathcal{T}_{s}(\bm{k},\bm{p},\bm{q},\omega,\omega_{\bm{p}},\omega_{\bm{q}})$,
are the rate of energy injection into wavevector-frequency mode
$(\bm{k},\omega)$ due to nonlinear interaction with modes
$(\bm{p},\omega_{\bm{p}})$ and $(\bm{q},\omega_{\bm{q}})$. Likewise, symmetric
temporal diagrams, $\mathcal{T}_{s}(\omega,\omega_{\bm{p}},\omega_{\bm{q}})$,
give the rate of energy injection into frequency mode $\omega$ due to
nonlinear interaction with modes $\omega_{\bm{p}}$ and $\omega_{\bm{q}}$. In
this paper, we do not focus on possible uses of asymmetric versus symmetric
diagrams. We instead use them to develop a potential method of applying tapers
to the nonlinear advection term. In equations where both $\mathcal{T}_{u}$ or
$\mathcal{T}_{s}$ can apply, we will leave off the subscript. We will show in
the following section that symmetric spatio-temporal diagrams can be
decomposed into the asymmetric spatio-temporal diagrams and symmetric temporal
diagrams can be decomposed into asymmetric temporal diagrams. Furthermore,
asymmetric and symmetric temporal diagrams can be derived from the asymmetric
and symmetric spatio-temporal diagrams, respectively.

If we use the STFT without a detrending operation, then we can apply the
convolution theorem to show that the spatio-temporal spectral transfer due to
the nonlinear advection term in Eq. 43 can be written as

$\displaystyle\mathcal{N}[\sigma]$
$\displaystyle(\bm{k},\omega,\tau;T):=\mbox{Re}[\widehat{\widetilde{\bm{u}}}^{*}[\sigma](\bm{k},\omega,\tau)\cdot(\widehat{\widetilde{\bm{\nabla}\Pi}}+\widehat{\widetilde{\bm{u}\cdot\bm{\nabla}\bm{u}}})[\sigma](\bm{k},\omega,\tau)]$
(56)
$\displaystyle=\sum_{\bm{p},\bm{q}}\int\>d\omega_{\bm{p}}d\omega_{\bm{q}}\mathcal{T}[\sigma_{1},\sigma_{2}](\bm{k},\bm{p},\bm{q},\omega,\omega_{\bm{p}},\omega_{\bm{q}},\tau;T)\delta_{\bm{k}-\bm{p}-\bm{q},\bm{0}}\delta(\omega-\omega_{\bm{p}}-\omega_{\bm{q}}).$
(57)

In this case, we choose to decompose the above equation in terms of diagrams,
setting $\mathcal{T}=\mathcal{T}_{u}$, defined explicitly by

$\displaystyle\mathcal{T}_{u}[\sigma_{1},\sigma_{2}](\bm{k},\bm{p},\bm{q},\omega,\omega_{\bm{p}},\omega_{\bm{q}},\tau):=\mbox{Im}\bigg{[}\Big{(}\widehat{\widetilde{\mathbf{u}}}^{*}[\sigma_{1}\sigma_{2}](\bm{k},\omega,\tau)\cdot\widehat{\widetilde{\mathbf{u}}}[\sigma_{1}](\bm{q},\omega_{\bm{q}},\tau)\Big{)}\Big{(}\mathbf{k}\cdot\widehat{\widetilde{\bm{u}}}[\sigma_{2}](\bm{p},\omega_{\bm{p}},\tau)\Big{)}\bigg{]}.\>\>$
(58)

We cannot get the above results with detrending because of the convolution
required to calculate Eq. 56. In the above definitions, $\sigma_{1}$ and
$\sigma_{2}$ are taper functions such that $\sigma_{1}\sigma_{2}=\sigma$,
where $\sigma$ is the taper function in Eq. 56. The constraint on taper
functions emerges due to the convolution necessary to calculate the STFT of
the nonlinear advection term. It is not a priori clear how to optimally choose
the tapers $\sigma_{1}$ and $\sigma_{2}$. The condition
$\sigma_{1}\sigma_{2}=\sigma$ implies that the only case where one can
trivially set $\sigma_{1}=\sigma_{2}$ is the case of applying a rectangular
function as a taper. If ringing is a potential problem, the two tapers must be
different. We suggest setting $\sigma_{2}$ as the rectangular taper function
(Eq. 20), which leaves $\sigma_{1}=\sigma$ such that $\sigma_{1}$ is the
user’s preferred taper function. This amounts to applying tapers to the two
modes whose energy exchange is described by the diagram. This choice is
appealing because it preserves the diagram’s antisymmetry in the
$(\mathbf{k},\omega)$ and $(\mathbf{q},\omega_{\mathbf{q}})$ modes. This
antisymmetry, identified and discussed in Appendix A of Skitka  _et al._ [19],
allows asymmetric diagrams to be interpreted as an energy exchange between
multiple modes. With the tapering in Eq. 58, one can compute such energy
exchanges with minimal error caused by ringing. Additionally, this definition
of tapered asymmetric diagrams allows us to define a tapered definition of
symmetric wavenumber-frequency diagrams, $\mathcal{T}=\mathcal{T}_{s}$, that
maintains their symmetry in $(\mathbf{p},\omega_{\mathbf{p}})$ and
$(\mathbf{q},\omega_{\mathbf{q}})$. Written out explicitly, we would then
define a symmetric wavenumber-frequency diagram as

$\displaystyle\mathcal{T}_{s}[\sigma_{1},\sigma_{2}](\bm{k},\bm{p},\bm{q},\omega,\omega_{\bm{p}},\omega_{\bm{q}},\tau):=\frac{1}{2}\Big{[}\mathcal{T}_{u}[\sigma_{1},\sigma_{2}](\bm{k},\bm{p},\bm{q},\omega,\omega_{\bm{p}},\omega_{\bm{q}},\tau)+\mathcal{T}_{u}[\sigma_{1},\sigma_{2}](\bm{k},\bm{q},\bm{p},$
$\displaystyle\omega,\omega_{\bm{q}},\omega_{\bm{p}},\tau)\Big{]}.$ (59)

Eq. 59, despite being defined in terms of diagrams, is equivalent to the
conventional definition of triads with tapered velocity fields [1, 63, 65, 38,
66].

The results of the rest of this paper pertaining to diagrams are independent
of $\tau$, $T$, and the choice of $\sigma$, so we leave these choices out of
the notation for triads.

Using Parseval’s theorem, we may obtain $\mathcal{T}(\bm{k},\bm{p},\bm{q})$
from
$\mathcal{T}(\bm{k},\bm{p},\bm{q},\omega,\omega_{\bm{p}},\omega_{\bm{q}})$ by
integrating over all frequencies $(\omega$, $\omega_{\bm{p}},$
$\omega_{\bm{q}})$. We may also obtain temporal diagrams by summing the
spatio-temporal diagrams over all wavevectors:

$\displaystyle\mathcal{T}(\omega,\omega_{\bm{p}},\omega_{\bm{q}}):=\sum_{\bm{k},\bm{p},\bm{q}}\mathcal{T}(\bm{k},\bm{p},\bm{q},\omega,\omega_{\bm{p}},\omega_{\bm{q}}).$
(60)

In two spatial dimensions we can similarly derive diagrams for enstrophy,
which are simply related to the energy diagrams:

$\displaystyle\mathcal{T}_{enstrophy}(\bm{k},\bm{p},\bm{q},\omega,\omega_{\bm{p}},\omega_{\bm{q}}):=k^{2}\mathcal{T}(\bm{k},\bm{p},\bm{q},\omega,\omega_{\bm{p}},\omega_{\bm{q}}).$
(61)

For two-dimensional turbulence, the spatial enstrophy transfers are
straightforwardly related to the spatial energy transfers
$\mathcal{T}(\bm{k},\bm{p},\bm{q})$:

$\displaystyle\mathcal{T}_{enstrophy}(\bm{k},\bm{p},\bm{q})$
$\displaystyle:=\sum_{\omega,\omega_{\bm{p}},\omega_{\bm{q}}}k^{2}\mathcal{T}(\bm{k},\bm{p},\bm{q},\omega,\omega_{\bm{p}},\omega_{\bm{q}})$
(62) $\displaystyle=k^{2}\mathcal{T}(\bm{k},\bm{p},\bm{q}).$ (63)

However, the temporal enstrophy transfers

$\displaystyle\mathcal{T}_{enstrophy}(\omega,\omega_{\bm{p}},\omega_{\bm{q}}):=\sum_{\bm{k},\bm{p},\bm{q}}k^{2}\mathcal{T}(\bm{k},\bm{p},\bm{q},\omega,\omega_{\bm{p}},\omega_{\bm{q}})$
(64)

have no clear relation to the temporal energy transfers
$\mathcal{T}(\omega,\omega_{\bm{p}},\omega_{\bm{q}})$. This lack of a
connection between $\mathcal{T}(\omega,\omega_{\bm{p}},\omega_{\bm{q}})$ and
$\mathcal{T}_{enstrophy}(\omega,\omega_{\bm{p}},\omega_{\bm{q}})$ suggests
that we cannot prove that energy and enstrophy must be transferred in opposite
directions in $\omega$-space in the same way that one proves that energy and
enstrophy must be transferred in opposite directions in $k$-space for two-
dimensional turbulence [1].

Symmetric temporal and spatio-temporal triads obey detailed conservation laws
similar to those of spatial triads [1]. In two and three dimensions, written
in terms of symmetric diagrams, detailed conservation of energy takes the form

$\displaystyle\mathcal{T}_{s}(\bm{k},\bm{p},\bm{q},\omega,\omega_{\bm{p}},\omega_{\bm{q}})+\mathcal{T}_{s}(\bm{p},\bm{q},\bm{k},\omega_{\bm{p}},\omega_{\bm{q}},\omega)+\mathcal{T}_{s}(\bm{q},\bm{k},\bm{p},\omega_{\bm{q}},\omega,\omega_{\bm{p}})$
$\displaystyle=0,$ (65)
$\displaystyle\mathcal{T}_{s}(\omega,\omega_{\bm{p}},\omega_{\bm{q}})+\mathcal{T}_{s}(\omega_{\bm{p}},\omega_{\bm{q}},\omega)+\mathcal{T}_{s}(\omega_{\bm{q}},\omega,\omega_{\bm{p}})$
$\displaystyle=0,$ (66)

for spatial triads $\bm{k}=\bm{p}+\bm{q}$ and temporal triads
$\omega=\omega_{\bm{p}}+\omega_{\bm{q}}$. Similarly, in two dimensions only,
written with symmetric diagrams, detailed conservation of enstrophy takes the
form

$\displaystyle
k^{2}\mathcal{T}_{s}(\bm{k},\bm{p},\bm{q},\omega,\omega_{\bm{p}},\omega_{\bm{q}})+p^{2}\mathcal{T}_{s}(\bm{p},\bm{q},\bm{k},\omega_{\bm{p}},\omega_{\bm{q}},\omega)+q^{2}\mathcal{T}_{s}(\bm{q},\bm{k},\bm{p},\omega_{\bm{q}},\omega,\omega_{\bm{p}})$
$\displaystyle=0,$ (67)
$\displaystyle\mathcal{T}_{s,enstrophy}(\omega,\omega_{\bm{p}},\omega_{\bm{q}})+\mathcal{T}_{s,enstrophy}(\omega_{\bm{p}},\omega_{\bm{q}},\omega)+\mathcal{T}_{s,enstrophy}(\omega_{\bm{q}},\omega,\omega_{\bm{p}})$
$\displaystyle=0,$ (68)

for symmetric spatial-temporal and temporal triads. These rules imply that
energy (and enstrophy if in two dimensions) is conserved within triads.

###  III.3 Locality of Triad interactions

In this section, we extend our discussion of mean flows and locality from Sec.
III.1 to triad interactions. The effect of a mean flow, approximated again by
a Galilean transformation, on the spatio-temporal diagrams is

$\displaystyle\mathcal{T}^{(2)}(\bm{k},\bm{p},\bm{q},\omega,\omega_{\bm{p}},\omega_{\bm{q}})=\mathcal{T}^{(1)}(\bm{k},\bm{p},\bm{q},\omega+\bm{V}\cdot\bm{k},\omega_{\bm{p}}+\bm{V}\cdot\bm{p},\omega_{\bm{q}}+\bm{V}\cdot\bm{q}),$
(69)

where $\bm{V}$ is the velocity of reference frame (2) in frame (1). The main
implication of Eq. 69 is that the imposition of a mean flow, or sweeping,
makes temporal triad interactions more nonlocal. This should be expected as it
has long been known [13] that in the presence of sweeping, Kolmogorov-type
arguments cannot be made for frequency spectra as they are made for wavenumber
spectra, at least at sufficiently high wavenumber.

The effect of isotropic sweeping on spectral transfers can be modeled by
considering a single sweeping velocity magnitude $V=|\bm{V}|$ and then
integrating Eq. 69 over all directions of $\bm{V}$. General radially symmetric
sweeping velocity distributions may then be handled by linear superposition.

Figure 2:  The effect of isotropic sweeping on spatio-temporal triad
interactions for three different spatial triads
$(\bm{K}_{1},\bm{K}_{2},\bm{K}_{3})$. Sweeping redistributes a triad
interaction at $(\Omega_{1},\Omega_{2},\Omega_{3})=(0,0,0)$ onto an ellipse in
the $\omega_{+}\omega_{-}$-plane, shown at the right. The projection of the
ellipse in the $\omega_{1}\omega_{2}$-plane is shown at the left. Arrows
indicate which spatial triad corresponds to which ellipse. Shaded regions
indicate which temporal triads are nonlocal, as defined by Eq. 74.

Suppose that the transfer in the absence of sweeping (i.e. in frame 1) takes
non-zero values only on the triad
$(\bm{K}_{1},\bm{K}_{2},\bm{K}_{3},\Omega_{1},\Omega_{2},\Omega_{3})$, and let
the sweeping velocity take the form

$\displaystyle\bm{V}=V\cos(\phi)\cos(\theta)\bm{\widehat{x}}+V\cos(\phi)\sin(\theta)\bm{\widehat{y}}+V\sin(\phi)\bm{\widehat{z}},$
(70)

where $\phi$ is the angle between $\bm{V}$ and the plane defined by the triad
$(\bm{K}_{1},\bm{K}_{2},\bm{K}_{3})$, and $\theta$ is the angle between
$\bm{K}_{1}$ and the projection of $\bm{V}$ onto that plane (making
$\bm{\widehat{x}}=\bm{K}_{1}/|\bm{K}_{1}|$). For consistency, we re-index the
frequencies in frame 2 with numbers, as in ($\omega_{1}$, $\omega_{2}$,
$\omega_{3}$). Restricting the analysis to two spatial dimensions, so that
$\phi=0$, the effect of integrating Eq. 69 over all $\theta$ (i.e. the effect
of isotropic sweeping by a single $V$) is to redistribute the transfer at
$(\Omega_{1},\Omega_{2},\Omega_{3})$ onto an ellipse in
$(\omega_{1},\omega_{2},\omega_{3})$-space centered at
$(\Omega_{1},\Omega_{2},\Omega_{3})$, with the ellipse lying in the plane
defined by $\omega_{\bm{k}}=\omega_{\bm{p}}+\omega_{\bm{q}}$. In that plane,
one may define an orthonormal coordinate system with components
$(\omega_{+},\omega_{-})$ given by
$\omega_{+}=\sqrt{3/2}(\omega_{2}+\omega_{3})=\sqrt{3/2}\omega_{1}$ and
$\omega_{-}=\sqrt{1/2}(\omega_{2}-\omega_{3})$.

Consider triads that satisfy $K_{2}=K_{3}$, where by definition
$K_{n}:=|\bm{K}_{n}|$. The more general case with arbitrary $\bm{K}_{2}$ and
$\bm{K}_{3}$ is not substantially more illuminating. Depending on the angle
between $\bm{K}_{2}$ and $\bm{K}_{3}$, the spatial triad
$(\bm{K}_{1},\bm{K}_{2},\bm{K}_{3})$ can be relatively local or nonlocal, with
maximal locality obtained when the wavevectors form an equilateral triangle
($K_{1}=K_{2}=K_{3}$).

The resulting ellipse in the $\omega_{+}\omega_{-}$-plane will take the form

$\frac{(\omega_{+}-\Omega_{+})^{2}}{(K_{+}V)^{2}}+\frac{(\omega_{-}-\Omega_{-})^{2}}{(K_{-}V)^{2}}=1,$
(71)

where $\Omega_{+}:=\sqrt{3/2}(\Omega_{2}+\Omega_{3})$ and
$\Omega_{-}:=\sqrt{1/2}(\Omega_{2}-\Omega_{3})$ and, similarly,
$\bm{K}_{+}:=\sqrt{3/2}(\bm{K}_{2}+\bm{K}_{3})$ and
$\bm{K}_{-}:=\sqrt{1/2}(\bm{K}_{2}-\bm{K}_{3})$. The axes of the ellipse will
be aligned with the $\omega_{+}$\- and $\omega_{-}$-axes, and the
corresponding major/minor axes will be $K_{+}V$ and $K_{-}V$. In terms of
$\bm{K}_{1}$ and $\bm{K}_{2}$, the major/minor axes can be written as

$\displaystyle K_{+}V$ $\displaystyle=$
$\displaystyle\sqrt{6}K_{3}|\cos(\theta_{23}/2)|V,$ (72) $\displaystyle
K_{-}V$ $\displaystyle=$ $\displaystyle\sqrt{2}K_{3}|\sin(\theta_{23}/2)|V,$
(73)

where $\theta_{23}$ is the angle between $\bm{K}_{2}$ and $\bm{K}_{3}$.

We use the locality measure originally defined by Lesieur [67] to quantify
locality of triads; any temporal triad that satisfies

$\gamma:=\frac{\max(|\omega_{1}|,|\omega_{2}|,|\omega_{3}|)}{\min(|\omega_{1}|,|\omega_{2}|,|\omega_{3}|)}<3$
(74)

is defined as local. Temporal triads are inherently more nonlocal than spatial
triads due to the fact that frequencies are scalars. While it is possible to
simultaneously satisfy $|\bm{k}|$=$|\bm{p}|$=$|\bm{q}|$ and
$\bm{k}=\bm{p}+\bm{q}$, it is not possible to satisfy both
$|\omega_{1}|=|\omega_{2}|=|\omega_{3}|$ and
$\omega_{1}=\omega_{2}+\omega_{3}$. Thus, the smallest possible value of
$\gamma$ for temporal triads is $\gamma=2$ rather than $\gamma=1$. Our choice
of threshold for nonlocality at $\gamma=3$ is natural because it divides the
$\omega_{+}\omega_{-}$-plane into equal areas (Fig. 2).

Fig. 2 demonstrates the effect of isotropic sweeping on three different
spatio-temporal triads with $(\Omega_{1},\Omega_{2},\Omega_{3})=(0,0,0)$. The
spatial triads for the three cases are $K_{1}=K_{2}=K_{3}$ (equilateral
triangle), $\bm{K}_{2}\approx\bm{K}_{3}$, and $\bm{K}_{2}\approx-\bm{K}_{3}$.
The figure shows the swept transfers in the $\omega_{+}\omega_{-}$-plane as
well as the projection of the swept transfers onto the
$\omega_{\bm{p}}\omega_{\bm{q}}$-plane. Shaded regions correspond to nonlocal
temporal triads, while unshaded regions correspond to local temporal triads,
according to Eq. 74.

According to Eq. 72 and Eq. 73, the major and minor axes will be equal (i.e.
the ellipse will be a circle) precisely when $\theta_{23}=\pm 2\pi/3$. In
other words, restricted to the case $K_{2}=K_{3}$, the ellipse will be a
circle if and only if the triad $(\bm{K}_{1},\bm{K}_{2},\bm{K}_{3})$ forms an
equilateral triangle. In contrast, the minor axis will be zero whenever
$\theta_{23}=0$ or $\theta_{23}=\pi$. When $\theta_{23}=0$ (i.e. $\bm{K}_{2}$
is parallel to $\bm{K}_{3}$) the major axis will lie along the
$\omega_{+}$-direction, which corresponds to local frequency triads. When
$\theta_{23}=\pi$ (i.e. $\bm{K}_{2}$ and $\bm{K}_{3}$ are antiparallel) the
major axis will lie along the $\omega_{-}$-direction, which corresponds to
nonlocal frequency triads.

Moreover, according to Eq. 72 and Eq. 73, the major/minor axes of the ellipses
in the $\omega_{+}\omega_{-}$-plane are directly proportional to the sweeping
velocity $V$ and the size $K$ of the spatial triad. Thus, for sufficiently
large $V$ or $K$ the ellipse will be large enough to encompass many regions of
locality and nonlocality in the $\omega_{+}\omega_{-}$-plane, causing any
originally local temporal triad to be distributed into both local and nonlocal
triads.

This simple example can be generalized to the case where $\Omega$’s are non-
zero if they are still approximately equivalent such that $\gamma\approx 1$
when defined on $\Omega$’s. If $\Omega$’s are non-zero, we can recast the
coordinate system of Fig. 2 as
$(\omega_{+}-\Omega_{+},\omega_{-}-\Omega_{-})$. This change in coordinates
would be accompanied subtracting $\Omega_{1}$, $\Omega_{2}$, and $\Omega_{3}$
from $\omega_{1}$, $\omega_{2}$, and $\omega_{3}$ respectively. If the
$\Omega$’s are nonlocal, shifting the $\omega_{+}\omega-$ coordinates system
could change the regions of locality and nonlocality relative to the origin of
such a coordinate plane. This would be the case where isotropic sweeping
localizes nonlocal spatio-temporal triads, as opposed to delocalizing local
spatio-temporal triads as shown in Fig. 2.

The three-dimensional case is a simple extension of the two-dimensional case
discussed above. In three dimensions the effect of integrating over all
$\theta$ is to form an ellipse in the $\omega_{+}\omega_{-}$-plane exactly as
in the two-dimensional case, except that $V$ is replaced by $V\cos(\phi)$. So,
upon integrating over all $\phi$, the interior of the circle is filled as
well, albeit non-uniformly with lowest density near the center of the circle.

##  IV Numerical investigation: 2D turbulence

In this section, we numerically simulate the two-dimensional Navier–Stokes
equations with a general dissipation and a forcing term. We use the spatio-
temporal spectral transfer diagnostic developed in Sec. II.2 to investigate
the dynamics of energy and enstrophy in statistically equilibrated simulations
forced at four different frequencies. To test the robustness of this
diagnostic, we also show the effects of varying the window size and sampling
rate on the numerically calculated spectral transfers, and we demonstrate the
effect of detrending when the simulation output is non-stationary.

###  IV.1 Numerical setup

In spectral space, assuming spatially periodic boundary conditions, the
simulated equation takes the following form:

$\frac{\partial}{\partial
t}k^{2}\widetilde{\psi}(\bm{k},t)=\widetilde{J(\psi,\nabla^{2}\psi)}(\bm{k},t)-\widetilde{D}(\bm{k},t)-\widetilde{F}(\bm{k},t),$
(75)

where the forcing term $\widetilde{F}$ and the dissipation term
$\widetilde{D}$ will be described below. The width of the square, periodic
domain is $L$.

Figure 3: Properties of the forcing, which is narrowband in both frequency and
wavenumber: (a) Temporal spectral densities $\mathcal{S}_{F}(\omega)$ of the
four different forcing frequencies (top) and area-preserving spectral
densities $\omega\mathcal{S}_{F}(\omega)$ (bottom). (b) Snapshot of the
forcing $F(\bm{x},t_{0})$ at initial time (top) and the corresponding area-
preserving spatial spectral density $k\mathcal{E}_{F}(k)$ (bottom), which was
the same for all four forcing frequencies. The scales of the vertical axes are
arbitrary.

The forcing term $\widetilde{F}(\bm{k},t)$ is narrowband in both wavenumber
and frequency. The forcing is statistically spatially isotropic with peak
spectral amplitude at the forcing wavenumber $k_{F}$ and is statistically
stationary with peak spectral amplitude at the forcing frequency $\omega_{F}$.
Fig. 3 shows the temporal and spatial spectral densities of the forcing for
four values of $\omega_{F}$. Fig. 3 also shows a snapshot of the forcing in
physical space for comparison with snapshots of stream function and vorticity
shown in Fig. 4. A full description of the forcing is provided in Appendix E.

For dissipation at the small and large scales, we apply two cutoff wavenumber
filters similar in form to hyperviscosity and inverse viscosity. Specifically,
we set

$\widetilde{D}(\bm{k},t):=\nu_{6}(k-k_{6})^{6}H(k-k_{6})k^{2}\widetilde{\psi}(\bm{k},t)+\nu_{-6}(k^{-1}-k_{-6}^{-1})^{6}H(k_{-6}-k)k^{2}\widetilde{\psi}(\bm{k},t),$
(76)

where $H(\cdot)$ is the Heaviside step function. The first term in Eq. 76 can
be thought of as a hyperviscosity applied to wavenumbers $k>k_{6}$, while the
second term in can be thought of as an inverse viscosity applied to
wavenumbers $k<k_{-6}$. This choice ensures a wide inertial range between
$k_{-6}<k<k_{6}$. We set $\nu_{6}=1.3\times 10^{-11}$, $k_{6}=366\approx
0.715k_{N}$, $\nu_{-6}=1.2\times 10^{9}$ , and $k_{-6}=7.9\approx
0.0154k_{N}$, where $k_{N}$ is the Nyquist wavenumber. The cut-off wavenumber
$k_{-6}$ for the inverse wavenumber filter was set sufficiently high so that
the stream function would not be significantly spatially correlated with
itself halfway across the domain. This ensures that the lowest frequencies
observed in the system are not appreciably affected by the spatial
periodicity.

To simulate Eq. 75 we use a parallelized pseudo-spectral method [68, 69], with
domain size $L=2\pi$ and spatial resolution of $1024^{2}$. The vorticity is
iterated using third-order Adams–Bashforth time-stepping, while the forcing
term is iterated separately using Euler’s method with a smaller time step.
While technically the inclusion of a stochastic forcing term changes the
strong (root-mean-square) convergence of the numerical scheme from third order
to first order, the error is better represented by $O(\Delta t^{3}+\Delta
t/M)$, where $M$ is the ratio of the vorticity-term time step to the forcing-
term time step. This effectively third-order numerical method (for
sufficiently large $M$) resembles numerical methods developed for stochastic
systems with small noise [70, 71]. We found that $M=1$ was sufficient for the
convergence of all statistics we present. The time-stepping increment was
either $\Delta t=2^{-12}$ or $2^{-13}$ for all runs.

Most runs were not fully de-aliased, but the wavenumber filter removed
practically all enstrophy at wavenumbers larger than $5/6$ times the Nyquist
wavenumber. A single fully de-aliased simulation was conducted for comparison,
and there was no significant change in the results.

We simulated the Eq. 75 for four different choices of the forcing frequency,
$\omega_{F}/2\pi\in\\{0.01,0.1,1,10\\}$, with corresponding integral time
scales of $\tau_{F}=5\times 2\pi/\omega_{F}$. For each of the four forcing
frequencies, we created three time series (each approximately one terabyte in
size) by sampling the simulation at three different sampling rates. Each time
series consisted of $32,768$ snapshots, and the three time series durations
were $t_{max}\in\\{256,1024,4096\\}$. For each time series the window size
used in calculating the spectral transfers was equal to the entire length of
the time series. Thus, the calculated spectral transfers took the form
$\mathcal{A}(\bm{k},\omega,t_{max}/2;t_{max})$, where the central time was set
equal to the center of the data set, $\tau=t_{max}/2$.

Figure 4:  Stream function and vorticity snapshots: (a) Snapshots of vorticity
$\nabla^{2}\psi(\bm{x},t_{0})$ after statistical equilibrium is reached for
three of the four forcing frequencies. Magnified insets show regions with and
without strong vortices. The vorticity snapshot for the missing forcing
frequency ($\omega_{F}/2\pi=0.1$) is visually and quantitatively similar to
the case $\omega_{F}/2\pi=1$. Only relative scale is important, so each plot
uses the same (but otherwise arbitrary) units. Note the decreased vorticity
range in the rightmost plots ($\omega_{F}/2\pi=10$). (b) Snapshots of stream
function $\psi(\bm{x},t_{0})$ with corresponding insets.  Figure 5:  Spatial
spectral fluxes of energy and enstrophy (top) and the corresponding spatial
spectral transfers (bottom) for the $\omega_{F}/2\pi=1$ simulation. A positive
transfer at wavenumber $k$ corresponds to the injection of energy (or
enstrophy) at that wavenumber. A positive flux corresponds to a downscale
transfer (direct cascade) to larger $k$, and a negative flux corresponds to an
upscale transfer (inverse cascade) to smaller $k$. A negatively sloping flux
corresponds to a positive transfer, while a positively sloping flux
corresponds to a negative transfer. The statistical noise in the fluxes and
transfers is very small because they have been averaged over $32,768$
snapshots.

###  IV.2 Spatial transfers and fluxes

In this section, we calculate the spectral transfers and fluxes. The spectral
fluxes capture all the contributions to the kinetic energy evolution at all
wavenumbers or frequencies above a threshold. Historically, this quantity has
been defined in wavenumber space either as [14]

$\Pi_{\mathcal{A}_{n}}^{>}(k,t):=\int_{k}^{\infty}dk^{\prime}\;\mathcal{A}_{n}(k^{\prime},t),$
(77)

or by the time average of the above equation,

$\Pi_{\mathcal{A}_{n}}^{>}(k):=\int_{k}^{\infty}dk^{\prime}\;\mathcal{A}_{n}(k^{\prime}),$
(78)

where $\mathcal{A}_{n}(k)$ is the time-average of $\mathcal{A}_{n}(k,t)$.
Alternatively, with the development of temporal spectral transfers, we can
define a temporal spectral flux:

$\Pi_{\mathcal{A}_{n}}^{>}(\omega,\tau;T):=\int_{\omega}^{\infty}d\omega^{\prime}\;\mathcal{A}_{n}(\omega^{\prime},\tau;T).$
(79)

While the theory presented earlier in this paper focuses primarily on spectral
transfers, we show both transfers and fluxes because some readers will be more
familiar with spectral fluxes. Fig. 5 shows the time-averaged spatial spectral
transfers and fluxes for the case $\omega_{F}/2\pi=1$. The spatial spectral
transfers and fluxes for the other three forcing frequencies are qualitatively
similar, with the main difference being that the run forced at the highest
frequency exhibited transfers and fluxes of smaller magnitude. Some readers
may be accustomed to significantly noisier spatial spectral transfers. Our
transfers are relatively smooth because we average over 32,768 snapshots in
time.

In the bottom left plot of Fig. 5, showing energy transfers, we see the
transfer by forcing (lavender-grey curve, labelled “F”) balanced out by a
corresponding negative transfer by nonlinear advection (dark purple curve,
labelled “N”) at the forcing wavenumber. We see a transfer of this energy to a
lower wavenumber, where the inverse wavenumber filter (yellow curve, labelled
“D2” on the plot and labelled “inv filter” in the legend) has a negative
transfer. The upper left plot, showing the corresponding energy fluxes
reflects the integrals of the transfer curves. The negative values of the
nonlinear advection curve reflects a transfer from high wavenumber to low
wavenumber. The forcing curve is positive since it represents an injection of
energy and the inverse filter is negative because it represents the
dissipation of energy from the system. The bottom right plot shows enstrophy
transfers, where we see similar behavior at the forcing wavenumber. However,
enstrophy is injected at a higher wavenumber by nonlinear advection. There,
the enstrophy is dissipated by the wavenumber filter (orange curve, labelled
“D1” on the plot and “filter” in the legend). This transfer at higher
wavenumber is captured in the top right plot of Fig. 5 by the positive
nonlinear advection curve. These results are expected for two-dimensional
turbulence. Energy is injected at a forcing wavenumber and is dissipated at
small wavenumber (large spatial scales). Enstrophy, on the other hand, is
injected at a forcing wavenumber and is dissipated at large wavenumber (small
spatial scales).

###  IV.3 Temporal spectral fluxes

We study temporal spectral fluxes of energy and enstrophy for all four forcing
frequencies in Fig. 6. We focus on spectral fluxes because the corresponding
spectral transfers are significantly noisier and potentially more difficult to
interpret. To illustrate this point, we also show the temporal spectral
transfers of energy for the highest two forcing frequencies in Fig. 6.

Fluxes are shown for each of three sampling rates, indicated by different line
styles. Solid lines correspond to the largest sampling rate , dashed lines
correspond to the intermediate sampling rate, and dotted lines correspond to
the lowest sampling rate. Sampling rates are defined by a fixed number of
snapshots (32,678) taken over a time duration $t_{max}$; thus, the largest
sampling rate corresponds to the smallest $t_{max}$ and vice-versa. There is
close agreement between the results with the different sampling rates, except
at the highest and lowest frequencies resolved by each time series. The close
agreement indicates that the temporal spectral transfers and fluxes are
reasonably robust diagnostics for time series that resolve most but not all
dynamical time scales.

In contrast with the spatial spectral fluxes, the temporal spectral fluxes
differ considerably between simulations with different forcing frequencies.
Just as with wavenumber fluxes, we use the nonlinear advection term’s (dark
purple curves, labelled “N”) fluxes as an indicator of the overall flow of
energy and enstrophy for the fluid. For temporal fluxes, this term reflects
the flow of a quantity from time scale to time scale [65, 19]. At the lowest
forcing frequency, $\omega_{F}/2\pi=0.01$, the energy flux is positive,
indicating a net transfer of energy to frequencies higher than the forcing
frequency. When $\omega_{F}/2\pi=1$, the energy flux is negative, indicating a
net transfer of energy to frequencies below the forcing frequency. In the
intermediate case, when $\omega_{F}/2\pi=0.1$, there are both upscale and
downscale transfers of energy. In all three of these cases, the energy is
injected at the forcing frequency and then transferred to the range of
frequencies over which large-scale dissipation occurs. We observe this in the
the energy and enstrophy flux plots for the three lowest forcing frequencies,
where the flux by forcing (lavender grey plot, labelled “F”), sharply
decreases at the forcing frequency. In the energy transfer plot (left most
column, top plot), this is reflected in the sharp peak in the transfers by
forcing. At the highest forcing frequency, $\omega_{F}/2\pi=10$, energy and
enstrophy are not injected at the forcing frequency. In the energy transfer
plot at the bottom left, we see both positive and negative peaks in the
transfers by forcing and then a surge of forcing amplitude at a range of lower
frequencies. In the middle plot of the bottom row, we observe that the energy
flux by forcing decreases to zero at frequencies almost an order of magnitude
lower than the forcing frequency. The same behavior is visible in the
enstrophy flux, in the bottom right plot. We interpret this as the forcing
frequency being too high for the fluid to respond at the forcing frequency.
Instead, energy and enstrophy are injected across a wider range of frequencies
below the forcing frequency. From there, energy is transferred to relatively
lower frequencies while enstrophy is transferred to relatively higher
frequencies. As with fluxes and transfers in wavenumber, the wavenumber filter
is responsible for dissipating the enstrophy (here, at high frequencies) while
the inverse wavenumber is responsible for dissipating the energy (here, at
frequencies higher or lower than the forcing frequency). Our results show that
while the direction of energy flux varies depending on the frequency of the
forcing, enstrophy appears to solely move to higher frequencies. In other
words, the dissipation of enstrophy always occurs at frequencies above the
injection frequency, which is evidenced by the positive nonlinear flux for all
cases. In the simulation with the highest forcing frequency, the fluid cannot
respond to the dominant time-scale of the forcing. Turbulence still develops
but at a less energetic level, as was previously found in other studies of
“modulated” turbulence in three-dimensions [47, 48, 49, 50, 51, 52]. Because
the fluid cannot respond to arbitrarily high frequencies (without changing the
forcing amplitude), there is an upper limit for the frequency at which
enstrophy can be injected. If this limiting injection frequency is lower than
the highest frequencies associated with the enstrophy dissipation range, then
there cannot be a transfer of enstrophy to frequencies below the injection
frequency. Thus, while energy can be transferred either to smaller or larger
timescales relative to the period of the forcing, enstrophy appears to only
move to small timescales in two-dimensional turbulence.

###  IV.4 Temporal spectral transfers

In Fig. 7, we show the temporal spectral transfers corresponding to the fluxes
in Fig. 6). We reduce the range of the y-axis, which corresponds to the scale
of the energy or enstrophy, to display fine-scale features of each transfer.

Figure 6:  Temporal spectral fluxes of energy and enstrophy for each of the
four forcing frequencies, and temporal spectral transfers of energy for the
two highest forcing frequencies. Each plot contains three curves, one solid,
one dashed, and one dotted. The curves vary by duration $t_{max}$, with longer
durations using lower sampling rates. The natural frequency $\omega_{nat}$ and
external frequency $\omega_{ext}$ are indicated by dashed lines. For each of
the four forcing frequencies (rows) temporal spectral _fluxes_ of energy
(middle column) and enstrophy (right column). For the two highest forcing
frequencies, the temporal spectral _transfers_ are shown as well (left
column), corresponding to the spectral fluxes directly to their right. The
upper legend defines each term in the spectral budget. The lower legend
indicates values of $t_{max}$ and thus the simulation sampling rate. The two
frequencies derived _post-priori_ (see Appendix G), $\omega_{diss}/2\pi$ and
$\omega_{nat}/2\pi$ are indicated by gray vertical dash-dotted lines. For all
plots, $\omega_{diss}/2\pi<\omega_{nat}/2\pi$.

We use the longest duration dataset for the lowest frequencies, the shortest
duration dataset for the highest frequencies, and the intermediate duration
dataset for the intermediate frequencies. This choice removes the noise due to
aliasing at the highest frequencies and removes incorrect transfers at the
lowest frequencies. When the fluid is forced at $\omega_{F}/2\pi=1$, energy
and enstrophy are injected within a narrow range around the forcing frequency.
When the fluid is forced at $\omega_{F}/2\pi=10$, there appear to be both
positive and negative transfers at the forcing frequency, but these average
roughly to zero as demonstrated by the nearly zero flux at the same frequency
in Fig. 6. Instead, both energy and enstrophy are injected over a relatively
wide range of frequencies centered at approximately $\omega/2\pi=0.4$.

Figure 7:  Temporal spectral transfers of energy and enstrophy for the
simulations forced at the two highest forcing frequencies, $\omega_{F}/2\pi=1$
and $\omega_{F}/2\pi=10$. Each plot is created using the same three sampling
rates as in Fig. 6, but the noisiest regions of each curve have been removed.
In the $\omega_{F}/2\pi=10$ simulation, the transfers near $\omega/2\pi=10$
are noisy but integrate roughly to zero, which is proved by the smoothness of
the fluxes near $\omega/2\pi=10$. The two frequencies derived _post-priori_
(see Appendix G), $\omega_{diss}/2\pi$ and $\omega_{nat}/2\pi$ are indicated
by gray vertical dash-dotted lines. The gaps in the bottom left plot are due
to differently sampled curves that do not overlap at that frequency. The curve
with $t_{max}=1024$ ends before the gap and the dashed plots with longer
$t_{max}=4096$ continues to lower frequencies.

###  IV.5 Spatio-temporal spectral transfers

Fig. 8 shows the full spatio-temporal spectral transfers for all four forcing
frequencies. The top plot combines the results for the three lowest forcing
frequencies. The combination of these three cases is possible because the
transfers due to dissipation remain largely unchanged while the transfers due
to the forcing do not overlap. The transfers due to the nonlinear term are not
shown, because they can be calculated using the fact that the sum of all the
transfers is zero everywhere in $(k,\omega)$-space (the residual is
effectively zero). The vertical smearing at fixed wavenumber in the bottom
plot of Fig. 8 again shows that energy injection takes place at a range of
frequencies well below the intended forcing frequency.

We show three colorbars, each for different sources of transfers. We have an
light orange colorbar for the inverse wavenumber filter that ranges from $-1$
to $0$, a lavender grey colorbar for the forcing that ranges from $0$ to $1$,
and an darker orange colorbar for the wavenumber filter ranging from $0$ to
$1$. All of the colorbars are white where their corresponding field is $0$ and
black where their corresponding field has a magnitude (absolute value) of $1$
or higher. We see black spots in the figure where that field is valued at $1$
or significantly higher. These are especially clear in the top and bottom
plots in at $0.1k_{F}$, where we observe the light orange streak denoting the
inverse wavenumber filter (labelled “energy dissipation” in the plots). The
spatio-temporal spectral transfers for these simulations are fairly
straightforward to understand. When the forcing frequency is not too high,
energy and enstrophy are injected at the forcing frequency, and the
dissipation ranges are not much affected by changes in the forcing frequency.
This allows energy to be dissipated at frequencies lower or higher than the
forcing frequency. When the forcing frequency is high, energy and enstrophy
are instead injected at relatively lower frequencies.

Figure 8: Spatio-temporal spectral transfers for the three lowest forcing
frequencies (top) and the highest forcing frequency (bottom). To reduce
statistical noise, the lower half of each plot is based on the simulation
output with the longest duration($t_{max}=4096$), while the upper half of each
plot is based on the simulation output with the shortest duration
($t_{max}=256$). Transfers of energy and enstrophy do not substantially
overlap in $(k,\omega)$-space except at the forcing frequencies, where we plot
energy transfers only. Dissipation transfers were quantitatively similar for
the three lowest forcing frequencies, so we simply show the result for
$\omega_{F}/2\pi=1$ in the upper plot. In the bottom plot, the symbol “x”
marks the location of $(k_{F},\omega_{F})$, which is notably outside the range
of frequencies in which energy (and enstrophy) is injected by the forcing. The
nonlinear term (not shown) is approximately the negative of the sum of the
other terms.

###  IV.6 Evolution of temporal transfers

Fig. 9 shows how the fluxes due to nonlinear advection of energy and enstrophy
evolve as the central time $\tau$ advances during spin-up of the simulations.
For each plot, each curve represents nonlinear advection at a distinct central
time $\tau$. Purple plots are earlier central times and orange plots are the
later central times. For the simulation with $\omega_{F}/2\pi=1$, energy and
enstrophy are injected at the forcing frequency. From there, energy is
transferred to lower frequencies while enstrophy is transferred to higher
frequencies. As the simulation spins up, the amount of energy transferred to
the lowest frequencies increases somewhat (and similarly for enstrophy at the
highest frequencies). For the simulation with $\omega_{F}/2\pi=10$, the
injection frequency range for both energy and enstrophy initially develops
over a range of frequencies well below the forcing frequency. As the
simulation evolves, the injection frequency range shifts to increasingly
higher frequencies. The dissipation ranges of energy and enstrophy also evolve
over time, with the energy dissipation range moving to lower frequencies and
the enstrophy dissipation range moving to higher frequencies. This evolution
of the dissipation ranges corresponds to the development and widening of the
spatial cascades of energy and enstrophy.

Figure 9:  Temporal spectral fluxes of the nonlinear advection of energy and
enstrophy during spin-up for the $\omega_{F}/2\pi=1$ and $\omega_{F}/2\pi=10$
simulations. The fluxes were calculated for sixteen central times
$\tau=\tau_{1}<\tau_{2}<\cdots<\tau_{16}$ using contiguous but non-overlapping
windows. The window durations were $T=64$ for the energy fluxes and $T=16$ for
the enstrophy fluxes. In both simulations, energy fluxes move to lower
frequencies while enstrophy fluxes move to higher frequencies as the system
evolves.

###  IV.7 Effects of varying the window size and sampling rate

In practice, the lowest and/or highest frequencies relevant to the dynamics of
a fluid might not be resolved by a data set, due to insufficient time series
duration or sampling rate. To aid in the interpretation of temporal fluxes
used as a diagnostic for real data, we present in Fig. 10 a systematic study
of the effect of temporal resolution on the resulting temporal fluxes for our
two-dimensional simulations. The investigation reveals that the temporal
fluxes typically give consistent results as the temporal resolution varies,
but there are specific cases where the diagnostic fails to give meaningful
results.

Figure 10:  The effect of window size ($T$) and sampling rate on temporal
spectral fluxes of energy and enstrophy. Each column corresponds to fluxes
calculated from the same simulation with the same forcing frequency. All
fluxes were calculated by sampling 512 snapshots in time from the same
simulation. For each column, the bottom set of fluxes samples the ensemble of
simulations the least frequently while the top set of fluxes samples the
ensemble the most. For each vertical axis, the lower-leftmost set of fluxes is
plotted correctly centered at zero, while all other fluxes are shifted
vertically for ease of comparing fluxes. The rightward tilt occurs because the
difference in sampling between plots causes a translation in the output
frequency scales of the STFT. The STFT of the simulation output with the
lowest sampling rate returns fluxes on the lowest range of frequencies, while
the STFT of the simulation output with the highest sampling rate returns
fluxes on the highest range of frequencies.

We show the temporal spectral fluxes $\Pi_{\mathcal{A}}(\bm{k},\omega,\tau;T)$
of energy and enstrophy for a wide range of window sizes $T$ and a wide range
of sampling rates. The window sizes were selected from
$T\in\\{2^{12},2^{10},2^{8},2^{6},2^{4}\\}$ and the sampling rate was selected
such that the total number of snapshots was $512$ for each data set. Each
column corresponds to one of the four simulations (indicated by its forcing
frequency). In each column, the curves in the lowest plot are sampled the
least frequently while the curves in the highest plot are sampled most
frequently. Each $y$-axis shows the scale for the lower-leftmost plot on that
axis. The remaining plots have been shifted vertically so that they do not all
overlap with each other. Each column “leans” to the right, because the
sampling rate increases from bottom to top while the number of samples remains
constant. To save space, we overlapped the graphs for the four forcing
frequencies, with the graphs for the lowest three forcing frequencies (the
left three columns) sharing the same vertical scale. Considering first the
energy dissipation by the inverse wavenumber filter (yellow curves, labelled
“D2”), we find that the temporal spectral fluxes are reasonably accurate down
to the lowest resolved frequencies. For simulation output that are too
infrequently sampled, the temporal spectral fluxes become inaccurate only at
the very highest resolved frequencies, due to aliasing of the unresolved
higher-frequency modes into lower modes. Fluxes due to the dissipation of
enstrophy by the wavenumber filter (orange curves, labelled “D1”) are also
accurate at low frequencies, but the aliasing effect at higher frequencies is
more apparent. Despite the aliasing, the total transfer (given by the flux at
the lowest resolved frequency) of energy and enstrophy remains roughly the
same throughout the simulation output, indicating that aliasing effects are
localized to the highest resolved frequencies. Energy dissipation by the
wavenumber filter and enstrophy dissipation by the inverse wavenumber filter
are negligible for all sampling rates and forcing frequencies. Considering the
forcing term (lavender grey curves, labelled “F”) for both energy and
enstrophy, we see that the calculated fluxes are correct if and only if the
forcing frequency is not too close to the edge of the interval of resolved
frequencies. For example, when the forcing frequency is far too low to be
resolved by the simulation output with the shortest duration, or when the
forcing frequency is far too high to be resolved by the sampling rate, the
calculated fluxes develop inaccuracies. When the forcing frequency is near the
edge of the interval of resolved frequencies, then the flux due to the forcing
tends to decrease to $0$ at too-low frequencies, particularly as especially
evident in the case for the energy and enstrophy flux when $\omega_{F}/2\pi=1$
and $T=2^{10}$ (second set of fluxes from the bottom). When the forcing
frequency is slightly lower than the smallest resolved frequency, especially
evident in the case $\omega_{F}/2\pi=0.01$ and $T=2^{8}$ (third from the
bottom), one sees the forcing flux becomes nonzero at higher than accurate
frequencies. For most of the simulations, energy and enstrophy budgets have
near zero residual flux. We can then determine the effect of the window size
and sampling rate on the temporal flux due to the nonlinear advection term
(dark purple curves, labelled ”nonlinear”) by combining the effects on the
forcing and dissipation terms. We note that the nonlinear flux term tends to
have good agreement between different sampling rates, _except_ when the
forcing frequency is barely _un_ resolved by the sampling rate, as
characterized by the energy fluxes when $\omega_{F}/2\pi=1$ and $T=2^{10}$ or
by the enstrophy fluxes when $\omega_{F}/2\pi=1$ and $T=2^{10}$ or $T=2^{8}$.
Those cases give particularly inaccurate temporal spectral fluxes because
neither the forcing frequency nor the enstrophy dissipation range is resolved
by the sampling rate. These problematic cases may be representative of real-
world time series (for large systems, like the ocean, where the shortest time
scales are not resolved by data). However, the fluxes are inaccurate only at
the highest resolved frequencies, while the fluxes at the remaining resolved
frequencies seem to be reliable. For both energy and enstrophy, we see
slightly non-zero residual fluxes across a wide range of frequencies when
$\omega_{F}/2\pi=10$ for each of the sampling rates. We also see this for low
frequencies in the simulations with the highest sampling rate where
$\omega_{F}/2\pi=0.01$ and $1$. Deviations in nonlinear advection between
ensembles in that column are explained by these slightly non-zero residual
fluxes. We expect the residuals to be non-zero for the instances of highest
sampling because these highly sampled simulation outputs do not capture lower
frequencies. As a result, aliasing in forcing could lead to a noticeable
residual, particularly in the cases where $\omega_{F}/2\pi=0.01$ and $1$ and
at the highest sampling rate (top of the first and third columns from the
left) because the range of frequencies is slightly higher than the forcing
frequencies. A similar thing appears to occur in the farthest right columns
where $\omega_{F}/2\pi=10$. For the bottom four plots where
$\omega_{F}/2\pi=10$, corresponding to the 4 lowest sampling rates, the
residual is likely non-zero where the range of frequencies does not include
either the frequency of forcing or the effective frequency of forcing observed
in Fig. 6 (bottom left corner plot), which is the frequency at which the
simulation responded to the forcing.

###  IV.8 Effects of varying the detrending method

For most of our simulations we calculate temporal spectral transfers both with
and without a temporal detrending operation (specifically, the removal of a
linear trend). In previous sections, we reported only on results that did not
include a detrending operation. For time series with sufficiently long
duration (Fig. 6, Fig. 7 and Fig. 8), because the simulations are
statistically stationary, the trends are effectively zero and the detrending
operation has on the temporal fluxes is only seen at the lowest two or three
resolved frequencies.

To study the effect of the detrending operation on the lowest resolved
frequencies, we specifically looked at time series that had durations
significantly shorter than the forcing period. The relatively low-frequency
forcing could be interpreted as the source for a long-term trend. In Fig. 11,
we analyze temporal spectral fluxes calculated with or without detrending for
eight independent time series of length $t_{max}=16$ that were forced at the
longest forcing period, $2\pi/\omega_{F}=100$. We retain the same color scheme
as previous plots: dark purple curves show nonlinear advection, lavender grey
curves show forcing, orange curves show wavenumber filter, yellow curves show
inverse wavenumber filter, and black dotted curves show residuals. For
comparison, the “true” values of the temporal spectral fluxes at the lowest
frequencies, determined from a much longer time series, are indicated by
arrows. Detrending brought temporal fluxes for the forcing and inverse
wavenumber filter into agreement with the true values for those fluxes. The
flux due to the nonlinear term, however, was not made more accurate by the
detrending operation. We note that the residual is significant here only at
lower frequencies because short duration datasets cannot capture the low
frequency behavior that the full duration dataset contains.

Figure 11:  Temporal spectral fluxes of energy with and without detrending.
Each plot shows temporal spectral fluxes for eight independent time series of
length $t_{max}=16$ obtained from the lowest frequency
($\omega_{F}/2\pi=0.01$) simulation. The low frequency forcing effectively
provides a trend for these short duration time series. The arrows indicate the
“true” values of the fluxes obtained from a much longer duration time series.
Agreements with the “true” fluxes, verified numerically, are indicated by
check marks ($\checkmark$).

##  V Conclusion

We have shown how temporal and spatio-temporal spectral transfers in fluid
turbulence may be defined quite generally in terms of time-frequency analysis
methods such as the short-time Fourier transform or the wavelet transform.
Moreover, these methods can be modified to include a fairly general detrending
operation in such a way that the temporal and spatio-temporal spectral budgets
remain exact even after detrending.

The interpretation of temporal and spatio-temporal spectral transfers is
similar to the interpretation of spatial-spectral transfers, with the main
difference being the replacement of time derivative with respect to $t$ by the
time derivative with respect to central time $\tau$. While the spatial
spectral transfer $\mathcal{A}(\bm{k},t)$ gives the $t$ rate of change of the
energy density $\mathcal{E}(\bm{k},t)$ in wavenumber mode $\bm{k}$, the
spatio-temporal spectral transfer $\mathcal{A}(\bm{k},\omega,\tau)$ gives the
$\tau$ rate of change of energy density $\mathcal{E}(\bm{k},\omega,\tau)$ in
wavenumber-frequency mode $(\bm{k},\omega)$. Similarly, the temporal spectral
transfer $\mathcal{A}(\omega,\tau)$ gives the $\tau$ rate of change of the
spectral density $\mathcal{E}(\omega,\tau)$ in frequency mode $\omega$.

We derived various theoretical properties of temporal and spatio-temporal
spectral transfers, including the effect of a mean flow. We showed that we can
decompose the Fourier transform of the nonlinear advection term into diagrams,
which contain information about the direction of energy transfer, and triads,
which fully capture the energetic interaction between any three scales. The
diagrams are introduced in this paper to provide a method of tapering the
components of the nonlinear advection term. They allow us to develop a
nonlinear turbulent triad definition where ringing can be mitigated without
changing the useful properties of the triad. Temporal and spatio-temporal
transfers of kinetic energy and enstrophy also satisfy detailed conservation
laws. We use metrics from turbulence studies of wavenumber triads to quantify
the effect of sweeping by a mean flow on the locality of the triad
interactions in the frequency domain [3, 67]. Much of the theory from this
section can be used outside of kinetic energy transfers and fluxes, as we show
in Appendix F. The issue of locality versus nonlocality of triad interactions
is particularly important to the interpretation of temporal and spatio-
temporal fluxes. One might be inclined to think of spectral fluxes as arising
from a local operator, but in wavenumber, frequency, or wavenumber-frequency
space, a flux is almost always a nonlocal function of the spectral quantity of
interest. In three dimensions, this operator is expected to be somewhat narrow
in $k$, and in two dimensions a bit wider [65]. However, we have shown how the
operator may become significantly wider in $\omega$ due to sweeping by a mean
flow. Thus, one should keep in mind the relatively greater degree of
nonlocality in the frequency domain when diagnosing spectral transfers and
fluxes. Studying the relationship between energy and enstrophy spatio-temporal
transfers, one cannot show that energy and enstrophy temporal fluxes
consistently flow in opposite directions (i.e. the dual cascade) as was shown
for spatial fluxes by Kraichnan [1]. Moreover, one can show that the presence
of a mean flow and high wavenumber modes can change the locality of triad
interactions. Local interactions can be delocalized due to a constant sweeping
velocity. For these reasons, the existence of a dual cascade in frequency
space seems less likely than its existence in wavenumber space. However, the
temporal and spatio-temporal triad interactions should not be discounted as
being fully determined by sweeping, because only the triad interactions
involving high wavenumbers are significantly affected by sweeping, and the
low-wavenumber, low-frequency modes of the flow are also of dynamical
interest. To illustrate the use of temporal and spatio-temporal spectral
transfers and fluxes, we applied the diagnostic to the output of numerical
simulations of two-dimensional turbulence. These simulations show how the
frequency affects the transfers of energy and enstrophy in both wavenumber and
frequency space. We also used the simulations to examine how the diagnostic is
affected by insufficient window size, insufficient sampling rate, and the use
of detrending, which are important for practical applications. The main
physical result of the numerical investigation is that temporal energy due to
the nonlinear advection term can move to larger or shorter timescales relative
to the scale of injection. In other words, it is possible to inject energy at
a wide range of frequencies, including frequencies above and below the energy
dissipation range. The lowest forcing frequency case ($\omega_{F}/2\pi=0.01$)
is particularly interesting, as it shows energy may be injected at very low
frequencies and then transferred to higher dissipation frequencies, indicated
by positive fluxes in Fig. 6. This result may be related to the positive
temporal energy flux calculated by Arbic  _et al._ [15, 16] for one region of
the ocean. In that one region, the time scale associated with energy injection
by baroclinic instability is much longer than the eddy time scale,
significantly more so than in the other regions considered. However, there may
be other explanations, such as inadequate spatio-temporal resolution of the
data, for the positive temporal energy flux observed in Arbic  _et al._ [15,
16]. We have shown here that a positive flux of energy in two dimensions is
possible when energy is injected at low frequencies. Enstrophy seems to only
transfer to shorter timescales.

We observed that if the forcing frequency is increased beyond some threshold,
the energy and enstrophy injection occurs at some lower frequency that remains
below the high frequency dissipation range. It is yet unclear what determines
this threshold for either the energy or the enstrophy.

Unlike the temporal fluxes, the spatial energy and enstrophy fluxes are
minimally affected by changing the forcing frequency, except for a decrease in
the magnitude of the flux for high forcing frequency. At the highest forcing
frequency, the decrease in flux magnitude is due to the decrease of the total
energy and enstrophy injected into the system.

The main practical result of the numerical investigation is that the temporal
and spatio-temporal spectral transfers and fluxes for finite time series are
quite robust to limitations typical of realistic data. Changes to window size,
sampling rate, and detrending method usually have little effect on the
temporal spectral transfers or fluxes except at marginal frequencies. For
short duration time series, temporal spectral fluxes tend to be accurate even
when the lowest frequencies in the system are not resolved by the data. When
the highest frequencies are not resolved, there is aliasing into the highest
resolved frequencies, but the temporal spectral fluxes at the remaining lower
frequencies remain accurate.

We found that one must be particularly careful in the special case where the
forcing frequency is close to – either above or below – the highest resolved
frequency. In this case, the temporal spectral fluxes become particularly
untrustworthy at nearby frequencies. The worst effect, fluxes having the wrong
sign, was seen when the forcing frequency was slightly too high to be resolved
and the dissipation range was also not fully resolved. In this situation, one
should be skeptical of the highest resolved frequencies in any case.

Our investigation also showed that when the forcing frequency is too low to be
resolved, a linear detrending operation may increase the accuracy of the
temporal spectral transfers at low frequencies. This effectively tells us that
detrending temporal transfers can be an effective method of dealing with
nonstationarity in time series. The accuracy was improved for all of the
temporal spectral transfers except for the transfer due to nonlinear
advection. The accuracy of the nonlinear transfer was possibly worsened by
detrending. For this reason, it would be reasonable to diagnose temporal
spectral transfers both with and without detrending. Future theoretical work
is needed to explain why the nonlinear transfers could be made less accurate
by detrending.

We find that our STFT-derived spatio-temporal transfers provide a relatively
robust diagnostic for studying spatial and temporal scales in fluid dynamics.
While the diagnostic has been used primarily in idealized simulations in the
past, we expect this method to be useful for studying more realistic, high-
resolution simulations and observational data. In particular, Skitka  _et al._
[19] demonstrates the usefulness of spatio-temporal diagrams in precisely
tracking the source and destination scales of energy transfer in simulations.
Given the broad applicability of spatio-temporal transfers and triads, we
expect that the Skitka  _et al._ [19]’s use of unsymmetric triads can be
generalized to studies of cascades outside of kinetic energy, such as the
temperature variance studied by Martin  _et al._ [24]. We discuss such
diagrams in Appendix F.

We briefly apply our diagnostic to studying the effect of sweeping by a
constant flow on spatio-temporal transfers. However, we expect the spatio-
temporal diagnostic to be useful in probing models of isotropic random
sweeping, discussed by Tennekes [13], Wilczek and Narita [62], Chen and
Kraichnan [72]. Future work can include theoretical and numerical research
that probes sweeping on kinetic energy transfers. Pending questions about the
method aside, this work has demonstrated the potential for broad application
of this diagnostic to fluids in a variety of contexts. We anticipate that this
method can be used for numerical studies of cascades in kinetic energy and
other quantities of interest in studies related to oceanography, meteorology,
and magnetohydrodynamics.

##  VI Acknowledgements

A.J.M. thanks Robert Deegan, John Boyd, and especially the late Charles
Doering for providing feedback on early versions of this research when that
material appeared in A.J.M.’s dissertation. A.J.M. and B.K.A. gratefully
acknowledge Mike Messina for support of the numerical simulations, which were
carried out on the supercomputer provided by the University of Michigan Office
of Research Cyberinfrastructure (ORCI). A.M., A.J.M. and B.K.A. acknowledge
support from National Science Foundation (NSF) Grant Nos. OCE-0960820,
OCE-1351837 and OCE-1851164 and from University of Michigan faculty startup
funds. G.R.F. received support from NSF Grant No. OCE-0960826.

##  VII Author Contributions

B.K.A., R.B.S., G.R.F., and A.J.M. conceived and designed the project. A.J.M.,
A.M., and J.S. developed the theory. A.M. prepared the final manuscript with
guidance from B.K.A. and A.J.M., based on an earlier manuscript that A.J.M.
prepared with guidance from B.K.A., G.R.F., and R.B.S. The numerical
experiments were carried out by A.J.M. with a codebase developed by G.R.F.
R.B.S. conducted this work within the France 2030 framework program, the
Centre Henri Lebesgue ANR-11-LABX-0020-01.

##  Appendix A Spatio-temporal Phase Evolution

While our focus is on spectral transfers and therefore on taking the real part
of the product of Eq. 22 and $\widehat{g}^{*}$, we provide insight into the
imaginary part as well. Multiplying Eq. 22 and taking the imaginary part gives

$-{\textstyle\frac{1}{2}}i\left(\widehat{g}^{*}\partial_{\tau}\widehat{g}-\widehat{g}\partial_{\tau}\widehat{g}^{*}\right)+(1-\alpha)\omega|\widehat{g}|^{2}=\mbox{Im}[\widehat{g}^{*}\widehat{f}].$
(80)

If we write $\widehat{g}$ in terms of its real amplitude and phase,

$\widehat{g}(\omega,\tau)=g_{0}(\omega,\tau)e^{i\phi(\omega,\tau)},$ (81)

then Eq. 80 can be rewritten as

$\partial_{\tau}\phi+(1-\alpha)\omega=\mbox{Im}[\widehat{g}^{*}\widehat{f}]/|\widehat{g}|^{2}$
(82)

when $|\widehat{g}|\neq 0$. Thus, taking the imaginary part tells us about the
$\tau$ rate of change of the phase of $\widehat{g}$. In the simpler situation
where only spatial spectra of $g$ are calculated, the imaginary part tells us
about the $t$ rate of change of the phase of $\widehat{g}$; this can be shown
with a near identical derivation to the one above.

Eq. 82 takes different forms depending on the choice of $\alpha$, because
different basis functions give different meaning to the phase. When $\alpha=0$
the phase is measured with respect to $t=0$, and when $\alpha=1$ the phase is
measured with respect to $t=\tau$.

Consider the case $\alpha=0$, which is the case where the basis functions
$e^{-i\omega t}$ do not move with the window. The evolution of the phase is

$\partial_{\tau}\phi=-\omega+\mbox{Im}[\widehat{g}^{*}\widehat{f}]/|\widehat{g}|^{2}.\;\;\;\;\;\;\;\;\;\;\text{(case
{\hbox{\alpha=0}})}$ (83)

If the window is very long, then the phase of $\widehat{g}$ measured with
respect to $t=0$ will change very little as the window advances. In other
words, $\partial_{\tau}\phi\approx 0$, so

$0\approx\partial_{\tau}\phi=-\omega+\mbox{Im}[\widehat{g}^{*}\widehat{f}]/|\widehat{g}|^{2}.\;\;\;\;\;\;\;\;\;\;\text{(case
{\hbox{\alpha=0}})}$ (84)

We get a balance of terms. The first term is simply $-\omega$, and the
remaining terms must approximately add up to $\omega$ so that the total is
approximately zero.

In the $\alpha=0$ case described above, the presence of the $\omega$ term may
seem somewhat mysterious. Why should the rate of change of the phase always
include an $\omega$ term that must be balanced by the remaining terms? The
situation is clarified by considering the $\alpha=1$ case, where the basis
functions $e^{-i\omega(t-\tau)}$ remain centered with the moving window. In
this case, the evolution of the phase is simply

$\partial_{\tau}\phi=\mbox{Im}[\widehat{g}^{*}\widehat{f}]/|\widehat{g}|^{2}.\;\;\;\;\;\;\;\;\;\;\text{(case
{\hbox{\alpha=1}})}$ (85)

If the window is very long, then the phase of $\widehat{g}$ measured with
respect to $t=\tau$ must change approximately at the rate $\omega$ as the
window advances in order to compensate for the moving window. In other words,
$\partial_{\tau}\phi\approx\omega$, so

$\omega\approx\partial_{\tau}\phi=\mbox{Im}[\widehat{g}^{*}\widehat{f}]/|\widehat{g}|^{2}.\;\;\;\;\;\;\;\;\;\;\text{(case
{\hbox{\alpha=1}})}$ (86)

By changing the basis functions, we have eliminated the $\omega$ term, instead
incorporating it into $\partial_{\tau}\phi$. As in the $\alpha=0$ case, the
remaining terms still approximately add up to $\omega$. The choice of $\alpha$
simply changes the meaning of $\partial_{\tau}\phi$ and helps to clarify the
presence of the $\omega$ term. The remaining terms are invariant concerning
the choice of $\alpha$.

To further illustrate the meaning of Eq. 82, let us consider the Rossby wave
example. We can begin by thinking purely about the spatial phase evolution
before we think about spatio-temporal phase evolution. The Rossby wave
equation is

$\partial_{t}\nabla^{2}\psi(\mathbf{x},t)=-\beta\partial_{x}\psi(\mathbf{x},t)$
(87)

where $\psi(\mathbf{x},t)$ is the stream function for a two-dimensional
incompressible fluid. Applying a spatial Fourier transform to Eq. 87 and
dividing by $-k$ gives

$k\partial_{t}\widetilde{\psi}(\mathbf{k},t)=i\beta\frac{k_{x}}{k}\widetilde{\psi}(\mathbf{k},t),$
(88)

where $k_{x}:=\hat{\mathbf{x}}\cdot\mathbf{k}$ is the component of
$\mathbf{k}$ in the x-direction. Eq. 88 can be written in the form of Eq. 2 by
letting ${\bf\it\chi}=\mathbf{k}$ and

$g({\bf\it\chi},t)=k\widetilde{\psi}(\mathbf{k},t).$ (89)

The general solution to Eq. 88 is a sum of linear Rossby waves. The Rossby
wave solutions take the form

$\widetilde{\psi}(\mathbf{k},t)=\Psi(\mathbf{k})e^{-i(\omega_{R}t+\theta_{0})}$
(90)

where $\psi(\mathbf{k})$ and $\theta_{0}$ and $\omega_{R}$ satisfies the
Rossby dispersion relation,

$\omega_{R}:=-\beta\frac{k_{x}}{k^{2}}.$ (91)

Each Rossby wave has constant amplitude $\Psi(\mathbf{k})$, so the energy in
each mode $\mathbf{k}$ does not change with time. Accordingly, the spectral
budget given by Eq. 13 is simply

$\partialderivative{t}\frac{1}{2}|g({\bf\it k},t)|^{2}=0$ (92)

The corresponding phase evolution would be given by

$\partialderivative{t}\phi(\mathbf{k},t)=\beta\frac{k_{x}}{k^{2}}$ (93)

Noting that the phase is $\phi(\mathbf{k},t)=-\omega_{R}t-\theta_{0}$, we see
that the phase evolution recovers the Rossby dispersion. Thus, the phase
evolution simply gives the contribution to the dispersion relation by each
term in the equation of motion.

Now, if we turn our attention to spatio-temporal phase evolution, we can take
Eq. 90 and apply an STFT. In the limit of an infinitely long time window, the
STFT applied to Eq. 90 gives

$\widehat{\widetilde{\psi}}(\bm{k},\omega,\tau)=2\pi\delta(\omega-\omega_{R})\Psi(\bm{k})e^{-i\theta_{0}}e^{-i\alpha\omega\tau},$
(94)

where $\delta(\cdot)$ is the Dirac delta function. With $g({\bf\it\chi},t)$
given by Eq. 89, this phase is

$\phi(\omega,\tau)=-\alpha\omega\tau-\theta_{0}$ (95)

when $\omega=\omega_{R}$ and is undefined when $\omega\neq\omega_{R}$. If
$\alpha=0$ then $\partial_{\tau}\phi=0$, whereas $\alpha=1$ then
$\partial_{\tau}\phi=\omega$. Substituting Eq. 95 into Eq. 94 gives the Rossby
dispersion relation $\omega=\omega_{R}$ for either choice of $\alpha$. Thus,
the contributions to the $\tau$ rate of change of the phase given by the
right-hand side of Eq. 82 once again tell us the contributions of the
dispersion relation by various terms in $f$.

##  Appendix B Wiener-Khinchin Theorem with STFT

In applied mathematics, the Wiener-Khinchin (or Wiener-Khintchine) theorem
states that the Fourier transform of the cross-correlation of two signals is
equivalent to their cross-spectrum if the signals are wide-sense stationary
[60, 61, 28]. In Chiu [25], frequency transfers for atmospheric kinetic energy
are calculated by applying a Fourier transform to the covariance of the
velocity and the time derivative of the velocity. By the Wiener-Khinchin
theorem, this quantity should be equivalent to $\partial_{t}\mathcal{E}$ in a
situation where a valid Fourier transform of the velocity exists. However, we
can use the STFT to formalize the method of Chiu [25] such that it can be
applied to practical applications. We start with the definition of the cross-
correlation,

$C(\theta)=\lim_{M\to\infty}\frac{1}{M}\int\limits_{0}^{M}dt\;\;\xi^{*}(t)\zeta(t+\theta),$
(96)

where we assume some arbitrary time signals $\xi$ and $\zeta$. $\theta$ here
is the time lag of a covariance. $M$ represents the time scale over which we
average the product of the time signals. However, this equation is not
appropriate in numerical applications because no real time series is
infinitely long. To make this theorem applicable to non-stationary data or
simulation output, we first need to redefine the cross-correlation:

$C[\sigma](\theta,\tau)=\lim_{M\to\infty}\frac{1}{M}\int\limits_{0}^{M}dt\;\;\xi^{*}(t,\tau)\zeta(t+\theta,\tau+\theta)\sigma(t-\tau)^{2}.$
(97)

Tapering allows us to take the integral from $0$ to $\infty$ even if we do not
realistically know what’s happening beyond a finite time limit in the data or
simulation output. We introduce the square bracket notation here for the
function of time since we are tapering the product of the two time series. The
taper for the $\theta$-shifted $\zeta$ is not shifted in time by $\theta$
because we also shift the central time $\tau$ by $\theta$. The shifts in both
$t$ and $\tau$ cancel out, allowing us to write $\sigma^{2}$. We shift both
$t$ and $\tau$ by $\theta$ so that both $\xi$ and $\zeta$ are the same size.
One can detrend both $\xi$ and $\zeta$ to deal with non-stationarity; we leave
this out of the notation for readability.

In the standard Wiener-Khinchin theorem proof, we would introduce the inverse
Fourier transform of $\xi$ and $\zeta$. Naively, one may expect to apply an
inverse STFT here to make a Wiener-Khinchin theorem for numerical
applications. However, an inverse STFT would remove the tapering introduced in
Eqn. 97. Instead, we suggest the following equation:

$\displaystyle\zeta[\sigma](t,\tau;T,\alpha)=\frac{1}{2\pi}\int\limits_{-\infty}^{\infty}d\omega\;\widehat{\zeta}[\sigma](\omega,\tau;T,\alpha)e^{i\omega(t-\alpha\tau)}.$
(98)

This inverse Fourier transform relates the tapered $\zeta$ in frequency space
to the tapered $\zeta$ in time. By maintaining the tapering, we can now
analytically evaluate integrals over the entire real domain, which is
necessary for calculating a numerically applicable Wiener-Khinchin theorem.

With Eq. 97 and Eq. 98, we can now calculate the relationship between tapered
transfers and cross-correlations. We can plug in the complex conjugate of Eq.
98 and the analogous inverse STFT of $\zeta$ into Eq. 96 and use the integral
definition of the Dirac delta function to derive an STFT version of the
Wiener-Khinchin theorem:

$\displaystyle C(\theta,\tau;T,\alpha)$
$\displaystyle=\lim_{M\to\infty}\frac{1}{M}\int\limits_{0}^{M}dt\;\;\xi^{*}[\sigma^{i}](t,\tau;T,\alpha)\;\zeta[\sigma^{i}](t+\theta,\tau+\theta;T,\alpha)$
(99)
$\displaystyle=\lim_{M\to\infty}\frac{1}{M}\int\limits_{0}^{M}dt\int\limits_{-\infty}^{\infty}d\omega\int\limits_{-\infty}^{\infty}d\omega_{2}\;\frac{\widehat{\xi^{*}}[\sigma](\omega,\tau;T,\alpha)}{2\pi}\;\frac{\widehat{\zeta}[\sigma](\omega_{2},\tau+\theta;T,\alpha)}{2\pi}e^{i\omega(t-\alpha\tau)}e^{-i\omega_{2}(t-\alpha\tau)}e^{-i\omega_{2}\theta(1-\alpha)}$
(100)
$\displaystyle=\int\limits_{-\infty}^{\infty}d\omega\int\limits_{-\infty}^{\infty}d\omega_{2}\Bigg{[}\lim_{M\to\infty}\frac{1}{M}\int\limits_{0}^{M}dt\;e^{i(\omega-\omega_{2})t}\Bigg{]}\frac{\widehat{\xi^{*}}[\sigma](\omega,\tau;T,\alpha)}{2\pi}\;\frac{\widehat{\zeta}[\sigma](\omega_{2},\tau+\theta;T,\alpha)}{2\pi}e^{-i\alpha\tau(\omega-\omega_{2})}e^{-i\omega_{2}\theta(1-\alpha)}$
(101)
$\displaystyle=\int\limits_{-\infty}^{\infty}d\omega\int\limits_{-\infty}^{\infty}d\omega_{2}\;2\pi\delta(\omega-\omega_{2})\frac{\widehat{\xi^{*}}[\sigma](\omega,\tau;T,\alpha)}{2\pi}\;\frac{\widehat{\zeta}[\sigma](\omega,\tau+\theta;T,\alpha)}{2\pi}e^{-i\alpha\tau(\omega-\omega_{2})}e^{-i\omega\theta(1-\alpha)}$
(102)
$\displaystyle=\int\limits_{-\infty}^{\infty}d\omega\frac{\widehat{\xi^{*}}[\sigma](\omega,\tau;T,\alpha)\;\widehat{\zeta}[\sigma](\omega,\tau+\theta;T,\alpha)}{2\pi}\;e^{-i\omega\theta(1-\alpha)}$
(103)

The result is the inverse Fourier transform of the cross-spectrum of tapered
time series $\xi$ and $\zeta$, where $\theta$ is the conjugate time variable
to $\omega$. Thus, we have proven that the STFT of the cross-correlation gives
you the cross-spectrum. We chose to start from a cross-correlation between
time series of the same size by shifting the central time of $\zeta$ to
$T+\theta$. We can demonstrate that this version of Wiener-Khinchin theorem is
necessary to use Chiu’s auto-correlation method of deriving kinetic energy
budgets on data or simulation output. Without the adjustments due to the STFT,
the budgets calculated by Chiu [25] suffer from the same problems as other
spectral budgets calculated in the frequency domain without adjusting for
convergence at $t=-\infty$ and $t=\infty$ and non-stationarity of time series.
We can calculate the cross-correlation of $g(t,\tau;T,\alpha)$ and
$\partial_{t}g(t,\tau;T,\alpha)$:

$\displaystyle\expectationvalue{g^{*}(t,\tau;T,\alpha)\partialderivative{g}{t}}$
$\displaystyle=\lim_{T\to\infty}\int\limits_{0}^{T}dt\;g^{*}(t,\tau;T,\alpha)\partialderivative{t}g(t,\tau;T,\alpha)$
(105)
$\displaystyle=\frac{1}{2\pi}\int\limits_{\omega=-\infty}^{\omega=\infty}d\omega\;\widehat{g}^{*}(\omega,\tau;T,\alpha)\bigg{[}\partialderivative{\tau}+i\omega(1-\alpha)\bigg{]}\widehat{g}(\omega,\tau;T,\alpha)e^{-i\omega\theta(1-\alpha)}.$
(106)

If we set $g=\mathbf{u}$, set $\theta=0$, and take the real part of the above
equation, we can show the derivation of the $\tau$-derivative of kinetic
energy:

$\displaystyle\mbox{Re}\bigg{[}\expectationvalue{\mathbf{u}^{*}(t,\tau;T,\alpha)\cdot\partialderivative{\mathbf{u}}{t}}\bigg{]}$
$\displaystyle=\frac{1}{2}\expectationvalue{\partialderivative{t}|\mathbf{u}(t,\tau;T,\alpha)|^{2}}$
(107)
$\displaystyle=\frac{1}{2\pi}\int\limits_{\omega=-\infty}^{\omega=\infty}d\omega$
$\displaystyle\;\partialderivative{\tau}\frac{1}{2}|\widehat{\mathbf{u}}(\omega,\tau;T,\alpha)|^{2}$
(108)
$\displaystyle=\frac{1}{2\pi}\int\limits_{\omega=-\infty}^{\omega=\infty}d\omega$
$\displaystyle\;\partialderivative{\tau}\mathcal{E}$ (109)

One can repeat this above calculation for every other term in kinetic energy
dynamics to calculate the full kinetic energy budget using this modified
Wiener-Khinchin theorem.

##  Appendix C Spectral budget for wavelets

Another commonly used bilinear time-frequency representation is the scalogram,
which is the modulus squared of the wavelet transform (WT). The continuous
wavelet transform of $f(t)$ is defined by

$\mbox{WT}_{f}(\omega,\tau):=\int_{-\infty}^{\infty}{f(t)|\omega|^{1/2}W((t-\tau)\omega)dt},$
(110)

where $W(\cdot)$ is the “mother wavelet,” which can be chosen according to
desired properties. While we focus mainly on the STFT, we show here that
analagous results hold for the wavelet transform.

We show that the form of the spectral budget (with and without detrending)
remains the same if we replace the short time Fourier transform with the
wavelet transform. In this appendix we redefine the “hat” operation to be the
wavelet transform:

$\widehat{f}(\omega,\tau):=WT_{f}(\omega,\tau),$ (111)

where $WT_{f}$ is defined by Eq. 110. Applying the WT to Eq. 12 gives

$\int_{-\infty}^{\infty}{(\partial_{t}g(t))|\omega|^{1/2}W((t-\tau)\omega)dt}=\widehat{f}(\omega,\tau).$
(112)

Integrating by parts and transforming $t$-derivatives into $\tau$-derivatives
gives

$\partial_{\tau}\widehat{g}(\omega,\tau)=\widehat{f}(\omega,\tau).$ (113)

Multiplying Eq. 113 by $\widehat{g}^{*}$ and then taking the real part gives

${\textstyle\frac{1}{2}}\partial_{\tau}|\widehat{g}(\omega,\tau)|^{2}=\mbox{Re}[\widehat{g}^{*}\widehat{f}_{i}],$
(114)

which has the same form as the spectral budget Eq. 24 based on the STFT with
$\alpha=1$.

In practice, time series are typically not detrended before applying a wavelet
transform. However, if we were to apply detrending then we would obtain the
same result as for the STFT with $\alpha=1$:

$\displaystyle\widehat{\left(\partial_{t}g\right)_{detrend}}(\omega,\tau)=\widehat{\left[(\partial_{\tau}+\partial_{t})g_{detrend}(t,\tau)\right]}(\omega,\tau)$
(115)
$\displaystyle=\int_{-\infty}^{\infty}dt\;|\omega|^{1/2}W((t-\tau)\omega)(\partial_{\tau}+\partial_{t})g_{detrend}(t,\tau)$
(116)
$\displaystyle=\partial_{\tau}\int_{-\infty}^{\infty}dt\;|\omega|^{1/2}W((t-\tau)\omega)g_{detrend}(t,\tau)$
(117) $\displaystyle=\partial_{\tau}\widehat{g_{detrend}}(\omega,\tau).$ (118)

The resulting spectral budget would take the same form as Eq. 41.

##  Appendix D Detrending

To perform a linear least squares fit we minimize the integral of the squared
residuals,

$S(\tau):=\int_{\tau-T/2}^{\tau+T/2}dt\left|f(t)-\sum_{n}c_{n}(\tau)\varphi_{n}(t-\tau)\right|^{2}.$
(119)

Changing the integration variable to $t^{\prime}:=t-\tau$ gives

$S(\tau)=\int_{-T/2}^{T/2}dt^{\prime}\left|f(t^{\prime}+\tau)-\sum_{n}c_{n}(\tau)\varphi_{n}(t^{\prime})\right|^{2}.$
(120)

To find $\\{c_{n}\\}$ that minimize $S$ we solve $\partial_{c_{n}}S=0$. The
result is

$\innerproduct{f(t^{\prime}+\tau)}{f(t^{\prime}+\tau)}{\varphi_{n}^{*}(t^{\prime})}=\sum_{m}c_{m}(\tau)\innerproduct{\varphi_{n}^{*}(t^{\prime})}{\varphi_{m}(t^{\prime})}_{t^{\prime}}.$
(121)

Define the matrix

$M_{nm}:=\innerproduct{\varphi_{n}^{*}(t^{\prime})}{\varphi_{m}(t^{\prime})}_{t^{\prime}},$
(122)

and note that it is independent of both $t$ and $\tau$ and only depends on the
form of the basis functions. Eq. 121 can then be solved for $c_{n}(\tau)$
using the matrix inverse of $M$, giving the result in the main text, Eq. 31.
Using Eq. 32 to calculate $(\dot{g})_{trend}$, $\partial_{t}g_{trend}$ and
$\partial_{\tau}g_{trend}$ we obtain

$(\dot{g})_{trend}(t,\tau)=\sum_{n}\varphi_{n}(t-\tau)\sum_{m}(M^{-1})_{nm}\innerproduct{\dot{g}(t^{\prime}+\tau)}{\dot{g}(t^{\prime}+\tau)}{\varphi_{m}^{*}(t^{\prime})}_{t^{\prime}}.$
(123)
$\partial_{t}g_{trend}(t,\tau)=\sum_{n}\dot{\varphi}_{n}(t-\tau)\sum_{m}(M^{-1})_{nm}\innerproduct{g(t^{\prime}+\tau)}{\varphi_{m}^{*}(t^{\prime})}_{t^{\prime}}.$
(124)

$\displaystyle\partial_{\tau}g_{trend}$
$\displaystyle(t,\tau)=-\sum_{n}\dot{\varphi}_{n}(t-\tau)\sum_{m}(M^{-1})_{nm}\innerproduct{g_{(}t^{\prime}+\tau)}{g_{(}t^{\prime}+\tau)}{\varphi_{m}^{*}(t^{\prime})}_{t^{\prime}}+\sum_{n}\varphi_{n}(t-\tau)\sum_{m}(M^{-1})_{nm}\innerproduct{\dot{g}(t^{\prime}+\tau)}{\varphi_{m}^{*}(t^{\prime})}_{t^{\prime}}.$
(125)

Combining these three results gives Eq. 34.

##  Appendix E Numerical forcing

We construct a statistically stationary, statistically isotropic forcing that
is narrowband in both wavenumber and frequency. We begin with defining a
deterministic forcing amplitude in wavenumber space:

$\displaystyle
f_{0}(\bm{k}):=\Bigg{\\{}\begin{array}[]{cc}\>\Big{[}(k_{F}+\Delta
k_{F})^{2}-|\bm{k}|^{2}\Big{]}\>\>\Big{[}|\bm{k}|^{2}-(k_{F}-\Delta