2.1. Linearised Navier–Stokes equations

The intermittent flow that we use as our example in this work is the incompressible 2-D Kolmogorov flow. The two dimensions are denoted by $x$ and $y$, respectively. The domain is doubly periodic, with size $L_x$ in the $x$-direction and $L_y$ in the $y$-direction. The flow is forced in the $x$-direction by a sinusoidal body forcing of the form $\zeta \sin (2{\rm \pi} n\grave {y}/L_y)$, where $\zeta$ is the amplitude of the forcing per unit mass of fluid, $n$ is the wavenumber of the forcing and $\grave {y}$ indicates the dimensional $y$-coordinate. The mean velocity is $U(y)\hat {\boldsymbol {i}}$, where $(\hat {\boldsymbol {i}},\hat {\boldsymbol {j}})$ are unit vectors in the $(x,y)$ directions. The mean is here defined over the $x$-direction and time. Additionally, it is symmetrised in the $y$-direction (why such a symmetrisation is required is discussed in § 2.3). The fluctuations of velocities around this mean are denoted by $u$ and $v$ in the $x$ and $y$ directions, respectively. Pressure is denoted by $p$ and time by $t$. The length scale $L_y/2{\rm \pi}$ and time scale $\sqrt {L_y/2{\rm \pi} \zeta }$ are used to non-dimensionalise the system. The non-dimensional number that characterises this system is the Reynolds number, defined as $Re:=(\sqrt {\zeta }/\nu )(L_y/2{\rm \pi} )^{(3/2)}$, where $\nu$ is the kinematic viscosity. Here, we consider a flow with $n=4$, $Re=40$ and $L_x=L_y=2{\rm \pi}$, which is a regime characterised by intermittent burst events as will be discussed in § 2.2.

The non-dimensional equations linearised around the mean state $(U(y),0)$ are

(2.1a,b)\begin{equation} \frac{\partial \boldsymbol{u}}{\partial t} + U\frac{\partial \boldsymbol{u}}{\partial x} + v \frac{\partial U}{\partial y} \hat{\boldsymbol{i}} + \boldsymbol{\nabla} p - \frac{1}{Re} \nabla^2 \boldsymbol{u} = \underbrace{\overline{\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u}} - \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u}}_{\boldsymbol{f}},\quad \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u} = 0, \end{equation}

where $\boldsymbol {u}=(u,v)$. The nonlinear terms of the equation $\overline {\boldsymbol {u}\boldsymbol {\cdot } \boldsymbol {\nabla } \boldsymbol {u}} - \boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {\nabla } \boldsymbol {u}$ are hereafter represented by a forcing term $\boldsymbol {f}=(f_x,f_y)$, where $f_x$ and $f_y$ represent the forcing to the $x$ and $y$ momentum equations, respectively.

Using the incompressibility condition and taking the curl of (2.1a), we obtain the vorticity ($\omega$) form of the equation as

(2.2)\begin{equation} \frac{\partial \omega}{\partial t} + \left[U\frac{\partial }{\partial x} -\frac{1}{Re} \varDelta - U'' \varDelta^{{-}1} \frac{\partial }{\partial x} \right] \omega= F, \end{equation}

where $\omega := \partial v/\partial x - \partial u/\partial y$ and $F:= \partial f_y/\partial x - \partial f_x/\partial y$. Let us consider a Fourier transform in the homogeneous $x$-direction. The one-dimensional (1-D) discrete Fourier transform of $\omega$ gives

(2.3)\begin{equation} \omega(x,y,t) = \sum_{k_x={-}(N_x/2)}^{(N_x/2)} \tilde{\omega}(y,t;k_x)\exp({{\rm i}k_x x}), \end{equation}

where $\tilde {\cdot }$ represents the 1-D Fourier transform, $k_x$ is the streamwise wavenumber non-dimensionalised by $L_y/2{\rm \pi}$ and $N_x$ is the number of grid points used to discretise the $x$-direction. Similarly, if $\tilde {F}$ is the 1-D Fourier transform of $F$, we obtain the following equation:

(2.4)\begin{equation} \dot{\tilde{\omega}}(y,t;k_x) + \boldsymbol{A}\tilde{\omega}(y,t;k_x) = \tilde{F}(y,t;k_x),\end{equation}

where $(\dot {\mbox {}})$ denotes the derivative in time. The matrix $\boldsymbol {A}$ contains the finite-dimensional discrete approximations of the linearised momentum equation from (2.2) in terms of the 1-D Fourier transforms with $\partial /\partial x$ replaced by $ik_x$.

2.2. Direct numerical simulation data

To obtain the DNS data for the 2-D Kolmogorov flow, we used the pseudo-spectral code as used in Chandler & Kerswell (Reference Chandler and Kerswell2013). We consider the flow at a Reynolds number $Re=40$ and forcing frequency $n=4$. The DNS was run on a grid with $N_x=256$ points in the $x$-direction and $N_y=257$ points in the $y$-direction. Figure 1(a) shows the time series of total dissipation $D(t):=\langle \omega, \omega \rangle _{x,y}$ normalised by the laminar dissipation rate $D_{lam}:=Re/(2n^2)$. Here, $\langle a,b \rangle _{x,y}$ denotes the inner product $\int _{0}^{L_x} \int _{0}^{L_y} b^*(x,y,t) a(x,y,t)\,{{\rm d}\kern0.7pt x}\,{{\rm d}y}$. The time series represented by the grey line is obtained from the full dataset. On the other hand, the time series in black is obtained by retaining only the smallest streamwise wavenumbers with $|k_x|\leq 3$. The black line follows the grey line reasonably well, thereby suggesting that the $|k_x|\leq 3$ modes are responsible for the dominant dynamics in the flow. Example snapshots of the flow at the five times denoted in figure 1(a) are shown in figure 1(c–g), where the red and blue contours represent the positive and negative vorticity fluctuations, respectively. Here again, we observe the presence of large-scale structures in the flow, further confirming that the smallest wavenumbers dominate the flow dynamics. This is consistent with the observations in Page et al. (Reference Page, Brenner and Kerswell2021). For the rest of this manuscript, we therefore project the full DNS data onto the $|k_x|\leq 3$ modes, and use the truncated data for the POD and resolvent-based analyses. This truncation using a Fourier transform is not essential for the work that follows, but just allows us to expedite the computations.

Figure 1.

(a) The time series of $D(t)/D_{lam}$ from the full data is shown in grey for a sample time interval. Additionally, the time series obtained from just the streamwise wavenumbers of $|k_x| \leq 3$ is shown in black. Burst events, defined as times when $D(t)/D_{lam} >0.15$, are shown as the red-shaded regions, with dashed black vertical lines indicating the beginnings of the burst events. (c–g) The vorticity field at the five time instances $t1$–$t5$ indicated in (a) are also shown, where the contours indicate positive (red) and negative (blue) vorticity, respectively. The colour limits are kept the same across the five time snapshots. In (b) the mean profile obtained from the data (black line) is compared with the symmetrised mean profile (red line).

Considering the burst events, if we define them as times when $D(t)/D_{lam} \geq 0.15$ (Page et al. Reference Page, Brenner and Kerswell2021), then they are indicated by the red-shaded regions in figure 1(a) (the vertical black dashed lines therefore indicate the beginnings of these burst events). Throughout the manuscript, we use this criterion of $D(t)/D_{lam} \geq 0.15$ to define burst events, and also $D(t)/D_{lam} \leq 0.1$ to identify regions that are assured to be quiet. Changing these numbers for $D(t)/D_{lam}$ will impact the quantitative values obtained. However, we are here interested only in the relative differences between the prediction methods, and we can expect these trends to remain insensitive to this change. Appendix C confirms that these relative trends are not impacted by this definition of the burst event.

2.3. Symmetrisation of the mean profile

Let us briefly consider the mean velocity profiles $U(y)$ obtained from this flow, where the mean is defined over the $x$-direction and time. The converged mean profile is expected to follow two symmetries: (i) a shift-and-reflect symmetry $\mathcal {S}:U(y) \to - U(y+{\rm \pi} /4)$ and (ii) a rotational symmetry $\mathcal {R}:U(y) \to - U(-y)$. Chandler & Kerswell (Reference Chandler and Kerswell2013) investigated the mean profiles obtained from $10^5$ time units of the flow, and observed that even this was not long enough to obtain a converged mean profile that shows the expected symmetries. The black line (which is asymmetric when investigated closely) in figure 1(b) represents the mean computed from 16 800 snapshots of data. This obtained mean does not follow the symmetries. Therefore, as done in Chandler & Kerswell (Reference Chandler and Kerswell2013), here, we use a symmetrised mean profile $U(y)$ which is extracted from the asymmetric mean obtained from the DNS $U_a(y)$ as

(2.5)\begin{equation} U(y) := \frac{1}{2n} \sum_{m=0}^{2n-1} \mathcal{S}^{{-}m}U^{\dagger}(\mathcal{S}^my),\quad \mbox{where}\ U^{\dagger}:=\frac{1}{2} [U_a(y) + \mathcal{R}^{{-}1}U_a(\mathcal{R}y)]. \end{equation}

Here, $n=4$. The definitions of the symmetry operations follows from before as $\mathcal {R}^{-1}:U(y) \to - U(-y)$, $\mathcal {S}^m:U(y) \to (-1)^m U(y+m{\rm \pi} /4)$ and $\mathcal {S}^{-m}:U(y) \to (-1)^m U(y-m{\rm \pi} /4)$, where $0\leq m \leq (2n-1)$. The red line in figure 1(b) shows $U(y)$, and this symmetrised mean will be used for the rest of this work.

3.1. Decomposing the flow: Fourier bases

In this section, we use the Fourier basis for decomposing the flow into the quiet region and the burst events. Looking ahead, isolating the time-localised burst events using the global Fourier basis may prove to be inefficient. Therefore, the aim of this section is to provide motivation for moving towards using a wavelet basis.

We are interested in analysing a burst event, and require multiple realisations of such events. Consider $N_e$ realisations of the burst event, each containing $N_t$ time snapshots of vorticity $\omega (x,y,t)$ spaced apart by time $dt$. Each of the $N_e$ realisations is therefore $T=N_t dt$ long. Here, to obtain one such realisation of the burst event, the start of a burst event $T_b$ is first identified using the condition $D(t)/D_{lam}\geq 0.15$. The realisation of the burst event is then taken to be $T_b-100\leq t \leq T_b+100$. (The reasons for choosing such time-windows, instead of the more generally used consecutive time blocks, will become apparent when considering the energy of the wavelet coefficients in figure 4 of the next section.) For the DNS dataset considered here, we split the data to obtain $N_e=50$ realisations of the burst event, with each realisation having $N_t=200$ and $dt=1$, and therefore $T=200$. It is ensured that no two realisations of the burst event have more than $50\,\%$ overlap.

Let us now consider the Fourier transform of the vorticity $\omega (x,y,t)$ in time

(3.1)\begin{equation} \omega(x,y,t) = \sum_{\varOmega=0}^{N_t-1} \breve{\omega}(x,y;\varOmega)\exp({{\rm i}(2{\rm \pi}/T)\varOmega t}), \end{equation}

