# Wavelet-based resolvent analysis for statistically-stationary and
temporally-evolving flows

Eric Ballouz 111Ph.D. Student, Mechanical and Civil Engineering, AIAA Student
Member California Institute of Technology, Pasadena, CA, 91125  Barbara Lopez-
Doriga 222Ph.D. Student, Mechanical, Materials, and Aerospace Engineering
Department, AIAA Student Member and Scott T. M. Dawson 333Assistant Professor,
Mechanical, Materials, and Aerospace Engineering Department, AIAA Senior
Member Illinois Institute of Technology, Chicago, IL, 60616  H. Jane Bae
444Assistant Professor, Graduate Aerospace Laboratories, AIAA Senior Member
California Institute of Technology, Pasadena, CA, 91125

###### Abstract

This work introduces a formulation of resolvent analysis that uses wavelet
transforms rather than Fourier transforms in time. This allows resolvent
analysis to be extended to turbulent flows with non-stationary means in
addition to statistically-stationary flows. The optimal resolvent modes for
this formulation correspond to the potentially time-transient structures that
are most amplified by the linearized Navier-Stokes operator. We validate this
methodology for turbulent channel flow and show that the wavelet-based and
Fourier-based resolvent analyses are equivalent for statistically-stationary
flows. We then apply the wavelet-based resolvent analysis to study the
transient growth mechanism in the buffer layer of a turbulent channel flow by
windowing the resolvent operator in time and frequency. The method is also
applied to temporally-evolving parallel shear flows such as an oscillating
boundary layer and three-dimensional channel flow, in which a lateral pressure
gradient perturbs a fully-developed turbulent flow in a channel.

##  1 Nomenclature

$x_{1}$ = | streamwise direction
---|---
$x_{2}$ = | wall-normal direction
$x_{3}$ = | spanwise direction
$\bar{u}_{i}$ = | total velocity in the $x_{i}$ direction
$U_{i}$ = | average velocity over ensembles and homogeneous directions in the $x_{i}$ direction
$u_{i}$ = | fluctuating velocity in the $x_{i}$ direction
$\bar{p}$ = | total pressure
$P$ = | average pressure over ensembles and homogeneous directions
$p$ = | fluctuating pressure
$\nu$ = | kinematic viscosity
$\rho$ = | density
$\langle\cdot\rangle$ = | average over ensembles and homogeneous directions
$F$, $(\hat{\cdot})$ = | Fourier transform in homogeneous directions and time
$W$, $(\tilde{\cdot})$ = | wavelet transform in time and Fourier transform in homogeneous directions
$L$ = | discretized Laplacian
$D_{t}$ = | discretized temporal derivative
$D_{i}$ = | discretized spatial derivative in the $x_{i}$ direction
$T$ = | time horizon
$N_{x_{i}}$ = | number of spatial points in the $x_{i}$ direction
$\Delta x_{i}$ = | grid size in the $x_{i}$ direction
$N_{t}$ = | number of discrete temporal points
$k_{i}$ = | wavenumber in the $x_{i}$ direction
$\omega$ = | temporal frequency
$\varphi$ = | wavelet function
$\alpha$ = | wavelet scale parameter
$\beta$ = | wavelet shift parameter
$\hat{\mathcal{H}},\tilde{\mathcal{H}}$ = | resolvent operator (Fourier- and wavelet-based resolvent, respectively)
$\hat{\phi},\tilde{\phi}$ = | principal resolvent forcing mode (for the Fourier- and wavelet-based resolvent, respectively)
$\hat{\psi},\tilde{\psi}$ = | principal resolvent response mode (for the Fourier- and wavelet-based resolvent, respectively)
$\psi_{i}$ = |  $i$-th component of the principal resolvent response mode
$\sigma$ = | principal resolvent singular value
$\mathcal{B}$ = | windowing matrix that restrict the forcing to a subset of the full space
$\mathcal{C}$ = | windowing matrix that restrict the response to a subset of the full space
$(\cdot)^{\dagger}$ = | Moore-Penrose pseudo-inverse
$\delta$ = | channel half-height
$Re_{\tau}$ = | friction Reynolds number
$u_{\tau}$ = | friction velocity
$(\cdot)^{+}$ = | viscous units
$\Omega$ = | channel wall oscillation frequency for the Stokes boundary layer
$Re_{\Omega}$ = | Reynolds number for the Stokes boundary layer
$\delta_{\Omega}$ = | Stokes laminar boundary layer thickness
$U_{\max}$ = | maximum streamwise wall velocity for the Stokes boundary layer
$U_{i,rms}$ = | root-mean-square velocities profiles in the $x_{i}$ direction
$\Pi$ = | ratio of the mean spanwise and streamwise pressure gradient
$u_{\tau,0}$ = | friction velocity at time $t=0$
$u_{\tau,T}$ = | friction velocity at time $t=T$
$\tau_{i}$ = | wall-shear stress in the $x_{i}$ direction
$\gamma$ = | wall-shear angle $\tan^{-1}(\tau_{3}/\tau_{1})$

##  2 Introduction

Fundamental studies of unsteady aerodynamics and turbulent flow have mostly
focused on statistically stationary configurations, where all statistics are
invariant under a shift in time. However, in real-world applications of
external aerodynamics, truly unsteady transient effects become important. Such
events include separation leading to stall, a sudden change in yaw angle, gust
encounters, and shockwave formation. The additional complexity of these
transient problems makes it harder to perform controlled experiments, both
numerical and in the laboratory. Our goal is to develop reduced-order models
that can be used to understand, predict, and control highly-unsteady transient
turbulent flow in various engineering systems.

Resolvent analysis has been a popular reduced-order model for understanding a
wide variety of turbulent flows. Resolvent analysis refers to the inspection
of the resolvent operator, a linear operator that consists of the linearized
Navier-Stokes and maps forcing inputs (which can be due to the nonlinear
advection terms, or other exogenous inputs) onto the flow states. This
operator governs how inputs are amplified by the linear dynamics of the
system. Its singular value decomposition (SVD) identifies the inputs to which
the linearized equations of motion are most receptive, their gains, and the
most amplified outputs [1, 2]. Often, the truncated SVD is enough to capture
most of the effect of the operator due to the fast decay of the singular
values. The resulting low-rank approximation of the forcing-response dynamics
of the full system is extremely valuable for modeling, controlling, and
understanding the physics of fluid flows [3, 4, 5, 6, 7, 8, 9], and has been
extended to compressible and stratified flows [10, 11, 12]. However, a Fourier
transform in time is traditionally used to formulate the resolvent operator,
which restricts its analysis and application to statistically-steady and
quasi-periodic flows [13]. Indeed, the resulting SVD modes will be Fourier
modes in time, and cannot represent temporally local effects.

To construct the resolvent operator, we instead propose using a wavelet
transform [14] in time. The most effective inputs and the most amplified
outputs given by the SVD of the newly-formulated resolvent operator will be
spatiotemporally localized. There have been numerous research efforts applying
wavelet methods in fluid mechanics, with early work predating the application
of resolvent analysis to turbulent flows and focusing on spatial
transformations [15, 16]. More recent work has focused on data-driven wavelet
transforms [17, 18]. The proposed research effort will provide an important
bridge connecting such data-driven analyses to the underlying physics,
analogous to recent work connecting and comparing resolvent analysis to data-
driven spectral proper orthogonal decomposition [19, 20, 21].

In the present work, we develop and validate wavelet-based resolvent analysis
for a variety of systems, ranging from quasi-parallel wall-bounded turbulent
flows to spatio-temporally evolving systems. In §3, we formulate the Fourier-
and wavelet-based resolvents; we highlight the properties of the wavelet
transform and justify the choice of wavelet basis. In §4 we first establish
the equivalence of Fourier- and wavelet-based resolvent analyses for the
statistically stationary turbulent channel flow and showcase the additional
capacity of the wavelet-based resolvent to capture the linear transient growth
of streaks under transient forcing. We also apply wavelet-based resolvent
analysis to statistically non-stationary flows in §5, notably the Stokes
boundary layer flow as well as a turbulent channel flow subjected to a sudden
lateral pressure gradient. Finally, the summary of the work and a discussion
of the results are given in §6.

##  3 Mathematical formulation

###  3.1 Fourier-based resolvent analysis

The incompressible Navier-Stokes equations in Einstein notation are given by

