We introduce new features of data-adaptive harmonic decomposition (DAHD) that are showcased to characterize spatiotemporal variability in high-dimensional datasets of complex and mutsicale oceanic flows, offering new perspectives and novel insights. First, we present a didactic example with synthetic data for identification of coherent oceanic waves embedded in high amplitude noise. Then, DAHD is applied to analyze turbulent oceanic flows simulated by the Regional Oceanic Modeling System and an eddy-resolving three-layer quasigeostrophic ocean model, where resulting spectra exhibit a thin line capturing nearly all the energy at a given temporal frequency and showing well-defined scaling behavior across frequencies. DAHD thus permits sparse representation of complex, multiscale, and chaotic dynamics by a relatively few data-inferred spatial patterns evolving with simple temporal dynamics, namely, oscillating harmonically in time at a given single frequency. The detection of this low-rank behavior is facilitated by an eigendecomposition of the Hermitian cross-spectral matrix and resulting eigenvectors that represent an orthonormal set of global spatiotemporal modes associated with a specific temporal frequency, which in turn allows to rank these modes by their captured energy and across frequencies, and allow accurate space-time reconstruction. Furthermore, by using a correlogram estimator of the Hermitian cross-spectral density matrix, DAHD is both closely related and distinctly different from the spectral proper orthogonal decomposition that relies on Welch’s periodogram as its estimator method.

The turbulent oceanic flows consist of ubiquitous complex motions—jets, vortices, and waves—that co-exist on very different spatiotemporal scales but also without a clear scale separation, and it brings natural challenge to characterize the whole complexity across the scales. In particular, the study of temporal scales has got less attention than of spatial scales. To that effect, we offer fresh perspectives and novel insights by introducing new features of data-adaptive harmonic decomposition (DAHD) that are applied to analyze complex high-dimensional spatiotemporal datasets of oceanic flows, including a synthetic example of identifying coherent oceanic waves embedded in high-amplitude noise and turbulent flows simulated by a hierarchy of oceanic models. DAHD results reveal striking low-rank behavior and a sparse representation of complex, multiscale, and chaotic flows by a relatively few data-inferred spatial patterns evolving with simple temporal dynamics, as well as well-defined scaling behavior across temporal frequencies, such as exponential-like shape and power law.

Over the past decade, approaches based on Dynamical Mode Decomposition (DMD) and Koopman analysis have quickly gained popularity to analyze datasets in the engineering fluids community (Schmid, 2010; Budišić et al., 2012; and Williams et al., 2015). DMD computes from the data of the eigenvalues and eigenvectors of a linear model that approximates the underlying nonlinear dynamics. Tu et al. (2014) have established a close connection of DMD with independently developed in climate science principal orthogonal patterns (Penland, 1989) and linear inverse modeling (Penland, 1996), which have been since generalized to include memory and nonlinear effects (Kondrashov et al., 2015; Mukhin et al., 2015; and Gavrilov et al., 2016), as well as state-dependent noise (Kravtsov et al., 2016; Martinez-Villalobos et al., 2018).

Spectral proper orthogonal decomposition (SPOD) (Towne et al., 2018) and spectral empirical orthogonal function analysis (SEOF) (Schmidt et al., 2019), are based on an eigenvalue decomposition of the estimated cross-spectral density matrix, yielding at each frequency a set of time-harmonic and orthogonal modes; these SPOD/SEOF modes can be interpreted, in turn, as ensemble DMD modes. Data-adaptive harmonic decomposition (DAHD) (Chekroun and Kondrashov, 2017; Kondrashov et al., 2018a) is based on spectral analysis of an integral shift operator with a two-point statistic kernel built from time-lagged cross correlations, and it is realized numerically by an eigendecomposition that is briefly reviewed in Sec. II A. Similar to SPOD/SEOF methods, DAHD is yielding orthogonal set of eigenvectors—data-adaptive harmonic modes (DAHMs) oscillating harmonically in time.

Furthermore, DAHD enables effective inverse modeling of the original dataset by a system of frequency-ranked nonlinear stochastic oscillators with memory effects, which has been successfully applied to challenging datasets across sciences, including oceanic turbulence (Kondrashov et al., 2018a), Arctic sea ice (Kondrashov et al., 2018b; 2018c), and space physics (Kondrashov and Chekroun, 2018).

This study introduces important DAHD modifications and new features, namely, an eigendecomposition of the Hermitian cross-spectral density matrix in the frequency-domain (Sec. II B) and the auxiliary energy spectrum (Sec. II C) obtained by using judicious projection of the data onto DAHMs. DAHD connection and difference with SPOD/SEOF is established (Sec. II B).