where $(\breve {\cdot })$ here represents the temporal Fourier transform. The corresponding Fourier spectrum is $E_{\omega }(\varOmega ) = [\int _0^{2{\rm \pi} } \int _0^{2{\rm \pi} } \breve {\omega } \breve {\omega }^* \,{{\rm d}\kern0.7pt x}\,{{\rm d}y}]_{N_e}$, with $({\cdot })^*$ representing complex conjugate. Here, $[{\cdot }]_{N_e}$ represents an averaging across the $N_e$ different realisations of the burst event. Figure 2(a) shows the obtained Fourier spectrum. To isolate the quiet region and the burst events, we use distinct sets of frequencies. The blue-shaded frequencies in figure 2(a) are used for the quiet region, and the red-shaded frequencies for the burst event. Let us denote the reconstruction of the data using this truncated Fourier basis as $\omega _{r}(x,y,t)$. Figure 2(b) shows the time series of the reconstructions $\langle \omega _r, \omega _r \rangle _{x,y}/D_{lam}$ of the quiet region in blue and the burst event in red. These reconstructions are compared with the full data $\langle \omega, \omega \rangle _{x,y}/D_{lam}$ in black.

Figure 2.

(a) The ensemble-averaged Fourier spectrum is shown in black, and the frequencies used to reconstruct the quiet region and the burst event, separately, are shaded in blue and red, respectively. (b) The full data (for $k_x={0,1,2,3}$) (black) are compared with the reconstructions of the quiet region (blue) and the burst event (red) using the shaded frequencies in (a). For comparison, the grey dashed-dot line shows the sum of the blue and the red lines.

The number of modes required for isolating the different regions is identified by defining an error $\epsilon (t)$ in the projection of $\omega _{r}(x,y,t)$ onto $\omega (x,y,t)$, where $\epsilon (t) := 1 - \langle \omega _r, \omega \rangle _{x,y}/\langle \omega, \omega \rangle _{x,y}$. Let $[{\cdot }]_{t}$ denote averaging in time. The blue-shaded frequencies in figure 2(a) ensure that $[ \epsilon (t_q)]_t < 0.2$ where $t_q$ is defined as the times $t$ when $D(t)/D_{lam} \leq 0.1$, i.e. the quiet region. Similarly, the red-shaded frequencies in figure 2(a) ensure that $[ \epsilon (t_b)]_t < 0.4$, where $t_b$ is defined as the times $t$ when $D(t)/D_{lam} \geq 0.15$, i.e. the burst event. The numerical values of $[ \epsilon (t_q)]_t < 0.2$ and $[ \epsilon (t_b)]_t < 0.4$ simply ensure that the dominant frequencies (or wavelets, as in the next section) are included, and changing the values does not significantly impact any of the discussions in this work (see Appendix C).

From figure 2(b), we find that the Fourier basis is not able to efficiently isolate the burst event, because the red line does not go to zero in the quiet regions of the flow. Additionally, we require a significant number of frequencies to reconstruct the burst event. This is expected since we are trying to represent a localised event using the Fourier basis that is global. In order to remedy these problems, in the next section we will explore using a localised basis – wavelets – for this purpose.

3.2. Daubechies 1 wavelet basis and wavelet transform

We now employ a wavelet basis to capture the burst events. The discrete Daubechies 1 (DB1) wavelet, which is also referred to as the Haar wavelet, will be used for a majority of this work. The only exception is Appendix B, where for a different choice of wavelets, that of Daubechies 2, we see that the discussions in this work remain similar.

To obtain a wavelet basis $\varTheta _{n_w}(t)$, we first make a choice of the wavelet corresponding to $n_w=1$, i.e. $\varTheta _{1}(t)$, which is the mother wavelet. For the Daubechies 1 wavelet basis used in this work, the mother wavelet is a step function, as shown in the first panel of figure 3(b). The mother wavelet covers the entire time domain, and here we denote this mother wavelet as belonging to level 1 of the wavelet basis. To get level 2 of this wavelet basis, the mother wavelet is compressed by half. Now we need two wavelets to cover the entire domain, and these are shown as $n_w=2$ and $n_w=3$ in figure 3(b). Level 3 consists of the mother wavelet compressed by a factor of four, and then repeated four times to cover the domain. This level is shown as $n_w=4\unicode{x2013}7$ in figure 3(b). The vertical black dashed-dot lines in figure 3(b) demarcate the different levels of the wavelet transform. Throughout this manuscript, we use such vertical black dashed-dot lines to demarcate the levels. (It should be noted that, in practice, there is a level $0$ corresponding to $n_w=0$ that contains the lowest frequencies of the data that are not included in the other levels. In more technical terms, the wavelet transform for $N=200$ includes the ‘detail coefficients’ from $7$ levels of the transform as well as the ‘approximation coefficient’ from level $1$ (see Daubechies Reference Daubechies1988).) We see that, when considering wavelets, there are two factors that are important: (i) the compression of the wavelet, which is related to the concept of frequency in the Fourier domain and (ii) the location of the wavelet in time.

Figure 3.

Schematic showing (a) the Fourier basis and (b) the discrete wavelet basis of Daubechies 1. The vertical dashed-dot lines in (b) demarcate the different levels of the wavelet basis, and will be used again in figure 5(a) when plotting the energy of the wavelet coefficients.

Similar to the Fourier transform in § 3.1, we obtain a wavelet transform of $\omega (x,y,t)$ as

(3.2)\begin{equation} \omega(x,y,t) = \sum_{n_w=0}^{N_w} \hat{\omega}(x,y;n_w)\varTheta_{n_w}(t), \end{equation}

where $(\hat {\cdot })$ represents the transform onto the basis $\varTheta _{n_w}$, which is here the wavelet basis. The energy of the wavelet coefficients can then be obtained as $E_{\omega }(n_\omega ) = [ \int _0^{2{\rm \pi} } \int _0^{2{\rm \pi} } \hat {\omega } \hat {\omega }^* \,{{\rm d}\kern0.7pt x}\,{{\rm d}y} ]_{N_e}$. It should be noted that the $\varTheta _{n_w}(t)$ used in this work is real valued. However, since it is a complete basis, it provides a basis for both real-valued and complex-valued data. This is relevant when we later consider the wavelet transform of the 1-D Fourier-transformed (in the $x$-direction) data.

A wavelet transform requires a choice of boundary conditions and throughout this manuscript, periodic boundary conditions have been used. The WPOD described later gives similar results for other physically relevant boundary conditions, such as symmetric and reflecting boundary conditions.

3.3. Burst centred windowing

Before looking at the energy of the wavelet coefficients, let us briefly discuss the choice of the $N_e$ realisations of the burst event. Generally, to obtain the different ensembles for ensemble averaging, we split the available data into several consecutive data blocks of equal length, and compute the spectrum for each of these ensembles. The spectra are then averaged across the ensembles. The energy of the wavelet coefficients obtained from such consecutive ensembles is shown in figure 4(a). Within each of the levels (indicated by the vertical dashed-dot lines), we see that there are no noticeable trends, and peaks are notably absent. This absence of peaks can be explained by recalling that the location in time is important when considering wavelets. In this case, for each ensemble, the burst events occur at different locations in time and therefore appear at different $n_w$.

Figure 4.

(a) The ensemble-averaged energy of the wavelet coefficients is shown for the case when ensembles are chosen as consecutive time blocks. (b) A different strategy for choosing ensembles is also illustrated, where the window is chosen such that a burst lies at the centre of it.

To remedy this, we here adopt the alternative windowing strategy that was briefly mentioned in § 3.1 and is depicted in figure 4(b). A realisation of the burst event is defined using the start of a burst event $T_b$ as $T_b-100\leq t \leq T_b+100$, where $T_b$ is identified using the condition $D(t)/D_{lam}\geq 0.15$ (note, time is here non-dimensionalised by $\sqrt {L_y/2{\rm \pi} \zeta }$). As in § 3.1, here, we consider $N_e=50$ realisations of the burst event, each containing $N_t=200$ time snapshots of vorticity $\omega (x,y,t)$ spaced apart by time $dt=1$. The spectrum obtained from this new windowing strategy is shown in figure 5(a). From the red-shaded regions of the spectrum in figure 5(a), we observe clear peaks that are present because of the burst events. In the next section, we will more clearly show that these peaks correspond to the burst events.

Figure 5.

(a) The ensemble-averaged energy of the wavelet coefficients is shown in black, and the wavelets used to reconstruct the separate regions in (b) are shaded in different colours: the $8$ wavelets for the quiet region are shaded in blue and the $39$ wavelets for the burst event are in red. The vertical dashed lines indicate the different levels of the wavelet transform (see figure 3). (b) The full data for $k_x={0,1,2,3}$ (black) are compared with the reconstructions (using the shaded frequencies in (a)) of the quiet region (blue) and the burst event (red). For comparison, the grey dashed-dot line shows the sum of the blue and the red lines.

3.4. Decomposing the flow: wavelet bases

The number of levels present in the spectrum is determined by the length of the time window (here $T=200$) and the type of wavelet used (here DB1). In figure 5(a) there are $8$ levels that are demarcated by the vertical dashed-dot lines. We here use distinct sets of wavelets to isolate the quiet region and the burst events. The first $3$ levels will be used to reconstruct the quiet region, and these are shaded in blue. Next, 20 % of the most energetic modes in levels $4$–$8$ will be used to reconstruct the burst event, and these are shaded in red (Appendix C shows the insensitivity of the results presented to this specific choice of $20\,\%$ of the wavelets). As in § 3.1, the number of wavelets used to reconstruct the flow is chosen such that the error $[ \epsilon (t_q)]_t < 0.2$ and $[ \epsilon (t_b)]_t < 0.4$ (see § 3.1 for the definition of the error metric). Let us denote the data reconstructed using a truncated wavelet basis as $\omega _{r}(x,y,t)$. Figure 5(b) shows the time series of the reconstructions $\langle \omega _{r},\omega _{r} \rangle _{x,y}/D_{lam}$ of the quiet region in blue and the burst event in red. The reconstructions are compared with the full data $\langle \omega,\omega \rangle _{x,y}/D_{lam}$ in black. We observe that isolating the quiet region and the burst events is approximately possible using wavelets.

In comparing the wavelet-based reconstruction in figure 5(b) with the Fourier-based reconstruction in figure 2(b), three observations stand out. Firstly, and most importantly, wavelets are able to isolate the burst events better than the Fourier basis. To understand why, let us turn our attention to the relatively small magnitude (relative to the burst events) oscillations in the quiet region of the flow. When considering the Fourier basis, the same Fourier frequencies contribute both to these quiet oscillations and the burst events. However, since these oscillations occur at different locations in time relative to the burst events, the wavelet basis is able to produce reconstructions of the burst event uncontaminated by the quiet oscillations. Secondly, to reconstruct the burst event, fewer wavelets are required in comparison with Fourier frequencies: $39$ wavelets as opposed to $108$ ($54$ positive) Fourier frequencies. Finally, from figure 5(b) we see that wavelets are able to isolate a single burst event, in contrast to the Fourier-based reconstruction in figure 2(b) where both the burst events are captured equally. Therefore, isolating single burst events, when there are multiple similar events present in a time window, is only possible using the wavelets bases.