$\frac{\partial\bar{u}_{i}}{\partial{t}}+\bar{u}_{j}\frac{\partial\bar{u}_{i}}{\partial{x_{j}}}=-\frac{1}{\rho}\frac{\partial\bar{p}}{\partial
x_{i}}+\nu\frac{\partial^{2}\bar{u}_{i}}{\partial{x}_{j}\partial{x}_{j}},\quad\frac{\partial\bar{u}_{i}}{\partial{x_{i}}}=0,$
(1)

where $\bar{u}_{i}$ is the total velocity (including the mean and the
fluctuating component), $\tilde{p}$ is the total pressure, $\rho$ is the
density, and $\nu$ is the kinematic viscosity. The total velocity can be
divided into $\bar{u}_{i}=U_{i}+u_{i}$, where $U_{i}$ is the average over
ensembles and homogeneous directions, and $u_{i}$ is the fluctuating
component. Similarly, pressure can be decomposed as $\tilde{p}=P+p$.

Equation (1) can then be divided into the mean and the fluctuating equations

$\displaystyle\frac{\partial{U_{i}}}{\partial
t}+\left\langle\bar{u}_{j}\frac{\partial\bar{u}_{i}}{\partial{x_{j}}}\right\rangle=-\frac{1}{\rho}\frac{\partial{P}}{\partial{x_{i}}}+\nu\frac{\partial^{2}U_{i}}{\partial{x_{j}}\partial{x_{j}}},\quad\frac{\partial{U_{i}}}{\partial{x_{i}}}=0,$
(2) $\displaystyle\frac{\partial{u_{i}}}{\partial
t}+U_{j}\frac{\partial{u_{i}}}{\partial{x_{j}}}+u_{j}\frac{\partial{U_{i}}}{\partial{x_{j}}}=-\frac{1}{\rho}\frac{\partial{p}}{\partial{x_{i}}}+\nu\frac{\partial^{2}u_{i}}{\partial{x_{j}}\partial{x_{j}}}+f_{i},\quad\frac{\partial{u_{i}}}{\partial{x_{i}}}=0,$
(3)

where $f_{i}$ is the remaining nonlinear terms in the fluctuating equations
and $\langle\cdot\rangle$ is the average over ensembles and in homogeneous
directions. Note that some of the terms in the fluctuating equations may be
zero depending on the flow configuration. The equations above do not have an
analytic solution unless in very particular situations and are most commonly
solved numerically. Discretizing the fluctuating equations, we get

$D_{t}{u_{i}}+U_{j}D_{j}{u_{i}}+u_{j}dU_{i,j}=-\frac{1}{\rho}D_{i}{p}+\nu
L{u_{i}}+f_{i},\quad D_{i}u_{i}=0,$ (4)

where $D_{t}$ is the discrete derivative in time, $D_{i}$ is the discrete
derivative in the $x_{i}$ direction, $L$ is the discrete Laplacian, and
$dU_{i,j}$ denotes the matrix form of $\partial U_{i}/\partial x_{j}$. Each
discretized equation is an $N_{t}\times N_{x_{1}}\times N_{x_{2}}\times
N_{x_{3}}$-dimensional system, where $N_{t}$ is the temporal resolution, and
$N_{x_{i}}$ are the spatial resolutions in the $x_{i}$, $i=1,2,3$ directions
respectively. The discretized velocity and velocity gradient, $U_{j}$ and
$dU_{i,j}$, are $(N_{t}\times N_{x_{1}}\times N_{x_{2}}\times N_{x_{3}})^{2}$
diagonal matrices. In the traditional resolvent analysis, we apply the Fourier
transform operator in homogeneous directions and time, $F$, to the left in Eq.
(4), leading to the equation

$(FD_{t}F^{-1})(F{u_{i}})+(FU_{j}D_{j}F^{-1})(F{u_{i}})+(Fu_{j})(FdU_{i,j}F^{-1})=\\\
-\frac{1}{\rho}(FD_{i}F^{-1})(F{p})+\nu(FLF^{-1})(F{u_{i}})+Ff_{i},\quad(FD_{i}F^{-1})(Fu_{i})=0,$
(5)

where $F^{-1}$ is the inverse transformation, or equivalently

$\widehat{D_{t}}\hat{u}_{i}+\widehat{U_{j}D_{j}}\hat{u}_{i}+\widehat{dU_{i,j}}\hat{u}_{j}=-\frac{1}{\rho}\widehat{D_{i}}\hat{p}+\nu\hat{L}\hat{u}_{i}+\hat{f}_{i},\quad\widehat{D_{i}}\hat{u}_{i}=0.$
(6)

Note that for an arbitrary matrix $M$ and vector $b$, $\hat{M}:=FMF^{-1}$ and
$\hat{b}:=Fb$. For temporally stationary systems, this equation can typically
be decoupled for each wavenumber and frequency combination. For example, in
the case of channel flow where the flow is homogeneous in the streamwise
($x_{1}$) and spanwise ($x_{3}$) directions, the linear operator can be cast
as

$\left[\begin{array}[]{c}\hat{u}_{1}(k_{1},k_{3},\omega)\\\
\hat{u}_{2}(k_{1},k_{3},\omega)\\\ \hat{u}_{3}(k_{1},k_{3},\omega)\\\
\hat{p}(k_{1},k_{3},\omega)\end{array}\right]=\hat{\mathcal{H}}(k_{1},k_{3},\omega)\left[\begin{array}[]{c}\hat{f}_{1}(k_{1},k_{3},\omega)\\\
\hat{f}_{2}(k_{1},k_{3},\omega)\\\ \hat{f}_{3}(k_{1},k_{3},\omega)\\\
0\end{array}\right]$ (7)

for a given $(k_{1},k_{3},\omega)$ triplet. Typically, the singular value
decomposition of the linear operator
$\hat{\mathcal{H}}\in\mathbb{C}^{4N_{x_{2}}}\times\mathbb{C}^{4N_{x_{2}}}$ is
taken to study the left and right singular vectors as response and forcing
modes, and the singular values as amplification factors or gains.

###  3.2 Wavelet-based resolvent analysis

####  3.2.1 Formulation

To account for transient behavior in the mean flow or the fluctuations, we
introduce the wavelet-based resolvent analysis. The benefit of the wavelet
transform in time is that it preserves both time and frequency information.
The wavelet transform projects a function onto a wavelet basis composed of
scaled and shifted versions of a mother function $\varphi(t)$. The transformed
function depends on the scale ($\alpha$) and shift ($\beta$) parameters
respectively linked to frequency and time information, whereas the Fourier
transform is a function of only frequency.

We propose using a wavelet transform in time while keeping the Fourier
transform in homogeneous directions. We denote the total transformation
operator (wavelet in time and Fourier in homogeneous directions) as $W$ and
the left inverse operator as $W^{-1}$, which is also the right inverse for
unitary transforms. The inverse operator is well-defined and unique for
orthogonal wavelet bases. We can then apply $W$ on the left in Eq. (4), which
gives

$(WD_{t}W^{-1})(W{u_{i}})+(WU_{j}D_{j}W^{-1})(W{u_{i}})+(WdU_{i,j}W^{-1})(Wu_{j})=\\\
-\frac{1}{\rho}(WD_{i}W^{-1})(W{p})+\nu(WLW^{-1})(W{u_{i}})+Wf_{i},\quad(WD_{i}W^{-1})(Wu_{i})=0,$
(8)

or

$\widetilde{D_{t}}\tilde{u}_{i}+\widetilde{U_{j}D_{j}}\tilde{u}_{i}+\widetilde{dU_{i,j}}\tilde{u}_{j}=-\frac{1}{\rho}\widetilde{D_{i}}\tilde{p}+\nu\tilde{L}\tilde{u}_{i}+\tilde{f}_{i},\quad\widetilde{D_{i}}\tilde{u}_{i}=0.$
(9)

Note that for an arbitrary matrix $M$ and vector $b$, $\tilde{M}:=WMW^{-1}$
and $\tilde{b}:=Wb$. These equations can be separated for each spatial
wavenumber in the homogeneous direction, and thus the dimension of each linear
equation is smaller than the full Navier-Stokes equations. Note that if the
transformation in time is given by the Fourier transform rather than a wavelet
transform, this would recover the traditional Fourier-based resolvent analysis
[2] (if the flow is temporally stationary) or the harmonic resolvent [13]
analysis (if the flow is periodic in time). Similar to the Fourier-based
resolvent analysis, for flows that are homogeneous in the $x_{1}$ and $x_{3}$
directions, this can be written in matrix form as

