Sparse space–time resolvent analysis for statistically stationary and time-varying flows

Abstract

Resolvent analysis provides a framework to predict coherent spatio-temporal structures of the largest linear energy amplification, through a singular value decomposition (SVD) of the resolvent operator, obtained by linearising the Navier–Stokes equations about a known turbulent mean velocity profile. Resolvent analysis utilizes a Fourier decomposition in time, which has thus far limited its application to statistically stationary or time-periodic flows. This work develops a variant of resolvent analysis applicable to time-evolving flows, and proposes a variant that identifies spatio-temporally sparse structures, applicable to either stationary or time-varying mean velocity profiles. Spatio-temporal resolvent analysis is formulated through the incorporation of the temporal dimension to the numerical domain via a discrete time-differentiation operator. Sparsity (which manifests in localisation) is achieved through the addition of an $l_1$-norm penalisation term to the optimisation associated with the SVD. This modified optimisation problem can be formulated as a nonlinear eigenproblem and solved via an inverse power method. We first showcase the implementation of the sparse analysis on a statistically stationary turbulent channel flow, and demonstrate that the sparse variant can identify aspects of the physics not directly evident from standard resolvent analysis. This is followed by applying the sparse space–time formulation on systems that are time varying: a time-periodic turbulent Stokes boundary layer and then a turbulent channel flow with a sudden imposition of a lateral pressure gradient, with the original streamwise pressure gradient unchanged. We present results demonstrating how the sparsity-promoting variant can either change the quantitative structure of the leading space–time modes to increase their sparsity, or identify entirely different linear amplification mechanisms compared with non-sparse resolvent analysis.

1. Introduction

Despite the highly nonlinear nature of turbulent fluid flows, linearised analyses of the governing Navier–Stokes equations have proven to be effective at capturing several pertinent properties of such systems. Resolvent analysis, first applied as a model for turbulent pipe flows (McKeon & Sharma 2010), has been informative in other problems that involve wall-bounded turbulence, spatio-temporal flow statistics (Towne, Lozano-Durán & Yang 2020) and the identification of coherent structures in turbulent jets (Lesshafft et al. 2019; Pickering et al. 2021), airfoils (Yeh & Taira 2019) supersonic boundary layers (Bae, Dawson & McKeon 2020a,b) and turbulent rectangular duct flows (Lopez-Doriga, Dawson & Vinuesa 2022). Combining a set of triadically consistent resolvent modes has also been used for the representation of hairpin structures in Sharma & McKeon (2013) and, more generally, for the reconstruction of phenomena observed in wall-bounded turbulence (McKeon 2017). This framework has further been applied for the estimation of flow states (Gómez et al. 2016; Beneddine et al. 2017; Illingworth, Monty & Marusic 2018; Symon et al. 2020), the prediction of coherent structures (Abreu et al. 2020; Tissot, Cavalieri & Mémin 2021) and statistical quantities and scalings (Hwang & Cossu 2010; Zare, Jovanović & Georgiou 2017; Towne et al. 2020), designing control strategies for drag reduction (Luhar, Sharma & McKeon 2014b; Toedtli, Luhar & McKeon 2019) and modelling the effect of complex surfaces (Luhar, Sharma & McKeon 2015; Chavarin & Luhar 2019). The broad applicability of such linearised analysis relies on (and can be seen as evidence to infer) the importance of linear amplification mechanisms in the generation and evolution of empirically observed coherent structures within turbulent flows, such as near-wall streaks (Kline et al. 1967), hairpin vortices (Theodorsen 1952; Head & Bandyopadhyay 1981), superstructures (Kim & Adrian 1999) and a range of other coherent features described in Jiménez (2018).

The linearised analyses discussed thus far assume that the linear system under investigation is time invariant, as is the case when the underlying flow is statistically stationary. The recently developed harmonic resolvent analysis (Padovan, Otto & Rowley 2020; Padovan & Rowley 2022) enables the resolvent framework to extend to statistically time-periodic flows, enabling cross-frequency analysis capturing triadic interactions between a time-periodic base flow and fluctuations about this mean state at other frequencies. In the context of flow control over a periodically plunging cylinder, Lin, Tsai & Tsai (2023) utilizes a Lyapunouv–Floquet transformation to map the corresponding linear time-periodic system to a time-invariant equivalent, enabling the application of standard resolvent analysis methods. Other noteworthy contributions to the study of time-varying linear systems include the linear stability analyses compiled in Blennerhassett & Bassom (2002, 2006, 2007), for flat and high-frequency oscillatory Stokes boundary layers, as well as an oscillating cylinder.

Whether considering a statistically stationary or a time-periodic mean state, the methods discussed thus far consider a Fourier decomposition in time. Typically, this involves identifying the forcing (input) and response (output) structures corresponding to the largest energy amplification by the linearised system (represented by the resolvent operator). While this decomposition arises naturally for such methods, it can potentially obscure the intermittent nature of velocity fluctuations present in turbulent flows. Alternative linear analyses methods can be similarly restrictive, with asymptotic stability analysis also identifying eigenmodes each associated with a single (possibly complex) frequency. Conversely, transient growth analysis (Böberg & Brösa 1988; Butler & Farrell 1992; Reddy & Henningson 1993; Schmid 2007) considers the unforced response to a specific initial condition, corresponding to maximal energy growth over a specified time horizon. This again is unrealistic for systems subject to continuous perturbations (Jovanović & Bamieh 2005), though such analysis has been used in turbulent flows, such as to predict the emergence of near-wall streamwise streaks (Del Alamo & Jimenez 2006) and vortices (Schoppa & Hussain 2002) in wall-bounded turbulence. To overcome these limitations, here we introduce a space–time formulation of the resolvent operator that is firstly applicable to non-statistically stationary systems with an arbitrarily time-varying mean profile, and secondly allows for the identification of optimal input and output trajectories that can have arbitrary time dependence. This builds upon preliminary work first reported in Lopez-Doriga et al. (2023). While not explored here, related work also considers explicitly replacing the Fourier transform used in standard resolvent analysis with a wavelet transform (Ballouz et al. 2023, 2024).

This generalization of operator-based decompositions to enable non-Fourier temporal modes is somewhat analogous to efforts to similarly generalize data-driven proper orthogonal decomposition (POD) methodology to identify intermittent behaviour in turbulent flows, such as the conditional POD formulated in Schmidt & Schmid (2019) and time-windowed space–time POD described in Frame & Towne (2022). Note that spectral POD (Towne, Schmidt & Colonius 2018) has also been recently generalized for time-periodic systems using cyclostationary analysis (Heidt & Colonius 2024).

Methods to identify non-modal linear energy amplification such as resolvent or transient growth analysis involve computing the leading singular values and vectors of an appropriately defined linear operator. The singular value decomposition (SVD), by design, is defined as an optimisation problem that involves an $l_2$-energy norm. In the context of resolvent analysis, this optimisation problem relates to the energy ratio between input and output flow states, and naturally yields spatio-temporal structures that are Fourier modes in time. Here, we consider modifications to the standard optimisation problem that yield alternative temporal functions, which are inclined to be localised in time. This is achieved by incorporating an $l_1$-norm term into the optimisation problem. The use of $l_1$ norms to promote localisation and/or sparsity has origins in compressive sensing (Candès & Wakin 2008).

In the context of fluid mechanics, sparsity-promoting methods have been utilized for developing reduced-complexity models across a number of contexts. These include the identification of sparse nonlinear reduced-order models (Brunton, Proctor & Kutz 2016; Loiseau & Brunton 2018; Rubini, Lasagna & Da Ronch 2020), the selection of a sparse set of active dynamic modes (Jovanović, Schmid & Nichols 2014) and in the reconstruction of temporal spectral content from data that is under-resolved in time (Tu et al. 2014). Recently, sparsity-promoting methods have also been incorporated in the resolvent analysis framework in Skene et al. (2022), where they are used to identify spatially localised forcing modes, which can be more directly useful for actuator placement in flow control applications. In Skene et al. (2022) a Riemannian optimisation process is used to solve an $l_1$-based optimisation problem, following a similar approach used by Foures, Caulfield & Schmid (2013) to identify spatially localised structures in transient growth analysis. The present work is similarly motivated, though we focus here on achieving localisation in time as well as space. We also use a different formulation of the optimisation problem, which allows for a balance between $l_1$- and $l_2$-norm contributions.

The structure of the paper is as follows. A discussion of the fundamentals of pseudospectral analysis and the wall-normal derivation of the governing equations, along with the space–time form of the resolvent operator, and a description of the algorithm that promotes sparsity on the resolvent modes are presented in § 2. The main results of our investigation are discussed in § 3: the sparse formulation of the standard resolvent operator is applied in the streamwise and spanwise directions of a turbulent channel flow in § 3.2; and the space–time and sparse space–time formulations of the resolvent operator are applied on a turbulent channel flow in § 3.3; a turbulent Stokes boundary layer in § 3.4; and a channel flow with a sudden lateral pressure gradient in § 3.5. Finally, we discuss the main findings and future prospects of our investigation in § 4.

2. Methodology

This section begins with a brief overview of the fundamentals of pseudospectral analysis of linear operators in § 2.1. This is followed by a derivation of the resolvent formulation of the incompressible Navier–Stokes equations in wall-normal velocity and vorticity variables in § 2.2, assuming homogeneity in both the spatial and temporal dimensions. This is followed by the development of a space–time resolvent operator where homogeneity is not assumed in the temporal dimension in § 2.3, also in wall-normal velocity and vorticity variables. Following this, § 2.4 introduces a formulation of resolvent analysis that promotes sparsity on the optimal resolvent modes.

2.1. Pseudospectral analysis of linear operators

Let us consider a dynamical system governed by

(2.1)\begin{equation} \dot{\boldsymbol{q}} (\boldsymbol{x},t)+ \mathcal{L} \boldsymbol{q}(\boldsymbol{x},t) = \boldsymbol{f}(\boldsymbol{x},t), \end{equation}

where $\boldsymbol {q}$ denotes the state of the system with respect to a reference state $\boldsymbol {q}_0$, $\mathcal {L}$ is a linear operator and $\boldsymbol {f}$ represents an exogenous input or forcing. The space and time dimensions are denoted by $\boldsymbol {x}$ and $t$, respectively. Assuming that the system is homogeneous in the temporal dimension, we propose solutions of the form (Schmid & Henningson 2001) $\boldsymbol {q}(\boldsymbol {x},t)=\hat {\boldsymbol {q}}(\boldsymbol {x})\exp {(-{\rm i}\omega t)}$ with $\omega \in \mathbb {C}$, and substituting in (2.1) gives

(2.2)\begin{equation} (-{\rm i}\omega \boldsymbol{\mathsf{I}} + \mathcal{L})\hat{\boldsymbol{q}}(\boldsymbol{x}) = \hat{\boldsymbol{f}}(\boldsymbol{x}). \end{equation}

In the case the forcing term is non-zero, the elements can be rearranged so that the governing equation represents the following system:

(2.3)\begin{equation} \hat{\boldsymbol{q}}(\boldsymbol{x}) = (-{\rm i}\omega \boldsymbol{\mathsf{I}} + \mathcal{L})^{{-}1 }\hat{\boldsymbol{f}}(\boldsymbol{x}):=\mathcal{H}_\omega \hat{\boldsymbol{f}}(\boldsymbol{x}). \end{equation}

We refer to $\mathcal {H}_\omega$ as the resolvent operator. Note that the subscript $\omega$ is retained to highlight the dependence on the temporal frequency. The original dynamical system in (2.1) has been recast as a linear mapping between a forcing $\hat {\boldsymbol {f}}$ and the state $\hat {\boldsymbol {q}}$.

According to (2.3), the properties of the state $\hat {\boldsymbol {q}}$ will be affected by both the nature of the forcing $\hat {\boldsymbol {f}}$ and the properties of the resolvent $\mathcal {H}_\omega$. In this work we focus in particular on the pairs of forcing and response modes that produce the largest amplification through the action of $\mathcal {H}_\omega$. That is, a forcing of small magnitude yields a response of large magnitude. Such structures can be identified via a SVD of the resolvent operator $\mathcal {H}_\omega$ as

(2.4)\begin{equation} \mathcal{H}_\omega = \sum_{j=1}^N \boldsymbol{\psi}_j \sigma_j \boldsymbol{\phi}_j^*, \end{equation}

where $\sigma _j \geqslant \sigma _{j+1} \geqslant 0$ for all $j$, and $({\cdot }^*)$ denotes the adjoint. Note that here the resolvent operator will take the form of a discretised operator, therefore, the summation in (2.4) is truncated to $N$ terms.

In particular, we seek $\hat {\boldsymbol {f}} = \boldsymbol {\phi }_1$ that maximises the largest singular value $\sigma _1$, where

(2.5)\begin{equation} \sigma_1 = \max_{\boldsymbol{\phi}} \frac{\|\mathcal{H}_{\omega} \boldsymbol{\phi}\|_2}{\|\boldsymbol{\phi}\|_2}=\min_{\boldsymbol{\psi}}\frac{\|\boldsymbol{\psi}\|_2}{\|\mathcal{H}_{\omega}^*\boldsymbol{\psi}\|_2}, \end{equation}

where the $l_2$ norm is taken over the spatial domain $\varOmega _{\boldsymbol {x}}$, so that, for example,

(2.6)\begin{equation} \|\boldsymbol{\psi}(\boldsymbol{x}) \|_2 = \left(\int_{\boldsymbol{x}\in \varOmega_{\boldsymbol{x}}} \left|\boldsymbol{\psi}(\boldsymbol{x})\right|^2\text{d}\kern 0.7pt\boldsymbol{x}\right)^{1/2}. \end{equation}

Alternatively, we can write this optimisation problem in terms of the leading forcing mode $\boldsymbol {\phi }_1$ as

(2.7)\begin{equation} \boldsymbol{\phi}_1 = \mathop{\rm arg\,max}\limits_{\boldsymbol{\phi}} \frac{\|\mathcal{H}_{\omega} \boldsymbol{\phi}\|_2}{\|\boldsymbol{\phi}\|_2}, \end{equation}

or the leading response mode $\boldsymbol {\psi }_1$,

(2.8)\begin{equation} \psi_1 = \mathop{\rm arg\,min}\limits_\psi \frac{\|\psi\|_2}{\|\mathcal{H}_{\omega}^*\psi\|_2}. \end{equation}

2.2. Resolvent formulation of the mean-linearised incompressible Navier–Stokes equations

The incompressible Navier–Stokes equations enforce conservation of momentum and mass, respectively, and are written in a Cartesian coordinate system as follows:

(2.9)$$\begin{gather} \partial_t \boldsymbol{u} ={-}\boldsymbol{u} \boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}-\boldsymbol{\nabla} p+\frac{1}{Re}\Delta\boldsymbol{u}, \end{gather}$$

(2.10)$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=0. \end{gather}$$

Here, the instantaneous velocity field has three components: $\boldsymbol {u}=[u(\boldsymbol {x},t),v(\boldsymbol {x},t),$ $w(\boldsymbol {x},t)]^{\rm T}$ with $\boldsymbol {x}=[x,y,z]^{\rm T}$ and $p=p(\boldsymbol {x},t)$ representing the instantaneous pressure field. In this reference frame, $x$ and $z$ correspond to the streamwise and spanwise directions, respectively, and are nominally considered to be infinite in extent. The other variable, $y$, represents the wall-normal dimension. Here $\partial _t$ denotes a time (partial) derivative, the spatial gradient operator is given by $\boldsymbol {\nabla } = [\partial _x,\partial _y,\partial _z]^{\rm T}$, and the Laplacian operator is defined as $\varDelta = \boldsymbol {\nabla }\boldsymbol {\cdot}\boldsymbol {\nabla }$.

We can write a given instantaneous velocity state $\boldsymbol {u}$ as the sum of the temporal mean $\boldsymbol {U}$ and a fluctuating component $\boldsymbol {u}'$, such that

(2.11)\begin{equation} \boldsymbol{u}(\boldsymbol{x},t) = \boldsymbol{U}(\boldsymbol{x})+\boldsymbol{u}'(\boldsymbol{x},t). \end{equation}

Applying this decomposition in (2.9)–(2.10) and subtracting the temporal average gives the governing equations used in this work,

(2.12)$$\begin{gather} \partial_t \boldsymbol{u}' +\boldsymbol{U}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}' +\boldsymbol{u}'\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{U} + \boldsymbol{\nabla}p'-\frac{1}{Re} \Delta\boldsymbol{u}' ={-}\boldsymbol{u}'\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}'-\overline{\boldsymbol{u}'\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}'}=\boldsymbol{f}', \end{gather}$$

(2.13)$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}' = 0, \end{gather}$$

that is, conservation of momentum and continuity of the fluctuating components. Note that here the right-hand side of (2.12) has been condensed into a forcing term $\boldsymbol {f}'$ that represents the effect of the fluctuations about the mean state of the nonlinear terms, and can be regarded in this context as an exogenous input to a linear system comprising of the remaining terms. Wall-bounded parallel flows with a mean/base flow in the streamwise and spanwise dimensions $\boldsymbol {U}(y)=[U(y),0,W(y)]^{\rm T}$, admit a transformation of variables from a primitive reference $\{u',v',w',p'\}$ towards a reference in terms of the wall-normal velocity $v'$ and vorticity $\eta '$ (where $\eta ' = \partial u'/\partial z-\partial w'/\partial x$), without loss of generality (i.e. resulting in the Orr–Sommerfeld and Squire equations). This formulation is proven to be equally informative (Moarref et al. 2013; Rosenberg & McKeon 2019; McMullen, Rosenberg & McKeon 2020) for resolvent analysis of planar flows. In this reference, the no-slip and no-penetration conditions translate into $v'(y=0)=v'(y=2h)=0$, $\partial _y v'(y=0)=\partial _y v'(y=2h)=0$ and $\eta '(y=0)=\eta '(y=2h)=0$, where $h$ represents the semi-height of the domain in the wall-normal dimension. Throughout this paper, the location of the no-slip walls will coincide with $y=0$ and $y=2h$. This transformation is achieved according to the process described in Schmid & Henningson (2001), while including a spanwise component of the mean/base flow, $W$. Note that while this spanwise mean component is included in the derivation, it will only be non-zero for the configuration presented in § 3.5. The resulting wall-normal formulation of the conservation laws shown in (2.12)–(2.13) is formed by the following two equations:

(2.14)$$\begin{gather} \left[\left(\partial_t+U\partial_x+W\partial_z\right) \nabla^2-\partial_y^2 U\partial_x-\partial_y^2 W\partial_z-\frac{1}{\textit{Re}}\nabla^4\right]v'=f'_v, \end{gather}$$

(2.15)$$\begin{gather}\left(\partial_t+U\partial_x+W\partial_z-\frac{1}{\textit{Re}}\nabla^2\right)\eta'+\left(\partial_y U\partial_z -\partial_y W\partial_x\right)v'=f'_\eta. \end{gather}$$

Assuming that the system is homogeneous in the temporal dimension and the streamwise and spanwise directions, we introduce assumed solutions of the form

(2.16a)$$\begin{gather} v'(\boldsymbol{x},t)=\hat{v} (y)\exp{[{\rm i}(k_x x+k_z z-\omega t)]}, \end{gather}$$

(2.16b)$$\begin{gather}\eta'(\boldsymbol{x},t)=\hat{\eta} (y)\exp{[{\rm i}(k_x x+k_z z-\omega t)]}, \end{gather}$$

(2.16c)$$\begin{gather}\boldsymbol{f}'(\boldsymbol{x},t)=\widehat{\boldsymbol{f}} (y)\exp{[{\rm i}(k_x x+k_z z-\omega t)]}, \end{gather}$$

where $k_x$ and $k_z$ denote the streamwise and spanwise wavenumbers, respectively, and $\omega$ denotes temporal frequency. Substituting these assumed solutions in (2.14)–(2.15) gives the following system of equations:

(2.17)\begin{equation} (-{\rm i}\omega \boldsymbol{\mathsf{M}}+ \boldsymbol{\mathsf{L}}) \begin{pmatrix} \hat{v} \\ \hat{\eta} \end{pmatrix} = \begin{pmatrix} \,\hat{f}_{v} \\ \hat{f}_{\eta} \end{pmatrix}. \end{equation}

Here

(2.18a,b)\begin{equation} \boldsymbol{\mathsf{M}}= \begin{pmatrix} \hat{\varDelta} & 0\\ 0 & \boldsymbol{\mathsf{I}} \end{pmatrix} \quad\text{and}\quad \boldsymbol{\mathsf{L}}= \begin{pmatrix} \mathcal{L}_{os} & 0 \\ {\rm i}k_z \partial_yU-{\rm i}k_x \partial_yW & \mathcal{L}_{sq} \end{pmatrix}. \end{equation}

The modified Laplacian operator $\hat {\varDelta }=\partial _{yy}-(k_x^2+k_z^2)$ is introduced for simplicity, and

(2.19)$$\begin{gather} \mathcal{L}_{os}=\left({\rm i}k_xU+{\rm i}k_zW\right) \hat{\varDelta}-\frac{1}{\textit{Re}}\hat{\varDelta}^2-{\rm i}k_x\partial_y^2U-{\rm i}k_z\partial_y^2W, \end{gather}$$

(2.20)$$\begin{gather}\mathcal{L}_{sq}={\rm i}k_x U+{\rm i}k_zW-\frac{1}{\textit{Re}} \hat{\varDelta} \end{gather}$$

represent the Orr–Sommerfeld and Squire operators, respectively. Premultiplying both sides of (2.17) by $\boldsymbol{\mathsf{M}}^{-1}$ and solving for the state $[\hat {v}(y),\hat {\eta }(y)]^{\rm T}$ gives