We can therefore conclude that wavelets are better at isolating burst events. Using wavelets, we are able to obtain a signal of the burst event that is uncontaminated by the oscillations in the quiet region. Additionally, wavelets enable a lower-order representation of the burst event by requiring fewer wavelets (than Fourier frequencies) to reconstruct the burst event. In the next section, we will therefore concentrate on using the wavelet basis to probe the flow patterns active in the quiet region and the burst events.

4.1. A description of WPOD

Many adaptations of POD, tailored for various classes of problems, are found in the literature (for instance, see reviews by Taira et al. (Reference Taira, Brunton, Dawson, Rowley, Colonius, McKeon, Schmidt, Gordeyev, Theofilis and Ukeiley2017) and Rowley & Dawson (Reference Rowley and Dawson2017) and references therein). We are here interested in finding structures that are coherent in space and time, and a recently introduced POD-based technique for this purpose is spectral-POD (SPOD) (Towne, Schmidt & Colonius Reference Towne, Schmidt and Colonius2018). In this method, the flow is first projected onto a Fourier basis. Thereafter, a POD is performed at each Fourier frequency that gives coherent structures at that particular frequency. These structures are coherent in space due to the properties of POD and coherent in time since they are Fourier modes. However, in § 3.4, we saw that a wavelet basis is better at characterising intermittent flows. Therefore, instead of a SPOD, here we use a WPOD. Using this method we obtain structures that are coherent in space, and both coherent and localised in time since they are wavelet coefficients.

Before describing WPOD, to provide context, let us briefly consider the regular POD. In this case, we have velocity fields $\boldsymbol {q}(x,y,t)=(u(x,y,t),v(x,y,t))$ for a range of time. (For POD, we are required to choose a norm that is maximised to obtain the POD modes. Kinetic energy is a physically relevant norm, and choosing velocity as data, instead of vorticity, ensures that the kinetic energy is used as the POD norm.) The POD modes $\boldsymbol {\phi }_i(x,y)$ give an orthogonal basis for $\boldsymbol {q}(x,y,t)$, such that the $1$st POD mode captures the largest variance of the data, the $i$th POD mode captures the $i$th largest variance, and so on. In other words, if $\mathcal {E}^{pod}_i$ is the projection of the data $\boldsymbol {q}(x,y,t)$ onto the $i$th POD mode, $\boldsymbol {\phi }_i(x,y)$

(4.1)\begin{equation} \mathcal{E}^{pod}_i = \frac{[|\langle \boldsymbol{q}(x,y,t),\boldsymbol{\phi}_i(x,y)\rangle_{x,y}|^2]_t}{\langle \boldsymbol{\phi}_i(x,y),\boldsymbol{\phi}_i(x,y)\rangle_{x,y}}, \end{equation}

then $\mathcal {E}^{pod}_i>\mathcal {E}^{pod}_{i+1}$. Here, as before, $[{\cdot }]_t$ denotes averaging in time, $\langle a,b \rangle _{x,y}$ denotes the inner product $\int _{0}^{L_y} \int _{0}^{L_x} b^*(x,y,t) a(x,y,t)\,{{\rm d}\kern0.7pt x}\,{{\rm d}y}$ and $|{\cdot }|$ represents absolute value. The problem of finding POD modes reduces to finding the eigenvalues of the correlation matrix $QQ^*$ where $Q$ is a matrix with columns containing $\boldsymbol {q}(x,y,t)$, i.e. the $1$st column of $Q$ is $\boldsymbol {q}(x,y,t_1)$ and the $i$th column is $\boldsymbol {q}(x,y,t_i)$. The eigenvectors of $QQ^*$ are equivalently the left singular vectors of the matrix $Q$, and therefore the problem can be solved by performing a singular value decomposition (SVD) of the data matrix $Q$.

Rather than finding the coherent structures in the entire flow using POD, we are interested in finding the coherent structures at particular wavelets. Let us therefore consider wavelet-POD. For this, first consider the $N_e$ different realisations of the burst event (see § 3.3 for how each realisation is chosen from the data). From these $N_e$ realisations, consider one realisation $n_e$, where now $1 \leq n_e \leq N_e$. At this $n_e$, we have the data $\boldsymbol {q}(x,y,t;n_e) = (u(x,y,t;n_e), v(x,y,t;n_e))$. A wavelet transform of this $n_e$th realisation of the burst event will give us $\hat {\boldsymbol {q}}(x,y;n_w,n_e)=(\hat {u}(x,y;n_w,n_e),\hat {v}(x,y;n_w,n_e))$, where $\hat {u}$ and $\hat {v}$ represent the wavelet transform of $u$ and $v$. Therefore, for each wavelet $n_w$, and from each realisation of the burst event $n_e$, we have the state vector $\hat {\boldsymbol {q}}(x,y;n_w,n_e)$. If we now concentrate on just one wavelet $n_w$, we have $N_e$ data vectors $\hat {\boldsymbol {q}}(x,y;n_w,1:N_e)$.

The objective now is to find basis vectors $\hat {\boldsymbol {\phi }}_i(x,y;n_w)$ such that

(4.2)\begin{equation} \hat{\boldsymbol{q}}(x,y;n_w,n_e) = \sum_{i=0}^{N_e} a_i(n_w,n_e) \hat{\boldsymbol{\phi}}_i(x,y;n_w). \end{equation}

For each $n_w$, the equivalent of the POD energy in (4.1) now becomes

(4.3)\begin{equation} \mathcal{E}^{wpod}_i(n_w) = \frac{[|\langle \hat{\boldsymbol{q}}(x,y;n_w,n_e), \hat{\boldsymbol{\phi}}_i(x,y;n_w)\rangle_{x,y}|^2]_{n_e}}{\langle \hat{\boldsymbol{\phi}}_i(x,y;n_w), \hat{\boldsymbol{\phi}}_i(x,y;n_w)\rangle_{x,y}}, \end{equation}

where $\mathcal {E}^{wpod}_i(n_w)\geq \mathcal {E}^{wpod}_{i+1}(n_w)$. It should be noted that, unlike in (4.1) where the averaging is over time, in (4.3) the averaging is across the $N_e$ different realisations of the burst event (denoted by $[{\cdot }]_{n_e}$). If we consider the WPOD modes across all $n_w$ and pick the dominant mode, i.e. the mode with maximum energy across $n_w$, then the projection that is maximised is

(4.4)\begin{equation} \max_{n_w}(\mathcal{E}^{wpod}_1(n_w)) = \frac{[|\langle \boldsymbol{q}(x,y,t;n_e), \boldsymbol{\phi}_1(x,y,t)\rangle_{x,y,t}|^2]_{n_e}}{\langle \boldsymbol{\phi}_1(x,y,t), \boldsymbol{\phi}_1(x,y,t)\rangle_{x,y,t}}, \end{equation}

where $\langle a,b \rangle _{x,y,t}$ denotes the inner product $\int _{0}^{T} \int _{0}^{L_y} \int _{0}^{L_x} b^*(x,y,t) a(x,y,t)\,{{\rm d}\kern0.7pt x}\,{{\rm d}y}\,{\rm d} t$. Here, $\boldsymbol {\phi }_i(x,y,t)$ represents the WPOD mode $\hat {\boldsymbol {\phi }}_i(x,y;n_w)$ in time, i.e. $\boldsymbol {\phi }_i(x,y,t)$ is an inverse wavelet transform of $\hat {\boldsymbol {\phi }}_i(x,y;n_w)$. Similar to POD, for WPOD we perform a SVD of the data matrix $\hat {Q}_{n_w}$, where columns of $\hat {Q}_{n_w}$ are the data $\hat {\boldsymbol {q}}(x,y;n_w,1:N_e)$. The $i$th left singular vector is the $i$th WPOD mode, and the square of the $i$th singular value is the corresponding WPOD energy $\mathcal {E}^{wpod}_i(n_w)$ (4.3).

4.2. The WPOD-based analysis of the 2-D Kolmogorov flow

Figure 6(a) shows the WPOD energies $\mathcal {E}^{wpod}_i(n_w)$ (4.3) as a function of wavelet $n_w$. In other words, each vertical line in figure 6(a) corresponds to the WPOD obtained at that particular $n_w$. Let us first focus on the wavelets responsible for the quiet region shaded in blue, i.e. the first three levels of the wavelet transform. Notably, only these lower levels show any appreciable low-rank behaviour, i.e. within these levels the first few WPOD modes capture significantly more energy than the higher modes. Low rankness in the WPOD spectrum suggests that there is an energetically dominant mechanism that is responsible for the quiet region of the flow. To probe this further, let us now look at the WPOD modes at $n_w=2$ in figure 6(b–f). The WPOD modes $1$ and $2$, as well as modes $3$ and $4$, are shifted versions of the same mode, a result of the inherent symmetries in the flow. (While symmetries can be incorporated into the POD modes, the instantaneous fluctuations of the flow do not adhere to these symmetries, and so we present modes without symmetrisation). Together, modes $1$ and $2$ account for 77 % of the energy at this $n_w$. Crucially, these first two modes closely resemble the unstable eigenfunction obtained from the Navier–Stokes equations linearised around the mean flow (see Appendix A). This indicates that this unstable eigenfunction is mainly responsible for the quiet region of this flow. This is consistent with the observations in Farazmand & Sapsis (Reference Farazmand and Sapsis2017), where they employed an alternative strategy to capture this unstable eigenfunction, and used it to predict the burst events. Moving on to modes $3$ and $4$, while these modes are energetically significant, they capture less energy compared with modes $1$ and $2$. Together they capture $19\,\%$ of the energy at this $n_w$. Structurally, modes $3$ and $4$ capture shearing motions.

Figure 6.

(a) The WPOD spectrum is shown as a function of wavelet $n_w$. The energies of the first $10$ WPOD modes are shown. The vertical dashed-dot lines demarcate the levels of the wavelet transform. The shaded regions represent the wavelets that are responsible for the quiet region in blue and the burst event in red. The first $5$ WPOD modes are shown for two different wavelets: (b–f) $n_w=2$ responsible for the quiet region marked by the blue vertical line in (a) and (g–k) $n_w=22$ responsible for the burst event marked by the red vertical line in (a).

Let us now turn our attention to the wavelets responsible for the burst event. These wavelets are shaded in red in figure 6(a). One initial observation is that there is no significant low rankness at a majority of these wavelets. This lack of low rankness generally suggests that there is no one dominant mechanism that is responsible for the burst event. To investigate these modes further, consider the WPOD modes at $n_w=22$ in figure 6(g–k). Notably, the first $4$ WPOD modes in the burst event are structurally similar to the modes in the quiet region. Additionally, we observe that modes $2$–$4$, while being similar to the unstable eigenfunction, are fragmented versions of this eigenfunction. It should be noted that, this fragmentation of the modes is more apparent for certain $n_w$ in the burst region (for instance, it is not as clearly apparent for $n_w=12$ in level 4). Therefore, for a more conclusive illustration of these fragmented modes, we will consider composite-WPOD modes in the next section.

From these observations, we can hypothesise that the shearing motions are responsible for disrupting the flow due to the unstable eigenfunction. Intermittently, this shearing disrupts the flow enough to cause a burst event. While a WPOD analysis can suggest that such a mechanism causes the burst events, to conclusively show this, we need to do a dynamic mode decomposition (Schmid Reference Schmid2022) analysis that identifies the relevant instabilities of the flow generated by the unstable eigenfunction. However, this falls beyond the scope of the current manuscript.