$\left[\begin{array}[]{c}\tilde{u}_{1}(k_{1},k_{3})\\\
\tilde{u}_{2}(k_{1},k_{3})\\\ \tilde{u}_{3}(k_{1},k_{3})\\\
\tilde{p}(k_{1},k_{3})\end{array}\right]=\tilde{\mathcal{H}}(k_{1},k_{3})\left[\begin{array}[]{c}\tilde{f}_{1}(k_{1},k_{3})\\\
\tilde{f}_{2}(k_{1},k_{3})\\\ \tilde{f}_{3}(k_{1},k_{3})\\\
0\end{array}\right],$ (10)

where the wavelet-based resolvent operator
$\tilde{\mathcal{H}}\in\mathbb{C}^{4N_{t}\times
N_{x_{2}}}\times\mathbb{C}^{4N_{t}\times N_{x_{2}}}$ is defined as

$\tilde{\mathcal{H}}=\left[\left(\tilde{D}_{t}-\nu\tilde{L}+\widetilde{U_{i}D_{i}}\right)\left(\begin{array}[]{cccc}1&0&0&0\\\
0&1&0&0\\\ 0&0&1&0\\\
0&0&0&0\end{array}\right)\right.+\left.\left(\begin{array}[]{cccc}\widetilde{dU}_{1,1}&\widetilde{dU}_{1,2}&\widetilde{dU}_{1,3}&\frac{1}{\rho}\tilde{D}_{1}\\\
\widetilde{dU}_{2,1}&\widetilde{dU}_{2,2}&\widetilde{dU}_{2,3}&\frac{1}{\rho}\tilde{D}_{2}\\\
\widetilde{dU}_{3,1}&\widetilde{dU}_{3,2}&\widetilde{dU}_{3,3}&\frac{1}{\rho}\tilde{D}_{3}\\\
\tilde{D}_{1}&\tilde{D}_{2}&\tilde{D_{3}}&0\end{array}\right)\right]^{-1}.$
(11)

This formulation allows us to study transient flows using resolvent analysis.

####  3.2.2 Wavelet-based resolvent analysis with windowing

We can reformulate a resolvent map between forcing and response at specific
time shifts and scales by defining a windowed resolvent operator

$\left[\begin{array}[]{c}\tilde{u}_{1}(k_{1},k_{3})\\\
\tilde{u}_{2}(k_{1},k_{3})\\\ \tilde{u}_{3}(k_{1},k_{3})\\\
\tilde{p}(k_{1},k_{3})\end{array}\right]=\mathcal{C}\tilde{\mathcal{H}}(k_{1},k_{3})\mathcal{B}\left[\begin{array}[]{c}\tilde{f}_{1}(k_{1},k_{3})\\\
\tilde{f}_{2}(k_{1},k_{3})\\\ \tilde{f}_{3}(k_{1},k_{3})\\\
0\end{array}\right],$ (12)

where $\mathcal{B}$ and $\mathcal{C}$ are windowing matrices on the forcing
and response modes, respectively [22, 23]. The windowing matrices select a
subset of the full forcing and response states. For example, to select a
particular scale and shift parameter $(\alpha_{s},\beta_{s})$ for the forcing
mode, we set

$\mathcal{B}=\text{diag}\big{(}\mathbbm{1}(\alpha=\alpha_{s})\mathbbm{1}(\beta=\beta_{s})\big{)},$
(13)

where $\mathbbm{1}(\cdot)$ is an indicator function. The SVD of the windowed
resolvent operator, $\mathcal{C}\tilde{\mathcal{H}}(k_{1},k_{3})\mathcal{B}$,
allows us to identify forcing and response modes restricted to a limited
frequency and time interval.

The windowed resolvent operator, however, is equivalent to taking the Moore-
Penrose pseudo-inverse of
$\mathcal{B}^{\dagger}\tilde{\mathcal{H}}(k_{1},k_{3})^{-1}\mathcal{C}^{\dagger}$,
where the superscript $\dagger$ indicates the pseudo-inverse. In this case,
$\mathcal{B}^{\dagger}=\mathcal{B}$ and $\mathcal{C}^{\dagger}=\mathcal{C}$.
By applying the pseudo-inverse to the linearized Navier-Stokes operator prior
to inverting it to compute $\tilde{\mathcal{H}}(k_{1},k_{3})$, we can reduce
the computational cost. The matrix inversion and SVD will apply to a matrix
with its size defined by the nonzero block of $\mathcal{B}\mathcal{C}$ rather
than the full system.

####  3.2.3 Choice of wavelet basis

(a)

(b)

Figure 1: Shannon wavelet in (a) time and (b) frequency domains.

Wavelet transforms are not unique and are defined by the choice of the mother
wavelet $\varphi(t)$. Consider an arbitrary function $f$ in
$L_{2}(\mathbb{R})$. Its Fourier and wavelet transforms are

$\hat{f}(\omega)=\int_{-\infty}^{+\infty}f(t)e^{-\mathrm{i}\omega t}dt,$ (14)
$\tilde{f}(\alpha,\beta)=\int_{-\infty}^{+\infty}f(t)\varphi^{*}_{\alpha,\beta}(t)dt,\quad\varphi_{\alpha,\beta}(t)\equiv\frac{1}{\sqrt{\alpha}}\varphi\left(\frac{t-\beta}{\alpha}\right).$
(15)

Note that there is a trade-off between precision in frequency and precision in
time, i.e., one cannot find a function $\varphi(t)$ that is well localized in
both time and frequency [24]. As two extreme examples, consider the Dirac
delta centered at $t=1$, which is perfectly localized in space but with an
infinite spread in frequency space, and the Fourier mode $e^{\mathrm{i}t}$,
$\mathrm{i}=\sqrt{-1}$, which is perfectly localized in frequency space at
$\omega=1$ but has infinite spread in time.

For this study, we choose the Shannon wavelet for its frequency properties.
Each Shannon wavelet acts as a perfect band-pass filter and covers a band
$2^{-i}/T([-2\pi,-\pi]\cup[\pi,2\pi])$, $i\in\mathbb{Z}$ (Fig. 1), where $T$
denotes the time horizon considered. Though the Shannon wavelet does not have
perfect frequency localization provided by the Fourier transform, it allows
the separation of the frequency content into distinct non-overlapping bands
for different scales. In addition, the discrete transform matrix $W$ [25, 24]
for the Shannon wavelet is an orthonormal basis for $L_{2}(\mathbb{R})$,
making $W$ unitary [25]. However, we note that the results in the current
study do not change significantly with different choices of the wavelet basis.

####  3.2.4 Computational cost

The construction of $\tilde{\mathcal{H}}$ requires the inversion of a
$4N_{y}N_{t}\times 4N_{y}N_{t}$ matrix, a computation that costs
$O(64N_{y}^{3}N_{t}^{3})$ operations. The full SVD of $\tilde{\mathcal{H}}$
would also require $O(64N_{y}^{3}N_{t}^{3})$ operations. With a direct solve,
the wavelet-based resolvent analysis would cost $O(N_{t}^{2})$-times more than
performing $N_{t}$ separate Fourier-based resolvent for each temporal scale,
though the latter would fail to capture the interactions between the different
time scales. This penalty of $O(N_{t}^{2})$ is the nominal cost of
constructing time-localized resolvent modes.

One method for reducing the memory and computational cost of wavelet-based
resolvent analysis is to use sparse finite difference operators and wavelet
transforms when constructing $\tilde{\mathcal{H}}^{-1}$. The resulting sparse
matrix can be factored with specialized packages like PARDISO [26, 27, 28],
which we utilize in the current study, in order to efficiently solve linear
equations of the form $\tilde{\mathcal{H}}^{-1}v=w$, where $v$ and $w$ are
arbitrary vectors. To further take advantage of the sparsity of
$\tilde{\mathcal{H}}$, we opt for an iterative method to perform the SVD. In
this work, we use a one-sided Lanczos bidiagonalization [29], an iterative
algorithm that allows us to compute a truncated SVD and accurately estimate a
number $q<4N_{y}N_{t}$ of the most significant singular input and output
modes. Other efficient SVD algorithms rely on randomized approaches, in
particular by sub-sampling the high-dimensional matrix and performing the SVD
on the lower-dimensional approximation [30, 31, 32]. A randomized SVD of a
high-dimensional discrete resolvent operator was used in [33, 34].