(2.21)\begin{equation} \begin{pmatrix} \hat{v} \\ \hat{\eta} \end{pmatrix} =(-{\rm i}\omega \boldsymbol{\mathsf{I}} + \mathcal{L})^{{-}1} \begin{pmatrix} \hat{g}_{v} \\ \hat{g}_{\eta} \end{pmatrix}:=\mathcal{H}_{\boldsymbol{\omega}} \begin{pmatrix} \hat{g}_{v} \\ \hat{g}_{\eta} \end{pmatrix}, \end{equation}

where

(2.22a,b)\begin{equation} \mathcal{L}= \boldsymbol{\mathsf{M}}^{{-}1}\boldsymbol{\mathsf{L}} \quad\text{and}\quad \begin{pmatrix} \hat{g}_{v} \\ \hat{g}_{\eta} \end{pmatrix}=\boldsymbol{\mathsf{M}}^{{-}1} \begin{pmatrix} \hat{f}_{v} \\ \,\hat{f}_{\eta} \end{pmatrix}. \end{equation}

This transfer function $\mathcal {H}_{\boldsymbol {\omega }}$ is denoted as the resolvent operator, in analogy to the definition introduced for a general linear system in (2.3) when $z$ becomes $z=-{\rm i}\omega$. Note that the resolvent operator is dependent on the triad $\{\omega,k_x,k_z\}$ but, for the sake of readability, this dependence is indicated by the subscript $\boldsymbol {\omega }$. Expanding the terms in accordance with the derivation described in Rosenberg & McKeon (2019) gives

(2.23)\begin{equation} \begin{pmatrix} \hat{v} \\ \hat{\eta} \end{pmatrix}= \mathcal{H}_{\boldsymbol{\omega}} \begin{pmatrix} \hat{g}_v \\ \hat{g}_\eta \end{pmatrix}:= \begin{pmatrix} \mathcal{H}_{vv} & 0 \\ -\mathcal{H}_{\eta\eta} ({\rm i}k_z \partial_yU-{\rm i}k_x\partial_yW)\mathcal{H}_{vv} & \mathcal{H}_{\eta\eta} \end{pmatrix} \begin{pmatrix} \hat{g}_v \\ \hat{g}_\eta \end{pmatrix}, \end{equation}

where the scalar operators are defined as

(2.24)$$\begin{gather} \mathcal{H}_{vv}=(-{\rm i}\omega+\hat{\varDelta}^{{-}1}\mathcal{L}_{os})^{{-}1}, \end{gather}$$

(2.25)$$\begin{gather}\mathcal{H}_{\eta\eta}=(-{\rm i}\omega+\mathcal{L}_{sq})^{{-}1}. \end{gather}$$

Note that the formulation in wall-normal variables enables the study of the dynamical properties of each of the variables independently. Nevertheless, it is possible to define a direct transformation of the resolvent operator, as well as the resolvent modes, from wall-normal velocity and vorticity $\{v,\eta \}$ to primitive velocity variables $\{u,v,w\}$ (and vice versa), according to the mapping that was introduced in Meseguer & Trefethen (2003) and further developed and applied in Jovanović & Bamieh (2005), McKeon & Sharma (2010), Moarref et al. (2013) and Sharma, Moarref & McKeon (2017). The cited transformation recasts the response and forcing modes in primitive variables as

(2.26)\begin{equation} \boldsymbol{\varPsi}_{(u,v,w)} = \boldsymbol{\mathsf{C}} \boldsymbol{\varPsi}_{(v,\eta)} \end{equation}

and

(2.27)\begin{equation} \boldsymbol{\varPhi}_{(u,v,w)} = (\boldsymbol{\mathsf{M}}^{{-}1}\boldsymbol{\mathsf{B}})^{{-}1} \boldsymbol{\varPhi}_{(v,\eta)}. \end{equation}

Here

(2.28)\begin{equation} \boldsymbol{\mathsf{C}} = \frac{1}{k_\perp^2} \begin{pmatrix} {\rm i}k_x\partial_y & -{\rm i}k_z \\ k_\perp^2 & 0 \\ {\rm i}k_z\partial_y & {\rm i}k_x \end{pmatrix} \end{equation}

and

(2.29)\begin{equation} \boldsymbol{\mathsf{B}} = \begin{pmatrix} -{\rm i}k_x\partial_y & -(k_x^2+k_z^2) & -{\rm i}k_z\partial_y \\ {\rm i}k_z & 0 & -{\rm i}k_x \end{pmatrix}, \end{equation}

represent the input and output matrices, respectively.

2.3. Space–time resolvent analysis

Here we present a form of resolvent analysis that is applicable to time-varying systems. This generalization is achieved by limiting the assumed directions of homogeneity to the streamwise and spanwise spatial dimensions. Thus, in this formulation, both components of the mean state $\boldsymbol {U}=[U,0,W]^{\rm T}$ are also assumed to be temporally dependent, and we write a generalized instantaneous state $\boldsymbol {u}$ as

(2.30)\begin{equation} \boldsymbol{u}(\boldsymbol{x},t) = \boldsymbol{U}(\boldsymbol{x},t)+\boldsymbol{u}'(\boldsymbol{x},t). \end{equation}

In analogy to the trajectories presented in (2.16) we let the solutions take the following form:

(2.31a)$$\begin{gather} v'(\boldsymbol{x},t)=\hat{v} (y,t)\exp{[{\rm i}(k_x x+k_z z)]}, \end{gather}$$

(2.31b)$$\begin{gather}\eta'(\boldsymbol{x},t)=\hat{\eta} (y,t)\exp{[{\rm i}(k_x x+k_z z)]}, \end{gather}$$

(2.31c)$$\begin{gather}\boldsymbol{f}'(\boldsymbol{x},t)=\hat{\boldsymbol{f}} (y,t)\exp{[{\rm i}(k_x x+k_z z)]}. \end{gather}$$

Note the more general dependence of these trajectories on both $y$ and $t$, allowing the solutions to adopt any sort of temporal function. Substituting these spatio-temporal solutions in the governing equations in wall-normal formulation in (2.14)–(2.15), and solving for the current state $[\tilde {v}(y,t),\tilde {\eta }(y,t)]^{\rm T}$ provides the following definition of the space–time resolvent operator $\mathcal {H}_{\boldsymbol {t}}$:

(2.32)\begin{equation} \begin{pmatrix} \hat{v} \\ \hat{\eta} \end{pmatrix}= \mathcal{H}_{\boldsymbol{t}} \begin{pmatrix} \hat{g}_v \\ \hat{g}_\eta \end{pmatrix}:= \begin{pmatrix} \tilde{\mathcal{H}}_{vv} & 0 \\ -\tilde{\mathcal{H}}_{\eta\eta} ({\rm i}k_z \partial_yU-{\rm i}k_x\partial_yW)\tilde{\mathcal{H}}_{vv} & \tilde{\mathcal{H}}_{\eta\eta} \end{pmatrix} \begin{pmatrix} \hat{g}_v \\ \hat{g}_\eta \end{pmatrix}. \end{equation}

The modified scalar operators are given by

(2.33)$$\begin{gather} \tilde{\mathcal{H}}_{vv}=(\boldsymbol{\mathsf{D}}_t+\hat{\varDelta}^{{-}1}\mathcal{L}_{os})^{{-}1}, \end{gather}$$

(2.34)$$\begin{gather}\tilde{\mathcal{H}}_{\eta\eta}=(\boldsymbol{\mathsf{D}}_t+\mathcal{L}_{sq})^{{-}1}. \end{gather}$$

Here, $\boldsymbol{\mathsf{D}}_t$ represents a generalized discrete time-differentiation operator, and the subscript in $\mathcal {H}_{\boldsymbol {t}}$ represents the triad $\boldsymbol {t}=\lbrace t,k_x,k_z\rbrace$. Note that definitions (2.33) and (2.34) have the symbol $(\tilde {\cdot})$ to emphasise on the temporal dependence and disambiguate from (2.24) and (2.25). As before, resolvent analysis proceeds by taking an SVD of the associated resolvent operator, $\mathcal {H}_t$. The leading resolvent forcing and response modes satisfy the same optimisation problems described in (2.5), (2.7)–(2.8), though now the norm is computed over both space and time, so that

(2.35)\begin{equation} \|\boldsymbol{\psi}(\boldsymbol{x},t) \|_2 = \left(\int_{t\in \varOmega_t}\int_{\boldsymbol{x}\in \varOmega_{\boldsymbol{x}}} \left|\boldsymbol{\psi}(\boldsymbol{x},t)\right|^2 \text{d}\kern 0.7pt\boldsymbol{x}\,\text{d}t\right)^{1/2}, \end{equation}

where $\varOmega _t$ denotes the temporal domain under consideration.

While the theory of this generalization of resolvent analysis is straightforward, it does come with a potential increase in the computational cost. Upon discretization, the size of the matrix representation of the resolvent operator is increased by a factor of the number of time steps, $N_t$, in both the row and column dimensions with respect to the space-only resolvent operator defined in (2.21). For a case where one spatial dimension ($y$) is discretised, this means that the total size is $(2N_y N_t \times 2N_y N_t)$ and each of the block elements is $(N_y N_t \times N_y N_t)$, where $N_y$ and $N_t$ are the number of discretisation points in the space and time dimensions, respectively. Note that the space-only formulation is constituted by a resolvent operator of total size $(2N_y \times 2N_y)$ and with block elements of size $(N_y \times N_y)$. This increase is due to the fact that each of the entries of the operator $\boldsymbol{\mathsf{D}}_t$ corresponds to a temporal instance of a given spatial location in the wall-normal axis. For the purposes of this study, however, this computational cost remains feasible.

Note that in order to disambiguate between the space-only modes and the space–time modes, the symbol $(\hat {\cdot})$ will be used to denote the space-only modes in § 3.

2.4. Sparse resolvent analysis

The theory presented in § 2.1 formulates finding the leading resolvent modes and corresponding gain as an optimisation problem (e.g. (2.5)) in terms of the spatial $l_2$ norm of forcing and response modes (defined in (2.6)). Similarly, the leading resolvent modes for the space–time resolvent formulation described in § 2.3 is formulated using the norms computed over both the spatial and temporal domains (2.35).

Such optimisation problems involving the $l_2$ norm are ubiquitous across a broad range of methods, and arise naturally for methods based on the SVD. However, it is possible to modify such optimisation problems such that their solution has different characteristic features.

Here, we introduce a variant of resolvent analysis that seeks to achieve localisation or sparsity, while also desiring the large energy-amplification levels that are obtained in the standard resolvent formulation. This is achieved by incorporation of the variant of sparse principal component analysis (PCA) described in Hein & Bühler (2010). Similar approaches are discussed in Jolliffe, Trendafilov & Uddin (2003), Zou, Hastie & Tibshirani (2006), Sigg & Buhmann (2008), Journée et al. (2010) and Zou & Xue (2018).