4.3. Coherent structures in quiet and burst regions – a composite-WPOD analysis

In this section, we introduce composite-WPOD. Instead of identifying structures at specific wavelets $n_w$ as done for WPOD in the previous section, in composite-WPOD we identify the coherent structures for a set of wavelets $\{n_{w1},n_{w2},\ldots \}$. This is important since, from § 3.4 we know that, rather than individual wavelets, sets of wavelets contribute to the quiet region and the burst events. For example, to find modes for the quiet region, we need modes that represent the set of wavelets $\boldsymbol {n}^q=\{n^q_{w1},n^q_{w2},\ldots n^q_{wN_q}\}$ that are responsible for the quiet region (represented by the blue-shaded values of $n_w$ in figure 6a). We therefore use a composite version of the WPOD used in § 4.1.

In § 4.1, to obtain the WPOD modes, we computed the SVD of a data matrix $\hat {Q}_{n_w}$. The columns of $\hat {Q}_{n_w}$ were taken to be the data $\hat {\boldsymbol {q}}(x,y;n_w,1:N_e)$, such that the $1$st column corresponds to $\hat {\boldsymbol {q}}(x,y;n_w,1)$, the $i^{th}$ column to $\hat {\boldsymbol {q}}(x,y;n_w,i)$, and so on. Here, we instead consider the data matrix $\hat {P}$ that contains coefficients at all the $n_w$ that is contained in $\boldsymbol {n}^q$ such that

(4.5)\begin{align} \hat{P} = [\hat{\boldsymbol{q}}(x,y;n^q_{w1},1), \ldots, \hat{\boldsymbol{q}}(x,y;n^q_{w1},N_e), \ldots\hat{\boldsymbol{q}}(x,y;n^q_{wi},1), \ldots, \hat{\boldsymbol{q}}(x,y;n^q_{wi},N_e), \ldots]. \end{align}

In effect, to obtain $\hat {P}$, we horizontally stack all the $N_q$ number of $\hat {Q}_{n_w}$ that correspond to the $n_w$ in $\boldsymbol {n}^q$, i.e. $\hat {P} = [\hat {Q}_{n_{w1}^q}, \hat {Q}_{n_{w2}^q}, \ldots, \hat {Q}_{n_{wi}^q}, \ldots ]$. The SVD of $\hat {P}$ will give us the dominant structures that are responsible for the set of wavelets $\boldsymbol {n}^q$.

We denote the composite-WPOD energies for wavelets $\boldsymbol {n}^q$, obtained as the square of the singular values of $\hat {P}$, as $\mathcal {E}_i^q$. In other words, $\mathcal {E}_1^q$ is the energy of the dominant composite-WPOD mode that captures the dynamics across the set of wavelets $\boldsymbol {n}^q$. (To clarify the notation for the POD norms used in this section, we used $\mathcal {E}^{pod}$ for the POD norm, $\mathcal {E}^{wpod}$ for the WPOD norm and $\mathcal {E}$ for the composite-WPOD norm.) The first few such composite-WPOD modes for the set of wavelets $\boldsymbol {n}^q$ are shown in figure 7(a–d). The percentage energy captured by the mode is also shown in the figure. Similarly, composite-WPOD modes can be obtained for the wavelets $\boldsymbol {n}^b$ that contribute to the burst events, and here we denote the WPOD energies of these modes as $\mathcal {E}_i^b$. The first $10$ among these WPOD modes for the burst event are shown in figure 7(e–n). Figure 7(o) shows the cumulative contribution of the first $i$ modes to the total energy. In other words, the blue line in figure 7(o) shows $(\sum _{j=0}^{i}\mathcal {E}_j^q)/(\sum _{j=0}^{M_q}\mathcal {E}_j^q)$, where $M_q=N_e\times N_q$ is the total number of composite-WPOD modes obtained for the set of wavelets $\boldsymbol {n}^q$ for the quiet region. The red line shows the same quantity for the $M_b=N_e\times N_b$ composite-WPOD modes, corresponding to the set of wavelets $\boldsymbol {n}^b$ for the burst event. In this study $N_e=50$, $N_q=8$ and $N_b=39$ (see figure 5a), and therefore $M_q=400$ and $M_b=1950$. The horizontal dashed-dot line in figure 7(o) indicates $95\,\%$.

Figure 7.

The first few composite-WPOD modes are shown for two different sets of wavelets: (a–d) wavelets responsible for the quiet region marked by the blue-shaded region in figures 6(a) and 6(e–n) wavelets responsible for the burst event marked by the red-shaded regions in figure 6(a). The titles of these plots show the percentage energy captured by the mode computed as $100\times \mathcal {E}_i^q/\sum _{j=0}^M \mathcal {E}_j^q$ for the quiet region, and similarly for the burst events. (o) Also shown is the cumulative contribution of the first $i$ modes to the total energy for the quiet regions computed as $(\sum_{j=0}^{i}\mathcal{E}_j^q)/(\sum_{j=0}^{M_q}\mathcal{E}_j^q)$ (blue line, circle markers) and the burst events computed as $(\sum_{j=0}^{i}\mathcal{E}_j^b)/(\sum_{j=0}^{M_b}\mathcal{E}_j^b)$ (red line, square markers).

Let us first focus on the quiet region. From the blue line in figure 7(o) we see that the curve very quickly approaches the dashed-dot line indicating $95\,\%$. The first four modes capture most of the energy of the flow. These four modes are shown in figure 7(a–d), and they can be compared with their counterparts in figure 6(b–e). We note that the leading modes obtained from both WPOD versions are similar. Consistent with the observations from figure 6, composite-WPOD modes $1$ and $2$ correspond to the unstable eigenfunction and modes $3$ and $4$ represent shearing motions in the $y$-direction. (The modes in figure 7 exhibit greater convergence compared with those in figure 6 due to the inclusion of data from a collection of $n_w$, thereby increasing the input data used.)

Now consider the burst region. From the red line in figure 7(o), we see that the increase in total energy is more gradual. The first four modes, shown in figure 7(e–h), still capture a considerable part of the energy. Structurally, these four modes resemble the modes in the quiet region. This is consistent with the observation in § 4.2, and shows that the dominant modes in the quiet region persist during the burst events as well. In effect, this flow cannot be precisely divided into the quiet region and the burst events. However, going back to figure 7(o), the red line shows a very slow increase, and many suboptimal modes are required to reach $95\,\%$ energy (the horizontal dashed-dot line). Therefore, although these suboptimal modes (i.e. mode 5 and onward) each contribute very little energy, large numbers of them together play a significant role in the burst events. A few of these modes are shown in figure 7(i–n). The leading among these modes appear to be modified, here sheared, versions of the unstable eigenfunction (see for example figures 7(i), 7(j) and 7(k)). This is also consistent with the observations in § 4.2. Notably, such fragmented and sheared versions of the unstable eigenfunction are modes that are typical to the burst events.

4.4. Tracking composite-WPOD modes in a time series

So far, we have identified coherent structures that exist in the quiet region and the burst events. In this section, we track these coherent structures in a time series obtained from the flow. In other words, we are interested in analysing how the contributions of these coherent structures to the flow evolve with time. This will pave the way for the discussion in the next section (§ 4.5), where we will use these coherent structures to introduce the WPOD-based method for predicting the burst events.

We have two sets of composite-WPOD modes: (i) let $L_q$ represent the set of $M_q = N_e \times N_q$ number of composite-WPOD modes corresponding to the quiet region with the $j$th mode in $L_q$ having WPOD energy $\mathcal {E}_j^q$ and similarly (ii) let $L_b$ represent the set of $M_b = N_e \times N_b$ modes for the burst event with energies $\mathcal {E}_j^b$. Consider the modes $\boldsymbol {\phi }^{q}_j(x,y)$ ($\kern 1.5pt j={1,\ldots,M_q}$) from $L_q$. Now consider a time series $\boldsymbol {q}(x,y,t)$ obtained from flow, where $\boldsymbol {q}=(u,v)$ is the state vector. Note that $\boldsymbol {q}(x,y,t)$ can be a time series of arbitrary length with any number of burst events occurring at any point in time. Additionally, $\boldsymbol {q}(x,y,t)$ could lie outside the time window used to obtain the $N_e$ realisations of the burst event for WPOD (as required for the problem of predicting the burst events in the next section). At each time $t$, we aim to assess the presence of the structures $\boldsymbol {\phi }^{q}_j(x,y)$ in the flow field $\boldsymbol {q}(x,y,t)$. For this, we first compute $E_{qj}(t)$, which is the magnitude of energy shared between a mode $\boldsymbol {\phi }^{q}_j(x,y)$ and $\boldsymbol {q}(x,y,t)$ as

(4.6) \begin{equation} E_{qj}(t) := \int_{k_y} \int_{k_x} |{FT}_{xy}(\boldsymbol{q})(k_x,k_y;t) {FT}_{xy}(\boldsymbol{\phi}^{q}_j)(k_x,k_y;t)^*|^2 \,{\rm d} k_x \,{\rm d} k_y. \end{equation}

Here, ${FT}_{xy}({\cdot })$ represents the 2-D Fourier transform in the spatial directions, and $k_x$ and $k_y$ are the wavenumbers in the $x$ and $y$-directions, respectively. We normalise $E_{qj}(t)$ using (i) $E_{jj} := \langle \boldsymbol {\phi }^{q}_j, \boldsymbol {\phi }^{q}_j \rangle _{x,y}$ which is the energy of $\boldsymbol {\phi }^{q}_j(x,y)$ and (ii) $E_{qq}(t) := \langle \boldsymbol{q}, \boldsymbol{q} \rangle _{x,y}$ which is the energy of $\boldsymbol {q}(x,y,t)$.

The coherence $\gamma ^q(t)$ between $\boldsymbol {q}(x,y,t)$ and the modes in $L_q$ can now be defined as

(4.7a,b)\begin{equation} \gamma^q_j(t) = \frac{E_{qj}(t)}{E_{qq}(t) E_{jj}},\quad \gamma^q(t) = \sum_{j=1}^{M_q}{\left(\frac{\mathcal{E}_i^q}{\sum_{j=0}^{M_q}\mathcal{E}_j^q}\right) \gamma^q_j(t)}. \end{equation}