In some cases such as §3.2.2, where only some wavelet scales or time shifts
are relevant, the windowing matrices $\mathcal{C}$ and $\mathcal{B}$ can be
chosen to select the significant regions of the time-frequency domain. We then
neglect the zero rows and columns in
$\mathcal{C}\tilde{\mathcal{H}}\mathcal{B}$ associated with non-significant
wavelet scales and shifts and perform the SVD and pseudo-inversion on the
lower-dimensional system.

##  4 Application to statistically-stationary flow

We first validate the wavelet-based resolvent analysis on a statistically-
stationary turbulent channel flow. We have shown that for unitary wavelet
transforms the Fourier-based and wavelet-based resolvent modes are equivalent.
Thus, for the channel flow case, we expect the two methods to produce
identical resolvent modes. We then utilize the temporally-local property of
wavelet-based resolvent analysis to study transient growth in a turbulent
channel flow.

###  4.1 Turbulent channel flow

(a)

(b)

Figure 2: (a) First ten singular values for the wavelet-based resolvent (red)
and the largest singular value for the Fourier-based resolvent (black)
computed for each $\omega_{i}$. The vertical gray lines indicate the temporal
resolution in frequency space, and the vertical black lines delimit the
frequency band covered by each of the chosen wavelet scales. (b) Magnitude of
the Fourier-transformed wavelet-based (red) and Fourier-based (black)
streamwise resolvent response mode for $\omega=17.14u_{\tau}/\delta$.

We first compare the results of both the traditional Fourier-based resolvent
analysis and wavelet-based resolvent analysis. For this, we consider a
turbulent channel flow. The mean profile of a turbulent channel flow at
friction Reynolds number $Re_{\tau}=186$ is obtained by a direct numerical
simulation (DNS) using a second-order staggered finite-difference [35] and a
fractional-step method [36] with a third-order Runge-Kutta time-advancing
scheme [37]. Periodic boundary conditions are imposed in the streamwise and
spanwise directions and the no-slip and no-penetration boundary conditions are
used at the top and bottom walls. The code has been validated in previous
studies in turbulent channel flows [38, 39, 40] and flat-plate boundary layers
[41]. The numerical domain is $8\pi\delta\times 2\delta\times 3\pi\delta$ and
is discretized using 768, 130, and 288 grid points in the streamwise, wall-
normal, and spanwise directions, respectively. We use uniform grid spacing in
the streamwise and spanwise directions of $\Delta x_{1}^{+}=6$ and $\Delta
x_{3}^{+}=6$, and a wall-normal grid stretched away from the wall using a
hyperbolic tangent with $\min(\Delta x_{2}^{+})=0.16$ and $\max(\Delta
x_{2}^{+})=7.2$. Here the superscript $+$ denotes wall units defined in terms
of $\nu$ and $u_{\tau}$, where $u_{\tau}$ is the friction velocity. The
simulations were run for 100 eddy turnover times (defined as
$\delta/u_{\tau}$) after transients to compute the mean quantities.

For the resolvent analysis, the wall-normal direction is discretized using a
Chebyshev collocation method using $N_{x_{2}}=128$, and the mean streamwise
velocity profile and its wall-normal derivative from the DNS are interpolated
to the Chebyshev collocation points. We use a periodic boundary condition for
the temporal domain, $T$, with a temporal resolution of $N_{t}=128$. We choose
$D_{t}$ to be a Fourier differentiation matrix, and $D_{x_{2}}$ a Chebyshev
differentiation matrix with a no-slip and no-penetration boundary condition at
the wall. Note that these matrices are not sparse; for this case, sparse
differentiation matrices were not needed, though finite difference matrices
may be used in higher-dimensional problems to improve efficiency. We choose
the spanwise and streamwise wavelengths of $\lambda_{1}^{+}\approx 1000$ and
$\lambda_{3}^{+}\approx 100$, in line with the most energetic structures close
to the wall.

Since the mean profiles are statistically steady, we have
$\widetilde{U_{j}D_{j}}=\widehat{U_{j}D_{j}}$,
$\widetilde{dU_{i,j}}=\widehat{dU_{i,j}}$, $\tilde{L}=\hat{L}$. The wavelet-
and Fourier-based cases thus only differ by the time differentiation matrix
such that $W^{-1}\widetilde{D_{t}}W=F^{-1}\hat{D}_{t}F=D_{t}$. Note that $W$
is unitary for the choice of wavelet. Since the singular value decomposition
is unique up to multiplication by a unitary matrix, we expect the singular
values of $\tilde{\mathcal{H}}$ to be the same as that of

$\hat{\mathcal{H}}(k_{1},k_{3})=\left(\begin{array}[]{cccc}\hat{\mathcal{H}}(k_{1},k_{3},\omega_{1})&&&\\\
&\hat{\mathcal{H}}(k_{1},k_{3},\omega_{2})&&\\\
&&\hat{\mathcal{H}}(k_{1},k_{3},\omega_{3})&\\\ &&&\ddots\\\
\end{array}\right),$ (16)

where $\omega_{i}=(2\pi i)/T$ for $i=-{N_{t}}/{2},\cdots,{N_{t}}/{2}-1$. Here,
we set $T=5.5\delta/u_{\tau}$. Moreover, we expect the response and forcing
modes of both systems to be related by the unitary transform given by the
Fourier and inverse-wavelet transform in time, $FW^{-1}$. In Fig. 2(a), we
show the singular values of the leading Fourier- and wavelet-based resolvent
response modes, $\hat{\psi}$ and $\tilde{\psi}$ respectively. The Fourier-
based resolvent modes were computed by applying the resolvent analysis at each
$\omega_{i}$ and taking the principal singular value. The wavelet-based
resolvent modes were computed by solving the full space-time system at once.
We also plot the most amplified streamwise resolvent response mode for the two
methods in Fig. 2(b). As expected, the singular values and the corresponding
modes are equivalent. The small deviations in the singular values for the two
cases are due to the truncation errors of the Shannon wavelet basis.

###  4.2 Transient growth mechanism of turbulent channel flow

Figure 3: Mean streamwise velocity profile for channel flow at
$Re_{\tau}=186$. The buffer layer is delimited by the dashed lines; the black
horizontal lines correspond to frequency bands covered by the chosen wavelets,
mapped to streamwise velocities using $U_{1}=\omega/k_{1}$.

(a)

(b)

Figure 4: (a) Magnitude of the principal response mode of
$\tilde{\mathcal{H}}$ for channel flow at $Re_{\tau}=186$ under transient
forcing, in the frequency-time plane. The forcing window is highlighted in
red. (b) Magnitude of the inverse wavelet-transformed principal response mode
as a function of time.

The added advantage of the wavelet-based method lies in its ability to
preserve temporal localization, which lets us formulate a time-scale sensitive
resolvent for transient flows even when the mean profiles are statistically
stationary. For this, we use the windowed wavelet-based resolvent analysis
framework from §3.2.2, but setting the windowing of the response modes,
$\mathcal{C}$, to the identity matrix. We then calculate the resolvent modes
of $\tilde{\mathcal{H}}\mathcal{B}$. The principal forcing mode will lie in
the span of a subset of the complete wavelet basis, as selected by
$\mathcal{B}$. By picking a particular wavelet scale and shift via
$\mathcal{B}$, we can restrict the forcing to a chosen frequency range, mostly
localized in a particular time range.

To demonstrate this capability, we study the effect of a time-localized
forcing on the buffer layer ($x_{2}^{+}\in[15,40]$). We choose $\mathcal{B}$
such that the forcing modes are limited to a single wavelet with the frequency
range containing the critical layer corresponding to the buffer layer, i.e.,
$[U(x_{2}^{+}=15),U(x_{2}^{+}=40)]=[\omega_{\min},\omega_{\max}]/k_{1}$ (Fig.
3). For the current case, we have $\omega_{\min}=12.52{u_{\tau}}/{\delta}$ and
$\omega_{\max}=17.17{u_{\tau}}/{\delta}$. We pick a temporal resolution of
$N_{t}=64$, and a time horizon of $T=5.5{\delta}/{u_{\tau}}$. We project the
flow variables onto Shannon wavelets that cover frequency intervals
$(2^{-i-1}{N_{t}}/{T})([-2\pi,-\pi]\cup[\pi,2\pi])$, for $i=0,1,\cdots 5$. The
windowing matrix $\mathcal{B}$ is chosen to restrict the forcing to the
unshifted wavelet ($\tau=0$, i.e., centered at $t=0$) that has been scaled to
cover $[9.14,18.28]{u_{\tau}}/{\delta}$ (i.e. $i=1$).