Sparsity is promoted through the incorporation of an additional term in the form of the $l_1$ norm, which produces the following minimization problem:

(2.36)\begin{equation} \boldsymbol{\psi}_1 = \mathop{\rm arg\,min}\limits_{\boldsymbol{\psi}} \frac{(1-\alpha)\|\boldsymbol{\psi}\|_2+ \alpha\|\boldsymbol{\psi}\|_1}{\|\mathcal{H}^*\boldsymbol{\psi}\|_2}. \end{equation}

Here, the sparsity parameter $\alpha \in [0,1]$ determines the number of non-zero elements in $\boldsymbol {\psi }_1$. Moreover, the sparsest solution will be retrieved when $\alpha =1$; while the case where $\alpha =0$ will give the least sparse outcome and will in fact match the result given by (2.5).

Note that the numerator in (2.36) is a convex function. The solution to this optimisation problem is achieved by reformulating it as a nonlinear eigenproblem that enables the use of an inverse power method to find its optimum (Hein & Bühler 2010). In practice, this method often produces solutions with sharp gradients, where some of the entries change drastically from zero to non-zero values. In order to produce coherent structures that resemble observable mechanisms, these solutions are regularised using the resolvent operator while maintaining sparsity. Here, we refer to the response modes computed on the first step as ‘raw’ modes, and use a superscript to disambiguate from the regularised or updated modes. The full collection of steps that produces the sparse leading forcing $\boldsymbol {\phi }_1$ and response $\boldsymbol {\psi }_1$ resolvent modes are described in algorithm .

Algorithm 1 Simplified sparse resolvent analysis

The forcing modes could potentially be updated by substituting the updated response modes in (2.37), although in practice there is not a significant difference between the updated and the first forcing modes. In addition, it is possible to compute higher-order resolvent modes within this framework using the deflation scheme described in Bühler (2014). According to this, the components that have already been identified are removed from the optimisation space before following the steps presented above. Moreover, observe that the method described here promotes sparsity on the response modes $\boldsymbol {\psi }$, although exchanging $\boldsymbol {\psi }$ with $\boldsymbol {\phi }$ and $\mathcal {H}$ with $\mathcal {H}^*$ would yield sparsity promotion in the forcing modes. Lastly, the subscript in the resolvent operator was removed in this section in order to indicate that this methodology can be applied to both the spatial and space–time resolvent operators. For an explicit and expanded form of this algorithm describing how the solution to (2.36) is obtained, refer to Appendix A.

Note that to preserve consistency in the cases presented here, we prescribe the value of the sparsity ratio $\gamma$, instead of the sparsity parameter $\alpha$, as an input to the algorithm. This value is defined as the number of non-zero entries divided by the number of total entries in $\boldsymbol {\psi }$. The relationship between the sparsity parameter $\alpha$ and the sparsity ratio $\gamma$ is further described in Appendix A. The advantage of prescribing an input value of the sparsity ratio $\gamma$ instead of the sparsity parameter $\alpha$ directly is the fact that the algorithm gives control to the user in terms of the desired non-zero elements contained in the computed sparse modes. To some extent, it also allows the user to have control over the number of key features of the flow that are captured by the sparse mode. Note prescribing a set value of $\alpha$ instead requires careful adjustment to achieve the desired sparsity, since it could yield modes with a varying number of non-zero elements.

3. Results

In this section we present the application of the proposed framework on four different systems. After describing numerical details associated with the resolvent calculations in § 3.1, in § 3.2 we showcase the implementation of sparse resolvent analysis on a statistically stationary turbulent channel flow. Here, we consider conventional resolvent analysis with a Fourier transform in time, but enable sparsity promotion in the spanwise direction. We then consider the space–time resolvent operator for this problem in § 3.3, and show that the proposed method can produce temporally localised modes for a statistically stationary flow. We next apply sparse and non-sparse space–time resolvent analysis to two non-statistically stationary systems: a periodic turbulent Stokes boundary layer in § 3.4 and a turbulent channel flow with a sudden lateral pressure gradient in § 3.5.

3.1. Numerical methods

The mean velocity profiles used for resolvent analysis are all obtained from direct numerical simulations (DNS). These simulations use a staggered second-order finite difference scheme (Orlandi 2000), with a fractional step method (Kim & Moin 1985) and third-order Runge–Kutta time-advancing scheme (Wray 1990). Further details regarding the use and validations of these methods, and their application to the specific cases considered here, can be found in Bae et al. (2018, 2019) and Lozano-Durán & Bae (2019).

For resolvent analysis, the wall-normal direction is discretised using a Chebyshev collocation method. In the cases where the spanwise dimension is explicitly discretised, we use a Fourier discretisation scheme with periodic boundary conditions along the spanwise domain. Moreover, if homogeneity is not assumed in the temporal dimension, we adopt a Fourier discretisation scheme when the system is assumed to be statistically stationary or time periodic (i.e. §§ 3.2–3.4), and an explicit Euler finite-differentiation scheme in the temporal dimension with Neumann boundary conditions at the boundaries in § 3.5. The corresponding differentiation operators for both Chebyshev and Fourier discretisations are defined according to the specifications given in Weideman & Reddy (2000). The number of collocation points used for each of the examples considered in this work will be indicated in the corresponding section. In each case, we verify that the numerical resolution gives converged results, with details of these convergence studies given in Appendix B.

In the space–time implementations of the analysis showcased in this paper (i.e. §§ 3.3–3.5), no-slip and no-penetration conditions are enforced at $y=0$ (lower wall), while free-slip and no-penetration conditions are enforced at $y=h$ (the channel centreline). The space–time analysis identifies modes that are sufficiently localised away from the mid-plane using the chosen set of parameters, allowing us to reduce the size of the numerical domain to the area below $y/h=1$.

3.2. Spatially sparse resolvent analysis of turbulent channel flow

Here, we apply the sparse resolvent analysis methodology on a fully developed turbulent channel flow, where we consider spatial (rather than temporal) sparsity. The mean velocity profile is obtained from DNS at a friction Reynolds number of $\textit {Re}_\tau =186$, defined as $\textit {Re}_\tau =u_\tau h/\nu$ and the friction velocity, expressed as $u_\tau =\sqrt {\tau _w/\rho }$. The variables $\nu$ and $\rho$ represent the kinematic viscosity and density of the flow, respectively, the shear stress at the wall boundary is denoted as $\tau _w$, and $h$ is the channel half-height. Hereafter, the superscript $({\cdot })^+$ denotes viscous (inner) units. In this reference, the velocities are scaled by the friction velocity $u_\tau$ and the length scale is given by $\nu /u_\tau$. In order to showcase the application of the sparse formulation of resolvent analysis introduced in § 2.4, we consider two configurations, both of which assume homogeneity in the streamwise direction and the temporal dimension: a first case where the spanwise dimension is assumed to have an infinite extent, followed by a second implementation in which the spanwise dimension is limited to a finite periodic domain.

The first configuration represents a one-dimensional analysis, where the chosen wavelengths in the streamwise $\lambda _x^+=1000$ and spanwise $\lambda _z^+=100$ directions correspond to the average size of the streaks and vortices that arise in the near-wall cycle (Jiménez & Pinelli 1999). The temporal frequency is fixed at $\omega =17.14$, where the temporal scale is normalised by the friction velocity ($u_\tau$) and the half-height of the channel ($h$). The frequency is chosen to be that which yields the maximum resolvent gain at these wavelengths, which will best enable direct comparison with the following space–time resolvent results. Note that this frequency gives a critical layer at a location $y^+\approx 40$ inner units from the wall, which is further from the wall than typical near-wall streaks (which are characteristically found at $y^+\approx 15$). To avoid any misinterpretation, henceforth we will generally refer to the resulting structures coming from resolvent response modes (which tend to have amplitude peaks near the critical layer) simply as streamwise (rather than near-wall) streaks and vortices. Note also that the location of these structures could be moved closer to the wall either by reducing the choice of $\omega$ (which we find yields qualitatively similar results) or by introducing an eddy viscosity term to the resolvent formulation (Symon et al. 2023). In this application, the number of collocation points in the wall-normal direction is $N_y=101$. No-slip and no-penetration conditions are enforced for the velocity at the upper and lower boundaries.

Results obtained from applying both standard and sparsity-promoting (in the wall-normal direction) resolvent analysis are shown in figure 1, where we show the amplitude of the streamwise velocity component of the leading two resolvent response modes. For standard resolvent analysis, the chosen wavenumbers and frequency produces modes with two peaks, each centred near one of the two critical layer locations (where $U(y) = \omega /k_x$). The leading two singular values are essentially identical and represent the fact that the two peaks are separated from each other, and can each represent structures that have an arbitrary phase shift between them. Mathematically, these modes are basis vectors for a subspace of dimension two, and thus, an equally valid choice of basis for this subspace would be given by localised modes with one peak each.

Figure 1. Streamwise components ($u$) of the first (a,b) and second (c,d) standard and sparse resolvent response modes computed for a turbulent channel flow with $\textit {Re}_\tau =186$, $\lambda _x^+=1000$, $\lambda _z^+=100$ and temporal frequency $\omega =17.14$. Note that the standard modes (black) have been computed using (2.8), while (2.36) produced the raw sparse modes $\hat {\boldsymbol {\psi }}^{raw}$ (a,c) from which we extract the adjusted modes $\hat {\boldsymbol {\psi }}$ (b,d). The streamwise mean flow $U$ is added to all plots for reference. The black dashed lines denote critical layer locations.

To compare to the results of sparse resolvent analysis, several solutions are computed for different sparsity ratios $\gamma$ via the optimisation problem posed in (2.36). The obtained raw sparse resolvent modes shown in figures 1(a) and 1(c) seem to be highly dependent on the corresponding value of $\gamma$, although they are again all located near one of the two critical layer locations. This dependence on the sparsity parameter subsides after regularising the modes following algorithm , as shown in figures 1(b) and 1(d). We observe that these regularised modes each recover one of the two peaks identified by regular resolvent analysis, indicating that the sparse variant is finding basis elements for the leading resolvent subspace that are spatially sparse. The sparse singular values (computed via (2.39)) are also consistent with the standard resolvent case, with $\sigma _1$ being about 1 % lower for the sparse case with $\gamma = 0.1$ ($\alpha = 0.0486$) in comparison with the non-sparse equivalent.

As the regularisation step appears to both remove the dependence on the sparsity parameter and recover the results for regular resolvent analysis for this example, we will exclusively use this regularised version for the remaining cases.

We now consider the same turbulent channel flow, but instead of assuming a Fourier decomposition in the spanwise direction, we instead explicitly discretise this dimension, applying periodic boundary conditions with a spanwise extent $L_z$ twice the channel height (that is, $L_z/h=4$). We use a Fourier basis in this spanwise dimension, with $N_z=92$ collocation points (with $N_y=101$ as before). We keep the same frequency and streamwise wavelength that was used in the one-dimensional analysis. This configuration is motivated by the fact that structures and correlations that are observed in wall-bounded turbulent flows typically exist only over a finite spanwise extent (Kim, Moin & Moser 1987; Hutchins & Marusic 2007; Dennis & Nickels 2011; Sillero, Jimenez & Moser 2014; Jiménez 2018). While such localised structures could be represented by an appropriate combination of spanwise Fourier modes, this would be less efficient than a single-mode representation.