Here, $\gamma ^q_j(t)$ is the coherence between a point in the time series from the flow $\boldsymbol {q}(x,y,t)$ and the $j$th composite-WPOD mode in $L_q$. The value of $\gamma ^q_j(t)$ is bounded in the range $0 \leq \gamma ^q_j(t)\leq 1$, with $0$ indicating no coherence between $\boldsymbol {\phi }^{q}_j(x,y)$ and $\boldsymbol {q}(x,y,t)$ and $1$ indicating perfect coherence. Hence, $\gamma ^q(t)$ is the weighted average of the coherence across the $M_q$ different modes in $L_q$, weighted by the fraction of energy that each mode contributes to the energy of the quiet region. The value of $\gamma ^q(t)$ therefore also lies between $0$ and $1$, where $0$ indicates that the modes in $L_q$ are not present in $\boldsymbol {q}(x,y,t)$ at that time $t$ and $1$ indicates that the modes in $L_q$ are the only structures present in $\boldsymbol {q}(x,y,t)$. Similarly, we can also define $\gamma ^{b}(t)$ as the average coherence of $\boldsymbol {q}(x,y,t)$ with structures in $L_{b}$. It should be noted that, the actual value of $\gamma ^q(t)$ and $\gamma ^b(t)$ do not hold physical significance. Instead, our interest is in the trends of $\gamma ^q(t)$ and $\gamma ^b(t)$ over time. Let us first consider $\gamma ^q(t)$ for a sample time series. Figure 8(a) displays $\gamma ^q(t)$ in blue, along with the time series of $D(t)/D_{lam}$ in black. To obtain a smoother curve, the average of $\gamma ^{q}(t-2)$, $\gamma ^{q}(t-1)$ and $\gamma ^{q}(t)$ is computed at each time $t$. Similar averaging is also later carried out for $\gamma ^b(t)$ (at each time $t_1$, only values $t\leq t_1$ are used for this averaging). Looking at $\gamma ^q(t)$ we observe that: (i) within the quiet region, $\gamma ^q(t)$ remains relatively high and (ii) $\gamma ^q(t)$ plummets down when burst events occur. It is evident that $\gamma ^q(t)$ exhibits distinctive changes in its trends during a burst event. This is consistent with the observations in the literature that shows that the occurrence of a burst event is correlated with how near or far the flow is from equilibrium solutions of the flow (Farazmand Reference Farazmand2016; Page et al. Reference Page, Brenner and Kerswell2021). Looking ahead to the problem of predicting the burst events, $\gamma ^q(t)$ should, therefore, become a valuable tool.

Figure 8.

(a) Tracking coherent structures for the quiet region using $\gamma ^q$ (blue line). Two components of $\gamma ^q$ are also shown: $\gamma ^q_{1:2}(t)$ ($\triangle$) and $\gamma ^q_{3:4}(t)$ ($\triangledown$). (b) Also shown are the mean-removed and normalised profiles of $D(t)/D_{lam}$, $\gamma ^q_{1:2}(t)$ and $\gamma ^q_{3:4}(t)$ for the time window indicated by the grey-shaded box in (a).

To further probe this, figure 8(a) also shows two components of $\gamma ^q(t)$: (i) $\gamma ^q_{1:2}(t)$, which shows the weighted average of the coherence of just the first two modes in $L_q$, i.e. the unstable eigenfunction (modes in figure 7a,b) and (ii) $\gamma ^q_{3:4}(t)$, which is the weighted average coherence of modes $3$ and $4$, i.e. the shearing motions. (modes in figure 7c,d). The first significant observation is that $\gamma ^q_{1:2}(t)$ closely follows the trends of $\gamma ^q(t)$. This is not surprising given that, together, modes $1$ and $2$ capture more than $70\,\%$ of the energy of the quiet region (as seen in figure 7). Therefore, we can say that the trends of $\gamma ^q(t)$ are most significantly impacted by the flow due to modes $1$ and $2$, i.e. the unstable eigenfunction. The second noteworthy observation is that $\gamma ^q_{3:4}(t)$ tends to move out of phase with $\gamma ^q_{1:2}(t)$. This becomes more evident in figure 8(b) where mean-removed and normalised profiles of $D(t)/D_{lam}$, $\gamma ^q_{1:2}(t)$ and $\gamma ^q_{3:4}(t)$ are shown. The curves are vertically shifted in order to make the trends clearer. From this figure, we see that, when $D(t)/D_{lam}$ increases, the presence of the flow due to the unstable eigenfunction ($\gamma ^q_{1:2}(t)$) decreases and that of the shearing structure ($\gamma ^q_{3:4}(t)$) increases. This observation provides additional support to our initial hypothesis that the shearing motions disrupt the flow generated by the unstable eigenfunction, thereby increasing dissipation.

Let us now examine $\gamma ^b(t)$ for the same time series. Figure 9 shows $\gamma ^b(t)$ in red alongside the time-series of $D(t)/D_{lam}$ in black. Upon comparing $\gamma ^b(t)$ with $\gamma ^q(t)$ in figure 8(a), the first apparent observation is that the trends between them are strikingly similar. This observation is consistent with the earlier finding from figure 7, that nearly $60\,\%$ of the energy of the burst events is captured by modes that are present in the quiet region. Looking ahead to the problem of predicting the burst events, this similarity between $\gamma ^q(t)$ and $\gamma ^b(t)$ brings the unfortunate conclusion that the trends of the full $\gamma ^b(t)$ may not be useful for prediction. However, our earlier observations showed the existence of suboptimal structures that are unique to the burst events (see figure 7). Presumably, these structures exhibit trends dissimilar to $\gamma ^q(t)$. The question, therefore, is whether we can identify a component of $\gamma ^b(t)$ that has predictive value?

Figure 9.

Tracking coherent structures for the burst events using $\gamma ^b$ (red line). Two components of $\gamma ^b$ are also shown: $\gamma ^b_{1:4}(t)$ ($\triangle$) and $\gamma ^b_{5:N_b}(t)$ ($\triangledown$).

To explore this further, figure 9 also shows two components of $\gamma ^b(t)$: (i) $\gamma ^b_{1:4}(t)$ computed using just the first four modes in $L_b$, i.e. the modes that resemble those from the quiet region (modes in figure 7e–h) and (ii) $\gamma ^b_{5:N_b}(t)$ computed using modes $5$ and onward (the first few of these modes are those in figure 7i–n). Two observations become immediately apparent. Firstly, as expected, $\gamma ^b_{1:4}(t)$ follow trends similar to the full $\gamma ^b(t)$ as well as $\gamma ^q(t)$. Secondly, and also more interestingly, $\gamma ^b_{5:N_b}(t)$ increases during the burst event. This again shows that these suboptimal modes are unique to the burst events. Again looking ahead to the problem of predicting the burst events, these observations show that, although the full $\gamma ^b(t)$ may not be helpful in identifying the burst events, the distinctive trends of $\gamma ^b_{5:N_b}(t)$ during the burst event will be. (Note, Appendix C shows that it is possible to predict the burst events using just $\gamma ^q(t)$, ignoring the trends of $\gamma ^b(t)$. However, incorporating the trends of $\gamma ^b_{5:N_b}(t)$ does improve the obtained predictions.)

4.5. Obtaining WPOD-based predictions of burst events

In this section, to predict burst events, we use the distinctive trends of $\gamma ^q(t)$ and $\gamma ^b_{5:N_b}(t)$ identified in the previous section. For this purpose, a predictor $\lambda =\gamma ^{q}(t)-\gamma _{5:N_b}^{b}(t)$ is defined. The average of $\lambda (t-2)$, $\lambda (t-1)$ and $\lambda (t)$ is computed at each time $t$. Figure 10 shows $D(t)/D_{lam}$ along with the evolution of $\lambda$ for three different time windows. The time window in figure 10(a) was included among the $N_e=50$ realisations of the burst events that were used for computing the WPOD modes, while the time windows in figures 10(b) and 10(c) fall outside these $N_e$ realisations. In other words, while the data in figure 10(a) fall within the ‘training dataset’, figures 10(b) and 10(c) fall outside this training data, and will therefore serve to illustrate the generalisability of the prediction method discussed here. Since we are interested in the trends of $\lambda$, and not the magnitudes, we here consider normalised $\lambda$ (normalised by the minimum and maximum values obtained within a time window $t=0\unicode{x2013}15\,000$).

Figure 10.

Predictions of the burst events obtained using the WPOD-based method for three separate time series from DNS. The green-shaded regions indicate the identified burst regions, where the predictor $\lambda$ (grey line) goes below the defined threshold value, i.e. $\lambda <0.95\lambda _t$. The vertical green lines mark the onset of these burst regions, and therefore represent the predictions of the burst event obtained from the WPOD-based method. Predicted times are compared with the black dashed vertical lines that indicate the predictions obtained from the energy-based method.

From figure 10, we see that $\lambda$ decreases during burst events. Consequently, we designate a threshold $\lambda _t$, and values of $\lambda$ below $\lambda _t$ will be identified as burst events. To calculate $\lambda _t$, we begin by computing the mean of $\lambda$ minus the variance, where both these statistics are calculated for time instances that correspond to the quiet regions (i.e. when $D(t)/D_{lam} \leq 0.1$) in the range $t=0\unicode{x2013}15\,000$. The threshold $\lambda _t$ is set to be $0.95$ times this mean minus the variance. (The impact of varying this coefficient $0.95$ to other values will be discussed in § 6.1.) The green-shaded regions in figure 10 mark the times when $\lambda <\lambda _t$, and the solid green vertical lines mark the beginnings of these regions. A majority of these green-shaded regions correspond to burst events. The green vertical lines are therefore the predictions of the burst events obtained using this WPOD-based method.

To assess the prediction performance, the obtained prediction times are here compared with those obtained from a more straightforward strategy of tracking the energy of the flow. In this case, the predictor becomes $\lambda = -D(t)/D_{lam}$ (where the negative sign makes comparison with the wavelet-based predictors more straightforward). Again, the average of $\lambda (t-2)$, $\lambda (t-1)$ and $\lambda (t)$ is computed at each time $t$. Using the same procedure as for the WPOD-based method, a threshold is computed as $0.95$ times the mean of $\lambda$ minus the variance computed for the quiet region in the range $t=0\unicode{x2013}15\,000$. The vertical black dashed lines in figure 10 indicate the beginnings of the time windows where the predictor goes below the threshold value. These lines therefore represent the predictions obtained from tracking the energy of the flow, which can directly be compared with the predictions obtained from the WPOD-based method in green.

From figure 10, we see that the WPOD-based method is able to predict the burst events well. However, for this chosen threshold, we obtain false positives, i.e. predictions of oncoming burst events in the absence of such events (e.g. $t \approx 18\,325$ in figure 10c). When compared with the energy-based method, the WPOD-based method does seem to give improved prediction times (see, for instance, $t \approx 17\,235$ and $t \approx 18\,355$ in figure 10). However, we also seem to obtain more false positives from the WPOD-based method.

To quantify the prediction performance, we should consider predictions over long time windows. Here we take a time window of length $T=20\,000$ and $dt=1$ outside the $N_e$ realisations used for WPOD (i.e. outside the ‘training dataset’), which contains $133$ burst events. We compute three quantities. First, the average prediction time $\tau$, which is computed as the difference between the obtained prediction time and the time when the burst begins, i.e. when $D(t)/D_{lam}\geq 0.15$. (Note, $\tau$ is shown in units of $\sqrt {L_y/2{\rm \pi} \zeta }$.) Second, the percentage of predictions that are false positives, denoted by FP%. A false positive is here defined as a prediction where $D(t)/D_{lam}$ of the flow does not go above $0.15$. Third, the percentage of false negatives, denoted by FN%. These are instances when the method does not identify a region with $D(t)/D_{lam}\geq 0.15$ as a burst event. Both FP% and FN% are represented as a percentage of the total number of predictions of the burst event obtained from the method. For the WPOD-based method $\tau =1.06$, ${\rm FP}\%=20$ and ${\rm FN}\%=2$, which can be compared with the values for the energy-based method $\tau =0.36$, ${\rm FP}\%=13$ and ${\rm FN}\%=0$. The WPOD-based method has difficulty in capturing burst events that persist for less than 1–2 time units, and these contribute to the very small number of false negatives obtained. From these numbers, we can conclude that the WPOD-based method is capable of predicting burst events with improved prediction times compared with the more straightforward method of tracking energy, albeit with the possibility of slightly increased number of false positives and false negatives.