The resulting principal response mode is confined to the frequency band
determined by the forcing, as shown in Fig. 4(a), which is expected since the
time scales are decoupled for statistically stationary flows. The response
mode additionally varies in time; the magnitude of the mode growth as a
function of time is shown in Fig. 4(b). The principal response mode peaks at a
time $tu_{\tau}/\delta=0.79$ before decaying. This transient growth can be
explained by the non-normality of the linearized system [42]. The response
modes at three different time shifts are shown in Fig. 5. The modes are
predominantly in the streamwise direction, forming alternating low- and high-
speed streamwise streaks. Although not shown, the principal forcing mode only
has a small contribution from the streamwise direction and is in the form of
streamwise rolls, in line with the linear theory [43, 44, 45] and the self-
sustaining process of wall turbulence [46, 47, 48, 49, 50, 51, 10].

(a)

(b)

(c) Figure 5: Principal response mode of the wavelet-based resolvent for the
channel flow at $Re_{\tau}=186$ and under transient forcing, shown at (a)
$tu_{\tau}/\delta=0.26$, (b) $tu_{\tau}/\delta=0.79$, and (c)
$tu_{\tau}/\delta=1.92$. The arrows in the right panels are colored according
to the sign of
$\sqrt{|\tilde{\psi_{2}}|^{2}+|\tilde{\psi_{3}}|^{2}}sgn(\tilde{\psi_{1}})$ to
show the intensity of the streamwise rolls relative to the streamwise streaks.

##  5 Application to non-stationary flow

We now apply wavelet-based resolvent analysis to problems with a time-varying
mean flow. In particular, we study the turbulent Stokes boundary layer and a
turbulent channel flow with a sudden lateral pressure gradient. The Stokes
boundary layer is a purely oscillatory flow in time, and thus, Fourier-based
resolvent analysis still may be used [13]. However, in the case of the
temporally-changing channel flow, the flow is truly unsteady, and a Fourier
transform in time will not be applicable.

###  5.1 Turbulent Stokes boundary layer

The Stokes boundary layer is simulated through a channel flow with the lower
and upper walls oscillating in tandem at a velocity of $U_{w}(t)=U_{\max}$
cos$(\Omega t)$ with no imposed pressure gradient. The relevant nondimensional
number is the Reynolds number $Re_{\Omega}=U_{\max}\delta_{\Omega}/\nu$, where
$\delta_{\Omega}=\sqrt{2\nu/\Omega}$ denotes the Stokes boundary layer
thickness. For the current case, we consider $Re_{\Omega}=1500$, which lies
within the intermittently turbulent regime [52, 53, 54, 55, 56]. This problem
has been well-studied numerically and experimentally in the literature [52,
57, 58, 53, 54, 55, 56, 59, 60, 61, 62, 63]. The same numerical solver used to
simulate the turbulent channel flow is used to generate the statistics, with
modifications in the boundary condition to accommodate the oscillating walls.
The domain size of the channel for the DNS is given by
$6\pi\delta_{\Omega}\times 80\delta_{\Omega}\times 3\pi\delta_{\Omega}$ and
discretized using $64$, $385$ and $64$ points in each direction. We compute
the mean velocity profiles by averaging in homogeneous directions and phase.
Fig. 6 shows the mean and the root-mean-square (rms). velocity profiles at
three different temporal locations. Note that $U_{1}(t\Omega+\pi)=-U(t\Omega)$
and $U_{i,rms}(t\Omega+\pi)=U_{i,rms}(t\Omega)$. We observe that the turbulent
energy peaks near the wall at $t\Omega=2.65$ and propagates away from the wall
in time.

(a)

(b)

Figure 6: (a) Mean streamwise velocity profile and (b) streamwise r.m.s.
velocity at from $t\Omega=0$ (blue) to $t\Omega=\pi$ (red). The profiles shown
are at
$t\Omega=0,\,\pi/8,\,\pi/4,\,3\pi/8,\,\pi/2,\,5\pi/8,\,3\pi/4,\,7\pi/8,\pi$.

To construct the resolvent operator, we first choose the spatial scales for
the homogeneous directions. Using the DNS data, we calculate the streamwise
energy spectrum at $x_{2}/\delta_{\Omega}=1.43$ and $t\Omega=2.51$, which
correspond to the peak spatio-temporal location of the $U_{1,rms}$. The most
energetic streamwise and spanwise scales at that location are given by
$k_{1}\delta_{\Omega}=0.67$ and $k_{3}\delta_{\Omega}=2.67$, which we choose
as the streamwise and spanwise scales for the resolvent operator. We use
$N_{x_{2}}=140$ and $N_{t}=780$ to solve the discrete system. We construct
$D_{t}$ and $D_{x_{2}}$ as sparse finite difference matrices. We additionally
choose the time-derivative operator $D_{t}$ to be circulant to enforce
periodicity in time. We compute the modes for the half-channel and modify
$D_{x_{2}}$ to enforce a no-slip and no-penetration boundary condition at the
wall, and a free-slip and no-penetration boundary condition at the centerline.

(a)

(b)

Figure 7: (a) Streamwise velocity component of the principal resolvent input
(top) and output (bottom) modes of the Stokes boundary layer. The black
contour lines are $U_{1,rms}$ with the levels indicating $30\%,50\%,75\%,95\%$
of its maximum value. The vertical dashed lines show the location of the
amplitude peak for the input mode ($t\Omega=2.54$) and output mode
($t\Omega=2.67$). (b) Magnitude of the principal response mode of $\tilde{H}$
for the Stokes boundary layer in the frequency-time plane.

(a)

(b)

Figure 8: Mean streamwise velocity profiles normalized by $u_{\tau}(t\Omega)$
at (a) $t\Omega=2.35$ and (b) $t\Omega=2.67$, which correspond to 10% and 100%
of the maximum intensity of the principal response mode. The vertical red
lines show the wall-normal location of the instantaneous peak of the
streamwise component of the output mode, the dashed black lines show the
viscous and logarithmic layers, and the horizontal blue line marks where the
velocity is zero.

The principal input and output modes corresponding to the chosen spatial
scales and boundary conditions in time are shown in Fig. 7(a). We observe that
the principal input and output modes are synchronized with the peaks in
$U_{1,rms}$. The energy propagation towards the center occurs at a similar
rate for the resolvent modes as for the DNS results. This suggests that the
energy amplification in the Stokes boundary layer can partially be explained
by the optimal linear mechanism, as for turbulent channel flow [45]. Moreover,
we observe that the principal input modes precede the principal output modes
with time delay, and thus a transient growth mechanism can be explained
through the wavelet-based resolvent modes. This is in line with a physical
interpretation of the modes in which the input modes ‘cause’ the output modes.
In Fig. 7(b), we show the magnitude of the response modes in the frequency-
time plane, which highlights that the bulk of the response is limited to a
single wavelet scale and a few time shifts. Using the windowed resolvent
operator to restrict our forcing and our response to the first wavelet scale,
we indeed obtain nearly identical modes as in 7 while reducing the effective
dimension of our resolvent operator by a factor of four.

Fig. 8 shows that the peak amplification occurs at the location of zero mean
streamwise velocity. This is in line with the fact that Fourier-based
resolvent modes are often centered around a critical layer [2], where the
critical layer $x_{2}=x_{2}^{c}$ is defined as $U_{1}(x_{2}^{c})=c$ and
$k_{x}c=\omega$. For the Stokes boundary layer with the chosen length scales,
the critical layer seems to occur at $U_{1}(x_{2}^{c})=0$, which corresponds
to $\omega=0$. For this time-periodic problem, we note that wavelet resolvent
analysis should be equivalent to a harmonic resolvent analysis [13] that
includes the interactions between all the resolved time scales. The wavelet
resolvent modes should indeed map to the harmonic resolvent modes via an
inverse wavelet transform and a Fourier transform in time. We expect that
harmonic resolvent analysis will also reveal a peak amplification for waves at
$\omega=0$; we will conduct a harmonic resolvent analysis on this system and
compare it with the wavelet-based resolvent modes in future works.

###  5.2 Channel flow with sudden lateral pressure gradient

(a)

(b)

Figure 9: Mean (a) streamwise and (b) spanwise velocity profile from
$tu_{\tau,0}/\delta=0$ (blue) to $tu_{\tau,0}/\delta=2.34$ (red). The times
shown are $tu_{\tau,0}/\delta=0,\,0.58,\,1.17,\,1.76,\,2.34$. Data taken from
[64].