Figure 2 shows the leading resolvent forcing and response modes obtained from standard and sparse resolvent analyses, visualised in the $y\unicode{x2013}z$ plane. Note that the contours represent the streamwise component ($u$) of the modes, which correspond to streamwise streaks of fast- and slow-moving regions. The vector fields represent the wall-normal ($v$) and spanwise ($w$) velocity components of the modes, which here form streamwise-aligned vortical structures. Since the system is homogeneous in the spanwise dimension, the standard resolvent modes (figure 2a,c) give Fourier modes in this direction. The response mode consists of alternating slow- and fast-moving streamwise streaks, with streamwise vortices located between each streak. The wall-normal profile of these modes is close to matching those from the one-dimensional analysis shown in figure 1, where the identified spanwise wavelength for the leading mode is slightly different from that selected in the one-dimensional analysis shown in figure 1. The configuration of these streamwise streaks is consistent with the lift-up mechanism (Landahl 1975, 1980), through which streamwise vortices lead to the formation of streamwise streaks by transporting slow-moving fluid away from the wall and vice-versa.

Figure 2. Real component of leading response modes $\hat {\boldsymbol {\psi }}_1$ (a,b) and forcing modes $\hat {\boldsymbol {\phi }}_1$ (c,d), for standard (a,c) and sparse (b,d) resolvent analysis with $\gamma = 0.001$ ($\alpha = 0.953$) resolvent analysis, applied to a turbulent channel flow at $\textit {Re}_\tau = 186$, with a finite periodic domain in the spanwise direction defined by $z/h \in \lbrace -2,2\rbrace$ and a frequency $\omega =17.14$. Contours of streamwise velocity are shown, with arrows indicating the velocity in the spanwise ($z$) and wall-normal ($y$) directions. Green dashed lines indicate critical layer locations. The streamwise wavelength is $\lambda _x^+$ = 1000.

The sparse resolvent modes shown in figures 2(b) and 2(d) contain a spatially localised unit of the periodic structures identified by standard resolvent analysis in figures 2(a) and 2(c). In particular, the sparse response mode consists of a primary central vortex surrounded by a fast and slow streak, flanked by lower-amplitude secondary vortices and streaks. Note that while we are only showing the lower half of the domain, on the upper half the same structure is present for standard resolvent analysis, but not for the sparse equivalent (again consistent with figure 1) with $\gamma = 0.001$ ($\alpha = 0.953$). While not shown here, higher-order sparse resolvent modes consist of repetitions of this localised structure both on the upper wall and translated in the spanwise direction (with the spanwise location of the leading mode being arbitrary). Note also that the relative phase of the forcing and response modes is consistent between the regular and sparse modes. In terms of the energy content of these structures, the first non-sparse singular value is about 1.178 times larger than its sparse counterpart.

To further explore the relationship between the sparse and non-sparse modes for this configuration, we show in figure 3 modes computed using the standard version of resolvent analysis, but the region where amplification is measured is restricted in the spatial domain along the $z$ axis (while maintaining a domain of infinite extent in the streamwise direction). This is achieved by modifying the weight function associated with the resolvent operator to only consider amplification at spatial locations within the indicated regions. Note that in this case, the spanwise domain is defined over a larger region $-4 \leqslant z/h \leqslant +4$ to further emphasise this localisation. As in the previous results for both sparse and non-sparse analyses, the identified modes consist of alternating slow- and fast-moving streamwise streaks, with streamwise vortices between each pair of streaks. The size of these structures is approximately the same once the window reaches a width $L_z \gtrapprox h$, with the number of such structures present dependent on the width of the window. Restricting the spanwise domain also reduces energy amplification associated with these structures, in comparison with the Fourier mode extending across the entire domain width.

Figure 3. Real component of leading response modes $\hat {\boldsymbol {\psi }}_1$ for standard resolvent analysis, applied to a turbulent channel flow at $\textit {Re}_\tau = 186$, with a finite periodic domain in the spanwise direction restricted to $L_z/h \in \{0.42,0.83,1.17,1.58,2 \}$ with $\omega =17.14$ and $\lambda _x^+$ = 1000. Contours of streamwise velocity are shown, with arrows indicating the velocity in the spanwise ($z$) and wall-normal ($y$) directions.

Figure 4 shows the leading singular value of the resolvent modes as a function of the mode width $(L_z/h)$, along with the singular values associated with the sparse and non-sparse modes depicted in figure 2 (blue dotted line). As the mode width increases, we observe that the singular value increases rapidly for small $L_z$, and then becomes more gradual as the singular value approaches that obtained using the full domain. When $L_z = 0.83h$ (figure 3b), we observe a mode structure similar to that found for the sparse mode, with a pair of central positive and negative streamwise streaks and a central streamwise vortex sitting between the streak pair. In figure 4 it is observed that this case is located near the elbow of the amplification against the mode width curve, with a singular value approximately matching that obtained for the sparse mode (blue dotted line). This suggests that the sparse mode is identifying a structure that finds an appropriate trade-off between amplification and sparsity, by identifying a structure that is much more localised that the original spanwise Fourier mode, while still having a comparable singular value. To show the effect that the choice of sparsity parameter has on this trade-off between amplification and localisation, we also show in figure 4 (green dotted line) the sparse singular value for a larger sparsity ratio of $\gamma =0.0078$ ($\alpha =0.625$). This amplification level intersects the $\sigma _1$ vs mode width curve at a larger mode width. The associated sparse resolvent response mode, while not plotted, was found to be similar to the configuration of streamwise streaks and vortices identified in figure 3(d), which is near this intersection point.

Figure 4. Leading singular values of the resolvent operator for a restricted spanwise domain of length $0.08 \leqslant L_z/h \leqslant 8$. The values corresponding to the modes compiled in figure 3 are indicated with red markers. Note that the leading singular value of the non-sparse and sparse resolvent modes depicted in figure 2 are also indicated with dashed black and blue dotted lines, respectively, for reference. In addition, the leading singular value of the sparse mode retrieved for $\gamma =0.0078$ ($\alpha =0.625$) is indicated with a green dotted line.

3.3. Space–time resolvent analysis of a turbulent channel flow

In this section we demonstrate the implementation of the non-sparse and sparse space–time variants of resolvent analysis on the statistically stationary turbulent channel flow considered in § 3.2. Here, we do not perform a Fourier transform in the temporal dimension but instead implement the space–time formulation of resolvent analysis introduced in (2.32). All the parameters adopted in the analysis conducted in § 3.2 are also adopted here (except for the frequency, which is no longer specified), and time is non-dimensionalised by the friction velocity, $u_\tau$, and the half-height of the channel, $h$.

First, the streamwise ($u$) and wall-normal velocity ($v$) components of the resolvent modes obtained from the non-sparse space–time analysis are shown in figure 5 with a time horizon $\tau _{tot}=20\tau$ that spans over 20 cycles of length $\tau =2{\rm \pi} /\omega =0.3659$, where $\omega = 17.14$ is the same frequency used in §§ 3.2–3.3. The number of collocation points is $N_y=101$ in the wall-normal direction and $N_t=551$ in the temporal dimension. The location of maximum mode amplitude along the $y$ axis is indicated with a horizontal dotted line, and the value is indicated for each subplot with the symbol with $y_0/h$ and its equivalent in inner units $y_0^+$ for reference.

Figure 5. Real component of space–time leading response $\boldsymbol {\psi }_1$ (a,b) and forcing $\boldsymbol {\phi }_1$ (c,d) modes of the resolvent operator in a turbulent channel flow with $\textit {Re}_\tau =186$, and $\lambda _x^+=1000$, $\lambda _z^+=100$ for a time horizon $\tau _{tot}=20\tau =20(2{\rm \pi} /\omega )$. Horizontal dotted lines in (a–d) indicate the $y$ location where each mode component achieves its maximum amplitude ($y_0/h$). Panels (a,c) show the streamwise velocity component ($u$), while (b,d) give the wall-normal component ($v$).

As expected, we observe that the leading space–time resolvent modes are Fourier modes in time, with a frequency that matches that used in the previous analyses (i.e. the mode exhibits twenty periods over the time domain of length $20\tau$).

To further demonstrate that the space–time variant gives equivalent results to standard resolvent analysis, in figure 6 we compare the leading 30 singular values from both versions, where in the standard resolvent analysis we compile the results across all frequencies that are permissible when using this temporal domain (i.e. $\omega =2{\rm \pi} n/\tau _{tot}$ with $n \in \mathbb {Z}$).

Figure 6. Leading $30$ singular values of the space–time resolvent (black) and singular values of the space (standard) resolvent computed at $\omega _n=n(2{\rm \pi} /\tau _{tot})$ with $-100\leqslant n \leqslant 100$ sorted in descending order (red) for a turbulent channel flow with $\textit {Re}_\tau =186$, $\lambda _x^+=1000$, $\lambda _z^+=100$ with a time horizon (for the space–time variant) $\tau _{tot}=20\tau =20(2{\rm \pi} /\omega )$, where $\omega =17.14$.

We now consider the spatio-temporally sparse variant of this analysis. The total time horizon is again set to $\tau _{tot}=20\tau$, which allows for the potential growth and decay of temporally localised modes without the influence of the periodic boundary conditions. Here $N_y=101$ collocation points are used in the wall-normal direction and $N_t=501$ in the time dimension. Figure 7 contains the streamwise and wall-normal components of the updated response and forcing modes with a sparsity ratio $\gamma =0.001$ ($\alpha = 0.883$). The location of the maximum mode amplitude along the $y$ axis is again indicated with a horizontal dotted line. Notice how the analysis identifies structures that are sparse both in the spatial and temporal dimensions. Observe that the forcing modes in figures 7(c) and 7(d) precede the response in figures 7(a) and 7(b). The phase variation in time (corresponding to the temporal frequency) and wall-normal location of the modes approximately match those identified from the leading non-sparse space–time resolvent mode shown in figure 5. The sparse mode, however, identifies aspects of the physics that cannot be directly observed without sparsity promotion (and thus localisation), such as the change in inclination angle and wall-normal location of the mode components with time. In particular, the inclination angle of structures tend to lean further backwards in the $y\unicode{x2013}t$ plane as time increases, which would correspond to an increased downstream tilt over time in the $x\unicode{x2013}y$ plane, consistent with the Orr mechanism (Orr 1907).