5.1. Resolvent operator

To define the resolvent operator, we start with the Fourier-transformed linearised Navier–Stokes equations that were laid out in (2.4) and reproduced here

(5.1)\begin{equation} \dot{\tilde{\omega}}(y,t;k_x) + \boldsymbol{A}\tilde{\omega}(y,t;k_x) = \tilde{F}(y,t;k_x). \end{equation}

Let us consider the basis vectors $\varTheta _{n_w}(t)$. We can write $\tilde {\omega }$ in terms of its projection onto this basis as (Farge Reference Farge1992)

(5.2)\begin{equation} \tilde{\omega}(y,t;k_x) = \sum_{n_w=1}^{N_{w}} \hat{\tilde{\omega}}(y;k_x,n_w) \varTheta_{n_w}(t), \end{equation}

where

(5.3)\begin{equation} \hat{\tilde{\omega}}(y;k_x,n_w) = \int_0^{T} \tilde{\omega}(y,t;k_x) \varTheta_{n_w}^*(t)\,{\rm d} t. \end{equation}

We can similarly obtain the coefficients $\hat {\tilde {F}}(y;k_x,n_w)$ of $\tilde {F}(y,t;k_x)$. Here, $T$ is the length of the time window used. For the case considered here, $T$ (here $T=200$) is the length of a realisation of the burst event that we used for WPOD (see § 3.3). Please note that, for the rest of this manuscript, just to simplify notation, we replace the notation $(\widehat {\tilde {\cdot }})$ with $(\hat {\cdot })$ and also denote $n_w$ as a subscript. In other words, we use $\hat {\omega }_{n_w}(y;k_x)$ and $\hat {F}_{n_w}(y;k_x)$ to represent $\hat {\tilde {\omega }}(y;k_x,n_w)$ and $\hat {\tilde {F}}(y;k_x,n_w)$, respectively.

The basis $\varTheta _{n_w}(t)$ is chosen such that it is complete and orthonormal, i.e. $\int _0^T \varTheta ^*_{m_w} \varTheta _{n_w} \,{\rm d} t = \delta _{n_wm_w}$, where $\delta _{n_wm_w}$ is the Dirac delta function. The basis could, for instance, be the Fourier basis, where $\varTheta _{n_w}(t)= \exp {{\rm i}(2{\rm \pi} /T) n_w t}$, or as will be considered here, the wavelet basis discussed in § 3.2. In terms of these transforms, the linearised equations in (5.1) become

(5.4)\begin{align} & \sum_{n_w=1}^{N_w} \hat{\omega}_{n_w}(y;k_x) \left(\frac{{\rm d}\varTheta_{n_w}(t)}{{\rm d} t}\right) + \boldsymbol{A} \sum_{n_w=1}^{N_w} \hat{\omega}_{n_w}(y;k_x)\varTheta_{n_w}(t) \nonumber\\ &\quad =\sum_{n_w=1}^{N_w} \hat{F}_{n_w}(y;k_x)\varTheta_{n_w}(t). \end{align}

Here, $N_w$ is the number of wavelet coefficients obtained from the wavelet transform (determined by the type of wavelet chosen and the number of time instances $N_t$ chosen for the time window). The bases vectors $\varTheta _{n_w}(t)$ are not guaranteed to be continuous or differentiable. For instance, the discrete Daubechies $1$ wavelet that we use (see § 3.2) is discontinuous and not differentiable. To consider the formulation in terms of the most general bases, including discontinuous ones, we re-write (5.4) in terms of the integrals of the bases $\varTheta _{n_w}(t)$, therefore giving us