The Hermitian form of DAHD is applied first for identification of coherent oceanic waves in noisy synthetic data (Sec. III) and then to analyze several high-dimensional datasets of complex turbulent geophysical flows, namely, Regional Oceanic Modeling System (ROMS) simulation of the equatorial region (Sec. IV) and a wind-driven gyres circulation by a eddy-resolving three-layer quasigeostrophic ocean model (Sec. V). Discussion and conclusions follow in Sec. VI.

We consider a multivariate time series X(n)=(X1(n),,Xd(n)) formed with d spatial channels and n=1,,N time points (sampled evenly). Double-sided (unbiased) cross-correlation coefficients ρ(p,q)(m) are estimated for all pairs of channels p and q and the time lag m up to a maximum M1,

ρ(p,q)(m)={1Nmn=1NmXp(n+m)Xq(n),0mM1,ρ(q,p)(m),m<0.
(1)

Here, M is the embedding window and its size should be larger than typical decorrelation times in the data; it is the slowest temporal scale captured by DAHD. To reduce the biases of cross-correlation estimates, common lag windowing can be applied.

Combining all the lagged cross-correlation coefficients between a given pair of channels, p and q lead to a cross-correlation Hankel matrix H(p,q), which is symmetric and obtained by a left shift of the row vector (ρM+1(p,q),,ρ0(p,q),,ρM1(p,q)); therefore, every anti-diagonal consists of the same elements. The block-matrix C of size d(2M1)d(2M1) is made symmetric by construction and obtained by arranging all the Hankel matrices H(p,q) for each pair of channels (p,q),

C(p,q)=H(p,q),pq,C(p,q)=H(q,p),elsewhere.
(2)

An important property of C is that its eigenvalues λ come in pairs of opposite values, while the eigenvectors Wj=(E1j,,Edj) represent a collection of global space-time patterns, namely, data-adaptive harmonic modes (DAHMs) oscillating at a single temporal frequency. In particular, Ekj is M-long time series that corresponds to a harmonic oscillation with a frequency f,

Ekj(s)=Bkjcos(2πfs+θkj),1sM,1kd,
(3)

where the amplitudes Bkj and phases θkj are data-adaptive, while the frequency f is equally spaced in the Nyquist interval [00.5] with M values,

f=(1)M1,=1,,M+12.
(4)

In total, j=1,,d(2M1) spectral DAHD eigenelements (|λj|,Wj) are computed, and the DAHD spectrum is obtained by plotting eigenvalues |λ| according to their frequency f.

Several DAHD studies (Chekroun and Kondrashov, 2017; Kondrashov et al., 2018a) have relied on eigendecomposition of C [Eq. (2)] to perform the analysis on relatively low-dimensional datasets, i.e., d not exceeding 40. High-dimensional datasets can be compressed first by PCA, aiming to retain nearly all of the variance (98%–99%) in leading PCs, which are then analyzed. In practice, DAHD results are invariant to the orthogonal rotation of the data and to Principal Component Analysis (PCA). in particular. Still, even after PCA compression, the matrix C can be very large and its direct eigendecomposition is computationally prohibitive. For such cases, Ryzhov et al. (2019; 2020) have utilized frequency-domain DAHD formulation, described next.

Theorem V.1 of Chekroun and Kondrashov (2017) has established that DAHD eigenvalues are related to singular values of a symmetrized complexd×d cross-spectral matrix S(f) whose elements are given by

Sp,q(f)={ρp,q^(f)ifqp,ρq,p^(f)ifq<p,
(5)

where ρp,q^(f) is the Fourier transform at the frequency f of the cross-correlation sequence ρp,q(m),

ρp,q^(f)=m=M+1M1ρp,q(m)e2πif.
(6)

In particular, for each singular value σk(f) of S(f), there exists, when f0, a pair of negative-positive DAHD eigenvalues (λk+(f),λk(f)) of C such that

λk+(f)=λk(f)=σk(f),1kd;
(7)

i.e., 2d eigenvalues are associated with each Fourier frequency f0, while respective DAHMs are shifted by the quarter of the period; i.e., θk+=θk+π/2. There are only d (not paired) eigenvalues for the frequency f=0.

As a point of departure from the DAHD formulation presented above, it is important to note that the matrix S(f) in Eq. (5) is symmetrized by the intent for practical numerical reasons because ρp,q^(f) and ρq,p^(f) are generally different when pq. However, ρp,q^(f) and ρq,p^(f) are necessarily related to each other due to the reversal property of estimated time-lagged cross-correlation sequences, namely,

ρp,q(m)=ρq,p(m),
(8)

where M+1mM1. This is due to the fact that for ρp,q(m), channel p is leading channel q by m lags, and for ρq,p(m), it is the opposite; see Eq. (1). Equivalently, Eq. (8) can be written also as

ρp,q(2Mm)=ρq,p(m),
(9)

where 1m2M1. It is easy to show that Fourier transforms of ρq,p(m) and ρp,q(m) sequences are then related by the frequency-dependent phase shift,

ρp,q^(f)=eiϕρq,p^¯(f),
(10)

where ρq,p^¯ is a complex conjugate, f=0.5(l1)M1, ϕ=2π(l1)2M1, 1lM. Furthermore, since Eqs. (8)–(10) hold when p=q, we obtain for the diagonal elements of S(f),

Sp,p(f)=eiϕ/2|ρp,p^(f)|.
(11)

We can thus define the cross-spectral density matrix S(f),

S(f)=eiϕ/2S(f),
(12)

where S(f) is non-symmetrized; i.e., Sp,q(f)=ρp,q^(f) for 1p,qd. The matrix S is Hermitian, and S=S=ST¯ by the way of Eqs. (10) and (11), and it is interpreted as a correlogram estimate of a multidimensional Hermitian cross-spectral density E{X(f)^X(f)^}, where an expectation operator E is an ensemble average over different realizations of a multidimensional random process X. The phase-shift factor eiϕ/2 in Eq. (11) accounts for discrete case estimation. For the univariate case in particular (i.e., p=q=1), S(f)=|ρ(f)^| represents a well established correlogram method based on the Wiener–Khinchin theorem, which states that the power spectral density is equal to the Fourier transform of its autocorrelation function (Percival and Walden, 1993). Furthermore, since S(f) is a Hermitian matrix, its eigendecomposition,

S(f)=UΛU,
(13)

yields a set of d real eigenvalues Λ={diag}(λ1,,λd) (their absolute values are used for the analysis and results presented), while associated d eigenvectors U form an orthonormal set at each frequency f in the Nyquist interval [see Eq. (10)].

The Hermitian DAHD formulation [Eqs. (11)–(13)] thus closely resembles a SPOD/SEOF approach (Towne et al., 2018; Schmidt et al., 2019), with the key difference being that the latter considers Welch’s overlapped averaging periodogram estimate of the cross-spectral density matrix, while the former is based on a correlogram estimate [Eq. (11)]. However, the same as in SPOD, Hermitian DAHD eigenvectors Uk are associated with the Fourier transform of a spatial wave, which is harmonic in time at a single frequency f; i.e., the time-domain DAHMs Wk [see Eq. (3)] have the following representation:

Wk^(f)=Uk(f).
(14)

For the purpose of reconstruction in the time-domain and the computing energy spectrum (see Sec. II C), for each λk, pairs (Wk+^(f)=Uk(f),Wk^(f)=eiπ/2Uk(f)) are formed. In turn, Wk+(f) and Wk(f) are obtained by an inverse Fourier transform to yield space-time patterns at a single temporal frequency f and shifted by a quarter of a period (aka sin and cos) and thus provide an orthonormal basis set of modes W in the time-domain across all frequencies; i.e., WTW=I. By analogy with SPOD/SEOF modes, DAHD modes can be interpreted then as optimal response modes of the forced linear system if the forcing is white in space and time (Schmidt et al., 2019; Towne et al., 2018).

The frequency-domain approach is fully parallelizable and thus is very computationally efficient since the eigendecomposition of each frequency can be performed in parallel.

The maximum resolution in the frequency [see Eq. (4)] is achieved when M=(N+1)/2, and in this case, DAHMs attain the size Nxd; i.e., they have exactly the same dimensions as the original data X in time and space. Since the set of DAHMs is an orthonormal basis, i.e., WTW=I, it can be used to perform the following expansion of X:

X=j=1dNμjWj.
(15)

Here, μj’s are scalar weights—Data-Adaptive Harmonic coefficients (DAHCs) that are uniquely defined by projecting X onto DAHMs, namely,

μj=WjTX=n=1Nk=1dEkj(n)Xk(n).
(16)

The original data can then be partially (or fully) reconstructed by selecting a subset (whole set) of μ’s and associated DAHMs in Eq. (15). The physical meaning of DAHCs readily follows; namely, when squared, they establish relative contributions of respective DAHMs to L2 energy of X, i.e., with the above definitions,

n=1Nk=1dXk2(n)=j=1Ndμj2.
(17)

Because of pairing in DAHD eigenvectors and their association with a single temporal frequency f [see Sec. II B and Eq. (14)], it is useful to isolate in the sum on the right hand side of Eq. (17) the energy contribution for any given pair,

αk2(f)=(μk+(f))2+(μk(f))2,1kd.
(18)

The DAHD energy spectrum is then obtained by plotting a collection of α2(f) as a function of their frequency, similar to the plot of λ(f)’s.

We consider here a synthetic example of several propagating 2D waves,

un(x,y,t)=Ancos(knx/Lx+lny/Ly+ωnt),
(19)

where An are the weights, while frequency ωn and wave numbers (kn,ln) obey the Rossby dispersion relation,

ωn=βknkn2+ln2+R2.
(20)

The full dataset v(x,y,t) consists of a coherent component s(x,y,t)=n=1Kun(x,y,t) embedded in a temporally correlated and a spatially uncorrelated noise r(x,y,t) represented locally by independent AR(1) processes and a randomly chosen AR(1) coefficient,

v(x,y,t)=(1ν)1/2s(x,y,t)+ν1/2r(x,y,t).
(21)

For this example, we choose Lx=Ly=64 and uniform spacing in x and y with Nx=Ny=64 of grid points; R=5, β=2107, and K=5 waves with the spatial wave numbers k1=10, k2=20, k3=5, k4=4, k5=2, l1=10, l2=20, l3=5, l4=3, l5=1; and the wave weights A1=1, A2=0.5, A3=0.75, A4=0.6, and A5=0.9. With an appropriately chosen time step, the wave periodicities |2π/ωn| in sampling units are T163, T2126, T332, T420, and T58, and we generate the dataset with N=999 points in time. By using the coefficient ν=0.9 in Eq. (21), we consider a case of a low signal-to-noise ratio; i.e., the wavy signal accounts only for 14% of the variance in the full data v(x,y,t); see the comparison of “clean” 2D wave snapshots with their noisy sum (Fig. 1, Multimedia view) and the time series at the (x,y) point (Fig. 2).

FIG. 1.

Temporal snapshots of 2D reference waves and their sum contaminated by noise. Multimedia view: https://doi.org/10.1063/5.0012077.1

FIG. 1.

Temporal snapshots of 2D reference waves and their sum contaminated by noise. Multimedia view: https://doi.org/10.1063/5.0012077.1

Close modal
FIG. 2.

Time series of the wavy signal and full data with imposed red noise at the selected (x,y) point.

FIG. 2.

Time series of the wavy signal and full data with imposed red noise at the selected (x,y) point.

Close modal

Then, we apply the frequency-domain DAHD Hermitian algorithm (Sec. II B) to the noisy data v(x,y,t) with the goal to diagnose and reconstruct reference waves un(x,y,t) in space and time. We use the embedding window M=(N+1)/2=500 in sampling units to obtain a maximum spectral resolution.

Figure 3 shows DAHD eigenvalues [λ(f), see Eq. (13)] and the energy spectrum α2(f) [see Eq. (18)], which is evenly spaced with M=500 bins in the Nyquist range [0 0.5]. Both λ(f) and α2(f) spectra are characterized by sharp peaks at the frequencies of reference waves, located high above noise background and associated with a pair of DAHMs [W’s, Eq. (3)]. Furthermore, the energy spectra α2(f) has an additional advantage by providing cleaner identification of coherent waves from noise, i.e., without multiple background lines present in λ(f) spectra.

FIG. 3.

DAHD spectra of eigenvalues λ (left panel) and energy α2 (right panel) for the 2D synthetic dataset: red circles—pairs of modes associated with the largest values at reference wave frequencies.

FIG. 3.

DAHD spectra of eigenvalues λ (left panel) and energy α2 (right panel) for the 2D synthetic dataset: red circles—pairs of modes associated with the largest values at reference wave frequencies.

Close modal

Figure 4 (Multimedia view) shows fairly accurate reconstruction [see Eq. (15)] of reference wave patterns by using a pair of DAHMs [W’s, Eq. (3)] and DAHCs [μ’s, Eq. (16)] associated with the spectral peaks in Fig. 3.

FIG. 4.

Temporal snapshot for DAHD reconstruction of waves by using DAHMs associated with spectral peaks in Fig. 3 and a comparison with respective reference patterns (see Fig. 1, Multimedia view). Multimedia view: https://doi.org/10.1063/5.0012077.2

FIG. 4.

Temporal snapshot for DAHD reconstruction of waves by using DAHMs associated with spectral peaks in Fig. 3 and a comparison with respective reference patterns (see Fig. 1, Multimedia view). Multimedia view: https://doi.org/10.1063/5.0012077.2

Close modal

The dataset consists of N=1803 hourly snapshots of the surface vorticity field in the square domain (with a 1002×1002 spatial resolution) of the equatorial region simulated by the Regional Oceanic Modeling System (Srinivasan et al., 2017)—a primitive equation ocean model. To enable efficient DAHD analysis for such a high-dimensional dataset, it was first compressed by PCA. The leading d=1000 empirical orthogonal function (EOF) modes and principal components (PCs) (from a total of 1803 modes) capture 99% of variance. Next, we apply the Hermitian DAHD algorithm (Sec. II B) to the dataset composed by the time series of d=1000 PCs and with a window M=(N+1)/2=902 (in sampling units) to obtain maximum spectral resolution.

Figure 5 shows DAHD eigenvalues [λ(f), see Eq. (13)] and the energy spectrum α2(f) [see Eq. (18)] that is evenly spaced with M=902 bins in the Nyquist range [00.5]h1. At each frequency, 100 largest values of α2 are shown from the maximum possible d=1000. The α2(f) spectrum is characterized by a thin spectral line formed by a pair of DAHMs capturing the largest energy per frequency and located high above diffuse background. In comparison, the eigenspectrum λ(f) is more compact and with a smaller gap from the background, similar to the results of Sec. III.

FIG. 5.

DAHD spectra of eigenvalues λ(f) [panels (a) and (b)] and energy α2(f) [panels (c) and (d)] of the surface vorticity field in ROMS simulation: blue—modes with the largest λ at each frequency; green—power law fit f0.4; magenta—exponential law fit e2τf with the characteristic time scale τ10 h; black dots—modes at selected frequencies shown in Fig. 6.

FIG. 5.

DAHD spectra of eigenvalues λ(f) [panels (a) and (b)] and energy α2(f) [panels (c) and (d)] of the surface vorticity field in ROMS simulation: blue—modes with the largest λ at each frequency; green—power law fit f0.4; magenta—exponential law fit e2τf with the characteristic time scale τ10 h; black dots—modes at selected frequencies shown in Fig. 6.

Close modal

Comparison of the top spectral α2 line in log–log and log-linear plots in Fig. 5 reveals power-like dependence at very low frequencies f<0.02h1 (green line) and exponential decay for f[0.020.3]h1 (magenta line). Several spectral structures with a much smaller magnitude and varying scaling laws are also revealed in the diffuse background for f[0.150.35]h1.

It should be noted that power spectra with exponential law have not been widely reported or discussed in the context of turbulent or chaotic data. Maggs and Morales (2012) have observed exponential spectra in UCLA plasma physics experiments and have interpreted it as nonlinear signatures of chaotic dynamics in Lorentzian pulses, given by g(t)=τ2(tt0)2+τ2 with power spectra e2ffs, where fs=1τ is the scaling frequency. By fitting values in the top spectral line, we obtain τ10 h for exponential spectra and f0.4 for the power-law spectra at very low frequencies.

Figure 6 shows the magnitude |B| of the DAHM after transformation into a physical space from the EOF space and associated with the largest value of α2(f) at selected frequencies—black circles in Fig. 5. While each of the DAHMs oscillates harmonically in time [see Eq. (3)], it is inherently multiscale with large and small spatial features. Their dynamical interpretation as optimal response modes of the linear system (Sec. II B) hinges on the validity of approximating effective forcing of the vorticity field, by the white noise in space and time, and it is clearly beyond the scope of this study.

FIG. 6.

Magnitude |B| of DAHM [see Eq. (3)] associated with the largest α2(f) and λ(f) at selected frequencies (black dots in Fig. 5).

FIG. 6.

Magnitude |B| of DAHM [see Eq. (3)] associated with the largest α2(f) and λ(f) at selected frequencies (black dots in Fig. 5).

Close modal

Figure 7 (Multimedia view) demonstrates that the instantaneous full vorticity flow field can be reconstructed with high accuracy [see Eq. (15)] by using a small subset of DAHD eigenelements in the top spectral energy line across frequencies, i.e., a pair of DAHMs [Ws, Eq. (3)] and associated DAHCs [μ’s, Eq. (16)], which yield the largest energy contribution (α2(f)) at each frequency f[0.00.35]h1 and in total accounting for 97% of the variance in the vorticity field. The correlation between reconstruction and reference patterns is very high over whole time; see Fig. 8.

FIG. 7.

An instantaneous snapshot of surface vorticity in the original data ROMS (left) vs reconstruction (right) obtained by using DAHCs and DAHMs associated with the top spectral energy line, i.e., the largest α2(f) at a given frequency; see the text for details, units (s1). Multimedia view: https://doi.org/10.1063/5.0012077.3

FIG. 7.

An instantaneous snapshot of surface vorticity in the original data ROMS (left) vs reconstruction (right) obtained by using DAHCs and DAHMs associated with the top spectral energy line, i.e., the largest α2(f) at a given frequency; see the text for details, units (s1). Multimedia view: https://doi.org/10.1063/5.0012077.3

Close modal
FIG. 8.

Correlation of instantaneous spatial patterns of reconstruction and reference ROMS vorticity fields (see Fig. 7, Multimedia view).

FIG. 8.

Correlation of instantaneous spatial patterns of reconstruction and reference ROMS vorticity fields (see Fig. 7, Multimedia view).

Close modal

Such high-degree of information compression in the reconstruction is achieved by 2900=1800 DAHD-inferred spatial patterns (i.e., DAHMs and taking into account the phase-quadrature property of the latter), which evolve with trivial temporal dynamics, namely, harmonic oscillations. In contrast, PCA compression yields 1000 spatial EOFs but with very complex temporal dynamics in associated PCs.

The dataset is obtained by the eddy-resolving oceanic quasigeostrophic (QG) model (Shevchenko and Berloff, 2015; Kondrashov and Berloff, 2015; and Ryzhov et al., 2019; 2020) that produces a double-gyre flow pattern, characterized by a well-developed and turbulent eastward jet extension of the western boundary currents with its adjacent recirculation zones. Dataset consists of N=5999 5-day snapshots of potential vorticity in the upper layer (Fig. 9) solution with a 513×513 spatial resolution. The leading d=2000 EOFs and PCs from PCA capture 99% of the total variance. Next, we apply the Hermitian DAHD algorithm (Sec. II B) to the dataset composed by the time series of d=2000 PCs and with a window M=(N+1)/2=3000 (in sampling units).

FIG. 9.

Snapshot of an upper-layer potential vorticity (PV) field in the three-layer quasigeostrophic oceanic model simulation of wind-driven gyres (Ryzhov et al., 2019). Nondimensional color scale units, the same for Fig. 12 (Multimedia view).

FIG. 9.

Snapshot of an upper-layer potential vorticity (PV) field in the three-layer quasigeostrophic oceanic model simulation of wind-driven gyres (Ryzhov et al., 2019). Nondimensional color scale units, the same for Fig. 12 (Multimedia view).

Close modal

Figure 10 shows DAHD eigenvalues [λ(f), see Eq. (13)] and the energy spectrum α2(f) [see Eq. (18)]. For energy spectra, only the largest 30 α2(f) values are shown from the total number of d=2000 at each frequency. The shape of the spectra is characterized by a narrow spectral line of one α2(f) value per frequency with a wide gap above diffuse background of lower values. Furthermore, the top spectral line exhibits a sharp spectral peak at 20 years, followed by an exponential law e2τf with the characteristic time scale τ25 days, while Kolmogorov’s (53) power law is diagnosed in the high-frequency range (green).

FIG. 10.

DAHD spectra of eigenvalues λ(f) [panels (a) and (b)] and energy α2(f) [panels (c) and (d)] of the PV field (Fig. 9); blue—modes with the largest λ at each frequency; green—the power-law fit f5/3; magenta—the exponential law fit e2tsf with the characteristic time scale τs25 days; black dots—DAHD modes shown in Fig. 11.

FIG. 10.

DAHD spectra of eigenvalues λ(f) [panels (a) and (b)] and energy α2(f) [panels (c) and (d)] of the PV field (Fig. 9); blue—modes with the largest λ at each frequency; green—the power-law fit f5/3; magenta—the exponential law fit e2tsf with the characteristic time scale τs25 days; black dots—DAHD modes shown in Fig. 11.

Close modal

Figure 11 shows the magnitude |B| of the DAHM associated with the largest value of α2(f) at selected frequencies (see black circles in Fig. 10). Notably, the DAHM associated with the decadal spectral peak represents a coherent spatial pattern along the jet, while the spatial patterns in exponential and Kolmogorov ranges are patchy and inherently multiscale.

FIG. 11.

Magnitude |B| of DAHM [see Eq. (3)] associated with the largest α2(f) at selected frequencies f (cycle/year); see black dots in the top spectral line in Fig. 10.

FIG. 11.

Magnitude |B| of DAHM [see Eq. (3)] associated with the largest α2(f) at selected frequencies f (cycle/year); see black dots in the top spectral line in Fig. 10.

Close modal

Figure 12 (Multimedia view) shows that the instantaneous PV anomaly field can be reconstructed very accurately [see Eq. (15)] by using DAHD eigenelements in the top spectral line (i.e., blue dots in Fig. 10), accounting for 96% of variance in the reference field, while the correlation between reconstruction and reference PV patterns is high over the whole time interval; see Fig. 13.

FIG. 12.

An instantaneous snapshot of a PV anomaly in the original QG data (left) vs reconstruction (right) obtained by using DAHD elements in the top spectral line (blue dots in Fig. 10). Multimedia view: https://doi.org/10.1063/5.0012077.4

FIG. 12.

An instantaneous snapshot of a PV anomaly in the original QG data (left) vs reconstruction (right) obtained by using DAHD elements in the top spectral line (blue dots in Fig. 10). Multimedia view: https://doi.org/10.1063/5.0012077.4

Close modal
FIG. 13.

Correlation of instantaneous spatial patterns of reconstruction and reference PV anomaly fields (see Fig. 12, Multimedia view).

FIG. 13.

Correlation of instantaneous spatial patterns of reconstruction and reference PV anomaly fields (see Fig. 12, Multimedia view).

Close modal

The energy spectra and scaling laws of turbulent atmospheric and oceanic flows (Delsole, 2004; Callies and Ferrari, 2013; McWilliams, 2016; and Chapman, 2017) are commonly analyzed in terms of spatial scales (wavenumber), implicitly implying ergodic assumption, i.e., larger/smaller spatial scales corresponding to their slower/faster evolution. In comparison, DAHD explicitly focuses on the temporal scales but without any built-in assumptions on their spatial content and thus can provide new insights, such as an exponential shape of energy spectra not widely reported before. Furthermore, DAHD yields a sparse representation of complex, multiscale, and chaotic dynamics by relatively few data-inferred spatial patterns oscillating harmonically in time. The observed low-rank behavior can be interpreted as dominance of a given physical mechanism of the energy distribution and transfer across temporal frequencies and revealed by scaling laws of DAHD spectra. By utilizing the correlogram estimator of the Hermitian cross-spectral density matrix, DAHD is closely related but distinctly different from SPOD that relies on Welch’s periodogram.

The results of this study pave the way to several promising DAHD applications, such as (i) obtaining a highly accurate spectral fingerprint of scaling laws in geophysical and astrophysical turbulent fluid flows, (ii) global spectral diagnostics for comparison and characterization of model simulations and observations, (iii) parsimonious extraction and reconstruction of mutiscale features of interest enabled by DAHD compression of energetic content, (iv) uncovering invariant laws across different spatial and temporal scales and phenomena, (v) predictive capability enabled by simple temporal dynamics contained in DAHMs, and (vi) identification of coherent structures in noisy data.

The authors would like to thank Professor James McWilliams and Dr. Kaushik Srinivasan for many useful discussions and providing a ROMS dataset. Two anonymous reviewers provided valuable inputs, which helped to improve the manuscript. This research was supported by the National Science Foundation (NSF) grant (No. OCE-1658357) and the NERC grant (No. NE/R011567/1). P.B. also gratefully acknowledges funding by the NERC (Grant No. NE/T002220/1) and Leverhulme (Grant No. RPG-2019-024). We would also like to acknowledge high-performance computing support from Cheyenne (https://doi.org/10.5065/D6RX99HX) provided by NCAR’s Computational and Information Systems Laboratory, sponsored by NSF. DAHD Toolbox is available at http://research.atmos.ucla.edu/tcd/dkondras/Software.html.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1
Budišić
,
M.
,
Mohr
,
R.
, and
Mezić
,
I.
, “
Applied Koopmanism
,”
Chaos
22
,
047510
(
2012
).
2
Callies
,
J.
and
Ferrari
,
R.
, “
Interpreting energy and tracer spectra of upper-ocean turbulence in the submesoscale range (1–200 km)
,”
J. Phys. Oceanogr.
43
,
2456
2474
(
2013
).
3
Chapman
,
C. C.
, “
New perspectives on frontal variability in the Southern Ocean
,”
J. Phys. Oceanogr.
47
,
1151
1168
(
2017
).
4
Chekroun
,
M. D.
and
Kondrashov
,
D.
, “
Data-adaptive harmonic spectra and multilayer Stuart-Landau models
,”
Chaos
27
,
093110
(
2017
).
5
Delsole
,
T.
, “
Stochastic models of quasigeostrophic turbulence
,”
Surv. Geophys.
25
,
107
149
(
2004
).
6
Gavrilov
,
A.
,
Mukhin
,
D.
,
Loskutov
,
E.
,
Volodin
,
E.
,
Feigin
,
A.
, and
Kurths
,
J.
, “
Method for reconstructing nonlinear modes with adaptive structure from multidimensional data
,”
Chaos
26
,
123101
(
2016
).
7
Kondrashov
,
D.
and
Berloff
,
P.
, “
Stochastic modeling of decadal variability in ocean gyres
,”
Geophys. Res. Lett.
42
,
1543
1553
, https://doi.org/10.1002/ 2014GL062871 (
2015
).
8
Kondrashov
,
D.
and
Chekroun
,
M. D.
, “
Data-adaptive harmonic analysis and modeling of solar wind-magnetosphere coupling
,”
J. Atmos. Solar-Terrestrial Phys.
177
,
179
189
(
2018
).
9
Kondrashov
,
D.
,
Chekroun
,
M. D.
, and
Berloff
,
P.
, “
Multiscale Stuart-Landau emulators: Application to wind-driven ocean gyres
,”
Fluids
3
,
21
(
2018a
).
10
Kondrashov
,
D.
,
Chekroun
,
M. D.
, and
Ghil
,
M.
, “
Data-driven non-Markovian closure models
,”
Physica D
297
,
33
55
(
2015
).
11
Kondrashov
,
D.
,
Chekroun
,
M. D.
, and
Ghil
,
M.
, “
Data-adaptive harmonic decomposition and prediction of Arctic sea ice extent
,”
Dyn. Stat. Clim. Syst.
3
,
dzy001
(
2018b
).
12
Kondrashov
,
D.
,
Chekroun
,
M. D.
,
Yuan
,
X.
, and
Ghil
,
M.
, “Data-adaptive harmonic decomposition and stochastic modeling of Arctic sea ice,” in Advances in Nonlinear Geosciences, edited by A. Tsonis (Springer, 2018c).
13
Kravtsov
,
S.
,
Tilinina
,
N.
,
Zyulyaeva
,
Y.
, and
Gulev
,
S. K.
, “
Empirical modeling and stochastic simulation of sea level pressure variability
,”
J. Appl. Meteorol. Climatol.
55
,
1197
1219
(
2016
).
14
Maggs
,
J. E.
and
Morales
,
G. J.
, “
Exponential power spectra, deterministic chaos and lorentzian pulses in plasma edge dynamics
,”
Plasma Phys. Controlled Fusion
54
,
124041
(
2012
).
15
Martinez-Villalobos
,
C.
,
Vimont
,
D. J.
,
Penland
,
C.
,
Newman
,
M.
, and
Neelin
,
J. D.
, “
Calculating state-dependent noise in a linear inverse model framework
,”
J. Atmos. Sci.
75
,
479
496
(
2018
).
16
McWilliams
,
J. C.
, “
Submesoscale currents in the ocean
,”
Proc. R. Soc. A: Math. Phys. Eng. Sci.
472
,
20160117
(
2016
).
17
Mukhin
,
D.
,
Gavrilov
,
A.
,
Feigin
,
A.
,
Loskutov
,
E.
, and
Kurths
,
J.
, “
Principal nonlinear dynamical modes of climate variability
,”
Sci. Rep.
5
,
15510
(
2015
).
18
Penland
,
C.
, “
Random forcing and forecasting using principal oscillation pattern analysis
,”
Mon. Weather Rev.
117
,
2165
2185
(
1989
).
19
Penland
,
C.
, “
A stochastic model of IndoPacific sea surface temperature anomalies
,”
Physica D
98
,
534
558
(
1996
).
20
Percival
,
D. B.
and
Walden
,
A. T.
,
Spectral Analysis for Physical Applications
(
Cambridge University Press
,
1993
).
21
Ryzhov
,
E.
,
Kondrashov
,
D.
,
Agarwal
,
N.
, and
Berloff
,
P.
, “
On data-driven augmentation of low-resolution ocean model dynamics
,”
Ocean Modell.
142
,
101464
(
2019
).
22
Ryzhov
,
E.
,
Kondrashov
,
D.
,
Agarwal
,
N.
, and
Berloff
,
P.
, “On data-driven induction of the low-frequency variability in a coarse-resolution ocean model,” Ocean Modell. (submitted) (2020).
23
Schmid
,
P. J.
, “
Dynamic mode decomposition of numerical and experimental data
,”
J. Fluid Mech.
656
,
5
28
(
2010
).
24
Schmidt
,
O. T.
,
Mengaldo
,
G.
,
Balsamo
,
G.
, and
Wedi
,
N. P.
, “
Spectral empirical orthogonal function analysis of weather and climate data
,”
Mon. Weather Rev.
147
,
2979
2995
(
2019
).
25
Shevchenko
,
P.
and
Berloff
,
P.
, “
Multi-layer quasi-geostrophic ocean dynamics in eddy-resolving regimes
,”
Ocean Modell.
84
,
1
14
(
2015
).
26
Srinivasan
,
K.
,
McWilliams
,
J. C.
,
Renault
,
L.
,
Hristova
,
H. G.
,
Molemaker
,
J.
, and
Kessler
,
W. S.
, “
Topographic and mixed layer submesoscale currents in the near-surface southwestern tropical pacific
,”
J. Phys. Oceanogr.
47
,
1221
1242
(
2017
).
27
Towne
,
A.
,
Schmidt
,
O. T.
, and
Colonius
,
T.
, “
Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis
,”
J. Fluid Mech.
847
,
821
867
(
2018
).
28
Tu
,
J. H.
,
Rowley
,
C. W.
,
Luchtenburg
,
D. M.
,
Brunton
,
S. L.
, and
Kutz
,
J. N.
, “
On dynamic mode decomposition: Theory and applications
,”
J. Comput. Dyn.
1
,
391
421
(
2014
).
29
Williams
,
M. O.
,
Kevrekidis
,
I. G.
, and
Rowley
,
C. W.
, “
A data-driven approximation of the Koopman operator: Extending dynamic mode decomposition
,”
J. Nonlin. Sci.
25
,
1307
1346
(
2015
).