(a)

(b)

Figure 10: (a) Friction velocity $u_{\tau}$ and (b) wall-shear stress angle
$\gamma=\tan^{-1}(\tau_{3}/\tau_{1})$ as a function of time. Data taken from
[64]. The vertical dashed lines are at $t\Omega=0,1.38,1.94$, and correspond
to the choice of $\lambda_{1}^{+}$ and $\lambda_{3}^{+}$ for the modes plotted
in Fig.11

Finally, we study a fully-developed turbulent channel flow at $Re_{\tau}=186$
that is subjected to a sudden lateral pressure gradient $dP/dx_{3}=\Pi
dP/dx_{1}$ at $t=0$ with $\Pi=30$ [65, 64]. This flow, commonly referred to as
a three-dimensional (3D) channel flow, has an initial transient period
dominated by 3D non-equilibrium effects. Eventually, the flow will reach a new
statistically steady state with the mean flow in the $(dP/dx_{1},dP/dx_{3})$
direction parallel to the wall. In the transient period, the tangential
Reynolds stress initially decreases before increasing linearly, with depletion
and increase rate that scales as $\Pi x_{2}/\delta$ [64].

The mean flow profiles are obtained from [64] and have nonzero streamwise and
spanwise components $U_{1}$ and $U_{3}$ (Fig. 9) as well as nonzero wall-
normal gradients of streamwise and spanwise components $dU_{1,2}$ and
$dU_{3,2}$. The time domain of the simulation is $Tu_{\tau,0}/\delta=2.34$,
where $u_{\tau,0}$ is the initial friction velocity. To construct the discrete
resolvent operator, we use a spatial resolution of $N_{x_{2}}=65$ for the
half-channel, and a finite difference matrix $D_{x_{2}}$ enforcing a no-slip
and no-penetration boundary condition at the wall and a free-slip and no-
penetration condition at the centerline. The boundary condition for the
temporal finite difference operator $D_{t}$ is chosen to enforce a Neumann-
type condition, $\partial_{t}(\cdot)|_{t=0}=\partial_{t}(\cdot)|_{t=T}=0$. To
reduce the impact of the boundary condition on the modes at $t=0$ we extend
$U_{1}$ and $U_{3}$ to the time interval $tu_{\tau,0}/\delta\in[-0.58,2.34]$
and assume $U_{i}(t\leq 0,y)=U_{i}(t=0,y)$ and $dP/dx_{3}(t<0)=0$. When the
modes are plotted, we only show the original time domain
$tu_{\tau,0}/\delta\in[0,2.34]$ and exclude the contribution from negative
times. We use a temporal resolution of $N_{t}=1000$ for the extended time
frame.

Regarding the spatial scales for the homogeneous directions, we choose
$(\lambda_{1}^{+},\,\lambda_{3}^{+})=(1000,\,100)$ as well as
($\lambda_{1}^{+},\,\lambda_{3}^{+})=(189,\,1890)$ and
$(\lambda_{1}^{+},\,\lambda_{3}^{+})=(297,\,2970)$. Here, $(\cdot)^{+}$
indicates the wall scaling with respect to $u_{\tau,0}$, before the lateral
pressure gradient is applied. The first combination of $\lambda_{1}^{+}$ and
$\lambda_{3}^{+}$ is the same as in §4.1 and corresponds to the spatial scales
preferred by the near-wall streaks at $Re_{\tau}=186$ prior to the lateral
pressure gradient. The resolvent modes for these scales are shown in Fig.
11(a,b). The magnitude of the modes in frequency-time space is also plotted in
Fig. 12(a). The resolvent modes are temporally centered around $t=0$ and
exhibit a predominant streamwise component. The modes are located in a region
$x_{2}/\delta<0.25$, which corresponds to $x_{2}^{+}<45$, i.e., the buffer
region. Thus, at $t=0$, the modes capture the highly energetic near-wall
streaks. The subsequent temporal decay of these modes can be explained by the
changing flow conditions, notably the growth of the spanwise wall-shear stress
$\tau_{3}$, and consequently $u_{\tau}$ (see Figure 10). Under these
conditions, the spatial scales preferred by the near-wall streaks stretch as
$u_{\tau}$ increases and the wall-shear stress tensor rotates toward the
$x_{3}$ direction.

The second and third pairs of spatial scales are chosen so that the resolvent
modes can capture the near-wall streaks under the new shear condition at times
$t\Omega=1.3$ and $t\Omega=1.94$ respectively. To take into account the
stronger mean shear in the spanwise direction, which increases proportionally
to $dP/dx_{3}$ [64], the quantities in wall units must be scaled by a factor
of $u_{\tau}(t)/u_{\tau,0}$. We also take into account the new orientation of
the streaks by applying a rotation by the wall-shear stress angle
$\gamma(t)=\tan^{-1}(\tau_{3}/\tau_{1})$, where $\tau_{i}$ is the
instantaneous wall-shear stress in the $x_{i}$ direction (see Fig. 10). The
new spatial scales are calculated by rotating a box of size
$\lambda_{1}^{+}=100(u_{\tau}(t)/u_{\tau,0})$ and
$\lambda_{3}^{+}=1000(u_{\tau}(t)/u_{\tau,0})$ with angle $\gamma(t)$ and
finding the length scales aligned with the $x_{1}$ and $x_{3}$ axis. The
response mode for the second pair of spatial scales, $\lambda_{1}^{+}=189$ and
$\lambda_{3}^{+}=1890$, correspond to $t\Omega=1.3$ and are plotted in
11(c,d). The frequency-time map of the modes is shown in Fig. 12(b). Similar
to the first case, the modes are centered around $t\Omega=1.3$, indicating
that the wavelet-based resolvent analysis is able to identify the
nonequilibrium effects of the non-stationary flow. We note that the spanwise
component of the response mode is much more dominant than the streamwise
component, which reflects the new wall-shear angle $\gamma=75.7^{\circ}$. We
also observe that the modes are closer to the wall, which depict the reduction
of the buffer layer as a result of the increase in $u_{\tau}$. Finally, for
the third case, we choose $\lambda_{1}^{+}=297$ and $\lambda_{3}^{+}=2970$. We
observe that these modes (Fig. 11e,f and Fig. 12c) are not centered around
$t\Omega=1.94$. We speculate that this is due to the temporal boundary
condition at $t\Omega=2.34$. As the flow is not at a statistically-steady
state at this time, a Neumann boundary condition may not be the most suitable
boundary condition. The modes cannot grow beyond the boundary due to the
boundary condition and are artificially damped near the end of the temporal
domain.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 11: Streamwise (left) and spanwise (right) velocity components of the
resolvent output modes, for spatial scales of (a, b) $\lambda_{1}^{+}=1000$,
$\lambda_{3}^{+}=100$, (c, d) $\lambda_{1}^{+}=189$, $\lambda_{3}^{+}=1890$,
and (e, f) $\lambda_{1}^{+}=297$, $\lambda_{3}^{+}=2970$. The vertical dashed
lines mark (a, b) $t\Omega=0$, (c, d) $t\Omega=1.38$, and (e, f)
$t\Omega=1.94$

(a)

(b)

(c)

Figure 12: Magnitude of the principal response mode in the frequency-time
plane for (a) $\lambda_{1}^{+}=1000$, $\lambda_{3}^{+}=100$, (b)
$\lambda_{1}^{+}=189$, $\lambda_{3}^{+}=1890$, and (c) $\lambda_{1}^{+}=297$,
$\lambda_{3}^{+}=2970$.

##  6 Summary

This work expands the resolvent analysis framework to non-stationary flow
problems. The resolvent operator is traditionally constructed for flow
quantities that are Fourier-transformed in the homogeneous spatial directions
and in time. Such a resolvent operator cannot be used to study time-localized
nonlinear forcing or a time-varying mean flow. Instead, we construct a
wavelet-based resolvent operator, applying a wavelet transform in time while
keeping the Fourier transform for the homogeneous spatial directions. This
resolvent operator, provided we use an orthonormal wavelet basis, is
equivalent to the Fourier-based resolvent analysis for statistically
stationary flows. Even in stationary cases, wavelet-based resolvent analysis
can be modified to explore the effects of transient forcing localized to time
scales of interest, such as those characterizing the buffer or logarithmic
layers.