(5.5)\begin{align} & \sum_{n_w=1}^{N_w} \hat{\omega}_{n_w}(y;k_x)\varTheta_{n_w}(t) + \boldsymbol{A} \sum_{n_w=1}^{N_w} \hat{\omega}_{n_w}(y;k_x) \left(\int_0^{t} \varTheta_{n_w}(t') \,{\rm d} t'\right) \nonumber\\ &\quad =\sum_{n_w=1}^{N_w} \hat{F}_{n_w}(y;k_x) \left(\int_0^{t} \varTheta_{n_w}(t') \,{\rm d} t'\right). \end{align}

Since the bases are complete, the integral of a particular basis vector can be written as a linear superposition of all basis vectors, i.e. $\int _0^{t} \varTheta _{n_w}(t') \,{\rm d} t' = \sum _{m_w=1}^{N_w} c_{n_wm_w}\varTheta _{m_w}(t)$, where $c_{n_wm_w}$ are scalar coefficients. This gives

(5.6)\begin{align} & \sum_{n_w=1}^{N_w} \hat{\omega}_{n_w}(y;k_x)\varTheta_{n_w}(t) + \boldsymbol{A} \sum_{n_w=1}^{N_w} \hat{\omega}_{n_w}(y;k_x) \left(\sum_{m_w=1}^{N_w} c_{n_wm_w}\varTheta_{m_w}(t)\right) \nonumber\\ &\quad =\sum_{n_w=1}^{N_w} \hat{F}_{n_w}(y;k_x) \left(\sum_{m_w=1}^{N_w} c_{n_wm_w}\varTheta_{m_w}(t)\right). \end{align}

Using the orthogonality of the basis vectors, the equation is now rewritten as

(5.7)\begin{equation} \hat{\omega}_{n_w}(y;k_x) + \boldsymbol{A} \sum_{m_w=1}^{N_w} c_{n_wm_w} \hat{\omega}_{m_w}(y;k_x) = \sum_{m_w=1}^{N_w} c_{n_wm_w} \hat{F}_{m_w}(y;k_x). \end{equation}

If we were dealing with the Fourier-based equations $c_{n_wm_w} = in_w\delta _{n_wm_w}$. The summations in (5.7) thereby drop off, giving independent equations for each $n_w$. This is not true for the case of all bases, such as a wavelet bases, where the response at one particular wavelet $\varTheta _{n_w}$ cannot be isolated from the responses at the other wavelets. In matrix form, this equation will therefore become

(5.8)\begin{align} & \begin{bmatrix} \hat{\omega}_1(y;k_x) \\ \hat{\omega}_2(y;k_x) \\ \hat{\omega}_3(y;k_x) \\ \vdots \end{bmatrix}+\underbrace{\begin{bmatrix} \boldsymbol{A} & 0 & 0 & \ldots \\ 0 & \boldsymbol{A} & 0 & \ldots \\ 0 & 0 & \boldsymbol{A} & \ldots \\ \vdots & \vdots & \vdots \end{bmatrix} \begin{bmatrix} c_{11}\boldsymbol{I} & c_{12}\boldsymbol{I} & c_{13}\boldsymbol{I} & \ldots \\ c_{21}\boldsymbol{I} & c_{22}\boldsymbol{I} & c_{23}\boldsymbol{I} & \ldots \\ c_{31}\boldsymbol{I} & c_{32}\boldsymbol{I} & c_{33}\boldsymbol{I} & \ldots \\ \vdots & \vdots & \vdots \end{bmatrix}}_{\bar{\boldsymbol{A}}(y,n_w;k_x)} \begin{bmatrix} \hat{\omega}_1(y;k_x) \\ \hat{\omega}_2(y;k_x) \\ \hat{\omega}_3(y;k_x) \\ \vdots \end{bmatrix} \nonumber\\ &\quad =\underbrace{\begin{bmatrix} c_{11}\boldsymbol{I} & c_{12}\boldsymbol{I} & c_{13}\boldsymbol{I} & \ldots \\ c_{21}\boldsymbol{I} & c_{22}\boldsymbol{I} & c_{23}\boldsymbol{I} & \ldots \\ c_{31}\boldsymbol{I} & c_{32}\boldsymbol{I} & c_{33}\boldsymbol{I} & \ldots \\ \vdots & \vdots & \vdots \end{bmatrix}}_{\bar{\boldsymbol{B}}(y,n_w;k_x)} \begin{bmatrix} \hat{F}_1(y;k_x) \\ \hat{F}_2(y;k_x) \\ \hat{F}_3(y;k_x) \\ \vdots \end{bmatrix} . \end{align}

We are considering the discretised equations with $N_y$ number of grid points in the $y$-direction and $N_w$ wavelets. The matrix $\boldsymbol {A}$ in (5.8) is therefore of size $N_y\times N_y$ and $\boldsymbol {I}$ represents an identity matrix of the same size. The projection coefficients $c_{n_wm_w}$ are scalars. The matrices $\bar {\boldsymbol {A}}(y,n_w;k_x)$ and $\bar {\boldsymbol {B}}(y,n_w;k_x)$ are therefore of sizes $(N_wN_y)\times (N_wN_y)$, and the vectors $\hat {\omega }=[\hat {\omega }_1,\hat {\omega }_2,\ldots \hat {\omega }_{N_w}]$ and $\hat {F}=[\hat {F}_1,\hat {F}_2,\ldots \hat {F}_{N_w}]$ are each of size $N_y N_w$. Using these definitions, we can write the wavelet-based linearised equation as

(5.9)\begin{equation} \hat{\omega}(y,n_w;k_x) + \bar{\boldsymbol{A}}(y,n_w;k_x) \hat{\omega}(y,n_w;k_x) = \bar{\boldsymbol{B}}(y,n_w;k_x) \hat{F}(y,n_w;k_x) .\end{equation}

To obtain the state-vector $\hat {\boldsymbol {q}}(y,n_w;k_x) = (\hat {u}(y,n_w;k_x), \hat {v}(y,n_w;k_x))$, we introduce an output matrix $\bar {\boldsymbol {C}}$ such that

(5.10)\begin{equation} \hat{\boldsymbol{q}}(y,n_w;k_x) = \bar{\boldsymbol{C}}(y,n_w;k_x) \hat{\omega}(y,n_w;k_x). \end{equation}

The matrix $\bar {\boldsymbol {C}}$ can also be used to ‘mask’ the resolvent. For instance, in the upcoming sections, we will compute the responses at specific values of $n_w$. By selectively including and excluding rows of the matrix $\bar {\boldsymbol {C}}$, we can mask the response such that we obtain responses exclusively at a particular $n_w$ as in § 5.2, or to obtain responses at specific sets of $n_w$ as in § 5.3. It should be noted that, similar to this, the matrix $\bar {\boldsymbol {B}}$ can be used to mask the forcing to the resolvent operator. However, in the present study, we do not mask the forcing. This is because, while we are interested in responses at specific values of $n_w$ that we can separately trace back to the quiet region or the burst events, we do not have a physical argument for restricting the forcing similarly. In the wavelet resolvent, the response at a particular $n_w$ can be forced by all other $n_w$.

From this point, the procedure to obtain the resolvent operator follows McKeon & Sharma (Reference McKeon and Sharma2010). Rearranging (5.9) gives us

(5.11)\begin{equation} \hat{\boldsymbol{q}}(y,n_w;k_x) = \underbrace{\bar{\boldsymbol{C}} [\boldsymbol{W}^{1/2}(I + \bar{\boldsymbol{A}})^{{-}1} \boldsymbol{W}^{{-}1/2}] \bar{\boldsymbol{B}}}_{\bar{\boldsymbol{H}}(y,n_w;k_x)} \hat{F}(y,n_w;k_x). \end{equation}

Here $\bar {\boldsymbol {H}}(y,n_w;k_x)$ represents the resolvent operator. If a non-uniform grid is used to discretise $y$, the weight matrix $\boldsymbol {W}$ will contain the weights corresponding to the grid. In this study, we employ a Fourier grid with $N=51$ grid points, resulting in $\boldsymbol {W}=\boldsymbol {I}$. As for the WPOD (see § 3.2), here the Daubechies 1 wavelets serve as the basis vectors $\varTheta _{n_w}$, with a time window of $T=200$ and a time step of $dt=1$. A similar wavelet-based operator as that obtained in (5.11) has recently also been used in Ballouz et al. (Reference Ballouz, Dawson and Bae2023a,Reference Ballouz, Lopez-Doriga, Dawson and Baeb).

5.2. The WRA response and forcing modes

Singular value decomposition is used to analyse the resolvent operator, $\bar {\boldsymbol {H}}(y,n_w;k_x)$

(5.12)\begin{equation} \bar{\boldsymbol{H}}(y,n_w;k_x) = \sum_{i=1}^{N_y N_w} \boldsymbol{\psi}_i(y,n_w) \sigma_i \phi_i(y,n_w). \end{equation}

The singular values are arranged in increasing order such that $\sigma _i \geq \sigma _{i+1}$. The left singular vectors $\boldsymbol {\psi }_i(y,n_w)$ represent the resolvent response modes and the right singular vectors $\phi _i(y,n_w)$ represent the resolvent forcing modes. Essentially, the forcing $\phi _i(y,n_w)$, when put through the resolvent operator, yields a response $\boldsymbol {\psi }_i(y,n_w)$ amplified by $\sigma _i$. Both $\boldsymbol {\psi }_i(y,n_w)$ and $\phi _i(y,n_w)$ are functions of $y$ and the wavelets $n_w$. The most amplified response is $\boldsymbol {\psi }_1(y,n_w)$ corresponding to $\sigma _1$, and the corresponding most sensitive forcing direction is $\phi _1(y,n_w)$.

First, let us consider resolvent response modes at specific values of $n_w$. To obtain these modes, we can selectively mask specific rows of the output matrix $\bar {\boldsymbol {C}}$ (see (5.10)). For the resolvent, we need to choose a wavenumber in the $x$-direction $k_x$. The dominant modes in this flow have $k_x=1$ (see figure 7). Therefore, in figure 11, we show resolvent modes at $k_x=1$. (The shearing motions in figure 7 can be obtained as $k_x\approx 0$ modes, which are not shown here for the sake of brevity.) Figure 11 shows the resolvent response modes for two values of $n_w$: (i) $n_w=2$ that corresponds to the quiet region and (ii) $n_w=22$ that corresponds to the burst event. The percentages in the title of the figure represent the fraction of energy (at that $n_w$) captured by the respective mode $\boldsymbol {\psi }_i(y,n_w)$ computed as $\sigma _i^2/\sum _j{\sigma _j^2}$.

Figure 11.

First five resolvent response modes (a–e) at $n_w=2$ that contribute to the quiet region, and (f–j) at $n_w=22$ that contribute to the burst region, are shown. The titles show the fraction of energy (at that $n_w$) captured by the respective mode, computed as $\sigma _i^2/\sum _j{\sigma _j^2}$ (rounded off to eight decimal places for a–e).

Consider the modes from the quiet region in the first row of figure 11. The first resolvent mode captures almost the full energy at this $n_w$ (the percentage is rounded off to eight decimal places). The mode corresponds to the unstable eigenfunction of the linearised Navier–Stokes equations (see Appendix A) (Chandler & Kerswell Reference Chandler and Kerswell2013). Crucially, when compared with the first row of figure 6, the first resolvent mode resembles the first WPOD mode. In § 4.4, we found that this first WPOD mode dominates the dynamics of the quiet region in DNS. Therefore, for the quiet region, we can conclude that the model is able to predict the dominant structure. Additionally, the model also recognises the significance of this structure, as indicated by the mode capturing almost all the energy at this $n_w$.

Now consider the modes for the burst event in the second row of figure 11. The first resolvent mode for this $n_w$ is the same as for the quiet region, representing the unstable eigenfunction and capturing the majority of the energy at this $n_w$. A similar trend was apparent when considering the WPOD modes in figure 6, where we found that the same modes as in the quiet region were also dominant for the burst event, and captured $60\,\%$ of the energy of the burst events. In the case of the resolvent modes, however, the quiet region mode captures almost the full energy. Examining the structure of the suboptimal resolvent modes, inclined and fragmented versions of the unstable eigenfunction are apparent, reminiscent of the WPOD modes. For instance, compare figure 7(j) with 11(g). Therefore, for the burst events, the suboptimal resolvent modes do capture the relevant structures. However, the model is not able to identify the significance of these suboptimal modes, as suggested by the percentage of energy captured by them.

5.3. Reconstructing DNS data using WRA modes

Our purpose in introducing the resolvent here is to find the forcing structures, i.e. the structures in the nonlinear terms, that produces the quiet region and burst events. Studies have used the Fourier-based resolvent operators to find the underling forcing structures that are important in a flow (e.g. Towne et al. Reference Towne, Colonius, Jordan, Cavalieri and Bres2015; Karban et al. Reference Karban, Martini, Cavalieri, Lesshafft and Jordan2022). The right and left singular vectors of the resolvent, $\phi _i$ and $\boldsymbol {\psi }_i$, form complete bases. Consequently, any response $\hat {\boldsymbol {q}}$, and any forcing $\hat {F}$, that is obtained from data (here from DNS) can be expressed in terms of these basis vectors as

(5.13)\begin{equation} \left.\begin{gathered} \hat{F} = \sum_{i=1}^{N_yN_w} \chi_i \phi_i,\quad \chi_i = \sum_{n_w=1}^{N_w} \int_0^{L_y} \hat{F}_{n_w}(y;k_x) (\phi_i)_{n_w}(y) \,{{\rm d}y} \\ \hat{\boldsymbol{q}} = \sum_{i=1}^{N_yN_w} \chi_i \sigma_i \boldsymbol{\psi}_i,\quad \chi_i \sigma_i = \sum_{n_w=1}^{N_w} \int_0^{L_y} \hat{\boldsymbol{q}}_{n_w}(y;k_x) (\boldsymbol{\psi}_i)_{n_w}(y)\,{{\rm d}y}, \end{gathered}\right\} \end{equation}

where $(\phi _i)_{n_w}(y)$ and $(\boldsymbol {\psi }_i)_{n_w}(y)$ are just representing the components of $\phi _i(y,n_w)$ and $\boldsymbol {\psi }_i(y,n_w)$ at a particular $n_w$. Therefore, given DNS data $\tilde {\boldsymbol {q}}(y,t;k_x)$, we first obtain the wavelet coefficients $\hat {\boldsymbol {q}}_{n_w}(y;k_x)$ through a wavelet transform, and thereafter compute $\chi _i \sigma _i$ by projecting the wavelet coefficients onto the resolvent response modes $\boldsymbol {\psi }_i(y,n_w)$. A reduced-order description of the flow can then be obtained by reconstructing the data using a truncated basis

(5.14)\begin{equation} \hat{\boldsymbol{q}}_{approx} = \sum_{i=1}^{N_1} \chi_i \sigma_i \boldsymbol{\psi}_i, \end{equation}

where $N_1$ is the number of modes that are retained in the truncated data. The approximate data $\boldsymbol {q}_{approx}$ are then obtained from an inverse wavelet transform in time and an inverse Fourier transform in the $x$-direction.

From the obtained product $\chi _i \sigma _i$, we then extract $\chi _i$. The $N_1$ is chosen such that we eliminate modes which contribute less than $0.1\,\%$ energy to any particular $k_x$ (i.e. the eliminated modes have $\chi _i\chi _i^* \sigma _i^2$ less than $0.1\,\%$ of the energy at that $k_x$). This $\chi _i$ is then used to compute the forcing $\hat {F}_{approx} = \sum _{i=1}^{N_1} \chi _i \phi _i$. The obtained forcing $\hat {F}_{approx}$ when put through the resolvent operator will give a response $\hat {\boldsymbol {q}}_{approx}$. In other words, we can find the component of the full forcing $\hat {F}$ from DNS (i.e. the nonlinear term) that is specifically responsible for amplifying the response $\hat {\boldsymbol {q}}_{approx}$. The forcing $F_{approx}$ (in time) can then be obtained from an inverse wavelet transform in time and an inverse Fourier transform in the $x$-direction. (It should be noted that, unlike in the case of the Fourier resolvent, the forcing obtained from the wavelet resolvent, $\hat {F}_{approx}$, does not represent a particular spatial mode shape. Instead, $\hat {F}_{approx}$ represents mode shapes at different wavelets $n_w$. This is because, unlike the Fourier resolvent, for the wavelet resolvent, response at a particular $n_w$ can be forced by multiple separate wavelets simultaneously.)

To obtain $F_{approx}$ specifically for the quiet region, we obtain resolvent modes for the set of wavelets $\boldsymbol {n}_w^q$ corresponding to the quiet region. The matrix $\bar {\boldsymbol {C}}$ can be used to mask the resolvent so as to obtain a response solely for this particular set of wavelets. The wavelet coefficients obtained from DNS for $\boldsymbol {n}_w^q$ are then projected onto the obtained resolvent modes. Thereafter, following the procedure above, we can obtain the forcing for the quiet region. Similarly, we can also obtain the forcing for the burst region. An illustration of this procedure is shown in figure 12, where figures 12(a) and 12(b) show the vorticity $\omega$ from DNS decomposed into the quiet (blue) and the burst (red) regions, respectively. The corresponding resolvent-based reconstructions $\omega _{approx}$ (black solid lines) are also shown, where $(\omega _q)_{approx}$ and $(\omega _b)_{approx}$ represent the reconstructions of the quiet region and the burst events, respectively. Additionally, the forcing $F_{approx}$ that generate these responses are shown (black dashed lines) in these figures. Using WRA, we have therefore constructed forcing $(F_q)_{approx}$ and $(F_b)_{approx}$ that separately generate the quiet region and the burst events, respectively.

Figure 12.

Values of $D(t)/D_{lam}$ computed from the full vorticity data for $k_x={0,1,2,3}$ (grey line) are shown alongside the wavelet-based decomposition of (a) the quiet region in blue and (b) the burst event in red. The decomposed data are compared with their respective resolvent-based reconstructions (solid black lines). The forcing $F_{approx}$ that generates these responses are also shown (dashed black lines). (The $y$-axis labels on the right in both (a) and (b) correspond to the forcing.)

5.4. Coherent structures of the forcing

Our aim is to identify the coherent structures in the forcing that generates the quiet region and the burst events. In order to find these coherent structures, we first start by using the resolvent to reconstruct the flow field as illustrated in figure 12, and thereby obtain the corresponding forcing. This reconstruction is performed for $M$ different realisations of the quiet region and the burst events (here $M=50$). Consequently, we obtain $M$ realisations of resolvent-based vorticity reconstructions $(\omega _{approx,1}(x,y,t),\omega _{approx,2}(x,y,t),\ldots,\omega _{approx,M}(x,y,t))$ and the corresponding forcing $(F_{approx,1}(x,y,t),F_{approx,2}(x,y,t),\ldots,F_{approx,M}(x,y,t))$. It is important to note that, for the reconstructions done here, $k_x=(0,1,2,3)$ wavenumbers in the $x$-direction are utilised. Subsequently, to extract the coherent structures, we perform a POD on the $M$ realisations of the forcing. This involves constructing a data matrix $P$

(5.15) \begin{align} P &= [F_{approx,1}(x,y,t_1), F_{approx,1}(x,y,t_2), \ldots, F_{approx,1}(x,y,T), \nonumber\\ &\quad\ \ \ F_{approx,2}(x,y,t_1), F_{approx,2}(x,y,t_2), \ldots, F_{approx,2}(x,y,T), \ldots]. \end{align}

A SVD of the data matrix $P$ gives us the POD modes. We obtain such POD modes for the quiet region and the burst events.

The top row of figure 13 shows the first few forcing POD modes obtained for the quiet region. The titles of these figures show the percentage of energy that is captured by the mode. Notably, the first two modes capture a majority of the energy, and the modes resemble the unstable eigenfunction itself. This again suggests that the dominant amplification mechanism is a normal-mode mechanism, where the forcing and response modes can be expected to be structurally similar. Interestingly, modes 3 and 4 resemble the dominant resolvent forcing to the unstable eigenfunction (see Appendix A), which suggests that there is a second route through which the unstable eigenfunction can be amplified. We can also obtain the forcing POD modes for the burst events, and the first five of these modes are shown in the bottom row of figure 13. Notably, inclined structures dominate the forcing to the burst events.

Figure 13.

The first few POD modes obtained from the forcing corresponding to the resolvent reconstruction of (a–d) the quiet region and (e–i) the burst events. The titles of these plots show the percentage energy captured by the mode computed as $100\times \mathcal {E}_i^q/\sum _{j=0}^M \mathcal {E}_j^q$ for the quiet region, and similarly for the burst events.

5.5. Obtaining WRA-based predictions of burst events

Now that we have identified the coherent structures in the forcing that generate the quiet regions and the burst events, we will track them in a time series obtained from the flow. More specifically, our interest is in analysing how the contributions of these structures to the nonlinear terms of the flow evolve with time. We have two sets of forcing coherent structures: (i) let $F_q$ represent the modes corresponding to the quiet region, with each mode $\phi _i^q(x,y)$ having POD energy $\mathcal {E}_i^q$, and (ii) $F_b$ the modes for the burst event, with each mode $\phi _i^b(x,y)$ having POD energy $\mathcal {E}_i^b$. Both $F_q$ and $F_b$ contain $M \times N_t$ modes, where $M$ is the number of resolvent reconstructions used (here $M=50$) and $N_t$ is the number of time snapshots of DNS data in each realisation (here, $N_t=200$). To simplify computation, here, we ignore modes that cumulatively contribute less than $1\,\%$ energy (i.e. POD modes $k$ with $\sum _{i=k}^M \mathcal {E}_i^q/\sum _{j=0}^M \mathcal {E}_j^q < 0.01$). At each time $t$, we aim to assess the presence of the structures in $F_q$ and $F_b$ in the nonlinear terms $F(x,y,t)$ obtained from the flow.

For this purpose, we follow the methodology laid out in § 4.4, with $\boldsymbol {q}(x,y,t)$ now taken to be the nonlinear terms $F(x,y,t)$ obtained from DNS. As explained in § 4.4, we obtain $\gamma ^q(t)$, which is the weighted average of the coherence of the modes in $F_q$ with each point in the time series $F(x,y,t)$. Equivalently, we also obtain $\gamma ^b(t)$ for the modes in $F_b$. (As in § 4.4, we are here interested in the trends of $\gamma ^q(t)$ and $\gamma ^b(t)$, and not the actual values.) Figure 14 shows the profiles of $\gamma ^q(t)$ and $\gamma ^b(t)$. We see that $\gamma ^q(t)$ and $\gamma ^b(t)$ follow distinctive trends during a burst event, with $\gamma ^q(t)$ (blue line) decreasing and $\gamma ^b(t)$ (red line) increasing. The trends of $\gamma ^q(t)$ and $\gamma ^b(t)$ can therefore be used for predicting the burst event.

Figure 14.

Tracking forcing structures in a time series. The evolution of $\gamma ^q(t)$ (blue) and $\gamma ^b(t)$ (red) obtained by tracking forcing structures for the quiet region and the burst events, respectively, are shown. Here, $D(t)/D_{lam}$ is also shown (in black).

Like in § 4.5, here, we define a predictor $\lambda = \gamma ^q(t) - \gamma ^b(t)$. Figure 15 shows the evolution of this predictor for three different time windows. (As in § 4.5, here, $\lambda$ is normalised by the minimum and the maximum values obtained within the time window $t=0\unicode{x2013} 15\,000$.) Using the same method as in § 4.5, we again choose a threshold $\lambda _t$. The times when $\lambda <\lambda _t$ are identified as predictions of the burst regions, and the red-shaded regions in figure 15 denote these predictions. These prediction times are compared with the green dashed lines, which denote the predictions obtained from the WPOD-based method in § 4.5. First, we note that, by tracking the forcing structures, we are indeed able to obtain predictions of the burst event. These prediction times are improvements over the WPOD method. We also note that false positives are obtained using this method as well (e.g. $t\approx 17\,195$ in figure 15). As in § 4.5, for a more accurate comparison, we compute the average prediction time $\tau$, the percentage of false positives FP% and the percentage of false negatives FN%. For the WRA-based method $\tau =1.88$, $\textrm {FP}\%=25$ and $\textrm {FN}\%=4$ (compared with $\tau =1.06$, $\textrm {FP}\%=20$ and $\textrm {FN}\%=2$ for the WPOD-based method and $\tau =0.36$, $\textrm {FP}\%=13$ and $\textrm {FN}\%=0$ for the energy-based method). Therefore, the WRA-based method gives improved prediction times in comparison with the WPOD-based method, albeit with slight increases in the false positives and false negatives obtained.

Figure 15.

Predictions of the burst events obtained using the WRA-based method for three separate time series from DNS. The red-shaded regions indicate the identified burst regions, where the predictor $\lambda$ (grey line) goes below the defined threshold value, i.e. $\lambda <0.95\lambda _t$. The vertical red lines mark the onset of these burst regions, and therefore represent the predictions of the burst event obtained from the WRA-based method. Predicted times are compared with the green dashed vertical lines that indicate the predictions obtained from the WPOD-based method.

6.1. A comparison of the three prediction methods

The prediction methods discussed in this manuscript rely on a predictor going below a threshold value. Therefore, the predictions will vary with this choice of the threshold. The thresholds were chosen as $0.95$ times the mean of the predictor $\lambda$ minus the variance, where these statistics were computed specifically for time instances corresponding to the quiet regions (see § 4.5). Figure 16 illustrates how varying this coefficient $0.95$ impacts predictions obtained from the energy-based method in grey, WPOD-based method in green and WRA-based method in red. The comparison between the methods is carried out using the three time averaged quantities defined in § 4.5: (i) average prediction time $\tau$ in figure 16(a), (ii) the percentage of false positives FP% and (iii) the percentage of false negatives FN% in figure 16(b).

Figure 16.

(a) The average prediction times $\tau$ are shown as a function of the threshold of the predictors for the energy-based method (grey), the WPOD-based method (green) and the WRA-based method (red). (b) The percentages of the obtained predictions that are false positives FP% (solid lines) and false negatives FN% (dashed-dot lines) are also shown for the three methods.

From figure 16(a), we see that, below a threshold value of approximately $1$, both the WPOD-based and WRA-based methods give better prediction times than the energy-based method. From figure 16(b), we see that this improvement comes at the cost of a slight increase in the number of false positives. Moving back to figure 16(a), we see that, above this threshold of approximately $1$, the energy-based method starts performing better than the WPOD-based method, but still worse than the WRA-based method. For very high values of the threshold, the energy-based method seems to perform better than both the wavelet-based methods. However, from figure 16(b) we note that this improvement in the performance of the energy-based method comes with a jump in the number of false positives obtained. For the very high thresholds, more than $50\,\%$ of the predictions given by the energy-based method are false, i.e. one in two predictions are false. We can therefore conclude that both the wavelet-based methods give a more robust prediction performance when compared with the energy-based method. Additionally, from figure 16(a), we see that the WRA-based method always outperforms the WPOD-based method. However, this does come with the cost of a slight increase in the false positives. Finally, considering the false negatives obtained, all three methods give a small number of false negatives. Both the wavelet-based techniques have difficulty in identifying short-lived burst events that exist for 1–2 time units. These events contribute to the false negatives obtained.

(Appendix B shows how this prediction performance changes for a different choice of wavelets, Daubechies 2 instead of the Daubechies 1 used here.)

6.2. Outlook

So far, we have looked at two wavelet-based methods to predict intermittent events in a flow: (i) tracking the coherent structures obtained from WPOD and (ii) tracking the forcing structures that generate these coherent structures through the resolvent operator. We illustrated these methods using the 2-D Kolmogorov flow, which is governed by the unstable eigenfunction of the linearised Navier–Stokes equations. The quiet region of this flow is dominated by this unstable eigenfunction, and the burst events seem to happen because of a disruption of this structure. In other words, there is a single eigenfunction of the linearised Navier–Stokes equations that seems to govern the flow. The energy amplification mechanism in this flow therefore is likely a normal-mode amplification mechanism. In this scenario, the resolvent analysis should not actually be expected to give a significant time delay between the forcing and the response.

Contrary to this, we could analyse a flow where the energy amplification occurs due to non-normal mechanisms, such as for instance the streak generation in channel flows. In this case, the forcing will have a time delay from the response. This is due to the transient growth mechanisms that cause energy amplification in non-normal flows (e.g. Trefethen et al. Reference Trefethen, Trefethen, Reddy and Driscoll1993; Schmid Reference Schmid2007). In this case, the WRA-based method will likely provide better prediction times in comparison with the WPOD-based method. Comparing these methods for non-normal flows is therefore an important future direction of work.

Applying these methods in practice in experiments would first involve finding the dominant energetic structures or forcing structures. For instance, a calibration experiment with good enough resolution to capture these structures or an large eddy simulation (LES) simulation that captures the relevant dynamics can aid in this purpose. The obtained structures can then be tracked in an experiment at a lower resolution. How sparse these measurements can be, and what spatial regions of the flow the measurements need to be concentrated in, are important questions to be addressed in future work. These factors will also closely depend on the flow considered. For instance, for the flow considered here, although $256$ wavenumbers are present in the $x$-direction, only $7$ ($k_x=[-3,-2,-1,0,1,2,3]$) are required for the POD-based prediction.