Figure 7. Real component of space–time sparse leading response $\boldsymbol {\psi }_1$ (a,b) and forcing $\boldsymbol {\phi }_1$ (c,d) modes of the resolvent operator in a turbulent channel flow with $\textit {Re}_\tau =186$, and $\lambda _x^+=1000$, $\lambda _z^+=100$, computed with a sparsity parameter $\gamma =0.001$ ($\alpha = 0.094$), for a time horizon $\tau _{tot}=20\tau =20(2{\rm \pi} /\omega )$. Horizontal dotted lines in (a–d) indicate the $y$ location where each mode component achieves its maximum amplitude ($y_0/h$). Panels (e,f) show contour levels of the absolute value of leading forcing (red) and response (blue) modes.

We also note that the temporal locations of these time-localised resolvent mode components are consistent with the inherent causal structure of the time domain, in which the forcing precedes the response in the same way that an input to a dynamical system precedes an output. This temporal structure is also consistent with known physical mechanisms as well, with forcing in $v$ proceeding the response in $v$ (figure 7f) that in turn proceeds the (energetically dominant) response in $u$ (figure 7e), consistent with the lift-up mechanism. In terms of the energy content, for a fixed time horizon of $20\tau$, the leading non-sparse singular value is larger than the leading sparse one by a factor of 5.60.

To explore the temporal profiles of the components of these sparse space–time modes, figure 8 shows the cross-sections along the $t$ axis of the sparse modes in figure 7 at the locations in $y$ where the amplitude of the signal finds its maximum. The temporal evolution of the forcing mode components seems to be more localised in these cross-sections, as they abruptly drop to zero near where the amplitude of the response modes finds its maximum value. The locations of the temporal peaks of these mode components is consistent with the discussion above related to figure 7, with the peak amplitude $\phi _v$ preceding that of $\psi _v$, which itself comes prior to the maximum in $\psi _u$. Quantitatively, the time duration between the peaks is of the same order of magnitude as $\tau$, with a time lag of approximately $0.41\tau$ between the peaks in the $v$ components of the forcing and response modes, and a further time lag of $0.68\tau$ between the peak amplitude of the responses in $v$ and $u$. This suggests that the energy transfer between these mode components occurs over a single oscillatory cycle. Note also that the phase lag in between the $\psi _v$ and $\psi _u$ peaks approximately matches the phase difference between the oscillating component of these signals, as both are entirely real near their peak amplitudes. The shape of the streamwise and wall-normal components of the modes (first row in figure 8) displays a Gaussian envelope with a nearly constant phase gradient and a uniform phase variation across all mode components. This could potentially suggest a natural choice of wavelet basis functions for a low-order temporally localised analysis of this flow (Ballouz et al. 2023, 2024; Madhusudanan & Kerswell 2024). Moreover, the particular shape of the response modes could enable the implementation of wavepacket pseudomode theory to approximate these response modes (Obrist & Schmid 2010, 2011; Mao & Sherwin 2011; Edstrand et al. 2018; Dawson & McKeon 2019, 2020).

Figure 8. Absolute (solid) and real (dashed) components of the temporal evolution of the sparse modes shown in figures 7(a)–7(d), at the $y$ locations of maximum amplitude ($y_0$).

3.4. Space-time resolvent analysis of a turbulent Stokes boundary layer

In this section we apply non-sparse and sparse space–time resolvent analysis to a turbulent Stokes boundary layer, with a time-periodic mean velocity field as shown in figure 9. This flow sits between two plates that oscillate synchronously with a velocity given by

(3.1)\begin{equation} U_w(t) = U_0\cos{(\varOmega t)}, \end{equation}

with $U_w$ denoting the velocity of the oscillating walls and $U_0$ representing a constant. There is no external pressure gradient, so the flow is entirely driven by the motion of the walls. The Reynolds number for this flow is defined as a function of the frequency of the oscillations $\varOmega$, such that $\textit {Re}_\varOmega =U_0\delta _\varOmega /\nu$, and $\delta _\varOmega =\sqrt {2\nu /\varOmega }$ is a length parameter denoting the boundary layer thickness. For the ensuing analysis, we choose $\textit {Re}_\varOmega =1500$ to ensure that the flow lies in the intermittently turbulent regime (Hino, Sawamoto & Takasu 1976; Akhavan, Kamm & Shapiro 1991; Verzicco & Vittori 1996; Vittori & Verzicco 1998; Costamagna, Vittori & Blondeaux 2003). Note that while there appears to be little prior work performing linear analyses of this configuration in the turbulent regime, linear analysis of laminar pulsatile channel flow has been considered using Floquet analysis (Pier & Schmid 2017) and optimally time-dependent modes (Kern et al. 2021).

Figure 9. (a) Streamwise mean flow $U(y,t)$ (colourmap) and streamwise turbulence intensity $u_{rms}$ (grey dashed) of a turbulent Stokes boundary layer, with levels indicating $\{20\,\%, 30\,\%, 50\,\%, 75\,\%, 95\,\%\}$ of its maximum value, over a full oscillating cycle; (b) time instances of the turbulent mean at $t\varOmega ={\rm \pi}$ (black) and $t\varOmega =3{\rm \pi} /2$ (red).

The mean velocity profile and the root mean square of the fluctuating turbulent velocity components are obtained through DNS as described in § 3.1. The size of the domain for the DNS is $6{\rm \pi} \delta _\varOmega \times 80\delta _\varOmega \times 3{\rm \pi} \delta _\varOmega$ (in the streamwise, wall-normal and spanwise directions, respectively) and contains 64, 385 and 64 grid points in each of the corresponding dimensions. To collect mean data, the DNS was run for 100 eddy turnover times after the decay of transient startup effects, where here we define this time scale as $\delta _\varOmega /u_\tau$, with $\delta _\varOmega =0.025h$. The time-periodic mean velocity profile is obtained through averaging in the dimensions of spatial homogeneity (i.e. the streamwise and spanwise directions) for each phase of wall motion. For resolvent analysis, the size of the numerical domain is $[0,h]\times [0,\tau _{tot}]$ with $h=1$ and $\tau _{tot}$ being the length of the time horizon.

For performing both sparse and non-sparse space–time resolvent analysis, we consider wavelengths corresponding to the extent of the DNS computational domain, i.e. $\lambda _x = 6{\rm \pi} \delta _\varOmega \approx 0.471 h$ and $\lambda _z = 3{\rm \pi} \delta _\varOmega \approx 0.236 h$. That is, we look at the largest three-dimensional structure that the DNS would be capable of resolving. The reason for focusing on such large structures in $x$ and $z$ is that we expect that they will correspond to space–time resolvent modes that also have a large extent in $t$ and $y$, allowing for the distinction between the sparse and non-sparse variants to be clearly observed.

For the non-sparse implementation of resolvent analysis, the domain is discretised using $N_y=105$ collocation points in the wall-normal axis and $N_t=651$ collocation points in the temporal dimension over a time window comprised of one oscillating cycle. The streamwise and wall-normal velocity components of the leading forcing and response modes obtained from this analysis are shown in figure 10, along with the corresponding streamwise and wall-normal turbulence intensity computed from the DNS. Observe that the identified coherent structures do not correspond to a single temporal Fourier mode. Instead, they are concentrated over approximately half of the full oscillation period. Although not shown here, the second leading resolvent mode depicts a coherent structure of almost identical energy content that is similar to the leading modes in figure 10, with a relative temporal shift of a half-period.

Figure 10. Real component of leading space–time resolvent (a,b) response $\boldsymbol {\psi }_1$ and (c,d) forcing $\boldsymbol {\phi }_1$ modes of the resolvent operator in a turbulent Stokes boundary layer with $\lambda _x/h=0.471$ and $\lambda _z/h=0.236$ during a full oscillating cycle. The grey dashed lines represent the (a,c) streamwise and (b,d) wall-normal turbulence intensity, with contours at $\{15\,\%,20\,\%30\,\%,40\,\%,60\,\%\}$ and $\{22\,\%,30\,\%,45\,\%,55\,\%,70\,\%\}$ of the maximum value of $u_{rms}$ and $v_{rms}$, respectively.

The location of the modes along the $y$ axis shifts away from the wall as the mean flow evolves over time. This angle in the $y\unicode{x2013}t$ plane appears to be similar to that observed for the turbulence intensities, suggesting that the mode depicts an energy transfer mechanism away from the wall that is present in the full nonlinear system. These modes are not concentrated in the near-wall region where the mean shear and turbulence intensity is largest. In this region the turbulent energy content is likely dominated by structures at smaller length scales than the wavelengths considered here. We do observe, however, in figure 10(a) that the streamwise velocity response is amplified as the mode enters a region of higher streamwise turbulence intensity at $(t\varOmega,y/h)\approx ({\rm \pi},0.2)$. In addition, the disturbance streamfunctions identified in Blennerhassett & Bassom (2002) for a flat Stokes boundary layer are similar to the behaviour of the space–time modes presented in figure 10 in two ways: first, the spatial shift of the structures along the wall-normal direction; second, the time-dependent oscillation frequency of these structures. The latter property will be addressed in more depth later in this section.

For the sparse implementation of space–time resolvent analysis at the same wavelengths, we adopt the same time horizon as in the non-sparse case, and a numerical domain with $N_y=101$ collocation points in the wall-normal axis and $N_t=631$ collocation points in the temporal dimension, with a sparsity ratio $\gamma =0.001$ ($\alpha = 0.625$). The location of maximum mode amplitude along the $y$ axis is indicated with a horizontal dotted line, and the value is indicated for each subplot with the symbol with $y_0/h$ and its equivalent in inner units $y_0^+$ for reference.

The leading sparse resolvent modes are shown in figure 11, again overlaid with contour levels of the streamwise and wall-normal turbulence intensity. We observe that at these streamwise and spanwise wavelengths, the structure of the leading resolvent modes differs substantially between the non-sparse and sparse analyses. In figure 11 the sparse analysis identifies localised structures for which the streamwise and wall-normal components are concentrated within one quarter of the total temporal domain, and within a small spatial region near $y/h=0.2$ (note that the vertical range shown in figure 11 is smaller than in figure 10). Moreover, both forcing mode components precede their corresponding response modes (shown directly in figure 11e,f), with a peak response in $v$ preceding that for $u$, which was also observed in the sparse modes computed for statistically stationary channel flow in § 3.3. Comparing the sparse modes with the non-sparse counterparts in figure 10, the sparse analysis favours structures that are closer to the peak of the turbulence intensity.

Figure 11. Real component of sparse space–time leading (a,b) response $\boldsymbol {\psi }_1$ and (c,d) forcing $\boldsymbol {\phi }_1$ modes of the resolvent operator for a turbulent Stokes boundary layer with $\lambda _x/h=0.471$ and $\lambda _z/h=0.236$, computed with a sparsity ratio $\gamma =0.001$, across a full oscillation cycle. Panels (e,f) show contour levels of the absolute value of leading forcing (red) and response (blue) modes. The grey dashed lines represent the corresponding streamwise and wall-normal turbulence intensity, with contour levels at $\{23\,\%,32\,\%40\,\%,50\,\%,70\,\%\}$ and $\{35\,\%,45\,\%,60\,\%,75\,\%\}$ of the maximum value of $u_{rms}$ and $v_{rms}$, respectively. Horizontal dotted lines in (a–d) indicate the $y$ location where each mode component achieves its maximum amplitude ($y_0/h$).

