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.

## I. INTRODUCTION

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.

## II. DATA-ADAPTIVE HARMONIC DECOMPOSITION

### A. Time-domain formulation

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

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 $(\rho \u2212M+1(p,q),\u2026,\rho 0(p,q),\u2026,\rho M\u22121(p,q))$; therefore, every anti-diagonal consists of the same elements. The block-matrix $C$ of size $d(2M\u22121)\u22c5d(2M\u22121)$ is made *symmetric by construction* and obtained by arranging all the Hankel matrices $H(p,q)$ for each pair of channels $(p,q)$,

An important property of $C$ is that its eigenvalues $\lambda $ come in pairs of opposite values, while the eigenvectors $Wj=(E1j,\u2026,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\u2032$-long time series that corresponds to a harmonic oscillation with a frequency $f$,

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

In total, $j=1,\u2026,d(2M\u22121)$ spectral DAHD eigenelements ($|\lambda j|,Wj$) are computed, and the DAHD spectrum is obtained by plotting eigenvalues $|\lambda |$ 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.

### B. Frequency-domain formulation

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

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

In particular, for each singular value $\sigma k(f)$ of $S(f)$, there exists, when $f\u22600$, a pair of negative-positive DAHD eigenvalues $(\lambda k+(f),\lambda k\u2212(f))$ of $C$ such that

i.e., $2d$ eigenvalues are associated with each Fourier frequency $f\u22600$, while respective DAHMs are shifted by the quarter of the period; i.e., $\theta k+=\theta k\u2212+\pi /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 $\rho p,q^(f)$ and $\rho q,p^(f)$ are generally different when $p\u2260q$. However, $\rho p,q^(f)$ and $\rho q,p^(f)$ are necessarily related to each other due to the reversal property of estimated time-lagged cross-correlation sequences, namely,

where $\u2212M+1\u2264m\u2264M\u22121$. This is due to the fact that for $\rho p,q(m)$, channel $p$ is leading channel $q$ by $m$ lags, and for $\rho q,p(m)$, it is the opposite; see Eq. (1). Equivalently, Eq. (8) can be written also as

where $1\u2264m\u22642M\u22121$. It is easy to show that Fourier transforms of $\rho q,p(m)$ and $\rho p,q(m)$ sequences are then related by the frequency-dependent phase shift,

where $\rho q,p^\xaf$ is a complex conjugate, $f=0.5(l\u22121)M\u22121$, $\varphi =2\pi (l\u22121)2M\u22121$, $1\u2264l\u2264M$. Furthermore, since Eqs. (8)–(10) hold when $p=q$, we obtain for the diagonal elements of $S(f)$,

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

where $S(f)$ is *non-symmetrized*; i.e., $Sp,q(f)=\rho p,q^(f)$ for $1\u2264p,q\u2264d$. The matrix $S$ is Hermitian, and $S=S\u2217=ST\xaf$ 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)^\u2217}$, where an expectation operator $E$ is an ensemble average over different realizations of a multidimensional random process $X$. The phase-shift factor $e\u2212i\varphi /2$ in Eq. (11) accounts for discrete case estimation. For the univariate case in particular (i.e., $p=q=1$), $S(f)=|\rho (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,

yields a set of $d$ real eigenvalues $\Lambda ={diag}(\lambda 1,\u2026,\lambda 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:

For the purpose of reconstruction in the time-domain and the computing energy spectrum (see Sec. II C), for each $\lambda k$, pairs ($Wk+^(f)=Uk(f),Wk\u2212^(f)=ei\pi /2Uk(f)$) are formed. In turn, $Wk+(f)$ and $Wk\u2212(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.

### C. Energy spectrum and space-time reconstruction

The maximum resolution in the frequency [see Eq. (4)] is achieved when $M=(N+1)/2$, and in this case, DAHMs attain the size $N$x$d$; 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$:

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

The original data can then be partially (or fully) reconstructed by selecting a subset (whole set) of $\mu $’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,

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,

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

## III. IDENTIFICATION OF ROSSBY WAVES

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

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

The full dataset $v(x,y,t)$ consists of a coherent component $s(x,y,t)=\u2211n=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,

For this example, we choose $Lx=Ly=64$ and uniform spacing in x and y with $Nx=Ny=64$ of grid points; $R=5$, $\beta =\u22122\u22c510\u22127$, and $K=5$ waves with the spatial wave numbers $k1=\u221210$, $k2=20$, $k3=\u22125$, $k4=4$, $k5=\u22122$, $l1=10$, $l2=20$, $l3=\u22125$, $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\pi /\omega n|$ in sampling units are $T1\u224863$, $T2\u2248126$, $T3\u224832$, $T4\u224820$, and $T5\u22488$, and we generate the dataset with $N=999$ points in time. By using the coefficient $\nu =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).

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 [$\lambda (f)$, see Eq. (13)] and the energy spectrum $\alpha 2(f)$ [see Eq. (18)], which is evenly spaced with $M=500$ bins in the Nyquist range [0 0.5]. Both $\lambda (f)$ and $\alpha 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 $\alpha 2(f)$ has an additional advantage by providing cleaner identification of coherent waves from noise, i.e., without multiple background lines present in $\lambda (f)$ spectra.

## IV. REGIONAL OCEANIC MODELING SYSTEM

The dataset consists of $N=1803$ hourly snapshots of the surface vorticity field in the square domain (with a $1002\xd71002$ 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 $\u224899%$ 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 [$\lambda (f)$, see Eq. (13)] and the energy spectrum $\alpha 2(f)$ [see Eq. (18)] that is evenly spaced with $M=902$ bins in the Nyquist range $[00.5]h\u22121$. At each frequency, $100$ largest values of $\alpha 2$ are shown from the maximum possible $d=1000$. The $\alpha 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 $\lambda (f)$ is more compact and with a smaller gap from the background, similar to the results of Sec. III.

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

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)=\tau 2(t\u2212t0)2+\tau 2$ with power spectra $\u2248e\u22122ffs$, where $fs=1\tau $ is the scaling frequency. By fitting values in the top spectral line, we obtain $\tau \u224810$ h for exponential spectra and $\u2248f\u22120.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 $\alpha 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.

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 [$W$s, Eq. (3)] and associated DAHCs [$\mu $’s, Eq. (16)], which yield the largest energy contribution ($\alpha 2(f)$) at each frequency $f\u2208[0.00.35]h\u22121$ and in total accounting for $\u224897%$ of the variance in the vorticity field. The correlation between reconstruction and reference patterns is very high over whole time; see Fig. 8.

Such high-degree of information compression in the reconstruction is achieved by $\u22482\u22c5900=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.

## V. EDDY-RESOLVING OCEANIC MODEL OF WIND-DRIVEN GYRES

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\xd7513$ spatial resolution. The leading $d=2000$ EOFs and PCs from PCA capture $\u224899%$ 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).

Figure 10 shows DAHD eigenvalues [$\lambda (f)$, see Eq. (13)] and the energy spectrum $\alpha 2(f)$ [see Eq. (18)]. For energy spectra, only the largest 30 $\alpha 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 $\alpha 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 $\u224820$ years, followed by an exponential law $e\u22122\tau f$ with the characteristic time scale $\tau \u224825$ days, while Kolmogorov’s $(\u221253)$ power law is diagnosed in the high-frequency range (green).

Figure 11 shows the magnitude $|B|$ of the DAHM associated with the largest value of $\alpha 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.

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 $\u224896%$ 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.

## VI. DISCUSSION AND CONCLUSIONS

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.

## ACKNOWLEDGMENTS

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.

## DATA AVAILABILITY

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

## REFERENCES

*Advances in Nonlinear Geosciences*, edited by A. Tsonis (Springer, 2018c).