In the channel flow case, the wavelet-based resolvent modes are able to
capture the transient growth expected for non-normal systems. We observe that
the input modes precede the output modes, opening the possibility to study
causality in turbulent flows using resolvent analysis. The wavelet-based
resolvent analysis is notable in its ability to reflect the effects of a non-
stationary mean flow. In the case of the Stokes oscillatory flow, the
resolvent modes show increased sensitivity to forcing and perturbation
amplification near the peaks of the streamwise root-mean-square velocity. The
wavelet-based resolvent modes allow us to track the spatial and temporal
location of the peak amplification alongside the varying mean flow. Finally,
for the 3D channel flow, the resolvent modes are able to identify the effect
of the varying flow conditions, mainly the increasing shear velocity and
rotating wall shear stress, on the principal resolvent modes. We compute the
resolvent modes using the length scales preferred by near-wall streaks for
flow conditions at three different times. The resulting resolvent response
modes peak around the chosen times, with the exception of the time close to
the end of the temporal domain. The predominant velocity component for the
resolvent modes also shifts from the streamwise component to the spanwise one,
mirroring the reorientation of the mean flow. Wavelet resolvent modes reflect
time-varying mean flow conditions and help locate energetic near-wall streaks
in space and time, and identify their preferred spatial scales. This can shed
light on the flow conditions that amplify these coherent structures. The cases
considered in this work thus show the value of wavelet-based resolvent
analysis as a new tool to study non-stationary turbulent flows.

## Acknowledgments

The authors acknowledge support from the Air Force Office of Scientific
Research under grant number FA9550-22-1-0109.