Figure 12 shows cross-sections of these sparse modes along the temporal axis at the $y$ locations where they achieve their maximum amplitude (see the dotted lines indicated in figure 11). Comparing figures 7 and 12, we observe greater variation in the phase both within and across each mode component for the time-varying Stokes boundary layer configuration. For example, the streamwise velocity response in figure 12(a) shows that the rate of phase variation increases over time. The phase variation in the streamwise component of the sparse forcing mode is much slower, which figure 11(c) shows is due to the almost horizontal inclination of this component. As was the case in figure 7, we observe that the forcing mode components have envelopes that are less symmetric in time, abruptly decaying to zero at a given cutoff time.

Figure 12. Absolute (solid) and real (dashed) components of the temporal evolution of the sparse modes shown in figures 11(a)–11(d), at the $y$ locations of maximum amplitude ($y_0$).

Figure 13 presents the spectral power of the non-sparse (a,b) and sparse (c,d) response modes shown in figures 10 and 11, respectively. Each plot is generated via a temporal Fourier transform at each spatial location along the wall-normal axis ($y$). On the one hand, the regions of higher spectral power corresponding to the sparse modes in figures 13(c) and 13(f) are localised in $y$, but consist of a relatively broad frequency range. Conversely, the regions corresponding to the location of maximum spectral power of the non-sparse modes do not show the same degree of localisation in $y$, and also exhibit several distinct peaks in the ($\omega, y/h$) plane. Furthermore, both the non-sparse and sparse spatio-temporal analyses identify continuous regions of active frequencies. Indeed, the identified modes in figures 10 and 11 depict structures with time-varying frequencies. This evolution of the frequency within a single space–time mode would not be directly captured via harmonic resolvent analysis, where each component mode is associated with a single (discrete) temporal frequency (though triadic interactions between different frequencies could be isolated with this approach). Instead, the superposition of a considerable amount of temporal frequencies would be required to capture the dominant spatio-temporal modes identified by the direct space–time analysis. This suggests that the present approach could be more insightful when applied to flows that do not necessarily possess sparse frequency content, such as the configuration considered in this section.

Figure 13. Spectral power of the streamwise ($u$) (a,c) and wall-normal ($v$) (b,d) components of the non-sparse (a,b) and sparse (c,d) spatio-temporal resolvent modes presented in figures 10 and 11, respectively. A discrete Fourier transform is computed during one oscillating period at each of the locations of the $y$ axis in the numerical domain.

3.5. Space–time resolvent analysis of a turbulent channel flow with sudden lateral pressure gradient

The last case considered in this work is a fully developed turbulent channel flow at $\textit {Re}_\tau =186$ that is subjected to a sudden lateral pressure gradient at $t=0$. The spanwise pressure gradient is related to the fixed streamwise gradient by $\partial p/\partial z=\varPi (\partial p/\partial x)$, where here we use $\varPi =30$. Additional information about this configuration can be found in Moin et al. (1990) and Lozano-Durán et al. (2020). The turbulent mean flow therefore contains both streamwise $U(y,t)$ and spanwise $W(y,t)$ non-zero components, which evolve over time until the system reaches its new statistically stationary state. These velocity components are shown in figure 14.

Figure 14. Temporal evolution of the streamwise $U(t,y)$ (black) and spanwise $W(t,y)$ (red) mean velocity profiles obtained by DNS of a turbulent channel flow at $\textit {Re}_\tau =186$ with a sudden lateral pressure gradient at $t=0$ with $\partial p/\partial z=\varPi (\partial p/\partial x)$ and $\varPi =30$ at $tu_\tau /h\in \lbrace 0,0.58,1.17,1.38,1.76,2.34\rbrace$ (Lozano-Durán et al. 2020).

To form the resolvent operator, the temporal domain is implicitly non-dimensionalised by the initial friction velocity $u_{\tau,0}$ at $t=0$ and $h$, and it is determined in two steps to achieve numerical accuracy. It is first defined as the interval $[-0.58,2.34]h/u_\tau$ to extend the numerical domain before the initial condition to include time prior to the application of the spanwise pressure gradient. Thus, for $t<0$, the base flow is constant, with $U(t<0,y)=U(T=0,y)$ and $W(t<0,y)=0$, and $\partial p/\partial z(t<0)=0$. After a first analysis is conducted, the time domain is restricted to a smaller temporal window near the observed temporal location of the leading modes that is large enough to produce unaltered results with a larger temporal resolution. Note that here we implement an explicit Euler finite-differentiation scheme in the temporal dimension. Thus, the final numerical domain used in this analysis is given by $[0,h]\times [0.4,2]h/u_\tau$.

As was the case for the statistically stationary channel flow studied in §§ 3.2–3.3, we concentrate our analysis on streamwise and spanwise wavenumbers corresponding to the typical size of near-wall streaks. After the pressure gradient is applied, the magnitude and direction of the mean flow changes, meaning that the wavelengths associated with near-wall streaks is also expected to change. Assuming that near-wall streaks adjust their size and orientation with the mean, we can select wavelengths corresponding to the typical size of near-wall streaks at a certain point in time (Ballouz et al. 2024) (though as before, we are not guaranteed that the wall-normal location of the resulting modes will match those typical of such streaks). Here, we select $\lambda ^+_x=189$ and $\lambda ^+_z=1890$, which correspond to the typical length of the near-wall streaks at a time instance $t u_\tau /h=1.38$.

For the non-sparse space–time resolvent analysis, the domain is discretised using $N_y=132$ and $N_t=500$ collocation points. The location of the maximum mode amplitude along the $y$ axis is indicated with a horizontal dotted line, and the value is indicated for each subplot with the symbol with $y_0/h$ and its equivalent in inner units $y_0^+$ for reference. Leading resolvent mode components for this analysis are shown in figure 15. We observe that the modes are localised in time, centred near to the nominal time ($t u_\tau /h=1.38$) where the streamwise and spanwise length scales of the modes coincide with near-wall streaks. The structure of the modes share some similarities with those observed for both the non-sparse and sparse analyses of stationary channel flow (figures 5 and 7). In all cases, the wall-normal response consists of vertically aligned structures, while the forcing is inclined upstream in the $y\unicode{x2013}t$ plane. The $u$ component of the response appears to incline further backwards as time progresses in both figures 7(a) and 15(a), though the latter does appear to move towards and move away from the wall as observed in the former. Note that the $u$ component plotted in figure 15 no longer entirely corresponds to the streamwise direction, due to the spanwise pressure gradient.

Figure 15. Real component of space–time $u$ and $v$ velocity components of the leading (a,b) response $\boldsymbol {\psi }_1$ and (c,d) forcing $\boldsymbol {\phi }_1$ modes of the space–time resolvent operator for a channel flow with a sudden lateral pressure gradient with $\textit {Re}_\tau =186$, and $\lambda ^+_x=189$, $\lambda ^+_z =1890$. (e,f) Contour levels of the absolute value of the leading forcing (red) and response (blue) modes. Horizontal dotted lines in (a–d) indicate the $y$ location where each mode component achieves its maximum amplitude ($y_0/h$). The grey dashed lines represent the boundaries of the temporal numerical domain.

Figures 15(e)–15(f) show that the forcing and response mode components are all located over approximately the same time interval for this case. This is different from the sparse analysis of statistically stationary channel flow (see figure 7e,f), where the forcing mode components decayed before their responses did, and the wall-normal forcing and response components started before the streamwise components. These differences can likely be attributed to the differences in how localised modes are obtained; for the three-dimensional channel configuration we are not explicitly promoting sparse modes, but rather the modes are aligning with a region of time corresponding to a mean flow with maximum linear energy amplification for structures at the specified spatial scales. This can be further explored by applying the sparse variant to this configuration.

For the sparse implementation of space–time resolvent analysis, we adopt the same wavelengths and a sparsity ratio $\gamma =0.001$ ($\alpha = 0.781$). The domain is discretised using $N_y=130$ and $N_t=500$ collocation points. The leading sparse modes are shown in figure 16, where both the $u$ and $v$ components display a similar structure, but with more temporal localisation, to those shown for the non-sparse analysis in figure 15. The location of the maximum mode amplitude along the $y$ axis is indicated with a horizontal dotted line, and the value is indicated for each subplot with the symbol with $y_0/h$ and its equivalent in inner units $y_0^+$ for reference.

Figure 16. Real component of sparse space–time $u$ and $v$ velocity components of the leading (a,b) response $\boldsymbol {\psi }_1$ and (c,d) forcing $\boldsymbol {\phi }_1$ modes of the space–time resolvent for a channel flow with a sudden lateral pressure gradient with $\textit {Re}_\tau =186$, and $\lambda ^+_x=189$, $\lambda ^+_z =1890$ with a sparsity ratio $\gamma =0.001$ ($\alpha = 0.781$). (e,f) Contour levels of the absolute value of the leading forcing (red) and response (blue) modes. Horizontal dotted lines in (a–d) indicate the $y$ location where each mode component achieves its maximum amplitude ($y_0/h$). The grey dashed lines represent the boundaries of the temporal numerical domain.

The sparse analysis thus appears to identify the same mechanism, localised around the same region in time. Unlike the non-sparse version, here we observe in figures 16(e)–16(f) differences in the temporal footprints of each mode component, with the decay in amplitude of both components of the forcing preceding the decay of the response.

The cross-sections along the time axis of the modes contained in figures 15 and 16 are shown in figure 17 in red and black, respectively. The amplitudes of both the non-sparse and sparse modes adopt qualitatively similar envelopes, with the sparse variant being narrower, and centred at a slightly later time. We emphasise that unlike the profiles shown in figure 8, here the time location of these modes is not arbitrary and corresponds to a region of time where the mean profile enables the largest linear amplification. The phase variation is similar across all mode components, indicating a characteristic frequency maximising amplification. Here, we thus find that both sparse and non-sparse space–time resolvent analyses appear to identify the same amplification mechanism, with the addition of the $l_1$-norm term in the objective function restricting the sparse modes to a smaller temporal window.

Figure 17. Absolute (solid) and real (dashed) components of the temporal evolution of the non-sparse modes shown in figures 15(a)–15(d) (red) and sparse modes shown in figures 16(a)–16(d) (black), at the $y$ locations of maximum amplitude.

