1. Introduction
The estimation of turbulent wakes behind airfoils has been a longstanding focus in fluid mechanics, with decades of experimental and numerical studies revealing their complex dynamics across a range of conditions (Gartshore Reference Gartshore1967; Reynolds Reference Reynolds1972; Ghaemi & Scarano Reference Ghaemi and Scarano2011; Shamsoddin & Porté-Agel Reference Shamsoddin and Porté-Agel2017; Gupta et al. Reference Gupta, Zhao, Sharma, Agrawal, Hourigan and Thompson2023; Schauerte & Schreyer Reference Schauerte and Schreyer2024). At moderate Reynolds numbers (e.g.
$ \textit{Re} \approx 20\,000$
) and mid-range angles of attack (
$6{-}10^{\circ }$
), a turbulent wake forms downstream of a laminar separation bubble and exhibits vortex shedding near a characteristic Strouhal number (
$St_{\alpha } \approx 0.2$
) based on airfoil thickness (Ducoin, Loiseau & Robinet Reference Ducoin, Loiseau and Robinet2016; Yeh & Taira Reference Yeh and Taira2019). In recent years, non-physics-based, data-driven models – such as neural-network surrogates, black-box system-identification methods and reduced-order models – have demonstrated that key wake features (e.g. dominant shedding frequency, lift and drag fluctuations and large-scale velocity structures) can be predicted directly from experimental measurements without explicitly incorporating the Navier–Stokes equations (Hočevar et al. Reference Hočevar, Širok and Grabec2005; Raissi, Perdikaris & Karniadakis Reference Raissi, Perdikaris and Karniadakis2019; Sekar et al. Reference Sekar, Jiang, Shu and Khoo2019). While these approaches can be remarkably accurate, they typically provide limited physical interpretability and are challenging to use for long-term, real-time estimation and control of high-dimensional turbulent flows. This gap motivates the present work. We seek estimation methods that remain faithful to the underlying flow physics, offer interpretable connections between measurements and target flow structures and are sufficiently low cost to be deployed in real time, while still achieving competitive accuracy relative to state-of-the-art data-driven models.
Among model-based approaches, resolvent analysis has emerged as a powerful tool for studying coherent structures in turbulent flows and for developing estimation and control strategies. By interpreting the linearised Navier–Stokes equations as a forced system, resolvent analysis identifies the input–output relationships associated with the most amplified flow responses in the frequency domain (Jovanović & Bamieh Reference Jovanović and Bamieh2005; McKeon & Sharma Reference McKeon and Sharma2010; Towne, Schmidt & Colonius Reference Towne, Schmidt and Colonius2018). These studies have demonstrated that resolvent modes provide a useful representation of coherent structures in turbulent flows, and have motivated subsequent developments in resolvent-based modelling and control. This framework has also been used to guide both open-loop and closed-loop control designs. For instance, Yeh & Taira (Reference Yeh and Taira2019) applied resolvent analysis to identify control parameters for separation control for a NACA airfoil at
$ \textit{Re}=23\,000$
, a Mach number of
$ \textit{Ma}=0.3$
and angles of attack of
$6^{\circ }$
and
$9^{\circ }$
. Similarly, Liu et al. (Reference Liu, Sun, Yeh, Ukeiley, Cattafesta and Taira2021) reduced pressure fluctuations over a cavity using resolvent-informed open-loop control, while Jin, Illingworth & Sandberg (Reference Jin, Illingworth and Sandberg2022) developed a resolvent-based iterative algorithm for closed-loop estimation and control in cylinder flow.
Building on these foundations, recent work has focused on resolvent-based estimation techniques that leverage the close connection between resolvent modes and spectral proper orthogonal decomposition modes (Towne et al. Reference Towne, Schmidt and Colonius2018) to construct data-informed models of coherent structures (Martini et al. Reference Martini, Cavalieri, Jordan, Towne and Lesshafft2020; Towne, Lozano-Durán & Yang Reference Towne, Lozano-Durán and Yang2020; Amaral et al. Reference Amaral, Cavalieri, Martini, Jordan and Towne2021; Amaral & Cavalieri Reference Amaral and Cavalieri2023; Ying et al. Reference Ying, Liang, Li and Fu2023). While these approaches have shown strong promise in capturing coherent structures from available measurements, most existing implementations have predominantly relied on non-causal formulations, which limits their direct applicability to real-time estimation and control. To address this gap, our recent work has focused on developing causal resolvent-based estimators using the Wiener–Hopf formalism (Noble Reference Noble1958), which enforces causality (Daniele & Lombardi Reference Daniele and Lombardi2007; Martinelli Reference Martinelli2009) and yields optimal estimation kernels suited for closed-loop control (Martini et al. Reference Martini, Jung, Cavalieri, Jordan and Towne2022) within the incompressible solver. Extending this approach, Jung, Bhagwat & Towne (Reference Jung, Bhagwat and Towne2025) developed an efficient computational implementation integrated into the compressible solver, successfully validating causal estimation and control methods for a laminar wake flow. Within the same framework, Towne et al. (Reference Towne, Bhagwat, Zhou, Jung, Martini, Jordan, Audiffred, Maia and Cavalieri2024) demonstrated the accuracy of linearised Navier–Stokes operators for various flows, leveraging methods proposed by Nielsen & Kleb (Reference Nielsen and Kleb2006) and Bhagwat (Reference Bhagwat2021).
While resolvent analysis is naturally suited to globally stable flows – where the linearised operator
$\unicode{x1D63C}$
has all eigenvalues in the stable half-plane and the resolvent
$(i\omega \unicode{x1D644} - \unicode{x1D63C})^{-1}$
(with
$\omega$
the temporal frequency and
$\unicode{x1D644}$
the identity) is well defined over the frequencies of interest – its extension to globally unstable systems is more challenging. In such cases, the usual interpretation of the resolvent as mapping a forcing to a bounded, steady time-periodic response breaks down. As a result, unstable eigenvalues produce unbounded growth rather than a steady state, and their contribution must therefore be treated or regularised with additional care (Colonius & Towne Reference Colonius and Towne2025). One strategy to mitigate this involves exponential discounting to shift unstable eigenvalues into the stable region of the spectrum (Jovanović & Bamieh Reference Jovanović and Bamieh2005; Yeh et al. Reference Yeh, Benton, Taira and Garmann2020). Alternatively, data-driven resolvent formulations have been proposed to bypass the reliance on a globally stable linear operator (Martini et al. Reference Martini, Jung, Cavalieri, Jordan and Towne2022; Jung et al. Reference Jung, Bhagwat and Towne2025). These methods maintain the benefits of resolvent-based estimation, specifically by remaining robust against global flow instability.
Linear stochastic estimation and related techniques have a long history in turbulent flow estimation (Cole, Glauser & Guezennec Reference Cole, Glauser and Guezennec1992; Adrian Reference Adrian1994; Tinney et al. Reference Tinney, Coiffet, Delville, Hall, Jordan and Glauser2006; Murray & Ukeiley Reference Murray and Ukeiley2007). These methods use second-order statistics to construct a linear map from measured quantities to an unobserved flow variable. Within the resolvent-based framework, the data-driven construction of the non-causal estimation kernel (Martini et al. Reference Martini, Cavalieri, Jordan, Towne and Lesshafft2020) is mathematically similar to frequency-domain stochastic estimators, in that it also yields a linear, statistically optimal mapping based on prior data. The present framework differs from this traditional approach in two key respects. First, we explicitly formulate a causal estimator that uses both past and current measurements to predict the current state, yielding a time-resolved input–output operator rather than a purely instantaneous mapping. Second, we embed the estimation problem within the resolvent formalism, so that the estimator is directly tied to the linear input–output structure of the Navier–Stokes equations. This connection provides a physics-based lens through which to understand the method and its results. In this sense, the non-causal estimator recovers a classical stochastic-estimation-like mapping from measurements to target quantities, while the resolvent-based causal formulation represents a systematic extension of these ideas.
Sensor placement is crucial for accurate flow estimation, influencing the fidelity of sensor-based measurements and subsequent control strategies. While formal optimisation of sensor placement can significantly enhance estimation accuracy (Manohar et al. Reference Manohar, Brunton, Kutz and Brunton2018; Sashittal & Bodony Reference Sashittal and Bodony2021; Jin et al. Reference Jin, Illingworth and Sandberg2022), such approaches often face practical limitations due to physical and operational constraints. Consequently, this study does not explicitly aim to optimise realistic sensor placement. Instead, effective sensor locations are determined empirically by assessing coherence between individual sensor–target pairs and by investigating single-sensor estimation errors, thereby providing practical insights relevant to realistic experimental scenarios. Within resolvent-based frameworks, spatial coherence has emerged as a particularly valuable metric for evaluating the suitability of sensor–target pairs. Recent experimental studies (Maia et al. Reference Maia, Jordan, Cavalieri, Martini, Sasaki and Silvestre2021; Audiffred et al. Reference Audiffred, Cavalieri, Brito and Martini2023, Reference Audiffred, Cavalieri, Maia, Martini and Jordan2024) and numerical analyses (Towne et al. Reference Towne, Bhagwat, Zhou, Jung, Martini, Jordan, Audiffred, Maia and Cavalieri2024; Jung et al. Reference Jung, Bhagwat and Towne2025) demonstrate that coherence-based approaches effectively identify regions that yield accurate estimation and control. Furthermore, analysing single-sensor estimation errors offers direct guidance for selecting sensor locations, particularly in high-coherence regions or along predictable disturbance propagation paths (Jung et al. Reference Jung, Bhagwat and Towne2025). Alternatively, streamline visualisations can effectively identify sensor placements dynamically connected to target regions in convectively dominated flows (Jung Reference Jung2024).
The objective of the present study is to estimate unsteady turbulent fluctuations in the wake of a NACA0012 airfoil at Re = 23 000 using limited measurements on the surface of the airfoil. Building on our prior work on a laminar wake (Jung et al. Reference Jung, Bhagwat and Towne2025), we employ a causal resolvent-based estimation approach and focus on estimating velocity fluctuations in the turbulent wake, for which the flow exhibits a broadband, multi-scale dynamics rather than periodic shedding. Unlike the laminar case, no upstream forcing is required to induce chaotic fluctuations; instead, we estimate naturally occurring turbulent fluctuations. The linearised Navier–Stokes operator in this regime is globally unstable, making traditional resolvent formulations poorly suited. To address this, we adopt a data-driven approach that constructs transfer functions from nonlinear simulation data, using empirical cross-spectral densities to approximate the resolvent operator (Martini et al. Reference Martini, Jung, Cavalieri, Jordan and Towne2022; Jung et al. Reference Jung, Bhagwat and Towne2025). These estimators leverage the statistical structure of coloured-in-time nonlinear forcing (Zare, Jovanović & Georgiou Reference Zare, Jovanović and Georgiou2017; Towne et al. Reference Towne, Lozano-Durán and Yang2020; Morra et al. Reference Morra, Nogueira, Cavalieri and Henningson2021) and can be constructed without a priori model reduction, maintaining computational efficiency. Furthermore, we enforce causality through a Wiener–Hopf-based formalism (Noble Reference Noble1958), enabling real-time estimation for control applications (Jung et al. Reference Jung, Martini, Cavalieri, Jordan, Lesshafft and Towne2020; Martini et al. Reference Martini, Jung, Cavalieri, Jordan and Towne2022). This work demonstrates that resolvent-based methods can provide accurate, interpretable and scalable estimators for spanwise-periodic turbulent wakes by exploiting coherent flow structures (Towne et al. Reference Towne, Schmidt and Colonius2018).
The remainder of this paper is organised as follows. In § 2, we describe the flow configuration, numerical simulation, validation procedures and the dataset used to construct data-driven estimation kernels. Based on this single simulation, we consider three estimation cases that target: (i) the spanwise-averaged flow, (ii) an individual spanwise-Fourier mode and (iii) the full three-dimensional flow field at the mid-span plane. In § 3, we briefly outline the resolvent-based estimation framework, including the theoretical formulation, the derivation of the causal resolvent-based estimation kernels via the Wiener–Hopf method and implementation of the resolvent-based estimation tools. Estimation results for the three cases are discussed in detail in § 4. Finally, concluding remarks and future research directions are presented in § 5.
2. Problem set-up
2.1. Problem description
We consider the turbulent wake of a NACA0012 airfoil at a moderate chord-based Reynolds number of
${Re}_{L_{c}} = 23{\,}000$
and an angle of attack of
$\alpha = 6^\circ$
, as studied in several previous works (Kojima et al. Reference Kojima, Nonomura, Oyama and Fujii2013; Munday & Taira Reference Munday and Taira2018; Yeh & Taira Reference Yeh and Taira2019; Towne et al. Reference Towne2023). The free-stream Mach number is set to
$Ma_{\infty }\equiv U_{\infty }/c_{\infty }=0.3$
, where
$U_{\infty }$
and
$c_{\infty }$
denote the free-stream velocity and sound speed, respectively. Shear-stress sensors are placed on the airfoil surface, while the targets are located downstream in the turbulent wake (see § 4 for details). In this paper, we present a resolvent-based estimation of the time-series fluctuations in the streamwise velocity
$u_{x}^{\prime }$
about the time-averaged mean flow
$\bar {u}_{x}$
.
2.2. Simulation
A large-eddy simulation (LES) utilising the high-fidelity compressible flow solver CharLES from Cadence Design Systems (Brès et al. Reference Brès, Ham, Nichols and Lele2017, Reference Brès, Jordan, Jaunet, Le Rallic, Cavalieri, Towne, Lele, Colonius and Schmidt2018, Reference Brès, Bose, Ivey, Emory and Ham2022; Yeh & Taira Reference Yeh and Taira2019) is conducted to simulate turbulent airfoil flow. The simulation is conducted on Koehr, the U.S. Department of Defense’s supercomputer, an HPE SGI 8600 system with a peak performance of 3.05 PFLOPS. For this study, 50 nodes of the supercomputer, each equipped with 48 cores, were utilised. A C-shaped mesh, generated using Pointwise, is depicted in figure 1(a). Figures 1(b) and 1(c) illustrate the mean and instantaneous streamwise velocity fields, respectively. The airfoil’s leading edge is positioned at the origin,
$x/L_{c} = y/L_{c} = 0$
. The computational domain spans
$x/L_{c} \in [-20, 21]$
,
$y/L_{c} \in [-20, 20]$
and
$z/L_{c} \in [-0.1, 0.1]$
. Compared with other studies, we do not employ adaptive mesh refinement. Instead, to reduce the total number of cells, we refine the mesh in the transverse and spanwise directions. A characteristic boundary condition
$[\rho , u_{x}, u_{y}, u_{z}, p] = [\rho _{\infty }, U_{\infty }, 0, 0, p_{\infty }]$
, where
$\rho$
is the fluid density,
$U_\infty$
is the freestream velocity, and
$p$
is the pressure, was applied at the domain edges, with a sponge layer implemented for
$x/L_{c} \in [11, 21]$
to mitigate spurious reflections at the outflow. Previous work has shown that airfoil wakes at moderate Reynolds numbers can be sensitive to numerous factors, including free-stream disturbances. No incoming perturbations are imposed in the present simulation. The grid contains approximately 39.5 million cells. Time integration is performed with a maximum Courant–Friedrichs–Lewy number of 0.84, consistent with previous studies (Yeh & Taira Reference Yeh and Taira2019; Towne et al. Reference Towne2023). To obtain the dataset used to construct the transfer function, flow-field snapshots are collected at a constant time step of
$\Delta t U_{\infty } / L_{c} = 4.32 \times 10^{-5}$
.
Computational mesh and simulation: (a) the complete C-shaped grid, with a zoomed-in view of the mesh near the wall shown in the inset; (b) the mean streamwise velocity
$u_{x}$
; and (c) a snapshot of the instantaneous streamwise velocity
$u_{x}$
obtained from the large-eddy simulation.
We validate the LES by comparing the pressure coefficient (
$C_p$
) and the aerodynamic force coefficients (
$C_D$
and
$C_L$
) with previous works in table 1 and figure 2. These are defined as
\begin{align} C_{p}=\frac {p-p_{\infty }}{\frac {1}{2} \rho _{\infty } U_{\infty }^{2}},\quad C_{D}=\frac {F_{D}}{\frac {1}{2} \rho _{\infty } U_{\infty }^{2} A},\quad \textrm {and} \quad C_{L}=\frac {F_{L}}{\frac {1}{2} \rho _{\infty } U_{\infty }^{2} A} \end{align}
where
$C_D$
is the drag coefficient,
$C_L$
is the lift coefficient,
$F_D$
and
$F_L$
are the drag and lift forces, respectively, and
$A$
is the reference area. Here, the overbar
$\bar {(\boldsymbol{\cdot })}$
is used to indicate a time average. Experimental data are provided by Kim et al. (Reference Kim, Yang, Chang and Chung2009). The closest analogue to our simulation is that of Yeh & Taira (Reference Yeh and Taira2019), which entails the same physical set-up and solver but a different grid. The simulations of Kojima et al. (Reference Kojima, Nonomura, Oyama and Fujii2013) and Munday & Taira (Reference Munday and Taira2018) differ due to their use of an incompressible flow solver. As summarised in table 1, the present time-averaged lift and drag coefficients differ from those of Yeh & Taira (Reference Yeh and Taira2019) by 4 % and 0.4 %, respectively. These discrepancies lie within the range typically reported between independent simulations of similar flows. In addition, the surface pressure distribution in figure 2 exhibits trends and overall behaviour that are consistent with both the experimental data of Kim et al. (Reference Kim, Yang, Chang and Chung2009) and the numerical results of Kojima et al. (Reference Kojima, Nonomura, Oyama and Fujii2013). These comparisons support the adequacy of the present grid resolution and LES set-up for the present study.
Comparison of the time-averaged drag and lift coefficients with the reference data for a NACA0012 airfoil at
${Re}_{L_{c}} = 23\,000$
and angle of attack
$6^{\circ }$
.
Comparison of the pressure coefficient
$C_{p}$
with previous works.
As shown in figure 3, the shear layer roll-up on the suction surface occurs at
$x/L_{c}\approx 0.4$
, leading to spanwise vortices. A laminar separation bubble is observed on the suction side of the airfoil in the region
$0.1 \lesssim x/L_{c} \lesssim 0.84$
, as indicated by the mean streamwise velocity shown in figure 1(b). Laminar–turbulent transition is detected at
$x/L_{c}\approx 0.6$
, as determined by the termination of the pressure coefficient plateau illustrated in figure 2. Finally, a turbulent wake develops behind the airfoil.
Visualisation of the
$Q$
-criterion showing the key flow features.
2.3. Dataset
We consider three distinct representations of the spanwise-periodic flow: spanwise-averaged, spanwise-Fourier and mid-span plane cases. The spanwise-averaged case is the simplest representation, defined mathematically as
\begin{align} {\boldsymbol{q}}_{\textit{span}\mbox{-}{avg}}(x,y,t) = \frac {1}{L_z} \int _{-L_{z}/2}^{L_{z}/2} \boldsymbol{q}(x,y,z,t) \, {\rm d}z, \end{align}
where
$L_z$
is the spanwise length. Here,
$\boldsymbol{q}_{{span\mbox{-}avg}}$
denotes a spatially averaged quantity in the spanwise direction. Spanwise-Fourier decomposition, a common technique in biglobal stability analyses of three-dimensional flows (Theofilis Reference Theofilis2003), has been previously utilised to investigate resolvent gains at various spanwise wavenumbers (Yeh & Taira Reference Yeh and Taira2019; Towne et al. Reference Towne2023). We similarly analyse spanwise-Fourier modes defined by
\begin{align} \hat {\boldsymbol{q}}(x,y,k_{z},t) = \frac {1}{L_z}\int _{-L_{z}/2}^{L_{z}/2} \boldsymbol{q}(x,y,z,t)e^{-ik_{z}z} \, {\rm d}z, \end{align}
where
$k_{z}$
denotes the spanwise wavenumber. This approach allows us to progressively incorporate more complex flow structures by examining the spanwise-Fourier modes. Finally, the mid-span plane case focuses on the flow at the central spanwise location, defined as
\begin{align} \boldsymbol{q}_{\textit{mid}}(x,y,t) = \boldsymbol{q}(x,y,z_{\textit{mid}},t), \end{align}
where
$z_{\textit{mid}} = 0$
is the mid-span location. This representation helps capture spanwise effects in the estimation process by examining the flow at a representative spanwise position and is consistent with a realistic use case where sensors are distributed along the airfoil at a single spanwise position. These representations enhance our understanding of spanwise effects on the flow and support investigations into sensor placement and flow estimation strategies for a spanwise-periodic flow.
The dataset comprises
$n_t = 75{\,}000$
time-resolved, three-dimensional snapshots of
$[\rho , \rho u_x, \rho u_y, \rho u_z, T]$
, where
$T$
is the temperature, spanning a non-dimensional time window of
$tU_\infty / L_c \in [0, 324]$
. These snapshots are sampled at a uniform interval of
$\Delta t U_\infty / L_c = 0.00432$
, which resolves the dominant vortex-shedding frequency at
$St_{\alpha } \equiv f L_c (\sin \alpha ) / U_\infty \approx 0.15$
. This finer temporal resolution and longer duration compared with previous datasets (Towne et al. Reference Towne2023) enable capturing of smaller-scale turbulent fluctuations in time and provide a larger number of snapshots, resulting in improved convergence of the estimation kernels. The time window captures approximately 466 shedding cycles, providing a sufficiently long time series to ensure statistical convergence of spectral quantities, including the cross-spectral densities used in constructing the estimation kernel. The spatial domain covers
$x / L_c \in [-0.5, 2.5]$
and
$y / L_c \in [-0.5, 0.5]$
, capturing the near-wake region with high spatial fidelity.
To construct resolvent-based estimation kernels, we follow the implementation of Jung et al. (Reference Jung, Bhagwat and Towne2025). Prior to estimation, the dataset was preprocessed by subtracting the temporal mean from each variable and applying a Hamming window (Hamming Reference Hamming1997) to reduce spectral leakage in the Fourier transform. The dataset is divided into 80 % for training and 20 % for testing, using a single fixed partition. This same training–testing split is used for all three estimation cases (spanwise-averaged, spanwise-Fourier and mid-span). Given that the dataset contains more than 400 shedding cycles, this fixed split is expected to provide statistically representative training and testing sets. For assessing generalisation performance, 80 % of the data are used to train the estimation model, while the remaining 20 % are reserved for testing.
3. Resolvent-based estimation framework
In this section, we briefly review the resolvent-based estimation framework developed in recent studies (Martini et al. Reference Martini, Cavalieri, Jordan, Towne and Lesshafft2020, Reference Martini, Jung, Cavalieri, Jordan and Towne2022; Towne et al. Reference Towne, Lozano-Durán and Yang2020; Jung et al. Reference Jung, Bhagwat and Towne2025).
3.1. Estimation system
We start with the compressible Navier–Stokes equations written as
\begin{align} \frac {\partial \boldsymbol{q}}{\partial t}= \mathcal{F} (\boldsymbol{q}), \end{align}
where
$\boldsymbol{q}$
is a state vector of flow variables
$[\rho ,\rho u_{x},\rho u_{y},\rho u_{z},\rho E]^{T}$
and
$\mathcal{F}$
is the nonlinear Navier–Stokes operator. The equations are linearised using a Reynolds decomposition, giving
\begin{align} \frac {\partial \boldsymbol{q}'}{\partial t}-\unicode{x1D63C}\boldsymbol{q'}=\boldsymbol{f}(\boldsymbol{\bar {q}},\boldsymbol{q'}), \end{align}
where
$\bar {\boldsymbol{q}}$
and
$\boldsymbol{q}'$
represent the mean and perturbation state vectors of the flow variables, respectively,
$\unicode{x1D63C} = ( {\partial \mathcal{F}}/{\partial \boldsymbol{q}}\big )|_{\bar {q}}$
is the linearised Navier–Stokes operator and
$\boldsymbol{f}$
represents the remaining nonlinear terms. For convenience, we drop the
$(\boldsymbol{\cdot })^{\prime }$
superscript for perturbations from this point on.
To derive the estimation method, we consider a generalisation of (3.1) in the form of the linear-time-invariant (LTI) system
\begin{align} \frac {{\rm d}\boldsymbol{q}}{{\rm d}t}(t)&= \unicode{x1D63C}\boldsymbol{q}(t)+\unicode{x1D63D}_{f}\boldsymbol{f}(t), \\[-12pt] \nonumber \end{align}
\begin{align} \boldsymbol{y}(t)&= \unicode{x1D63E}_{y}\boldsymbol{q}(t)+\boldsymbol{n}(t), \\[-12pt] \nonumber \end{align}
\begin{align} \boldsymbol{z}(t)&= \unicode{x1D63E}_{z}\boldsymbol{q}(t), \\[9pt] \nonumber \end{align}
where the forcing matrix
$\unicode{x1D63D}_{f} \in {\rm \mathbb{C}}^{n \times n_{f}}$
can be used to restrict the form of the forcing
$\boldsymbol{f}$
, and the measurement
$\boldsymbol{y}$
and target
$\boldsymbol{z}$
indicate readings of the state perturbation extracted by a measurement matrix
$\unicode{x1D63E}_{y}\in {\rm \mathbb{C}}^{n_{y} \times n}$
and a target matrix
$\unicode{x1D63E}_{z}\in {\rm \mathbb{C}}^{n_{z} \times n}$
. The number of sensors and targets is denoted
${n}_{y}$
and
$n_{z}$
, respectively. The vector
$\boldsymbol{n}$
indicates the sensor noise. The measurement
$\boldsymbol{y}$
is chosen to be the shear stress
$\boldsymbol{\tau }_{w}$
at one or more points on the surface of the airfoil. It is computed from the state
$\boldsymbol{q}$
using a suitably defined measurement matrix
$\unicode{x1D63E}_{y}$
, as in Jung et al. (Reference Jung, Bhagwat and Towne2025). The target
$\boldsymbol{z}$
represents one or more Gaussian sensors, each with spatial support given by
\begin{align} \alpha e^{-(x-x_{c})^{2}/ 2 \sigma _{x} ^{2} -(y-y_{c})^{2}/ 2 \sigma _{y} ^{2}-(z-z_{c})^{2}/ 2 \sigma _{z} ^{2}}, \end{align}
where
$\sigma _{x}$
,
$\sigma _{y}$
and
$\sigma _{z}$
define the width of the Gaussian distribution in each coordinate direction, and
$\alpha$
is a weighting factor representing the integrated contribution over the Gaussian support. In this study, we set all Gaussian widths to the same value,
$ \sigma _{x} = \sigma _{y} = \sigma _{z} = 0.2$
. These targets are extracted from the state using the target matrix
$\unicode{x1D63E}_{z}$
.
By Fourier transforming the linear system into the frequency domain and defining the resolvent operator as
$\unicode{x1D64D} = (i \omega \unicode{x1D644} - \unicode{x1D63C})^{-1}$
, where
$\unicode{x1D644}$
is the identity matrix and
$\omega$
is the temporal frequency, we derive modified resolvent operators used for constructing the resolvent-based kernels for estimation. The measurement and target vectors can be written as
\begin{align} \boldsymbol{\hat {y}}&=\unicode{x1D64D}_{yf}\boldsymbol{\hat {f}}+\boldsymbol{\hat {n}}, \\[-12pt] \nonumber \end{align}
\begin{align} \boldsymbol{\hat {z}}&= \unicode{x1D64D}_{zf}\boldsymbol{\hat {f}}, \\[9pt] \nonumber \end{align}
with
$\unicode{x1D64D}_{yf}=\unicode{x1D63E}_{y}\unicode{x1D64D}\unicode{x1D63D}_{f}$
, and
$\unicode{x1D64D}_{zf}=\unicode{x1D63E}_{z}\unicode{x1D64D}\unicode{x1D63D}_{f}$
. The notation (
$\hat {\boldsymbol{\cdot }}$
) indicates a quantity in the frequency domain.
3.2. Kernels for resolvent-based estimation
We perform estimation with three different kernels: non-causal, truncated non-causal (TNC) and causal (Martini et al. Reference Martini, Jung, Cavalieri, Jordan and Towne2022). In the context of this work, an estimator is called causal if it depends only on past and present sensor measurements, but not on future measurements. To make this precise, we introduce the time lag
$\tau$
, which appears in convolution as
$t-\tau$
. For a given estimation time
$t$
, values with
$\tau \gt 0$
correspond to past times,
$\tau =0$
corresponds to the present time, and
$\tau \lt 0$
corresponds to future times. A causal estimator therefore uses only
$\tau \geqslant 0$
, whereas an estimator that uses
$\tau \lt 0$
is non-causal.
First, we start with a fully non-causal estimator using a convolution function, written as
\begin{align} \tilde {\boldsymbol{z}}_{nc}(t) = \int _{-\infty }^{\infty } \unicode{x1D64F}_{nc}(t-\tau )\,\boldsymbol{y}(\tau )\,\mathrm{d}\tau , \end{align}
where
$\unicode{x1D64F}_{nc} \in \mathbb{R}^{n_{z} \times n_{y}}$
is a non-causal estimation kernel between the sensor measurements and the estimated target state. The objective is to select the kernel to minimise the expected value of the estimation error
$e(t)=\tilde {\boldsymbol{z}}(t)-\boldsymbol{z}(t)$
, formalised by the cost function
\begin{align} \unicode{x1D645}_{nc} = \int _{-\infty }^{\infty } \mathbb{E}\{ e(t)^{\dagger } e(t) \}\,\mathrm{d}t, \end{align}
where
$\mathbb{E}\{\boldsymbol{\cdot }\}$
is the expectation operator. The overall estimation framework is illustrated in the block diagram in figure 4(a). By minimising the cost function, we obtain a non-causal estimation kernel
\begin{align} \kern2.5pt \hat {\kern-2.5pt \unicode{x1D64F}}_{nc}(\omega ) = \bigl (\unicode{x1D64D}_{zf}\,\kern2.5pt \hat {\kern-2.5pt \unicode{x1D641}}\,\unicode{x1D64D}_{yf}^{\dagger }\bigr ) \bigl (\unicode{x1D64D}_{yf}\,\kern2.5pt \hat {\kern-2.5pt \unicode{x1D641}}\,\unicode{x1D64D}_{yf}^{\dagger } + \kern1.5pt \hat {\kern-1.5pt \unicode{x1D649}}\bigr )^{-1}, \end{align}
where
$\kern2.5pt \hat {\kern-2.5pt \unicode{x1D641}} = \mathbb{E}\{ \kern2.5pt \hat {\boldsymbol{\kern-2.5pt f}} \kern2.5pt \hat {\boldsymbol{\kern-2.5pt f}}^{\dagger } \}$
and
$\kern1.5pt \hat {\kern-1.5pt \unicode{x1D649}} = \mathbb{E}\{ \hat {\boldsymbol{n}} \hat {\boldsymbol{n}}^{\dagger } \}$
denote the forcing and sensor-noise cross-spectral density matrices, respectively. Note that the space–time colour of the forcing terms is explicitly included in (3.8), so that the nonlinear dynamics of the flow can be statistically accounted for in the transfer function.
Block diagram of resolvent-based estimation: (a) design of the resolvent-based estimator based on the linear system and (b) application to the nonlinear system. Here,
$(\boldsymbol{\cdot })'$
indicates fluctuations around the mean quantity,
$(\bar {\boldsymbol{\cdot }})$
denotes mean quantity and
$(\tilde {\boldsymbol{\cdot }})$
represents an estimated quantity.
Second, for real-time estimation, only current and past (but not future) sensor measurements are available. A common practical approach is therefore to truncate the non-causal kernel by setting the part associated with future measurements (
$\tau \lt 0$
) to zero. Using the time-lag variable explicitly, we define the TNC kernel as
\begin{eqnarray} \unicode{x1D64F}_{\textit{tnc}}(\tau ) = \begin{cases} \unicode{x1D64F}_{nc}(\tau ), & \tau \geqslant 0,\\ 0, & \tau \lt 0, \end{cases} \end{eqnarray}
so that the TNC estimates are computed as
\begin{align} \tilde {\boldsymbol{z}}_{\textit{tnc}}(t) = \int _{-\infty }^{0} \unicode{x1D64F}_{\textit{tnc}}(t-\tau )\,\boldsymbol{y}(\tau )\,\mathrm{d}\tau , \end{align}
which uses only past and present measurements relative to time
$t$
. The TNC estimator is causal in terms of the information it uses; the adjective ‘non-causal’ refers only to how the kernel is constructed (by truncation of a non-causal solution) rather than to its actual temporal support.
Lastly, we derive an optimal causal estimation kernel by enforcing the causality constraint directly in the cost function via a Wiener–Hopf formalism. A causal estimator is written as
\begin{align} \tilde {\boldsymbol{z}}_{c}(t)= \int _{-\infty }^{0} \unicode{x1D64F}_{c}(t-\tau )\,\boldsymbol{y}(\tau )\,\mathrm{d}\tau , \end{align}
where
$\unicode{x1D64F}_{c} \in \mathbb{R}^{n_{z} \times n_{y}}$
is the causal estimation kernel. To impose causality, we require that the kernel vanishes for negative time lags,
$\unicode{x1D64F}_{c}(\tau \lt 0)=0$
, so that the estimator depends only on past and present measurements. This constraint is enforced by introducing a Lagrange multiplier
${\boldsymbol{\varLambda}}_{-}(t)$
, which acts on the non-causal part of the kernel, i.e. the component whose inverse Fourier transform is supported on
$\tau \lt 0$
. The cost function in (3.7) is then modified to
\begin{align} \unicode{x1D645}_{c}= \int _{-\infty }^{\infty }\big ( \mathbb{E}\{ e(t)^{\dagger } e(t) \} + {\boldsymbol{\varLambda}}_{-}(t)\,\unicode{x1D64F}_{c}(t) + {\boldsymbol{\varLambda}}_{-}^{\dagger }(t)\,\unicode{x1D64F}_{c}^{\dagger }(t)\big )\,\mathrm{d}t. \end{align}
Here, the subscripts
$(\boldsymbol{\cdot })_{+}$
and
$(\boldsymbol{\cdot })_{-}$
denote, respectively, the components of a matrix or function whose inverse Fourier transforms are supported on
$\tau \geqslant 0$
(causal) and
$\tau \lt 0$
(non-causal), obtained via a Wiener–Hopf factorisation. The optimisation problem in (3.12) therefore seeks the kernel
$\unicode{x1D64F}_{c}(\tau )$
that minimises the mean-square estimation error subject to the causality constraint
$\unicode{x1D64F}_{c}(\tau \lt 0)=0$
.
In this framework, variation of
$\unicode{x1D645}_{c}$
with respect to
${\boldsymbol{\varLambda}}_{-}$
enforces the condition
$(\unicode{x1D64F}_{c})_{-}=0$
, i.e. that there is no support for
$\tau \lt 0$
, while variation with respect to
$\unicode{x1D64F}_{c}$
yields the optimal kernel under this constraint. The Lagrange multiplier
${\boldsymbol{\varLambda}}_{-}$
is thus associated specifically with the non-causal component
$(\boldsymbol{\cdot })_{-}$
, and its role is to drive that component of the kernel to zero. A brief summary of the Wiener–Hopf procedure is provided with further details in Appendix A. The resulting causal kernel in the frequency domain is
\begin{eqnarray} \kern2.5pt \hat {\kern-2.5pt \unicode{x1D64F}}_{c}(\omega )= \bigl [ \unicode{x1D64D}_{zf}\,\kern2.5pt \hat {\kern-2.5pt \unicode{x1D641}}\,\unicode{x1D64D}_{yf}^{\dagger }\, \bigl (\unicode{x1D64D}_{yf}\,\kern2.5pt \hat {\kern-2.5pt \unicode{x1D641}}\,\unicode{x1D64D}_{yf}^{\dagger } + \kern1.5pt \hat {\kern-1.5pt \unicode{x1D649}}\bigr )^{-1}_{-}\bigr ]_{+}\bigl (\unicode{x1D64D}_{yf}\,\kern2.5pt \hat {\kern-2.5pt \unicode{x1D641}}\,\unicode{x1D64D}_{yf}^{\dagger } + \kern1.5pt \hat {\kern-1.5pt \unicode{x1D649}}\bigr )^{-1}_{+}. \end{eqnarray}
All the transfer functions we consider involve the forcing cross-spectral density matrix that can account for the impact of the nonlinearity in a statistical sense. However, our method is more suitable for the real-time estimation approach since causality is constrained. In this paper, we present the results for both the causal and TNC transfer functions to assess how optimally enforcing causality impacts the estimation accuracy for turbulent flow.
The linearisation about the turbulent mean is globally unstable (Yeh & Taira Reference Yeh and Taira2019). While the resolvent operator
$\unicode{x1D64D}(i\omega )=(i\omega \unicode{x1D644}-\unicode{x1D63C})^{-1}$
is a standard frequency-domain object, the central difficulty in the present work is not its formal definition but the assumptions required to construct an optimal transfer function from it. Our derivation of the resolvent-based estimator (and its equivalent frequency-domain transfer function) assumes statistically stationary input/output processes, so that time-invariant correlations and cross-spectral densities exist. This assumption is satisfied by the statistically stationary nonlinear (LES) data used in this work. In contrast, trajectories generated by integrating an unforced unstable linear system grow exponentially and are not statistically stationary, so the correlations/spectra needed by the Wiener–Hopf method are not well defined in that setting. Therefore, even if one were to numerically regularise the unstable integration (e.g. via discounting), the resulting signals would correspond to a modified, non-physical process and would not provide an appropriate proxy for the stationary turbulent statistics for which the optimal kernels are derived. To circumvent this issue, we use a data-driven approach (Martini et al. Reference Martini, Jung, Cavalieri, Jordan and Towne2022) that estimates the required correlations/cross-spectra directly from statistically stationary flow data (Audiffred et al. Reference Audiffred, Cavalieri, Brito and Martini2023; Jung et al. Reference Jung, Bhagwat and Towne2025), thereby constructing the same resolvent-based kernels without explicitly integrating the unstable linearised operator.
Following the notation of the LTI system in § 3.1, the datasets of the sensor and target readings are Fourier transformed, expressed as
\begin{align} \begin{bmatrix} \kern2pt \hat {\kern-2pt \boldsymbol{y}}\\[5pt] \hat {\boldsymbol{z}} \end{bmatrix} = \begin{bmatrix} \unicode{x1D64D}_{yf} & 1 \\[5pt] \unicode{x1D64D}_{zf} & 0 \end{bmatrix} \begin{bmatrix} \kern2.5pt \hat {\boldsymbol{\kern-2.5pt f}} \\[5pt] \hat {\boldsymbol{n}} \end{bmatrix} . \end{align}
Assuming that the additive sensor noise is uncorrelated with the forcing, i.e.
$\mathbb{E}\{\hat{\boldsymbol n}\,\hat{\kern-3.5pt \boldsymbol f}^{\dagger}\}=\mathbb{E}\{ \ \hat{\kern-3.5pt \boldsymbol f}\hat{\boldsymbol n}^{\dagger}\}=0$
, and defining the noise-free measurement signal
$\hat{\boldsymbol y}_s=\unicode{x1D64D}_{yf}\,\hat{\kern-3.5pt \boldsymbol f}$
, computing the cross-spectral density of
$[\hat{\boldsymbol y}_s \quad \hat{\boldsymbol z}]^T$
gives
\begin{align}\begin{bmatrix}\unicode{x1D64E}_{yy} & \unicode{x1D64E}_{yz} \\[5pt] \unicode{x1D64E}_{zy} & \unicode{x1D64E}_{zz} \end{bmatrix} = \begin{bmatrix} \unicode{x1D64D}_{yf}\,\hat {\kern-2.5pt \unicode{x1D641}}\unicode{x1D64D}_{yf}^{\dagger} &\unicode{x1D64D}_{yf}\,\hat {\kern-2.5pt \unicode{x1D641}}\unicode{x1D64D}_{zf}^{\dagger} \\[5pt] \unicode{x1D64D}_{zf}\,\hat {\kern-2.5pt \unicode{x1D641}}\unicode{x1D64D}_{yf}^{\dagger} & \unicode{x1D64D}_{zf}\,\hat {\kern-2.5pt \unicode{x1D641}}\unicode{x1D64D}_{zf}^{\dagger} \end{bmatrix},\end{align}
with
$\unicode{x1D64E}_{yy}=\mathbb{E}\{\hat{\boldsymbol y}_s\hat{\boldsymbol y}_s^{\dagger}\}$
and
$\unicode{x1D64E}_{zy}=\mathbb{E}\{\hat{\boldsymbol z}\hat{\boldsymbol y}_s^{\dagger}\}$
. The sensor-noise cross-spectral density
$\kern1.5pt \hat {\kern-1.5pt \unicode{x1D649}}=\mathbb{E}\{\hat{\boldsymbol n}\hat{\boldsymbol n}^{\dagger}\}$
is incorporated separately when forming the estimation kernels. In practice, these cross-spectral density (CSD) matrices are obtained directly from the time-resolved data. Specifically, we compute the Fourier transform of the measured time series using a discrete Fourier transform, and then form the CSD by Welch’s method (Welch Reference Welch1967). Since the right-hand side of (3.15) contains the terms needed to build the estimation kernels, this shows that the correlations on the left-hand side can be used in their place in (3.8) and (3.13). The data-driven non-causal and causal estimation kernels in (3.8) and (3.13) are computed using the CSDs from (3.15), yielding
\begin{align} \kern2.5pt \hat {\kern-2.5pt \unicode{x1D64F}}_{nc} &= \unicode{x1D64E}_{zy} ( \unicode{x1D64E}_{yy} + \kern1.5pt \hat {\kern-1.5pt \unicode{x1D649}})^{-1}, \\[-12pt] \nonumber \end{align}
\begin{align} \kern2.5pt \hat {\kern-2.5pt \unicode{x1D64F}}_{c} &= \big[\unicode{x1D64E}_{zy} ( \unicode{x1D64E}_{yy} + \kern1.5pt \hat {\kern-1.5pt \unicode{x1D649}})^{-1}_{-} \big]_{+} (\unicode{x1D64E}_{yy} + \kern1.5pt \hat {\kern-1.5pt \unicode{x1D649}})_{+}^{-1}. \\[9pt] \nonumber \end{align}
Note that in the derivation, the forcing CSD
$\kern2.5pt \hat {\kern-2.5pt \unicode{x1D641}}$
appears explicitly in (3.16). As a result, the CSDs inherently contain statistical information about the nonlinearities of the flow through
$\kern2.5pt \hat {\kern-2.5pt \unicode{x1D641}}$
. In the data-driven approach,
$\kern2.5pt \hat {\kern-2.5pt \unicode{x1D641}}$
is accounted for implicitly in (3.15) through the empirically obtained CSDs
$\unicode{x1D64E}_{yy}$
and
$\unicode{x1D64E}_{zy}$
. The mean flow is subtracted from both the measurement readings and target values. Estimation kernels are computed offline for the linear system, as illustrated in figure 4(a). These kernels are then used for real-time (online) estimation of velocity fluctuations in the wake, applying the approach to the nonlinear system, as depicted in figure 4(b).
3.3. Implementation of resolvent-based estimation
We employ recently developed resolvent-based estimation and control tools integrated within the compressible flow solver CharLES (Towne et al. Reference Towne, Bhagwat, Zhou, Jung, Martini, Jordan, Audiffred, Maia and Cavalieri2024; Jung et al. Reference Jung, Bhagwat and Towne2025), which was carefully designed to scale to large-scale problems such as a turbulent wake. Originally inspired by the incompressible solver implementation by Martini et al. (Reference Martini, Jung, Cavalieri, Jordan and Towne2022), the approach was extended by Jung et al. (Reference Jung, Bhagwat and Towne2025) to compressible flows using parallel algorithms written in C++, built directly upon the existing solver architecture. These resolvent-based tools are tightly integrated with the solver, enabling in situ computations during simulations. Moreover, our modular software design simplifies integration with other relevant software packages and ensures straightforward deployment in high-performance computing environments.
Flowchart of the resolvent-based estimation tool within a compressible flow solver. Dashed-line boxes represent the resolvent-based tool integrated into the existing computational fluid dynamics (CFD) solver (solid-line box). Blue boxes, including those outlined in red (offline run) and orange (build kernels) dashed lines, indicate offline processes, while black boxes with green solid outlines represent online (estimation) processes. Here, WH decomposition refers to Wiener-Hopf decomposition.
An overview of the implementation is shown in figure 5. Our compressible resolvent-based tool integrates external libraries such as PETSc (Balay et al. Reference Balay2019) and FFTW (Frigo & Johnson Reference Frigo and Johnson2005), enabling scalable parallel computation through in situ matrix operations with flow quantities during the training process. The tool constructs measurement matrices (
$\unicode{x1D63E}_{y}$
and
$\unicode{x1D63E}_{z}$
) using control volume indices from the CFD mesh. For cases involving high-rank targets (
$n_{z} \gg n_{y}$
), an online discrete Fourier transform (DFT) is used to minimise memory consumption while preserving time-resolved information (Schmidt & Towne Reference Schmidt and Towne2019; Farghadan et al. Reference Farghadan, Jung, Bhagwat and Towne2024, Reference Farghadan, Martini and Towne2025). Training data are obtained from LES, and the corresponding CSDs are used to generate the estimation kernels. The DFT parameters used to obtain converged kernels are described in Appendix B.
Importantly, our framework uses only a small number of spatial locations for both sensors and targets, which makes the data requirements compatible with realistic experimental set-ups. In this sense, the extension to experimental data is straightforward. That is, one replaces the LES time series with experimentally measured time series, while the estimation procedure itself remains unchanged and continues to avoid any direct handling of the unstable resolvent operator.
4. Results
We demonstrate the effectiveness of the causal resolvent-based approach for predicting turbulent fluctuations over a NACA0012 airfoil at a
$6^{\circ }$
angle of attack. Unlike earlier resolvent-based studies that employed non-causal estimation methods (Martini et al. Reference Martini, Cavalieri, Jordan, Towne and Lesshafft2020; Towne et al. Reference Towne, Lozano-Durán and Yang2020; Amaral et al. Reference Amaral, Cavalieri, Martini, Jordan and Towne2021; Amaral & Cavalieri Reference Amaral and Cavalieri2023), the causal estimator uses only present and past measurements and kernels optimised for this purpose via the Wiener–Hopf method (Martini et al. Reference Martini, Jung, Cavalieri, Jordan and Towne2022). This causality enforcement makes our method suitable for real-time applications. This work is the first numerical investigation explicitly applying the Wiener–Hopf method to enforce causality in turbulent flow estimation, extending previous efforts primarily focused on experimental flow-control studies (Audiffred et al. Reference Audiffred, Cavalieri, Brito and Martini2023, Reference Audiffred, Cavalieri, Maia, Martini and Jordan2024) and laminar-flow scenarios (Martini et al. Reference Martini, Jung, Cavalieri, Jordan and Towne2022; Jung et al. Reference Jung, Bhagwat and Towne2025).
As previously mentioned, our analysis centres on the spanwise-averaged, spanwise-Fourier and mid-span plane cases, as depicted in figure 6. The specific target locations analysed in the wake region are also highlighted in figure 6(a). The spanwise-averaged case provides a concise yet informative view of spatial structure in the spanwise direction, enabling a focused analysis of sensor positioning in the streamwise and cross-stream planes. We progressively incorporate more complex flow structures by examining the spanwise-Fourier modes, with particular attention to the most energetic mode at
$k_{z}=10\pi$
, and the flow on the mid-span plane to better capture spanwise effects in the estimation process.
Streamwise velocity fluctuation fields (
$u_{x}'/U_{\infty }$
) for (a) spanwise-averaged, (b) spanwise-Fourier mode
$k_{z}=10\pi$
and (c) mid-span plane cases. Target locations (
$z_{1}, z_{2}, z_{3}, z_{4}$
) are indicated by blue circles and positioned downstream of the trailing edge at streamwise coordinates
$x/L_{c} = [1.2,1.5,2.0,2.5]$
for a constant transverse position
$y/L_c = -0.11$
.
In this study, we compare the causal estimation performance with that of a TNC estimator. The fully non-causal estimation kernel cannot be used directly for real-time applications because its time support extends into the future. Any practical real-time implementation must therefore truncate the kernel in time so that only past and present measurements are used. As shown in (3.9), this truncation yields the TNC kernel (which is causal but not optimal) used in practice, whose performance can be directly compared with that of the optimal causal estimator.
Given the large number of potential sensor positions available on the three-dimensional airfoil surface, our investigation begins with selecting the most effective sensor locations based on single-sensor estimation error metrics and coherence analyses. Subsequently, we examine the power spectral density (PSD) at these selected sensor and target locations to better understand their frequency-domain characteristics. We then construct estimation kernels specifically designed for multi-sensor estimation, analysing their characteristics in both the time and frequency domains to gain insight into their role in estimation performance. Finally, we present multi-sensor estimation results, quantifying errors at both single target points and over extended target regions, to demonstrate the accuracy of the resolvent-based method.
4.1. Investigation of sensor placement
Sensor placement is guided by evaluating single-sensor estimation error metrics, similar to the approach in Jung et al. (Reference Jung, Bhagwat and Towne2025), as well as assessing coherence between sensor and target locations (Maia et al. Reference Maia, Jordan, Cavalieri, Martini, Sasaki and Silvestre2021; Audiffred et al. Reference Audiffred, Cavalieri, Maia, Martini and Jordan2024). The estimation error used throughout this paper is quantified by the metric (Martini et al. Reference Martini, Jung, Cavalieri, Jordan and Towne2022; Jung & Towne Reference Jung and Towne2024; Jung et al. Reference Jung, Bhagwat and Towne2025)
\begin{align} E=\frac {\sum _{i} \int (\boldsymbol{\tilde {z}}_{i}(t)-\boldsymbol{z}_{i}(t))^{2} {\rm d}t}{\sum _{i} \int (\boldsymbol{z}_{i}(t))^{2} {\rm d}t}, \end{align}
where
$\boldsymbol{\tilde {z}}_{i}$
and
$\boldsymbol{z}_{i}$
represent the estimated and true values at the
$i$
th target location, respectively. Additionally, coherence provides valuable insights into frequency-dependent estimation performance and is defined as (Towne et al. Reference Towne, Lozano-Durán and Yang2020)
\begin{align} \hat {\gamma }_{yz}(\omega )=\frac {|\unicode{x1D64E}_{yz}(\omega )|}{\sqrt {\unicode{x1D64E}_{yy}(\omega )}\sqrt {\unicode{x1D64E}_{zz}(\omega )}}, \end{align}
where
$\unicode{x1D64E}_{yz}$
denotes the CSD between sensor and target, and
$\unicode{x1D64E}_{yy}$
and
$\unicode{x1D64E}_{zz}$
denote their respective PSDs. The CSDs and PSDs are computed as described in (3.15).
Estimation error maps for the streamwise velocity fluctuation
$u_{x}'$
in single sensor–target configurations: (a,b) spanwise-averaged flow, (c,d) spanwise-Fourier mode for
$k_z=10\pi$
and (e,f) mid-span plane flow. Results using the TNC approach are shown in panels (a,c,e), while the causal estimation results are shown in panels (b,d,f).
Estimation errors of streamwise (a) and cross-streamwise (b) velocities as a function of the sensor location on the airfoil surface. The blue, orange and purple lines represent spanwise-averaged (S-A), mid-span (M-S) and spanwise-Fourier (S-F) (
$k_{z}=10\pi$
) flows, respectively. The solid diamond (
$\diamond$
) and dashed square (
$\square$
) lines denote the causal and TNC approaches, respectively.
Figure 7 presents spatial maps of estimation errors for a single sensor estimating streamwise velocity fluctuations at a single target location. Results are compared between the TNC (panels a,c,e) and causal (panels b,d,f) estimation approaches, covering three different flow representations: spanwise-averaged case (panels a,b), spanwise-Fourier modes at wavenumbers
$k_{z}=10\pi$
(panels c,d) and mid-span plane case (panels e,f). For the spanwise-Fourier case, we focus specifically on the wavenumber
$k_{z}=10\pi$
, as it represents the lowest frequency and highest energy modes identified by previous resolvent gain studies (Yeh & Taira Reference Yeh and Taira2019; Towne et al. Reference Towne2023).
Sensor locations are characterised using the coordinate
$s$
, defined such that positive values correspond to the suction surface and negative values to the pressure surface of the airfoil, with the magnitude representing the normalised streamwise position (
$x/L_{c}$
). Overall, the causal estimator improves upon the TNC estimator across the tested configurations, with the most pronounced gains occurring when the non-causal kernel has significant energy for
$\tau \lt 0$
. For the spanwise-Fourier case at
$k_{z}= 10\pi$
, the improvement is small because the dominant kernel peak occurs for
$\tau \gt 0$
. This advantage is particularly noticeable near the trailing-edge region (
$1 \lt x/L_{c} \lt 1.5$
) and prominently so when sensors are positioned on the front surfaces of the airfoil (
$-0.5 \lt s \lt 0.5$
). The performance improvement is more pronounced for the spanwise-averaged flow compared with the mid-span flow. This difference arises primarily due to the characteristics of the laminar region spanning
$-0.5 \lt s \lt 0.5$
on the surface of the airfoil, where larger spanwise structures exist, resulting in less information about turbulent fluctuations in the wake, diminishing estimation accuracy. The detrimental impact of laminar regions on estimation accuracy is especially evident in the spanwise-Fourier modes. For the spanwise-Fourier modes, sensors located within the laminar region (
$-0.5 \lt s \lt 0.5$
), particularly on the pressure surface (
$-0.5 \lt s \lt 0$
), exhibit poor estimation performance. This effect is notably stronger at low wavenumbers, e.g.
$k_{z}=10\pi$
, likely due to its larger associated spanwise structures compared with higher wavenumbers, e.g.
$k_{z}=20\pi$
. However, downstream of the laminar-to-turbulent transition, estimation performance significantly improves due to the emergence of large-scale turbulent structures within the wake that are better correlated with the measurements.
Figure 8 presents the average estimation errors over the wake region as a function of sensor location
$s$
, comparing spanwise-averaged (S-A), spanwise-Fourier (S-F) and mid-span plane (M-S) flows. These results are shown separately for the streamwise velocity fluctuations (
$u_{x}'$
, panel a) and the cross-streamwise velocity fluctuations (
$u_{y}'$
, panel b). Here, we report the average error over a set of targets in the wake that are uniformly distributed along the line segment
$1 \lt x/L_{c} \lt 2.5$
and
$y/L_{c} = -0.11$
. The resolvent-based kernels, which inherently include the linearised Navier–Stokes operator, estimate both velocity components from shear-stress sensor measurements, although the errors can be substantial in some cases. This effectiveness implies potential applicability of our method for estimating other relevant flow quantities that are derivatives of the velocity field, such as the vorticity field. Sensors placed on the pressure surface generally yield higher accuracy in spanwise-averaged and mid-span flows. Interestingly, the laminar region on the suction surface also performs well for estimation in spanwise-averaged flows, which contrasts with the results shown previously in figure 7. This discrepancy arises because the targets analysed in figure 8 span both the near-wake and far-wake regions downstream of the trailing edge, thereby incorporating regions of different coherence characteristics. Additionally, we observe a deterioration in estimation accuracy for sensors within the laminar separation bubble region (
$0.4 \lt s \lt 0.7$
), although accuracy improves progressively as the sensor location approaches the trailing edge.
Overall, superior estimation performance is achieved for the spanwise-averaged flow compared with the mid-span flow, with the spanwise-Fourier modes exhibiting intermediate performance. This is primarily due to the limitations inherent in using only a single sensor. The S-F and M-S cases exhibit smaller-scale dynamics in the flow, as shown in figure 6. Using a single sensor, the dominant flow structures are already captured within the finite non-causal window, so enforcing causality provides only a limited additional benefit. In other settings where the relevant dynamics is more spatially extensive, the advantage of the causal estimator is more pronounced, as in the S-A case. This improvement becomes more pronounced as we use more sensors. These results are consistent with the multi-sensor estimation results we will investigate in § 4.4. Based on these observations, sensor positions near the trailing edge (
$0.5 \lt x/L_{c} \lt 1$
) and on the upstream suction (
$0 \lt x/L_{c} \lt 0.3$
) and pressure (
$-1 \lt x/L_{c} \lt -0.7$
) surfaces are identified as particularly effective locations for achieving high-performance resolvent-based estimation.
Coherence of the streamwise velocity with respect to the target
$z_{1} (x/L_{c}=1.2)$
for the spanwise-averaged, mid-span and spanwise-Fourier (
$k_{z}=10\pi$
) cases.
In figure 9, we examine coherence of the streamwise velocity fluctuations between the sensor and the wake regions for the spanwise-averaged, mid-span and spanwise-Fourier (
$k_{z}=10\pi$
) cases at the sensor locations where single-sensor estimation was effective, as demonstrated in figures 7 and 8. These locations include regions on the front suction side (
$0 \lt x/L_{c} \lt 0.3$
), the rear pressure side (
$-1 \lt x/L_{c} \lt -0.7$
) and the mid-chord area (
$0.5 \lt x/L_{c} \lt 1$
), as illustrated in figure 9. Since our goal in this coherence study is to gain physical insights into coherence between sensor locations and wake fluctuations, we focus our analysis on the near-wake target position at
$x/L_{c}=1.2$
. We first examine the spanwise-averaged and mid-span cases in figures 9(a) and 9(b), and then study the spanwise-Fourier case in figure 9(c). We confirmed that coherence trends at farther wake positions were similar but with reduced magnitudes. Specifically, the sensors representing the laminar regions were placed at 1 (suction side, laminar region
$0 \lt x/L_{c} \lt 0.4$
) and 6 (pressure side, laminar region
$-1 \lt x/L_{c} \lt 0$
), as shown in figure 9. Within each laminar region, the coherence patterns remain similar, differing only in magnitude. Notably, the spanwise-averaged case exhibits stronger coherence at frequencies of interest (
$St_{\alpha }\lt 0.2$
) compared with the mid-span case. This enhanced coherence results from spatially larger coherent fluctuations associated with spanwise averaging, resulting in stronger coupling between sensor and target fluctuations. Interestingly, laminar regions (sensors 1 and 6) exhibit significant coherence with the near-wake target, which is consistent with previous studies showing that upstream laminar disturbances can be strongly linked to downstream wake fluctuations (Zaman Reference Zaman1987; Dovgal, Kozlov & Michalke Reference Dovgal, Kozlov and Michalke1994; Hain, Kähler & Radespiel Reference Hain, Kähler and Radespiel2009; Yarusevych, Sullivan & Kawall Reference Yarusevych, Sullivan and Kawall2009; Yeh & Taira Reference Yeh and Taira2019; Jung et al. Reference Jung, Bhagwat and Towne2025). Furthermore, sensor 5, located in the turbulent region, demonstrates notable coherence with the target across a wide frequency range. Among middle locations (sensors 2, 3 and 4), each sensor displays pronounced coherence at the frequency bands of interest (
$St_{\alpha }\lt 0.2$
). These locations are empirically selected to capture dominant frequencies at the target, which are crucial for kernel-based amplification in the frequency domain.
We further investigated coherence between the sensor and the wake regions with respect to the streamwise velocity fluctuations for the spanwise-Fourier mode at
$k_{z}=10\pi$
, as shown in figure 9(c). In particular, we are interested in verifying whether coherence is indeed weak in the region
$-0.7 \lt x/L_{c} \lt 0$
, as shown in the single-sensor estimation result in figure 7. To assess this, we examined coherence at sensor location 7, placed at
$x/L_{c} = -0.3$
, as indicated by the yellow solid line in figure 9(c). Although the von Kármán vortex shedding and other dominant frequencies are evident for
$St_{\alpha }\lt 0.2$
, sensor 7 does not exhibit significant coherence in this frequency range. In contrast, sensor locations 4 (purple), 5 (blue) and 6 (red) display prominent coherence peaks, consistent with effective estimation performance. Figure 9(c) further indicates that sensors 5 and 6 exhibit increased coherence at a spanwise wavenumber of
$10\pi$
. Ultimately, we identify sensors 1, 5 and 6 as the most effective for estimation purposes, with an additional selection of sensor 4 from among sensors 2, 3 and 4.
4.2. Power spectral density of the sensor and target readings
Figure 10 presents the PSD of the sensor and target signals for the spanwise-averaged, spanwise-Fourier and mid-span plane flows. Panels (a) and (b) highlight three dominant frequencies evident in the signals at sensor positions
$y_{3}$
,
$y_{4}$
and
$y_{5}$
, as well as the target position
$z_{1}$
. The energy associated with these dominant frequencies decreases downstream, as demonstrated by the PSD at target positions
$z_2$
,
$z_3$
and
$z_4$
in panel (b). It is well established that laminar separation bubble (LSB) detachment over the suction surface promotes the Kelvin–Helmholtz instability, resulting in vortex shedding characterised by lower-frequency oscillations (
$St_{\alpha }\lt 0.12$
) (Ducoin et al. Reference Ducoin, Loiseau and Robinet2016). These lower-frequency peaks are clearly captured in panels (a) and (b) for sensor positions
$y_{3}$
,
$y_{4}$
,
$y_{5}$
and the near-wake target
$z_{1}$
. In contrast, the turbulent region on the airfoil surface and near wake exhibit dominant von Kármán vortex-shedding frequencies within the higher-frequency band (
$0.152 \leqslant St_{\alpha } \leqslant 0.211$
) (Yeh & Taira Reference Yeh and Taira2019).
The PSD of streamwise velocity fluctuations
$u_{x}'$
for sensor readings and target readings, obtained from the spanwise-averaged (a) and (b), spanwise-Fourier (
$k_z=10\pi$
) (c) and (d) and mid-span plane (e) and (f) cases presented in figure 6. In (a), the dashed lines indicate positions
$y_1$
and
$y_6$
, corresponding to the laminar-flow region.
Figures 10(c) and 10(d) respectively illustrate the PSD of sensor and target measurements for the spanwise-Fourier mode at the wavenumber
$k_{z}=10\pi$
. The PSD computations utilise the same parameters as those used for the spanwise-averaged case. Compared with the spanwise-averaged (a) and mid-span plane (e) cases, the PSD at the sensor and target locations (c,d) does not exhibit distinct dominant peaks, consistent with the coherence results shown in figure 9(c). As observed in Yeh & Taira (Reference Yeh and Taira2019) and Towne et al. (Reference Towne2023), the response modes for spanwise-Fourier modes at
$k_{z}=10\pi$
appear predominantly slightly away from the airfoil surface. Thus, the spanwise-Fourier modes at the sensor locations are relatively weak, leading to weaker PSD values, as illustrated in figure 10(c). Modes at lower wavenumbers correspond to larger, more energetic spanwise structures, whereas modes at higher wavenumbers represent smaller, less energetic structures. We confirmed that the primary difference between the wavenumbers
$k_z=10\pi$
and
$20\pi$
is the magnitude of their PSD values, while their PSD patterns remain similar.
Figures 10(e) and 10(f) display the PSD of sensor and target signals in the mid-span plane. Similar to the spanwise-averaged case, the laminar-flow region (
$y_{1}$
and
$y_{6}$
) exhibits lower energy fluctuations. Additionally, energy progressively dissipates as the target moves further downstream. Dominant frequencies in the approximate range
$0.1 \lt St_{\alpha } \lt 0.2$
are distinctly visible in the near-wake targets (
$z_{1}$
and
$z_{2}$
), as clearly illustrated in panel (f).
Overall, these PSD comparisons clarify why the estimator performance varies across sensor–target configurations. For the spanwise-averaged and mid-span cases, sensors located in regions where the dominant shedding frequencies are strong and well defined (e.g. near the LSB-induced and von Kármán peaks) provide rich spectral information that is well aligned with the energetic content of the targets, supporting more accurate estimates. In contrast, for the spanwise-Fourier mode at
$k_z = 10\pi$
, the weaker and less distinct spectral content at the sensor locations provides less useful information to reconstruct the target dynamics, which is consistent with the reduced coherence and lower estimation performance observed for these configurations.
4.3. Estimation kernels
Next, we examine the estimation kernels, building on our previous studies (Martini et al. Reference Martini, Cavalieri, Jordan, Towne and Lesshafft2020, Reference Martini, Jung, Cavalieri, Jordan and Towne2022; Jung et al. Reference Jung, Bhagwat and Towne2025), where we initially explored kernels within the resolvent-based framework. In this paper, we derive the kernels using a data-driven approach that accounts for the statistical effects of flow nonlinearity, implicitly incorporating the CSD matrix
$\kern2.5pt \hat {\kern-2.5pt \unicode{x1D641}}$
, which acts as a forcing in the linear dynamics. Note that when white noise forcing is applied, the estimation kernels correspond to the transfer functions of the Kalman filter (Kalman Reference Kalman1960), as shown by Martini et al. (Reference Martini, Jung, Cavalieri, Jordan and Towne2022). Incorporating the forcing matrix
$\kern2.5pt \hat {\kern-2.5pt \unicode{x1D641}}$
is suitable for this study, as our application involves strongly nonlinear turbulent flows. Modelling the nonlinear terms for estimation leads to better performance.
A key characteristic of the estimation kernels is the presence of ‘peak points’, defined as the times in the kernel where the amplitude attains its maximum value. These peak points indicate the time delay at which the sensor signal is most strongly correlated with the target, and thus provide an estimate of the travel time of the dominant hydrodynamic disturbances from the sensor location to the target. Since the estimation kernels are based on the compressible linearised Navier–Stokes operator, they also capture acoustic waves, resulting in secondary peaks in addition to the hydrodynamic waves (Jung, Bhagwat & Towne Reference Jung, Bhagwat and Towne2023; Towne et al. Reference Towne, Bhagwat, Zhou, Jung, Martini, Jordan, Audiffred, Maia and Cavalieri2024). In other words, the interaction between hydrodynamic and acoustic effects is already embedded in the structure of the kernels. As a result, the input–output statistics from which the kernels are obtained inherently capture the coupling between the two, ensuring that the estimated dynamics reflects both hydrodynamic and acoustic contributions and their mutual interactions. In an online implementation, however, the estimator tends to underestimate the target values during the initial transient – before a sufficiently long history of past measurements has been accumulated such that the time lag reaches the peak of the kernel. This underestimation occurs because part of the kernel’s support lies outside the available time window, so the corresponding information is effectively missing.
Another important characteristic is whether there are significant values in the non-causal part of the kernel (
$\tau \lt 0$
), including the peak point (Jung Reference Jung2024). This information is crucial for estimation, but it may be truncated when the estimator is applied to the convolution function if it resides in the non-causal part of the kernels. This truncation can lead to a significant difference in the performance between the causal and TNC approaches. From a numerical perspective, the estimation kernels function by amplifying the frequency signals of the sensor measurements through multiplication with the kernels to approximate the target signals, i.e.
$\hat {\tilde {\boldsymbol{z}}}(\omega )=\kern2.5pt \hat {\kern-2.5pt \unicode{x1D64F}}_{z}(\omega )\kern2pt \hat {\kern-2pt \boldsymbol{y}}(\omega )$
. Therefore, analysing the estimation kernels in the frequency domain alongside the PSD of the sensor and target readings in § 4.2 provides valuable numerical insights into accurate estimation.
We examine four kernels, specifically within the single-input and multiple-output configuration (sensor:
$y_{1}$
and targets:
$z_{1}$
,
$z_{2}$
,
$z_{3}$
,
$z_{4}$
), as illustrated in figure 11. This selection is based on the observation that when the sensor is positioned in the laminar region within the effective locations, the kernel’s peaks are distinct, with less noise in other frequencies, showing a clear correlation between the sensor and the target signals. In figure 11, prominent peak points are observed in the non-causal kernels (a,c,e,g) at [
$\tau _{{peak},z_{1}}^{},\tau _{{peak},z_{2}}^{},\tau _{{peak},z_{3}}^{},\tau _{{peak},z_{4}}^{}] = [-0.36, 0.17, 0.68, 1.25]$
, respectively. These points represent the primary hydrodynamic wave travel times between the sensor and the targets. The case in figure 11(a) serves as a good example of how a causal approach using Wiener–Hopf decomposition can enhance estimation accuracy. When this kernel is applied to real-time estimation for evaluating the convolution function, the peak at
$\tau =-0.36$
, which contains critical travel information of the fluctuations, is truncated. The TNC kernel in the frequency domain (magenta dashed line) is notably lower in magnitude compared with the causal kernel (blue line), as shown in figure 11(b). It is important to note that the kernels amplify the sensor measurements at dominant frequencies (
$0.1 \lt St_{\alpha } \lt 0.2$
).
Estimation kernels between sensor
$y_{1}$
and targets
$z_{1}$
(a,b),
$z_{2}$
(c,d),
$z_{3}$
(e,f) and
$z_{4}$
(g,h) for the spanwise-averaged case. Kernels are shown in the time domain (left column) and frequency domain (right column). The black line represents the non-causal kernel, the blue line represents the causal kernel obtained via Wiener–Hopf decomposition, and the magenta dashed line (when shown) indicates the TNC kernel.
Imposing sensor noise is necessary to ensure that the inversion of the PSD is well posed (Martini et al. Reference Martini, Cavalieri, Jordan, Towne and Lesshafft2020, Reference Martini, Jung, Cavalieri, Jordan and Towne2022). The sensor-noise covariance
$\hat {\kern-1.5pt \unicode{x1D649}}$
is assumed to be diagonal, i.e. the measurement noise is modelled as uncorrelated (white) with a specified variance, similar to other optimal estimation and control methods (Chevalier et al. Reference Chevalier, Hœpffner, Bewley and Henningson2006; Schmid & Sipp Reference Schmid and Sipp2016). This choice reflects the common situation in which only partial information on the noise statistics is available, and detailed noise correlations between sensors are not known. In this study, the noise level, which is used to regularise the kernels, is set at
$10^{-2}$
of the maximum PSD value of the sensor measurements, consistent with the approach used in previous work (Martini et al. Reference Martini, Jung, Cavalieri, Jordan and Towne2022; Jung et al. Reference Jung, Bhagwat and Towne2025). To assess the sensitivity of the estimation to the assumed noise level, we performed additional tests in which the noise variance was varied over a range wider than the value of
$10^{-2}$
. The estimation error increases with increasing noise level and decreases as the noise level is reduced. These observations suggest that, while the method is robust to moderate variations in noise level, improved characterisation of sensor noise, particularly correlated noise, could further enhance performance. Incorporating partially known or unknown noise correlations into
$\kern1.5pt \hat {\kern-1.5pt \unicode{x1D649}}$
is therefore an important topic for future work.
Estimation kernels between sensor
$y_{5}$
and target
$z_{1}$
for the spanwise-Fourier mode at
$k_{z}=10\pi$
, shown in the time (a) and frequency (b) domains. Sensor and target positions are indicated in figure 6. The black line represents the non-causal kernel, and the blue line represents the causal kernel obtained via Wiener–Hopf decomposition.
Figure 12 shows the kernels for the spanwise-Fourier modes at
$k_{z}=10\pi$
. Based on figure 7, sensor placement in the laminar region is less effective; hence, we focus on kernels corresponding to sensor
$y_{5}$
, rather than
$y_{1}$
, and target
$z_{1}$
, consistent with locations examined previously in figure 11. Because the spanwise-Fourier-transformed state is complex valued, the computational cost for estimation increases by a factor of four. Additionally, the kernels are no longer symmetric at
$St_{\alpha }=0$
, as shown in figure 12. The peak points are clearly visible (see figure 12
a) at
$\tau \gt 0$
, indicating the dominance of hydrodynamic waves. Consequently, performance differences between causal and TNC estimation methods are minimal. The values of the kernels are largest in the range of
$|St_{\alpha }|\lt 0.2$
, effectively amplifying sensor measurements to produce accurate estimates.
Estimation kernels between sensors
$y_{1}$
(a,b),
$y_{5}$
(c,d) and
$y_{6}$
(e,f) and target
$z_{1}$
for the mid-span plane case. Kernels are shown in the time domain (left column) and frequency domain (right column). The black line represents the non-causal kernel, and the blue line represents the causal kernel obtained via Wiener–Hopf decomposition.
Finally, we present the kernels for the mid-span plane flow in figure 13. Sensors positioned in both the laminar (
$y_{1}$
) and turbulent regions (suction side:
$y_{5}$
; pressure side:
$y_{6}$
) are shown. The prominent peak points in the kernels between sensors (
$y_{1}$
,
$y_{5}$
) and target
$z_{1}$
closely match those observed in the spanwise-averaged and spanwise-Fourier flows, respectively, as depicted in figures 11 and 12 for the kernels between
$y_1$
and
$z_1$
and
$y_5$
and
$z_1$
, respectively. The kernel between
$y_{6}$
and
$z_{1}$
also aligns well with the other cases, though it is not explicitly shown here. However, the kernels for the mid-span flow exhibit significant amplitudes over a larger range of
$\tau$
, indicating greater difficulty in achieving temporal convergence. A notable difference appears in figure 13(a), where the kernel peak lies in the positive domain (
$\tau \gt 0$
), unlike the corresponding kernel in the spanwise-averaged flow, which peaks in the negative domain. This discrepancy may be attributed to the fact that the spanwise-averaged case integrates the flow across multiple spanwise positions, capturing coherent large-scale structures that extend significantly in the spanwise direction. These large-scale structures are more capable of influencing upstream conditions (non-causal behaviour), appearing as peaks in the non-causal part (
$\tau \lt 0$
) of the kernel. In contrast, the mid-span plane case represents a single spanwise slice, emphasising smaller-scale, localised spatial structures. Such structures typically have limited upstream influence since their acoustic disturbances are weaker, more rapidly dissipated and less coherent over larger spanwise distances. Thus, such conditions enhance convective flow characteristics in the mid-span plane case, thereby improving the accuracy when applying TNC estimation to the mid-span plane flow. Additionally, we observe that the kernel between
$y_{4}$
and
$z_{1}$
is similar to that between
$y_{5}$
and
$z_{1}$
within the mid-span flow, so it is omitted from this figure. Both causal and TNC kernels for targets further downstream (
$z_{2}$
,
$z_{3}$
,
$z_{4}$
) closely resemble those in the spanwise-averaged flow (figure 11), with their peak points consistently located in the positive
$\tau$
domain. The mid-span case includes the spanwise-averaged mode (
$k_{z}=0$
), the specific harmonic
$k_{z}=10\pi$
and many additional spanwise scales. Consequently, the resulting kernel is broadband, and the associated estimation problem is more challenging.
4.4. Multi-sensor estimation
By selecting sensor locations
$y_{1}$
,
$y_{4}$
,
$y_{5}$
and
$y_{6}$
(notations from figures 9
a and 9
b) based on the investigation of the most effective locations using single-sensor estimation and coherence in § 4.1, we show the estimation errors for spanwise-averaged and mid-span flows along the trailing edge as additional sensors are introduced, as shown in figure 14. The use of sensor
$y_{5}$
(blue) alone does not yield significant improvement over TNC estimation. However, adding an additional sensor (
$y_{6}$
) markedly enhances estimation accuracy (yellow). With these two sensors, TNC estimation performs poorly near the trailing edge, whereas causal (C) estimation remains accurate. The inclusion of a third sensor
$y_{4}$
(purple) alongside
$y_{5}$
and
$y_{6}$
results in a slight improvement. The best estimation results are achieved using all four sensors
$y_{1}$
,
$y_{4}$
,
$y_{5}$
,
$y_{6}$
(green). These findings demonstrate greater accuracy than the use of six sensors uniformly distributed on the pressure side of the airfoil, as reported in our previous work (Jung & Towne Reference Jung and Towne2024). Additionally, the four-sensor TNC case (green dashed) can lie slightly above the two- and three-sensor cases (yellow and purple dashed). This behaviour reflects a limitation of the TNC estimation. Specifically, TNC applies a finite time-window truncation to the optimal non-causal multi-sensor kernel. Since each sensor has a distinct optimal correlation structure, a fixed truncation window (
$\tau \lt 0$
) can remove different portions of each sensor’s kernel unevenly, thereby distorting the joint multi-sensor weighting and producing a small increase in error. In contrast, the causal estimator does not exhibit this behaviour because it enforces a physically consistent causal kernel structure. In the mid-span case, the improvement from one to two sensors is larger than from two to three sensors. This is because the second sensor typically captures most of the additional independent information about the convecting structures (e.g. providing a second view that helps resolve lag ambiguity), whereas additional sensors tend to offer lower incremental coherence with the target for the present sensor spacing. Consequently, the benefit of adding sensors beyond two shows diminishing returns.
Estimation error along the trailing edge for spanwise-averaged and mid-span flows using one sensor (
$y_{5}$
), two sensors (
$y_{5}$
,
$y_{6}$
), three sensors (
$y_{4}$
,
$y_{5}$
,
$y_{6}$
) and four sensors (
$y_{1}$
,
$y_{4}$
,
$y_{5}$
,
$y_{6}$
). Results for the spanwise-Fourier case (
$k_{z}=10\pi$
) are not shown because the causal and TNC estimators yield nearly identical errors for this mode.
Estimation of
$u_x'$
as a function of time using four sensors (
$y_1$
,
$y_4$
,
$y_5$
,
$y_6$
): spanwise-averaged (panels a,d,g,j), mid-span plane (panels b,e,h,k) and spanwise-Fourier cases at
$k_z = 10\pi$
(panels c,f,i,l). Target locations are defined by the coordinates
$[x/L_c,\; y/L_c]$
: panels (a–c)
$z_{1}=[1.2, -0.11]$
, (d–f)
$z_{2}= [1.5, -0.11]$
, (g–i)
$z_{3} = [2.0, -0.11]$
, and (j–l)
$z_{4}= [2.5, -0.11]$
, as illustrated in the top figure. The true data from LES are shown in black, causal estimations in blue and TNC estimations as magenta dashed lines. The causal estimation error is indicated in red in the bottom right corner of each panel.
The LES and estimated snapshots of
$u_{x}'$
across the extended wake region for the spanwise-averaged (panels a,c,e,g,i,k) and mid-span plane (panels b,d,f,h,j,l) flows, utilising four sensors (
$y_{1},y_{4},y_{5},y_{6}$
). The times
$t_1$
,
$t_2$
and
$t_3$
were chosen to represent three distinct phases based on the dominant vortex-shedding frequency.
Figure 15 shows the time-series estimation of streamwise velocity fluctuations using four sensors (
$y_{1}, y_{4}, y_{5}, y_{6}$
) for spanwise-averaged (panels a,d,g,j), mid-span plane (panels b,e,h,k) and spanwise-Fourier (panels c,f,i,l) cases. The causal estimation error, computed using (4.1), is indicated in the bottom right corner of each panel. As expected, the estimation error increases with greater distance between the sensor and the target location. The Wiener–Hopf approach significantly improves estimation accuracy for targets located closer to the trailing edge, as evident in figure 15(a). The difference in accuracy between TNC and causal methods diminishes as the target moves further downstream. Among the investigated flow types, the spanwise-averaged flow, characterised by less chaotic spatial structures, is estimated most accurately. For the mid-span plane flow, presented in figure 15 (panels b,e,h,k), the fluctuations are more broadband and contain finer-scale variability than the spanwise-averaged case. However, despite this increased complexity, the broader trends in the velocity fluctuations over time are still effectively captured. Although the sensor may not fully capture the detailed footprint of turbulence, the estimation remains effective. These results are consistent with the trends shown in figure 14.
Figure 15 (panels c,f,i,l) displays time-series estimations of the real part of
$u_{x}'$
for the spanwise-Fourier modes with
$k_{z}=10\pi$
. As anticipated, estimating higher-
$k_{z}$
content is generally more challenging due to reduced spatial coherence and weaker wall footprints; however, a systematic study across spanwise wavenumbers is beyond the scope of the present work. Kernels derived from resolvent operators at frequencies associated with lower energetic gains provide lower estimation accuracy. Since the TNC and causal estimations for the spanwise-Fourier modes yield similar results, the TNC estimations are not shown explicitly in figure 15 for these cases.
Figure 16 shows three random snapshots of the causal estimation in the extended wake region for spanwise-averaged and mid-span plane cases, using the same sensor configuration as in figure 15. The estimation errors for the wake region, calculated using metric (4.1), are approximately 0.254 and 0.324 for causal estimation, and 0.352 and 0.389 for TNC estimation in the spanwise-averaged and mid-span plane cases, respectively. Although the small-scale structures are not fully captured due to the limitations of the sensors in detecting turbulence, the causal resolvent-based estimator effectively estimates the wake fluctuations.
Overall, the estimation framework shows good agreement with the reference fields across the tested configurations, particularly in capturing the dominant flow structures and overall trends. This performance is notable given that the approach relies solely on a limited number of real-time surface-based measurements. At the same time, the normalised errors remain of the order of 0.3 or higher for most mid-span flow cases, as shown in figures 14(b) and 15 (panels b,e,h,k), and the fully three-dimensional mid-span flow scenario exhibits some smoothing and a partial loss of small-scale features in figure 16 (panels d,h,l). These effects reflect the inherent limitations of reconstructing high-resolution, complex three-dimensional flow structures from sparse surface data. Despite these limitations, the results indicate that the method provides useful quantitative estimates and reliable qualitative reconstructions of the main flow features under practical sensing constraints.
5. Conclusions
In this study, we developed and demonstrated a causal resolvent-based framework for estimating turbulent velocity fluctuations in the wake of a spanwise-periodic NACA0012 airfoil at moderate Reynolds number
${Re}_{L{c}} = 23\,000$
,
$ \textit{Ma}_{\infty } \equiv U_{\infty }/c_{\infty } = 0.3$
and an angle of attack of
$\alpha = 6^\circ$
. Building on previous research focused primarily on a laminar and globally stable flow scenario (Jung et al. Reference Jung, Bhagwat and Towne2025), our approach extended resolvent-based estimation into the more complex regime of a turbulent, globally unstable wake flow.
To address the key challenges inherent in this turbulent regime, including global instabilities, broadband multiscale dynamics and large datasets, we employed a data-driven approach. This allowed construction of resolvent-based kernels from CSDs directly obtained via LESs, thereby implicitly accounting for the coloured-in-time statistical properties of the nonlinear forcing through the empirically determined correlations. A key improvement was the adoption of the Wiener–Hopf method to enforce causality (Martini et al. Reference Martini, Jung, Cavalieri, Jordan and Towne2022), which significantly enhanced real-time estimation accuracy compared with traditional TNC estimation, with the most substantial gains observed in the multi-sensor configurations. To handle computational complexity arising from high-dimensional datasets, our methodology utilised scalable parallel algorithms integrated within an existing compressible flow solver (Towne et al. Reference Towne, Bhagwat, Zhou, Jung, Martini, Jordan, Audiffred, Maia and Cavalieri2024; Jung et al. Reference Jung, Bhagwat and Towne2025).
Our results illustrated successful causal estimation across spanwise-averaged, spanwise-Fourier and mid-span plane representations of the turbulent wake, using only limited shear-stress measurements from the airfoil surface. While using a single sensor yielded only modest improvements over TNC estimation, adding more sensors substantially improved performance by better capturing the turbulent footprint through correlations in wall shear stress. By strategically analysing single-sensor estimation errors and coherence between sensor and target locations, we identified practical sensor-placement strategies that significantly enhanced estimation fidelity. Ultimately, our approach provided accurate, interpretable and scalable predictions of an extended portion of the wake, capturing dominant coherent structures crucial for both scientific understanding and potential flow-control applications.
Future research directions suggested by this study include developing formal optimisation methods for sensor placement and incorporating three-dimensional sensor configurations to further enhance estimation accuracy. The success of the resolvent-based estimation framework suggests that resolvent-based closed-loop flow-control strategies could help mitigate turbulent fluctuations, reduce drag, increase lift and minimise aerodynamic noise. The present data-driven formulation requires prior time-resolved data for each sensor–target configuration, which limits its immediate generalisability across different flow conditions or geometries. Extending the framework to more broadly applicable, parameterised estimators remains an important direction for future work. Integrating advanced machine-learning techniques, such as neural networks, into the resolvent-based framework would further improve real-time computational efficiency and effectively handle larger datasets. Finally, experimental validation of the proposed resolvent-based estimation approach in practical wind-tunnel tests or flight scenarios is recommended to confirm its real-world applicability and robustness. By pursuing these avenues, our resolvent-based estimation framework has the potential to significantly advance the predictive modelling, real-time monitoring and control of complex turbulent aerodynamic flows.
Funding
This research was supported by the Air Force Office of Scientific Research (AFOSR) under grant number FA9550-20-1-0214. Computational resources were provided by the U.S. Department of Defense High-Performance Computing Modernization Program.
Declaration of interests
The authors report no conflicts of interest.
Appendix A. Wiener–Hopf method
We briefly review the theoretical and numerical Wiener–Hopf decompositions employed for enforcing causality within the resolvent-based estimation and control framework described by Martini et al. (Reference Martini, Jung, Cavalieri, Jordan and Towne2022), Audiffred et al. (Reference Audiffred, Cavalieri, Brito and Martini2023), Jung et al. (Reference Jung, Bhagwat and Towne2025).
A.1. Theoretical Wiener–Hopf decomposition
The Wiener–Hopf method (Noble Reference Noble1958) decomposes frequency-domain functions into causal (analytic in lower half-plane,
$(+)$
) and non-causal (analytic in upper half-plane,
$(-)$
) components. This decomposition is employed here to enforce causality on estimation and control kernels, following Daniele & Lombardi (Reference Daniele and Lombardi2007) and Martinelli (Reference Martinelli2009).
We consider the Wiener–Hopf problem
\begin{align} \hat {\unicode{x1D643}}(\omega ){\hat {\varGamma }}_{+}(\omega ) = {\hat {\varLambda }}_{-}(\omega )+\hat {\unicode{x1D642}}(\omega ), \end{align}
where
$\hat {\unicode{x1D643}}$
and
$\hat {\unicode{x1D642}}$
are known matrices, while
${\hat {\varGamma }}_{+}$
and
${\hat {\varLambda }}_{-}$
are unknown.
The Wiener–Hopf method involves additive and multiplicative factorisations, defined as
\begin{align} {\hat {\varGamma }}(\omega )=[{\hat {\varGamma }}(\omega )]_{-}+[{\hat {\varGamma }}(\omega )]_{+}, \quad {\hat {\varGamma }}(\omega )=({\hat {\varGamma }}(\omega ))_{-}({\hat {\varGamma }}(\omega ))_{+}, \end{align}
respectively. The explicit solution for the above problem is
\begin{align} {\hat {\varGamma }}_{+}(\omega )=\big [\hat {\unicode{x1D642}}(\omega )(\hat {\unicode{x1D643}}(\omega ))_{-}^{-1}\big ]_{+}(\hat {\unicode{x1D643}}(\omega ))_{+}^{-1}. \end{align}
A.2. Numerical Wiener–Hopf decomposition
Numerically, the additive factorisation is straightforward, while the multiplicative factorisation requires solving an integral equation. Following Martini et al. (Reference Martini, Jung, Cavalieri, Jordan and Towne2022), this integral equation is given by
\begin{align} \hat {\boldsymbol{x}}_i(\omega )+\frac {1}{2 \mathrm{i}}\mathcal{H}(\hat {\boldsymbol{x}}_i)(\omega )-\frac {1}{2 \mathrm{i}}\hat {\unicode{x1D642}}^{-1}(\omega )\mathcal{H}(\hat {\unicode{x1D642}}\hat {\boldsymbol{x}}_i)(\omega )=\hat {\unicode{x1D642}}^{-1}(\omega )\frac {\hat {\boldsymbol{w}}_{i,-}(\omega _0)}{\omega -\omega _0}, \end{align}
where
$\mathcal{H}$
denotes the Hilbert transform defined by
\begin{align} \mathcal{H}(\hat {\boldsymbol{x}})=P.V.\frac {1}{\pi }\int _{-\infty }^{\infty }\frac {\hat {\boldsymbol{x}}(u)}{\omega -u}\,\mathrm{d}u. \end{align}
The reader is referred to Noble (Reference Noble1958), Daniele & Lombardi (Reference Daniele and Lombardi2007), Martinelli (Reference Martinelli2009), Martini et al. (Reference Martini, Jung, Cavalieri, Jordan and Towne2022) for a detailed discussion of the Wiener–Hopf method. We numerically solve this integral equation using the generalised minimal residual iterative method, implemented directly within the CFD solver. This approach ensures computational efficiency and minimises the need for external post-processing (Jung et al. Reference Jung, Bhagwat and Towne2025).
Appendix B. Convergence of the data-driven estimation kernels
Constructing estimation kernels using a data-driven approach requires careful selection of DFT parameters. Specifically, the streaming Fourier transform (Schmidt & Towne Reference Schmidt and Towne2019) depends on several key parameters, including the number of frequencies (
$N_{\textit{freq}}$
), the number of data blocks (
$N_{\textit{blk}}$
) and the number of overlapping windows (
$N_{ovlp}$
). To reduce spectral leakage between blocks (
$N_{\textit{blk}}$
), we apply a Hamming window with an overlap of
$50\,\%$
(Trethewey Reference Trethewey2000). Table 2 summarises a convergence study evaluating how these parameters influence estimation accuracy within the spanwise-averaged flow under best sensor placement conditions. Based on this analysis, parameter set 6 was selected, and these DFT parameters were subsequently applied to the spanwise-Fourier and mid-span plane cases examined in this study.
Convergence of the estimation kernels with respect to the parameters for the DFTs used to compute CSDs.
Open access