## References

  * Jovanović and Bamieh [2005] Jovanović, M. R., and Bamieh, B., “Componentwise energy amplification in channel flows,” _J. Fluid Mech._ , Vol. 534, 2005, pp. 145–183.
  * McKeon and Sharma [2010] McKeon, B. J., and Sharma, A. S., “A critical-layer framework for turbulent pipe flow,” _J. Fluid Mech._ , Vol. 658, 2010, pp. 336–382.
  * Sharma and McKeon [2013] Sharma, A. S., and McKeon, B. J., “On coherent structure in wall turbulence,” _J. Fluid Mech._ , Vol. 728, 2013, pp. 196–238.
  * Moarref et al. [2013] Moarref, R., Sharma, A. S., McKeon, B. J., and Tropp, J. A., “Model-based scaling and prediction of the streamwise energy intensity in high-Reynolds number turbulent channels,” _J. Fluid Mech._ , Vol. 734, 2013, pp. 275–316.
  * Luhar et al. [2014] Luhar, M., Sharma, A. S., and McKeon, B. J., “Opposition control within the resolvent analysis framework,” _J. Fluid Mech._ , Vol. 749, 2014, pp. 597–626.
  * Martini et al. [2020] Martini, E., Cavalieri, A. V. G., Jordan, P., Towne, A., and Lesshafft, L., “Resolvent-based optimal estimation of transitional and turbulent flows,” _J. Fluid Mech._ , Vol. 900, 2020, p. A2.
  * Towne et al. [2020] Towne, A., Lozano-Durán, A., and Yang, X. I. A., “Resolvent-based estimation of space–time flow statistics,” _J. Fluid Mech._ , Vol. 883, 2020, p. A17.
  * Bae et al. [2021] Bae, H. J., Lozano-Durán, A., and McKeon, B. J., “Nonlinear mechanism of the self-sustaining process in the buffer and logarithmic layer of wall-bounded flows,” _J. Fluid Mech._ , Vol. 914, 2021, p. A3.
  * Yeh and Taira [2019] Yeh, C.-A., and Taira, K., “Resolvent-analysis-based design of airfoil separation control,” _J. Fluid Mech._ , Vol. 867, 2019, pp. 572–610.
  * Bae et al. [2020a] Bae, H. J., Dawson, S. T. M., and McKeon, B. J., “Resolvent-based study of compressibility effects on supersonic turbulent boundary layers,” _J. Fluid Mech._ , Vol. 883, 2020a, p. A29.
  * Bae et al. [2020b] Bae, H. J., Dawson, S. T. M., and McKeon, B. J., “Studying the effect of wall cooling in supersonic boundary layer flow using resolvent analysis,” _AIAA Scitech 2020 Forum_ , 2020b, p. 0575.
  * Ahmed et al. [2021] Ahmed, M. A., Bae, H. J., Thompson, A. F., and McKeon, B. J., “Resolvent analysis of stratification effects on wall-bounded shear flows,” _Phys. Rev. Fluids_ , Vol. 6, 2021, p. 084804.
  * Padovan et al. [2020] Padovan, A., Otto, S. E., and Rowley, C. W., “Analysis of amplification mechanisms and cross-frequency interactions in nonlinear flows via the harmonic resolvent,” _J. Fluid Mech._ , Vol. 900, 2020, p. A14.
  * Meyer [1992] Meyer, Y., _Wavelets and Operators: Volume 1_ , Cambridge University Press, 1992\.
  * Meneveau [1991] Meneveau, C., “Analysis of turbulence in the orthonormal wavelet representation,” _J. Fluid Mech._ , Vol. 232, 1991, pp. 469–520.
  * Lewalle [1993] Lewalle, J., “Wavelet transforms of the Navier-Stokes equations and the generalized dimensions of turbulence,” _Appl. Sci. Res._ , Vol. 51, No. 1-2, 1993, pp. 109–113.
  * Ren et al. [2021] Ren, J., Mao, X., and Fu, S., “Image-based flow decomposition using empirical wavelet transform,” _J. Fluid Mech._ , Vol. 906, 2021, p. A22.
  * Floryan and Graham [2021] Floryan, D., and Graham, M. D., “Discovering multiscale and self-similar structure with data-driven wavelets,” _Proc. Natl. Acad. Sci. USA_ , Vol. 118, No. 1, 2021, p. e2021299118.
  * Towne et al. [2018] Towne, A., Schmidt, O. T., and Colonius, T., “Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis,” _J. Fluid Mech._ , Vol. 847, 2018, pp. 821–867.
  * Abreu et al. [2020] Abreu, L. I., Cavalieri, A. V. G., Schlatter, P., Vinuesa, R., and Henningson, D. S., “Spectral proper orthogonal decomposition and resolvent analysis of near-wall coherent structures in turbulent pipe flows,” _J. Fluid Mech._ , Vol. 900, 2020, p. A11.
  * Tissot et al. [2021] Tissot, G., Cavalieri, A. V. G., and Mémin, E., “Stochastic linear modes in a turbulent channel flow,” _J. Fluid Mech._ , Vol. 912, 2021, p. A51.
  * Jeun et al. [2016] Jeun, J., Nichols, J. W., and Jovanović, M. R., “Input-output analysis of high-speed axisymmetric isothermal jet noise,” _Phys. Fluids_ , Vol. 28.4, 2016, p. 047101.
  * Kojima et al. [2020] Kojima, Y., Yeh, C., Taira, K., and Kameda, M., “Resolvent analysis on the origin of two-dimensional transonic buffet,” _J. Fluid Mech._ , Vol. 885, 2020, p. R1.
  * Mallat [2001] Mallat, S., _A Wavelet Tour of Signal Processing_ , Academic Press, 2001.
  * Najmi [2012] Najmi, A.-H., _Wavelets: A Concise Guide_ , The Johns Hopkins University Press, 2012.
  * Alappat et al. [2020] Alappat, C., Basermann, A., Bishop, A. R., Fehske, H., Hager, G., Schenk, O., Thies, J., and Wellein, G., “A recursive algebraic coloring technique for hardware-efficient symmetric sparse matrix-vector multiplication,” _ACM Trans. Parallel Comput._ , Vol. 7, No. 3, 2020, pp. 1–37.
  * Bollhöfer et al. [2020] Bollhöfer, M., Schenk, O., Janalik, R., Hamm, S., and Gullapalli, K., “State-of-the-art sparse direct solvers,” _Parallel Algorithms in Computational Science and Engineering_ , Springer, 2020, pp. 3–33.
  * Bollhöfer et al. [2019] Bollhöfer, M., A., E., Scheidegger, S., and Schenk, O., “Large-Scale Sparse Inverse Covariance Matrix Estimation,” _SIAM J. Sci. Comput._ , Vol. 41, No. 1, 2019, pp. A380–A401.
  * Simon and Zha [2000] Simon, H. D., and Zha, H., “Low-rank matrix approximation using the Lanczos bidiagonalization process with applications,” _SIAM Journal on Scientific Computing_ , Vol. 21, No. 6, 2000, pp. 2257–2274.
  * Halko et al. [2011] Halko, N., Martinsson, P.-G., and Tropp, J. A., “Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions,” _SIAM review_ , Vol. 53, No. 2, 2011, pp. 217–288.
  * Drineas and Mahoney [2016] Drineas, P., and Mahoney, M. W., “RandNLA: randomized numerical linear algebra,” _Communications of the ACM_ , Vol. 59, No. 6, 2016, pp. 80–90.
  * Tropp et al. [2017] Tropp, J. A., Yurtsever, A., Udell, M., and Cevher, V., “Practical sketching algorithms for low-rank matrix approximation,” _SIAM Journal on Matrix Analysis and Applications_ , Vol. 38, No. 4, 2017, pp. 1454–1485.
  * Ribeiro et al. [2020] Ribeiro, J. H. M., Yeh, C.-A., and Taira, K., “Randomized resolvent analysis,” _Phys. Rev. Fluids_ , Vol. 5, No. 3, 2020, p. 033902.
  * Yeh et al. [2020] Yeh, C.-A., Benton, S. I., Taira, K., and Garmann, D. J., “Resolvent analysis of an airfoil laminar separation bubble at Re= 500 000,” _Physical Review Fluids_ , Vol. 5, No. 8, 2020, p. 083906.
  * Orlandi [2000] Orlandi, P., _Fluid Flow Phenomena: A Numerical Toolkit_ , Fluid Flow Phenomena: A Numerical Toolkit, Springer, 2000.
  * Kim and Moin [1985] Kim, J., and Moin, P., “Application of a fractional-step method to incompressible Navier-Stokes equations,” _J. Comp. Phys._ , Vol. 59, 1985, pp. 308–323.
  * Wray [1990] Wray, A. A., “Minimal-storage time advancement schemes for spectral methods,” Tech. rep., NASA Ames Research Center, 1990.
  * Bae et al. [2018] Bae, H. J., Lozano-Durán, A., Bose, S. T., and Moin, P., “Turbulence intensities in large-eddy simulation of wall-bounded flows,” _Phys. Rev. Fluids_ , Vol. 3, 2018, p. 014610.
  * Bae et al. [2019] Bae, H. J., Lozano-Durán, A., Bose, S. T., and Moin, P., “Dynamic slip wall model for large-eddy simulation,” _J. Fluid Mech._ , Vol. 859, 2019, pp. 400–432.
  * Lozano-Durán and Bae [2019] Lozano-Durán, A., and Bae, H. J., “Characteristic scales of Townsend’s wall-attached eddies.” _J. Fluid Mech._ , Vol. 868, 2019, p. 698\.
  * Lozano-Durán et al. [2018] Lozano-Durán, A., Hack, M. J. P., and Moin, P., “Modeling boundary-layer transition in direct and large-eddy simulations using parabolized stability equations,” _Phys. Rev. Fluids_ , Vol. 3, 2018, p. 023901\.
  * Schmid and Henningson [2000] Schmid, P. J., and Henningson, D. S., _Stability and Transition in Shear Flows_ , Vol. 142, Springer Science & Business Media, 2000.
  * Orr [1907] Orr, W., “The stability or instability of the steady motions of a perfect liquid and of a viscous liquid. Part I: A perfect liquid,” _Proc. R. Ir. Acad._ , Vol. A 27, 1907, pp. 9–68.
  * Landahl [1975] Landahl, M. T., “Wave breakdown and turbulence,” _SIAM J. Appl. Math_ , Vol. 28, No. 4, 1975, pp. 735–756.
  * Jiménez [2013] Jiménez, J., “How linear is wall-bounded turbulence?” _J. Fluid Mech._ , Vol. 25, 2013, p. 110814.
  * Jiménez and Moin [1991] Jiménez, J., and Moin, P., “The minimal flow unit in near-wall turbulence,” _J. Fluid Mech._ , Vol. 225, 1991, pp. 213–240.
  * Hamilton et al. [1995] Hamilton, J. M., Kim, J., and Waleffe, F., “Regeneration mechanisms of near-wall turbulence structures,” _J. Fluid Mech._ , Vol. 287, 1995, pp. 317–348.
  * Jiménez and Pinelli [1999] Jiménez, J., and Pinelli, A., “The autonomous cycle of near-wall turbulence,” _J. Fluid Mech._ , Vol. 389, 1999, pp. 335–359.
  * Waleffe [1997] Waleffe, F., “On a self-sustaining process in shear flows,” _Phys. Fluids_ , Vol. 9, No. 4, 1997, pp. 883–900.
  * Schoppa and Hussain [2002] Schoppa, W., and Hussain, F., “Coherent structure generation in near-wall turbulence,” _J. Fluid Mech._ , Vol. 453, 2002, pp. 57–108.
  * Farrell et al. [2017] Farrell, B. F., Gayme, D. F., and Ioannou, P. J., “A statistical state dynamics approach to wall turbulence,” _Phil. Trans. R. Soc. Lond. A_ , Vol. 375, No. 2089, 2017, p. 20160081.
  * Hino et al. [1976] Hino, M., Sawamoto, M., and Takasu, S., “Experiments on transition to turbulence in an oscillatory pipe flow,” _J. Fluid Mech._ , Vol. 75, No. 2, 1976, pp. 193–207.
  * Akhavan et al. [1991] Akhavan, R., Kamm, R. D., and Shapiro, A. H., “An investigation of transition to turbulence in bounded oscillatory Stokes flows Part 1. Experiments,” _J. Fluid Mech._ , Vol. 225, 1991, pp. 395–422.
  * Verzicco and Vittori [1996] Verzicco, R., and Vittori, G., “Direct simulation of transition in Stokes boundary layers,” _Phys. Fluids_ , Vol. 8, No. 6, 1996, pp. 1341–1343.
  * Vittori and Verzicco [1998] Vittori, G., and Verzicco, R., “Direct simulation of transition in an oscillatory boundary layer,” _J. Fluid Mech._ , Vol. 371, 1998, pp. 207–232.
  * Costamagna et al. [2003] Costamagna, P., Vittori, G., and Blondeaux, P., “Coherent structures in oscillatory boundary layers,” _J. Fluid Mech._ , Vol. 474, 2003, pp. 1–33.
  * Spalart and Baldwin [1989] Spalart, P. R., and Baldwin, B. S., “Direct simulation of a turbulent oscillating boundary layer,” _Turbulent shear flows 6_ , Springer, 1989, pp. 417–440.
  * Jensen et al. [1989] Jensen, B. L., Sumer, B. M., and Fredsøe, J., “Turbulent oscillatory boundary layers at high Reynolds numbers,” _J. Fluid Mech._ , Vol. 206, 1989, pp. 265–297.
  * Von Kerczek and Davis [1974] Von Kerczek, C., and Davis, S. H., “Linear stability theory of oscillatory Stokes layers,” _Journal of Fluid Mechanics_ , Vol. 62, No. 4, 1974, pp. 753–773.
  * Sarpkaya [1993] Sarpkaya, T., “Coherent structures in oscillatory boundary layers,” _Journal of Fluid Mechanics_ , Vol. 253, 1993, pp. 105–140.
  * Blondeaux and Vittori [1994] Blondeaux, P., and Vittori, G., “Wall imperfections as a triggering mechanism for Stokes-layer transition,” _Journal of Fluid Mechanics_ , Vol. 264, 1994, pp. 107–135.
  * Carstensen et al. [2010] Carstensen, S., Sumer, B. M., and Fredsøe, J., “Coherent structures in wave boundary layers. Part 1. Oscillatory motion,” _Journal of Fluid Mechanics_ , Vol. 646, 2010, pp. 169–206.
  * Ozdemir et al. [2014] Ozdemir, C. E., Hsu, T.-J., and Balachandar, S., “Direct numerical simulations of transition and turbulence in smooth-walled Stokes boundary layer,” _Physics of Fluids_ , Vol. 26, No. 4, 2014, p. 045108.
  * Lozano-Durán et al. [2021] Lozano-Durán, A., Giometto, M. G., Park, G. I., and Moin, P., “Non-equlibrium three-dimensional boundary layers at moderate Reynolds numbers,” _J. Fluid Mech._ , Vol. 883, 2021, p. A20.
  * Moin et al. [1990] Moin, P., Shih, T.-H., Driver, D. M., and Mansour, N. N., “Direct numerical simulation of a three-dimensional turbulent boundary layer,” _Phys. Fluids_ , Vol. 2, No. 10, 1990, pp. 1846–1853.