Lastly, to demonstrate the numerical convergence of the results presented in this section, figure 18 presents the first 10 singular values of the non-sparse and sparse space–time resolvent operators with $\gamma \in \lbrace 0.001,0.15 \rbrace$ space–time resolvent operators against the leading singular values obtained in the study described in Ballouz et al. (2024), where the numerical methods utilize a wavelet transform in time. The corresponding leading modes coincide with figure 15 and figure 16 in the non-sparse and sparse with $\gamma = 0.001$ ($\alpha = 0.781$) implementation of the space–time analysis, respectively. First, observe that the singular values of the non-sparse analysis closely match the reference values. Second, a sparsity ratio of $\gamma =0.15$ ($\alpha = 0.004$) is sufficiently large to obtain similar amplification levels to these reference values (note that the corresponding modes are not shown in this paper). This agreement also confirms the convergence of the sparse analysis. However, the singular values retrieved from the sparse analysis with $\gamma = 0.001$ ($\alpha = 0.781$) and $\gamma = 0.005$ ($\alpha = 0.258$) retrieve only a portion of the reference values. In particular, the leading reference value is larger than the leading sparse singular value with $\gamma = 0.001$ by a factor of 2.0958, and we observe that this ratio decreases for higher-order singular values. The fact that the sparse and non-sparse analyses yield similar results for this case is likely due to the highly transient nature of the flow, where the leading resolvent modes are naturally localised even without explicitly promoting sparsity.

Figure 18. Leading $10$ singular values of the non-sparse and sparse (with $\gamma \in \lbrace 0.001,0.15 \rbrace$) space–time resolvent operator for a channel flow with a sudden lateral pressure gradient with $\textit {Re}_\tau =186$, and $\lambda ^+_x=189$, $\lambda ^+_z =1890$. The reference values correspond to the analysis described in Ballouz et al. (2024).

4. Discussion and conclusions

In this work we have proposed a sparse space–time variant of resolvent analysis that can identify time-localised coherent structures. These spatio-temporal structures correspond to the inputs and outputs that optimise an objective function that promotes large linear energy amplification, while also promoting the localisation in space and time of these structures. Localisation is achieved through the addition of a sparsity-promoting $l_1$-norm term to the standard optimisation problem used for resolvent-type analyses. The new optimisation problem takes the form of a nonlinear eigenproblem, for which the optimal solution is achieved through an inverse power method.

We have demonstrated the implementation of this sparse variant of resolvent analysis on several different configurations, demonstrating that we obtain sparse modes first in the spatial domain, and also in the spatial and temporal domains when using a generalized space–time formulation of resolvent analysis. The first case studied was a statistically stationary turbulent channel flow, a canonical configuration that has received substantial prior study in the context of resolvent analysis (and for many other forms of analysis). When using the standard Fourier transform in time, the proposed sparse resolvent analysis identified spatially sparse modes, in either the wall-normal or wall-normal and spanwise dimensions, depending on the problem set-up. This provides a means for identifying localised spatial structures corresponding to similar amplification as the non-localised structures identified in standard resolvent analysis. For the one-dimensional channel flow analysis, this sparsity simply resulted in modes concentrated about one of the two critical layers. In this case, the symmetry of the configuration means that leading resolvent modes come with a multiplicity of two, and sparse resolvent analysis found a basis for this subspace where each basis vector was as sparse as possible. For the two-dimensional analysis where the spanwise direction was explicitly discretised, sparse resolvent analysis found a leading resolvent mode that consisted of a primary streamwise vortex, surrounded by a pair of streamwise streaks, with comparable total amplification to the Fourier mode identified from standard resolvent analysis. This suggests that the ordered, spanwise repetition of such structures are not essential for their appearance in turbulent flows (under the assumption that linear mechanisms are at least partially responsible for their emergence). Indeed, this is consistent with analysis from turbulent flow data, where spanwise correlations are observed to decay. Indeed, correlations at a wall-normal location corresponding to near-wall streaks are observed to decay (Kim et al. 1987) over a spanwise length scale similar to the streak width itself ($\lambda _z^+ \approx 100$), which is consistent with the spanwise localisation observed for the leading resolvent response mode in figure 2(b).

Next, a space–time generalization of resolvent analysis was applied to the same channel flow configuration. Being a statistically stationary system, each of the space–time modes converged to a single temporal Fourier mode in the absence of sparsity promotion. Therefore, the space–time resolvent modes are equivalent to the standard resolvent modes compiled across all frequencies permissible by the prescribed time horizon. The sparsity-promoting variant of this space–time resolvent analysis yielded structures that displayed localisation in both the spatial and temporal dimensions. These modes again reflect the fact that linear energy amplification can arise from forcing that is limited to a localised region in time. These time-localised structures exhibit features that cannot be seen in Fourier modes, such as the time evolution of inclination angles and wall-normal locations of structures within resolvent modes. Note that the change in the inclination angles of the modes can also be interpreted as regions of the modes with different wall-normal locations evolving with different wave speeds, consistent with the distribution of the mean velocity across the wall-normal extent of the modes. This phenomena is again not observable in standard Fourier-based resolvent analysis. It was noted in § 3.3 that when streamwise and spanwise wavelengths are chosen to be typical of near-wall streaks, the resulting resolvent modes are centred much further from the wall than where such structures are typically found. However, the corresponding sparse space–time mode components exhibited peak amplitudes much closer to the wall, with the streamwise component of the leading response mode having a peak at $y^+\approx 17.5$, which is more consistent with the location of near-wall streaks. Future work could explore in more detail the quantitative extent to which sparse modes predict phenomena present in turbulent channel flow data.

Beyond statistically stationary systems, we next explored the application of our proposed methodology to systems with a time-varying mean. The first such system that was considered was a time-periodic turbulent Stokes boundary layer. In this case, the non-sparse and sparse variants identified leading resolvent modes that were qualitatively different, with the sparse models being significantly more localised in both space and time. Our analysis focused on relatively large-scale structures by specifying streamwise and spanwise wavelengths corresponding to the size of the domain used for DNS. While not shown here, we found that this difference was less apparent for smaller wavelengths, where even the non-sparse variant identified structures that were localised in time, such as the case considered in Ballouz et al. (2024). For the case considered here, the sparse mode components identified were also localised near the peak in turbulence intensity at the wall-normal location where they were situated, suggesting that they identify a mechanism relevant to the generation and amplification of turbulent features at this temporal phase of the system.

The final set of results considered a non-periodic, time-evolving turbulent channel flow subjected to a sudden lateral pressure gradient. This flow was studied during a temporal window in which the flow transitions between two statistically stationary states. It was found that both the non-sparse and sparse modes showed temporal localisation, presumably due to the choice of wavelengths permitting largest energy amplification in a certain time window corresponding to a given orientation of the mean. Overall, we have demonstrated that spatio-temporally sparse resolvent analysis can identify features that can either be qualitatively similar or distinct from the non-sparse equivalent.

Across all of the examples considered, there are several common observations that can be made. First, across all configurations, the streamwise velocity component of the response modes are larger than the wall-normal component. Furthermore, for the time-localised modes, the peak in the amplitude of the streamwise velocity typically occurs slightly after the peak in the wall-normal component. This is consistent with the lift-up mechanism that transfers momentum between these components for cases with non-zero spanwise wavenumbers. For cases where sparsity is promoted in time, we observe that response mode shapes are typically localised within temporal envelopes that are approximately Gaussian in shape, with modes exhibiting several periods of oscillation within such envelopes. For all but the Stokes boundary layer case, the frequency of these oscillations remains approximately constant in time throughout the mode evolution. For the Stokes boundary layer, the sparse response mode oscillation frequency is found to increase continuously with time, approximately doubling throughout the time evolution of the mode. This suggests that for systems with large, rapid variations in the mean velocity field, analysis of linear amplification mechanisms could be more difficult to capture in the frequency (rather than time) domain, given the observed continual transfer of energy across frequencies.

To conclude, the contributions provided by both of the frameworks introduced in this paper, in terms of improved physical insight, deserve a separate discussion. Firstly, the main advantage of the space–time analysis with respect to traditional resolvent analysis resides in its applicability to not only statistically stationary systems, but also time-varying mean flows of any nature. This aspect enables the identification and characterisation of the relevant spatio-temporal coherent structures that are most highly amplified by the linearised governing equations. Secondly, the analysis also recognises the physically relevant temporal instance at which coherent structures of a targeted length scale naturally arises. This was observed in the cases presented in §§ 3.4 and 3.5. Thirdly, considering the application of the space–time framework to time-periodic systems, its capabilities surpass those of harmonic resolvent analysis, as it is able to infer triadic interactions between modes in which each of them can be associated with several temporal frequencies (see the cases presented in § 3.4 for a turbulent Stokes boundary layer). Lastly, even the application of the space–time analysis provides additional physical insight for time-invariant mean flows, as it highlights not only the coherent structures of maximal amplification, but also the associated temporal frequency within a prescribed temporal period (see figure 7 for a turbulent channel flow).

In addition, the incorporation of the sparsity-promoting term to the space–time analysis further expands its potential in a number of ways. For instance, by promoting localisation, it is possible to break down a given spatio-temporal coherent structure formed by an aggregate of temporal frequencies into smaller units (see case presented in figure 11 for a turbulent Stokes boundary layer). This presumably allows for the study of complex linear energy-amplification mechanisms through the examination and characterisation of localised coherent structures of a reduced frequency content. Moreover, the exploitation of the sparsity-promoting space–time analysis could assist in the design of control strategies of targeted length and temporal scales. Note that this analysis not only provides a localised structure that is optimal in terms of the energy gain (this was investigated in more depth in figures 3–4 in the study of the turbulent channel flow with discretisation in the spanwise dimension), but will also be identified in a physically relevant temporal instance. Thus, this technique could be advantageous in the determination of the most optimal location for actuator and/or sensor placement.

The focus of this work has been the development and demonstration of this method across a range of flow configurations. Additional investigations into each of these (and other) configurations are required to quantify the extent to which these methods can be used to predict and understand structures and statistics of such turbulent flows, over a wider range of spatial scales. There are several possible approaches towards assessing the capacity of sparse-resolvent-based models in predicting turbulent flow features. One approach involves a direct comparison between predicted mode structures from sparse resolvent analysis with a data-driven equivalent. This could involve modifying the optimisation problem associated with standard (space-only or space–time) POD methods to include a sparsity-promoting term, or through comparison with previously developed conditional (Schmidt & Schmid 2019) or windowed (Frame & Towne 2022) space–time POD. Beyond this, sparse resolvent modes could also be used to develop simplified ‘turbulence kernels’, which may facilitate modelling similar structures to those captured using sets of standard resolvent modes (Sharma & McKeon 2013; Luhar, Sharma & McKeon 2014a; McKeon 2017), but while also capturing the decay in two-point correlations with increasing spatio-temporal separation that are typically observed in turbulent flows, a property not readily captured with a small number of Fourier modes. Additionally, bases of sparse resolvent modes could also be utilized for real-time flow estimation tasks (Arun, Bae & McKeon 2023), which may be particularly useful if needing to estimate localised features within large domains. The numerical methods employed in the present work have required the explicit formation and decomposition of linear operators discretised in both space and time. Future work could also investigate time-stepping approaches to perform such analyses. Such methods, which have been applied for resolvent (Martini et al. 2020; Farghadan, Martini & Towne 2023) and other linear analyses (Barkley, Blackburn & Sherwin 2008), would likely be particularly computationally advantageous for our methodology.