1. Introduction
Turbulent flows are fundamental to a wide range of industrial applications, including jets, combustion systems and aerodynamic surfaces like aircraft wings and control surfaces, wind turbine rotor blades, as well as turbomachinery components like compressor and turbine blades. In wall-bounded flows, such as those over airfoils, flow separation is a critical phenomenon: when a boundary layer encounters a sufficiently strong adverse pressure gradient (APG) or a geometric discontinuity, it can detach from the surface and form a free shear layer. Under certain conditions, the flow may subsequently reattach downstream, forming a separation bubble. The onset of flow separation has significant implications for the aerodynamic performance of the flow, by increasing drag, decreasing lift and introducing an unsteady dynamics caused by instabilities of the free shear layer or the separation bubble (Simpson Reference Simpson1989).
1.1. Turbulent separation bubble dynamics
Turbulent separation bubbles (TSBs) exhibit strong unsteadiness across a broad frequency range. This leads to fluctuating structural and thermal loads, as well as noise, in many engineering applications. Turbulent separation bubbles have been the subject of extensive study over the past five decades, across a wide range of geometries and flow configurations. They are typically classified as either geometry-induced or APG-induced. Geometry-induced TSBs arise at features such as backward-facing steps (Eaton & Johnston Reference Eaton, Johnston, Bradbury, Durst, Launder, Schmidt and Whitelaw1982) and rectangular leading edges (Kiya & Sasaki Reference Kiya and Sasaki1983; Cherry, Hillier & Latour Reference Cherry, Hillier and Latour1984). Adverse pressure gradient-induced TSBs occur, for instance, on airfoils (Wang & Ghaemi Reference Wang and Ghaemi2022; Sarras et al. Reference Sarras, Tayeh, Mons and Marquet2024), backward-facing ramps (Kaltenbach et al. Reference Kaltenbach, Fatica, Mittal, Lund and Moin1999; Weiss et al. Reference Weiss, Steinfurth, Chamard, Giani and Combette2022) or flat plates (Patrick Reference Patrick1987; Na & Moin Reference Na and Moin1998; Mohammed-Taifour & Weiss Reference Mohammed-Taifour and Weiss2016; Abe Reference Abe2017; Wu et al. Reference Wu, Meneveau and Mittal2020; Cura et al. Reference Cura, Hanifi, Cavalieri and Weiss2024). A related phenomenon is the appearance of stall cells, which are three-dimensional recirculation regions that occur near the trailing edge of two-dimensional airfoils (Winkelman & Barlow Reference Winkelman and Barlow1980; Sarras et al. Reference Sarras, Tayeh, Mons and Marquet2024). Turbulent separation bubbles also play a key role in high-speed flows, particularly in shock–boundary layer interactions, where shock waves induce separation and reattachment (Delery Reference Delery1985; Dussauge, Dupont & Debiève Reference Dussauge, Dupont and Debiève2006; Poggie et al. Reference Poggie, Bisek, Kimmel and Stanfield2015; Hao Reference Hao2023).
Over the years, different phenomena in TSBs have been associated with distinct frequency bands, which are typically referred to as shedding, flapping or breathing. A dominant frequency was first identified by Mabey (Reference Mabey1972), who proposed the Strouhal scaling
${\textit{St}}_{\textit{sep}} = f L_{\textit{sep}} / U_{\infty }$
, where
$L_{\textit{sep}}$
is the separation length and
$U_\infty$
is the free-stream velocity. Shedding occurs at
${\textit{St}}_{\textit{sep}} = 0.35{-}0.8$
and is linked to vortex roll-up in the shear layer (Eaton & Johnston Reference Eaton, Johnston, Bradbury, Durst, Launder, Schmidt and Whitelaw1982; Kiya & Sasaki Reference Kiya and Sasaki1983; Cherry et al. Reference Cherry, Hillier and Latour1984; Weiss, Mohammed-Taifour & Schwaab Reference Weiss, Mohammed-Taifour and Schwaab2015), often attributed to Kelvin–Helmholtz instability (Tenaud et al. Reference Tenaud, Podvin, Fraigneau and Daru2016). Low-frequency dynamics (
${\textit{St}}_{\textit{sep}} \lt 0.02$
) with substantial energy content was already reported by Eaton & Johnston (Reference Eaton, Johnston, Bradbury, Durst, Launder, Schmidt and Whitelaw1982), and low-frequency trailing-edge oscillations were first observed by Zaman, Mckinzie & Rumsey (Reference Zaman, Mckinzie and Rumsey1989). These low-frequency motions are now commonly described as either flapping or breathing. Flapping refers to shear-layer oscillations at
${\textit{St}}_{\textit{sep}} \approx 0.08{-}0.18$
, typically observed in geometry-induced TSBs (Largeau & Moriniere Reference Largeau and Moriniere2006; Pearson, Goulart & Ganapathisubramani Reference Pearson, Goulart and Ganapathisubramani2013; Fang & Wang Reference Fang and Wang2024), whereas breathing occurs at
${\textit{St}}_{\textit{sep}} \approx 0.01$
or below, and is often interpreted as a global expansion and contraction of the separation bubble (Weiss et al. Reference Weiss, Mohammed-Taifour and Schwaab2015; Mohammed-Taifour & Weiss Reference Mohammed-Taifour and Weiss2016; Borgmann et al. Reference Borgmann, Cura, Weiss and Little2024).
The breathing motion has drawn particular attention due to its large-scale, low-frequency nature, which poses significant measurement and modelling challenges. Capturing this dynamics requires long observation times and sufficient spatial coverage. The origin of this dynamics is the subject of ongoing research. Three key insights into the breathing mechanism have recently emerged: first, resolvent analyses show maximal amplification at finite spanwise wavenumbers, indicating a strong spanwise dependence of the breathing mode (Cura et al. Reference Cura, Hanifi, Cavalieri and Weiss2024; Sarras et al. Reference Sarras, Tayeh, Mons and Marquet2024; Fuchs et al. Reference Fuchs, Steinfurth, von, Jakob, Weiss and Oberleithner2026). Second, global linear stability analysis reveals a stationary eigenmode with distinctly elevated growth rate at these spanwise wavenumbers (Cura et al. Reference Cura, Hanifi, Cavalieri and Weiss2024; Sarras et al. Reference Sarras, Tayeh, Mons and Marquet2024; Fuchs et al. Reference Fuchs, Steinfurth, von, Jakob, Weiss and Oberleithner2026) that is expected to be the origin of the low-frequency dynamics. Third, the eigenmode has been linked to a centrifugal instability (Barkley, Gomes & Henderson Reference Barkley, Gomes and Henderson2002; Rodríguez et al. Reference Rodríguez, Gennaro and Juniper2013; Savarino, Sipp & Rigas Reference Savarino, Sipp and Rigas2025).
1.2. The Gaussian bump benchmark experiment
APG-induced flow separation occurring on smooth surfaces is sometimes termed smooth-body separation (SBS). This process is particularly difficult to predict and model, often times requiring extensive numerical and experimental research. As part of the recent research efforts to improve the capability of computational fluid dynamics (CFD) to accurately capture SBS in turbulent flows, a new benchmark case, termed the Boeing Gaussian Bump was introduced (Williams et al. Reference Williams, Samuell, Sarwas, Robbins and Ferrante2020; Gray et al. Reference Gray, Gluzman, Thomas, Corke, Lakebrink and Mejia2021). Initial experiments using the geometry have been performed at the University of Washington by Sarwas (Reference Sarwas2019) and Williams et al. (Reference Williams, Samuell, Sarwas, Robbins and Ferrante2020), followed by an experimental campaign at the University of Notre Dame by Gray et al. (Reference Gray, Gluzman, Thomas, Corke, Lakebrink and Mejia2021, Reference Gray, Gluzman, Thomas and Corke2022a
,
Reference Gray, Gluzman, Thomas, Corke, Lakebrink and Mejiab
, Reference Gray, Corke, Thomas, Gluzman and Straccia2023a
,
Reference Gray, Lakebrink, Thomas, Corke, Gluzman and Stracciab
) and Gray (Reference Gray2023). The results of the latter are to a large extent archived on the ‘NASA Turbulence Modeling Resource’ website. In this study we work with this dataset and, unless otherwise noted, ‘experimental’ refers to the experiments performed at the University of Notre Dame. An overview of the three-dimensional geometry is provided in figure 1(a). Here, selected PIV interrogation windows are overlaid and coloured by the streamwise mean velocity component
$\bar {u}$
from the experiment. All data shown in figure 1 are for
$\textit{Re}=2\times 10^6$
, where the flow is fully separated, with a TSB on the downstream face of the bump. Figure 1(b) shows the streamwise and spanwise principal views, including the same PIV windows. Wall skin friction streamlines from a wall-modelled LES (Iyer & Malik Reference Iyer and Malik2023a
,
Reference Iyer and Malikb
) are shown in figure 1(c) to provide an overview of the flow’s three-dimensional structure, and figure 1(d) shows a schematic of the breathing and shedding dynamics in the flow.
Overview of the Gaussian bump test case. (a) Three-dimensional schematic of the bump mounted in the wind tunnel test section. Selected particle image velocimetry (PIV) interrogation windows are overlaid and coloured by the streamwise mean velocity component,
$\bar {u}$
. The wind tunnel sidewalls are located at
$z=-0.5$
and
$z=0.5$
. (b) Geometry of the bump and measurement locations shown in the streamwise (top) and spanwise (bottom) principal views. Mean pressure measurement positions are indicated by blue dots and instantaneous pressure measurement positions by purple crosses. The PIV measurement regions are coloured by
$\bar {u}$
using the same colour scale as in (a). Instantaneous velocity data are available within the spanwise stereo PIV window (green border) and the streamwise PIV windows labelled ‘FOV 1’ and ‘FOV 2’ (black borders). (c) Wall skin friction streamlines obtained from a wall-modelled large-eddy simulation (LES) of the flow; adapted from Iyer & Malik (Reference Iyer and Malik2023a
,
Reference Iyer and Malikb
) with permission. (d) Qualitative illustration of breathing and shedding dynamics. All data are shown for the separated flow at Reynolds number
$\textit{Re} = 2\times 10^6$
.
A number of computational studies have been performed to benchmark the performance of various CFD approaches on the test case. Williams et al. (Reference Williams, Samuell, Sarwas, Robbins and Ferrante2020) investigated the capability of two- and three-dimensional Reynolds-averaged Navier–Stokes (RANS) models to reproduce the results from their experiment. They found that all of their simulations failed to accurately reproduce the surface-pressure coefficient in the separated flow region, whereas the pressure outside the separated region was generally well matched. They also found that the RANS models did not reflect the Reynolds number invariance of the surface-pressure coefficient for
$\textit{Re}\geqslant 2\times 10^6$
(based on the wind tunnel width and free-stream velocity) that they observed in the experiment. Gray et al. (Reference Gray, Lakebrink, Thomas, Corke, Gluzman and Straccia2023b
) performed RANS and delayed detached eddy simulation of the case at
$\textit{Re}=4\times 10^6$
. They found that the RANS model failed to capture the pressure in the separated region by predicting almost no flow separation. Delayed detached eddy simulation qualitatively captures the flow separation but quantitative deviations from the experiment remain. Their comparison also includes pressure coefficients from the experiment at the University of Washington at
$\textit{Re}=3.4\times 10^6$
(Sarwas Reference Sarwas2019; Williams et al. Reference Williams, Samuell, Sarwas, Robbins and Ferrante2020), which are almost indiscernible from the University of Notre Dame experiment at
$\textit{Re} = 4\times 10^6$
, albeit there is a difference in Reynolds number, highlighting the Reynolds number invariance of the separated flow. Zhou & Bae (Reference Zhou and Bae2024) investigated the sensitivities of wall-modelled LES with respect to the specific modelling choices and mesh parameters. They found that especially the model for the subgrid-scale stresses significantly impacts the solution in the separated flow region, and that there is a significant dependence on the mesh resolution up to the point where the resolution approaches that of wall-resolved LES. Direct numerical simulations (DNS) of the flow have been performed at
$\textit{Re}=10^6$
(Balin & Jansen Reference Balin and Jansen2021; Uzun & Malik Reference Uzun and Malik2021),
$\textit{Re}=2\times 10^6$
(Uzun & Malik Reference Uzun and Malik2022) and
$\textit{Re}=4\times 10^6$
(Uzun & Malik Reference Uzun and Malik2025). For computational reasons, these studies do not take the three-dimensional geometry of the bump into account. Instead, the simulations employ a spanwise periodic domain where the profile of the bump corresponds to that at the spanwise centreline in the experimental configuration. They found evidence for relaminarisation of the boundary layer in the acceleration zone at the upstream face of the bump at
$\textit{Re}=10^6$
. At this Reynolds number, the flow undergoes only very weak separation in the deceleration region downstream of the bump apex. At
$\textit{Re}=2\times 10^6$
and
$4\times 10^6$
, the relaminarisation is suppressed and the flow undergoes much stronger separation. Comparison of the surface-pressure and skin-friction coefficients from the DNS at these Reynolds numbers with experimental values from Williams et al. (Reference Williams, Samuell, Sarwas, Robbins and Ferrante2020) shows relatively good agreement (Uzun & Malik Reference Uzun and Malik2022, Reference von Saldern, Schmidt, Jordan and Oberleithner2025). However, at
$\textit{Re}=4\times 10^6$
, the shear layer is tilted significantly more towards the wall in the DNS compared with the experiment, which the authors attribute to the spanwise periodic configuration of the DNS that neglects three-dimensionality and tunnel endwall effects (Uzun & Malik Reference Uzun and Malik2025). Iyer & Malik (Reference Iyer and Malik2023b
) performed wall-modelled LES of the case in both, a spanwise periodic and a fully three-dimensional configuration involving the wind tunnel sidewalls. They similarly found the shear layer to be tilted more towards the wall in the spanwise periodic simulation whereas the shear layer in the fully three-dimensional simulation shows good agreement with the experiment. This effect has been linked to the interaction of the shear layer with two counter-rotating vortices, which, at the spanwise centreline, helps to lift the shear layer away from the bump surface (Uzun & Malik Reference Uzun and Malik2025). The effect becomes apparent from the surface streamline pattern of their three-dimensional LES, which is reproduced in figure 1(c). These studies highlight the challenges in accurately modelling the separated flow region downstream of the bump with reduced-order models and even with high-fidelity CFD, if the three-dimensional structure of the flow is not accounted for.
1.3. Motivation and objectives
In evaluating the performance of CFD computations, most studies of the Gaussian bump flow focus on mean-flow quantities such as time-averaged velocity fields, Reynolds stresses or wall-pressure distributions, while investigations of dominant coherent structures are lacking. However, some of the challenges in accurately modelling the separated flow over the bump likely arise from the dynamics of the TSB that forms downstream of the bump. This dynamics generates large-scale, three-dimensional coherent structures that are difficult to capture in simulations (Mohammed-Taifour & Weiss Reference Mohammed-Taifour and Weiss2016; Manohar et al. Reference Manohar, Williams, Martinuzzi and Morton2023; Borgmann et al. Reference Borgmann, Cura, Weiss and Little2024; Cura et al. Reference Cura, Hanifi, Cavalieri and Weiss2024). In particular, the observation that unsteady simulations match experimental results significantly better when the full span of the wind tunnel is resolved motivates an investigation into the role of three-dimensional coherent structures in the flow.
We therefore address this gap by providing an extensive spectral characterisation of the coherent flow dynamics. The specific objectives of this study are to (i) identify the role of coherent structures in the broadband turbulent dynamics of the flow, (ii) compare their driving mechanisms between attached and fully separated flow conditions and (iii) assess the role of the finite span and tunnel sidewalls on the dominant flow structures. In light of the discrepancies between CFD and experiments for this flow, the analysis is based entirely on experimental data.
1.4. Structure
Section 2 provides a brief introduction to the Boeing Gaussian Bump benchmark case and the available experimental database. § 3 outlines the main methodologies used in the study, including spectral proper orthogonal decomposition (SPOD), linear stability analysis (LSA) and resolvent analysis (RA). In § 4, coherent structures are identified using SPOD, highlighting low-frequency streaky structures and medium-frequency vortex shedding in both flow configurations. In § 5, LSA and RA are employed to model this coherent dynamics. In the separated case, the streaky structures are found to be linked to a three-dimensional zero-frequency global mode, whereas no evidence for a modal origin is found in the attached case. The possible driving mechanisms are then discussed and compared between both cases. Spanwise-standing-wave dynamics resulting from the finite span of the wind tunnel is investigated in § 6, and its implications regarding the domain size and boundary conditions of numerical simulations are discussed. Finally, § 7 summarises the main findings of the study.
2. Database
The Boeing Gaussian Bump was developed as a new benchmark test case for high Reynolds number flows undergoing SBS (Williams et al. Reference Williams, Samuell, Sarwas, Robbins and Ferrante2020; Gray et al. Reference Gray, Gluzman, Thomas, Corke, Lakebrink and Mejia2021). In the streamwise direction, the bump follows a Gaussian profile, where favourable and adverse pressure gradients are induced at the up- and downstream faces of the bump, respectively. In the spanwise direction, the bump is tapered according to an error function to minimise sidewall interactions. The bump geometry is defined as
\begin{equation} y_\varGamma (x, z) = h\frac {1+\mathrm{erf}\left (({1}/{2}-2z_0 - |z|)/z_0\right )}{2}\exp \left (-\left (\frac {x}{x_0}\right )^2\right )\!, \end{equation}
where
$(x, y, z)$
are the streamwise, vertical and spanwise coordinates, respectively, and
$y_\varGamma (x,\,z)$
is the bump surface. All variables in this paper are expressed in their respective non-dimensional form with the (spanwise) wind tunnel width
$L$
and free-stream velocity
$U_\infty$
serving as integral reference scales. Here,
$h=0.085L$
is the bump height, erf is the (Gauss) error function,
$z_0=0.06L$
and
$x_0=0.195L$
. Throughout this paper, unless otherwise noted, the Reynolds number is defined as
$\textit{Re}=U_\infty L/\nu$
, where
$\nu$
is molecular viscosity.
2.1. Experimental database
The bump geometry and measurement locations are shown in figure 1(b). Mean and instantaneous pressure measurements are available on the bump surface. The positions of the 6 pressure sensors where simultaneous time series were recorded are indicated in figure 1(b) by purple crosses. The signals were sampled at 100 kHz for 20 s. In the streamwise plane, PIV measurements provide the streamwise and vertical mean velocity components, with instantaneous velocity data available in the black-bordered regions labelled ‘FOV 1’ and ‘FOV 2’. These measurements are located along the spanwise centreline at
$z=0$
. In the spanwise plane, downstream of the bump at
$x=0.361$
, mean and instantaneous velocity data for all three components are available from stereo PIV (SPIV) measurements. The pressure data and PIV mean fields are included in the aforementioned archive, and the PIV/SPIV snapshot data were additionally made available by P. Gray for use in this study. The PIV was recorded at 200 Hz for 3 intervals of 5 s each for the streamwise windows and for one continuous interval of 25 s for the spanwise SPIV window. It is noted that the PIV/SPIV measurements are not synchronised with the pressure measurements. The reader is referred to the ‘NASA Turbulence Modeling Resource’ and the associated publications (Gray et al. Reference Gray, Gluzman, Thomas, Corke, Lakebrink and Mejia2021, Reference Gray, Gluzman, Thomas and Corke2022a
,
Reference Gray, Gluzman, Thomas, Corke, Lakebrink and Mejiab
, Reference Gray, Corke, Thomas, Gluzman and Straccia2023a
,
Reference Gray, Lakebrink, Thomas, Corke, Gluzman and Stracciab
; Gray Reference Gray2023) for a more detailed description of the experimental set-up.
We consider two flow conditions with Reynolds numbers
$\textit{Re}=U_\infty L/\nu =10^6$
and
$2\times 10^6$
based on the free-stream velocity
$U_\infty$
and wind tunnel width
$L$
. In the experimental configuration (Gray Reference Gray2023), the free-stream velocities are
$U_\infty =17\,$
and
$34\,\mathrm{ms}^{-1}$
, respectively, and the wind tunnel width is
$L=0.91\,\mathrm{m}$
. The Reynolds numbers based on the bump height are
$8.5\times 10^4$
and
$1.7\times 10^5$
, and the Mach numbers are 0.05 and 0.1. The variation of the mean surface-pressure coefficient for both cases is shown in figure 4. The flow at
$\textit{Re}=10^6$
displays only intermittent or very weak separation with no reverse mean flow whereas the flow at
$\textit{Re}=2\times 10^6$
is fully separated. This distinction is based on experimental observations (Gray et al. Reference Gray, Lakebrink, Thomas, Corke, Gluzman and Straccia2023b
) and consistent with results from DNS (Uzun & Malik 2021, Reference Uzun and Malik2022). In the following, we refer to these flow conditions as the attached and separated cases, respectively. It is worth noting, however, that these labels are based on the mean state of the flow, and that the instantaneous state at any specific time might be different.
The difference in the flow dynamics between the two different cases is evident from figure 2, which shows power spectral density (PSD) spectra of surface-pressure measurements in the downstream region of the bump. The x-axis shows the Strouhal number
${\textit{St}}_h=\textit{fh}/U_\infty$
, where
$h=0.085L$
is the bump height. In the proposed scaling by Mabey (Reference Mabey1972), the Strouhal number is based on the separation length. For the fully separated flow over the bump, the separation length was reported as
$L_{\textit{sep}}\approx 0.3L$
(Uzun & Malik Reference Uzun and Malik2022; Gray Reference Gray2023). At
$\textit{Re}=10^6$
, however, there is no flow separation and the length scale would be undefined. Therefore, a simple geometric length scale is used for both cases throughout this paper. For the separated case, the reader may convert to a separation-length-based scaling using
$L_{\textit{sep}}=3.53h$
. From figure 2, two primary regions of elevated PSD are identified: There is a medium-frequency regime at
$0.1 \lessapprox {\textit{St}}_h \lessapprox 1$
, indicated by the blue area. This regime is evident at both Reynolds numbers. For the separated case at
$\textit{Re}=2\times 10^6$
, there is also a prominent low-frequency regime at
${\textit{St}}_h \lessapprox 0.05$
, indicated by the purple area. These spectra thus indicate a shedding dynamics at both Reynolds numbers and a pronounced low-frequency breathing dynamics in the separated case. In what follows, we identify coherent structures which govern the flow dynamics in these regimes.
Power spectral density of surface-pressure measurements in the downstream region of the bump for the attached (a) and separated (b) cases. Tick marks at the top denote the frequency resolution.
2.2. Data-assimilated mean flows
To examine the dominant coherent structures within the broadband turbulent flow over the bump, we employ LSA and RA. These approaches require the construction of a linearised operator, which depends on the mean-flow state, along with its spatial gradients, defined across the entire domain of interest. However, the experimentally obtained mean velocity fields measured via PIV are limited to disjoint regions and do not provide full spatial coverage. The discrete nature of these measurements also complicates the accurate computation of gradients. Moreover, an eddy-viscosity field, typically used in linearised mean field methods to represent turbulence effects on coherent structures (Reynolds & Hussain Reference Reynolds and Hussain1972), cannot be directly extracted from the sparse experimental data.
Data-assimilated mean flow used in the approximation of the linear operators for the attached (a) and separated cases (b), with the eddy-viscosity field shown as background contours and mean streamlines overlaid in blue. Based on data reported in (Klopsch et al. Reference Klopsch, Fuchs, Rigas, Oberleithner, Von and Jakob2025).
Validation of the streamwise variation of mean surface pressure at
$z=0$
for the attached (a) and separated cases (b), with the grey-shaded region indicating the 95 % confidence interval for the experimental data. The pressure coefficient is offset such that
$c_p=0$
at
$x=-0.4$
. Based on data reported in (Klopsch et al. Reference Klopsch, Fuchs, Rigas, Oberleithner, Von and Jakob2025).
To address these issues, the mean flows have been data assimilated using physics-informed neural networks (PINNs) in our previous study (Klopsch et al. Reference Klopsch, Fuchs, Rigas, Oberleithner, Von and Jakob2025). The data assimilation approach combines the measured velocity mean fields from PIV at the spanwise centreline of the bump (Gray et al. Reference Gray, Gluzman, Thomas, Corke, Lakebrink and Mejia2021, Reference Gray, Gluzman, Thomas and Corke2022a
,
Reference Gray, Gluzman, Thomas, Corke, Lakebrink and Mejiab
, Reference Gray, Corke, Thomas, Gluzman and Straccia2023a
,
Reference Gray, Lakebrink, Thomas, Corke, Gluzman and Stracciab
; Gray Reference Gray2023) with physical constraints in the form of the (two-dimensional) RANS and continuity equations and a no-slip-wall boundary condition. The PINN provides two-dimensional (automatically) differentiable velocity fields at the bump centreline, which cover the relevant region. Additionally, a corresponding mean-field-consistent eddy viscosity and pressure field are inferred in the process. The assimilated velocity and eddy-viscosity fields are shown in figure 3. In order to validate the data assimilation approach against unseen data, the assimilated pressure at the bump surface is compared with experimental measurements in the previous study. The results of this validation step are shown in figure 4. Here, the offset is adjusted so that the pressure coefficient
$c_p=0$
at
$x=-0.4$
. This is necessary because the RANS equations used in the data-assimilation procedure only include the pressure gradient and provide no information on the absolute value. The grey-shaded region in figure 4 indicates the
$95\,\%$
confidence interval for the experimental data. For details regarding the uncertainty analysis, we refer to the Appendix of Gray (Reference Gray2023). The veracity of the eddy-viscosity field was assessed through a breakdown of the RANS equations, which showed that the turbulent forces are represented with satisfactory accuracy (see figure 8 in Klopsch et al. (Reference Klopsch, Fuchs, Rigas, Oberleithner, Von and Jakob2025)). Relatively small residuals in the force balance remain, which likely result from three-dimensionality of the mean flow, which is not accounted for in the two-dimensional equations used for the PINN, and from limitations of the Boussinesq hypothesis. For a detailed description of the data assimilation methodology and validation of the assimilated mean fields, we refer to our previous studies (Klopsch et al. Reference Klopsch, Fuchs, Rigas, Oberleithner, Von and Jakob2025; von Saldern et al. Reference von Saldern, Reumschüssel, Kaiser, Sieber and Oberleithner2022). The data-assimilated mean fields serve as foundation for approximating the linear operators in this study.
3. Methodology
This section outlines the main methods used in this study. We begin with an introduction to the SPOD methodology, which, throughout this study, serves as the primary data-driven approach for identifying coherent structures, followed by an introduction of the linearised Navier–Stokes equations, LSA and the RA framework.
3.1. Spectral proper orthogonal decomposition
To investigate the dominant dynamics in the broadband turbulent flow, SPOD is applied (Lumley Reference Lumley1970; Towne, Schmidt & Colonius Reference Towne, Schmidt and Colonius2018). The method identifies structures of spatial and temporal coherence in a temporally and spatially resolved signal based on estimates of the corresponding cross-spectral density (CSD) matrix. In this work, SPOD is applied to both surface-pressure measurements and time-resolved PIV snapshot data.
The SPOD algorithm involves several steps. First, the time series of a signal
$\boldsymbol{q}'$
involving
$M$
degrees of freedom is divided into
$N$
(overlapping) blocks. Each block is then transformed into the frequency domain using a Fourier transform and the resulting coefficients of all blocks are partitioned into a data matrix
\begin{equation} \kern2.5pt \hat { \kern-2.5pt \unicode{x1D64C}}_\omega = \left ( \hat {\boldsymbol{q}}^1_\omega \quad \hat {\boldsymbol{q}}^2_\omega \quad \ldots \quad \hat {\boldsymbol{q}}^N_\omega \right ) \end{equation}
for each frequency
$\omega$
. Following Welch’s method (Welch Reference Welch1967), the CSD matrix at a given
$\omega$
is estimated based on the Fourier modes of all
$N$
blocks,
$\mathrm{CSD}\approx{\alpha }/{N}\kern2.5pt \hat { \kern-2.5pt \unicode{x1D64C}}_\omega \kern2.5pt \hat { \kern-2.5pt \unicode{x1D64C}}_\omega ^{\mathsf{H}}$
. The superscript
$\mathsf{H}$
denotes the complex-conjugate transpose, and the factor
$\alpha = 1/\Delta f$
ensures that the eigenvalues represent modal PSD. Here,
$\Delta f$
is the discrete Fourier transform bin spacing of a windowed segment,
$\Delta f = {1}/({N_w\,\Delta t})$
, with time step
$\Delta t$
, and window length
$N_w$
. This choice is consistent with a forward normalisation of the Fourier transform. If a taper is applied,
$\alpha$
must be adjusted to compensate for the power loss. An eigenvalue decomposition of the CSD matrix yields the SPOD eigenvalues and corresponding eigenmodes
\begin{equation} \lambda _\omega ^k\boldsymbol{\phi }_\omega ^k = \frac {\alpha }{N}\kern2.5pt \hat { \kern-2.5pt \unicode{x1D64C}}_\omega \kern2.5pt \hat { \kern-2.5pt \unicode{x1D64C}}_\omega ^{\mathsf{H}}\boldsymbol{\phi }_\omega ^k. \end{equation}
This process is repeated for each frequency. For processing the PIV snapshot data, where
$M \gt N$
, we make use of the ‘method of snapshots’, where
\begin{equation} \boldsymbol{\phi }^k_\omega = \sqrt {\frac {\alpha }{N \lambda _\omega ^k}}\kern2.5pt \hat { \kern-2.5pt \unicode{x1D64C}}_\omega \boldsymbol{\psi }_\omega ^k\quad \text{with}\quad \lambda ^k_\omega \boldsymbol{\psi }^k_\omega = \frac {\alpha }{N}\kern2.5pt \hat { \kern-2.5pt \unicode{x1D64C}}_\omega ^{\mathsf{H}}\kern2.5pt \hat { \kern-2.5pt \unicode{x1D64C}}_\omega \boldsymbol{\psi }_\omega ^k, \end{equation}
which reduces the eigenvalue problem to be of size
$N\times N$
. For a comprehensive overview of the method, we refer the reader to the SPOD guide by Schmidt & Colonius (Reference Schmidt and Colonius2020).
The method yields
$N$
modes for each frequency bin that can be ranked according to their PSD, which is contained in the corresponding eigenvalue. If over certain frequency bands large parts of the overall PSD are associated with a small number of modes, one refers to a low-rank dynamics. The leading modes are then considered as the dominant coherent structures of the flow in the respective frequency range.
For processing the PIV snapshot data, we make use of the PySPOD package (Mengaldo & Maulik Reference Mengaldo and Maulik2021).
3.2. Linear stability analysis
Linear stability analysis is based on a linearised formulation of the Navier–Stokes equations around the temporal mean flow. The governing equations can be derived by substituting the Reynolds decomposition
$(\boldsymbol{\cdot })=\bar {(\boldsymbol{\cdot })}+(\boldsymbol{\cdot })'$
into the Navier–Stokes equations and subtracting the mean of the equations
\begin{equation} \frac {\partial \boldsymbol{u}'}{\partial t} = -(\boldsymbol{\bar {u}}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{u}' - (\boldsymbol{u'}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{\bar {u}} - \boldsymbol{\nabla }\!p' + \frac {1}{\textit{Re}} \boldsymbol{\nabla} ^2\boldsymbol{u}' - \boldsymbol{\nabla }\boldsymbol{\cdot } \unicode{x1D64D}^{\kern1pt \prime}. \end{equation}
The equations are given in non-dimensional form;
$\boldsymbol{u}$
denotes the velocity,
$t$
time,
$p$
pressure and
$\textit{Re}=U_\infty L/\nu$
the Reynolds number;
$\bar {(\boldsymbol{\cdot })}$
indicates the (temporal) mean and
$(\boldsymbol{\cdot })'$
is the fluctuating component. Incompressibility and constant density are assumed, allowing the non-dimensional density to be omitted in the equation. Due to the turbulent nature of the flow, the linearised equations include an unknown turbulent term
$\unicode{x1D64D}^{\kern1pt \prime} = \boldsymbol{u'}\boldsymbol{u'} - \overline {\boldsymbol{u'}\boldsymbol{u'}}$
that is also referred to as the fluctuating Reynolds stress tensor (Reynolds & Hussain Reference Reynolds and Hussain1972). This term represents the nonlinear interaction between different scales and is essential for their energy transfer (Kuhn et al. Reference Kuhn, Müller, Knechtel, Soria and Oberleithner2022; von Saldern et al. Reference von Saldern, Reumschüssel, Kaiser, Sieber and Oberleithner2024).
A harmonic ansatz for all fluctuating quantities of the form
\begin{equation} \boldsymbol{q}'(x,y,z,t) = \int _\beta \int _\omega \hat {\boldsymbol{q}}(x,y) \exp (\mathrm{i} \beta z - \mathrm{i}\omega t) + \mathit{c.c.} \end{equation}
is chosen, where
$\beta$
denotes the spanwise wavenumber and
$\omega$
the temporal frequency. Through
$\beta$
, we introduce a harmonic ansatz in the spanwise direction that reduces the three-dimensional problem to a two-dimensional problem for each
$\beta$
. This treatment of the spanwise coordinate allows us to capture some three-dimensional effects while retaining the simplicity of two-dimensional computations. Implicitly, this approach assumes both the geometry and mean flow to be constant along the spanwise coordinate. Obviously, this is not the case in the configuration considered in this study. Nevertheless, it is shown in this paper, through extensive validation of the linear models against experimental data, that this modelling approach is sufficient to capture the key dynamics of the flow. By inserting this ansatz, the equation is transformed into frequency space
\begin{equation} \mathrm{i}\omega \boldsymbol{\hat {u}} = (\boldsymbol{\bar {u}}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{\hat {u}} + (\boldsymbol{\hat {u}}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{\bar {u}} + \boldsymbol{\nabla }\!\hat {p} - \frac {1}{\textit{Re}} \boldsymbol{\nabla} ^2\boldsymbol{\hat {u}} + \boldsymbol{\nabla }\boldsymbol{\cdot } \hat { \unicode{x1D64D}}, \end{equation}
where the coherent component of the Reynolds stress tensor
$\hat { \unicode{x1D64D}}$
represents the fluctuating Reynolds stresses in frequency space. In analogy to the Boussinesq model for the RANS equations, an eddy-viscosity model is employed to represent the deviatoric component of the coherent Reynolds stress tensor
\begin{equation} -\mathrm{i}\omega \boldsymbol{\hat {u}} + (\boldsymbol{\bar {u}}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{\hat {u}} + (\boldsymbol{\hat {u}}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{\bar {u}} + \boldsymbol{\nabla }\hat {q} - \boldsymbol{\nabla }\boldsymbol{\cdot } \left(\frac {1}{\textit{Re}} + \nu _{{t}}\right)\big[\big(\boldsymbol{\nabla }+ \boldsymbol{\nabla} ^{\mathsf{T}}\big)\boldsymbol{\hat {u}}\big] = \boldsymbol{\hat {f}}. \end{equation}
The remaining spherical share is absorbed into the pressure term forming the modified harmonic pressure
$\hat {q} = \hat {p} + 1/3 \text{Tr}(\hat { \unicode{x1D64D}})$
, where
$\text{Tr}$
is the trace operator. For the eddy viscosity, we follow the common approach and use the value consistent with the mean field equations (Rukes, Paschereit & Oberleithner Reference Rukes, Paschereit and Oberleithner2016; Tammisola & Juniper Reference Tammisola and Juniper2016; von Saldern et al. Reference von Saldern, Reumschüssel, Kaiser, Sieber and Oberleithner2024) that is available from the data assimilation. The vector
$\boldsymbol{\hat {f}}$
, also referred to as the nonlinear forcing vector, can be interpreted as the remaining share of coherent Reynolds stresses that is not captured by the Boussinesq eddy-viscosity model. Equation (3.7) is complemented by the continuity condition for the fluctuating velocity field
$\boldsymbol{\nabla }\boldsymbol{\cdot } \boldsymbol{\hat {u}} = 0$
. In a compact form, the two governing equations can be written as
\begin{equation} -\mathrm{i}\omega \boldsymbol{\hat {q}} = \mathcal{L}\boldsymbol{\hat {q}}+\mathcal{B}\boldsymbol{\hat {f}}, \end{equation}
where
$\boldsymbol{\hat {q}}$
is the coherent state vector, including the velocity components and pressure,
$\mathcal{L}$
is the linearised operator incorporating (3.7) and the continuity condition and
$\mathcal{B}$
is a restriction operator, constraining the forcing to the momentum equations. Considering
$\boldsymbol{\hat {f}}=\boldsymbol{0}$
, the eigenvalues and eigenvectors of this system are the linear stability modes. The real part of the eigenvalue
$\mathrm{Re}(\omega _{\textit{eig}})$
represents the frequency and the imaginary part
$\mathrm{Im}(\omega _{\textit{eig}})$
the growth rate of the mode. The spatial shape of the mode is given by the corresponding eigenvector.
3.3. Resolvent analysis
To additionally investigate coherent structures emerging from non-modal mechanisms, RA is applied (McKeon & Sharma Reference McKeon and Sharma2010). Resolvent analysis focuses on the forced system dynamics, (3.8) with
$\boldsymbol{\hat {f}}\neq \boldsymbol{0}$
that can be rearranged into an input–output transfer function
\begin{equation} \boldsymbol{\hat {q}} = \underbrace {(-\mathcal{L}-\mathrm{i}\omega \boldsymbol{I})^{-1}}_{\mathcal{R}}\mathcal{B}\boldsymbol{\hat {f}}, \end{equation}
where
$\mathcal{R}$
is the resolvent operator. The operator maps a given forcing to the corresponding response in terms of velocity and pressure fluctuations. Note that a separate resolvent operator is obtained for each
$\beta$
and
$\omega$
. Since the true nonlinear forcings
$\boldsymbol{\hat {f}}$
are generally not known, we analyse the operator in terms of its optimal input–output behaviour. Specifically, we formulate
\begin{equation} \sigma ^2 = \max _{\boldsymbol{\hat {f}}}\frac {\boldsymbol{\hat {q}}^{\mathsf{H}} \unicode{x1D652}_r\boldsymbol{\hat {q}}}{\boldsymbol{\hat {f}}^{\mathsf{H}} \unicode{x1D652}_{\!f}\boldsymbol{\hat {f}}}, \end{equation}
which seeks forcing–response pairs that yield the maximum amplification
$\sigma ^2$
. The matrices
$ \unicode{x1D652}_r$
and
$ \unicode{x1D652}_{\!f}$
define the discrete turbulent kinetic energy (TKE) norms in which the output and input are measured, respectively. This optimisation problem is solved through a singular value decomposition of the resolvent operator
\begin{equation} \mathcal{R} = \unicode{x1D651} \boldsymbol{\varSigma } \unicode{x1D641}^{\,\mathsf{H}}, \end{equation}
where
$ \unicode{x1D641}$
and
$ \unicode{x1D651}$
contain the forcing and response modes, respectively. The corresponding amplification factors are found as the singular values
$\sigma$
contained in the diagonal matrix
$\boldsymbol{\varSigma }$
. The leading (largest) singular value and its associated forcing and response modes represent the most amplified linear mechanism in the flow. Because of the strong linear amplification and the broadband excitation in turbulent flows, resolvent response modes associated with large gains have been shown to reliably model dominant coherent structures in turbulent shear flows (Pickering et al. Reference Pickering, Rigas, Nogueira, Cavalieri, Schmidt and Colonius2020; Müller et al. Reference Müller, Von Saldern, Kaiser and Oberleithner2024; Sarras et al. Reference Sarras, Tayeh, Mons and Marquet2024). In the theoretical case of fully uncorrelated true nonlinear forcings, it can even be shown that the resolvent response modes correspond exactly to the SPOD modes (Towne et al. Reference Towne, Schmidt and Colonius2018; Lesshafft et al. Reference Lesshafft, Semeraro, Jaunet, Cavalieri and Jordan2019).
Equation (3.9) becomes singular when
$\omega =\omega _{\textit{eig}}$
, that is when the resolvent is computed at a frequency corresponding to an eigenvalue of the system. In this case, the optimisation problem in (3.10) yields the eigenvector as response, the null vector as forcing and infinite gain. In conventional RA, this manifests as sharp peaks in the gain spectrum near marginally stable (i.e. zero-growth-rate) eigenmodes. This problem can be addressed by employing discounted RA (Jovanovic Reference Jovanovic2004; Rolandi et al. Reference Rolandi, Ribeiro, Yeh and Taira2024), which is evaluated at a complex-valued
$\omega$
. Typically, an imaginary offset larger than the maximum growth rate of all eigenvalues is added to
$\omega$
to ensure the resolvent operator is not evaluated at or very close to eigenvalues of the system. This way, an interpretable RA spectrum is obtained. This imaginary offset on the frequency acts as a temporal discounting parameter, enabling the system dynamics to be evaluated over a finite time horizon and separating shorter time scales from the exponential growth of the modal instability (Jovanovic Reference Jovanovic2004; Rolandi et al. Reference Rolandi, Ribeiro, Yeh and Taira2024).
Sketch of the computational domain. (a) Full extent of the computational domain with the blue area indicating the PINN domain and black lines show contour levels of the sponge: 0.1 (solid), 1 (dashed), 10 (dash-dotted) and 50 (dotted). Axis breaks are used to highlight the most relevant part of the domain. (b) Zoomed view on the bump, where the red area denotes regions with PIV data and the blue area the response domain for RA. (c) Visualisation of the mesh near the bump surface.
3.4. Computational domain
A sketch of the computational domain is shown in figure 5(a). Here, the extent of the axes indicates the size of the computational domain used for LSA and RA. The extent of the data-assimilated mean flows is indicated by the blue area. Outside the blue area in the far field, nearest-neighbour extrapolation is applied to the mean field quantities. In order to suppress spurious free-stream modes, the RA energy norm is constrained to the blue-coloured area in figure 5(b). This is realised by using a corresponding weight matrix
$ \unicode{x1D652}_r$
in (3.10). We do not apply a spatial restriction on the forcing. To ensure numerical stability, sponging is applied, which absorbs and minimises reflections from computational boundaries outside the region of interest (Bodony Reference Bodony2006). The sponge level is zero within the rectangle defined by
$-0.4\leqslant x\leqslant 3$
,
$0\leqslant y \leqslant 0.4$
and grows exponentially with the distance to the closest point on the rectangle. Contours of the sponge level are shown in figure 5(a). The linear analyses are performed using FELiCS (Kaiser et al. Reference Kaiser, Demange, Müller, Knechtel and Oberleithner2023), an open-source finite element solver for linearised mean field methods. Additional information regarding the implementation of the RA can be found in the documentation available at https://felics.eu/. The discretisation is performed on a triangular mesh. A close-up of the mesh near the bump surface is shown in figure 5(c).
3.5. Mode alignment
In addition to a qualitative comparison of the modes through visualisation, we consider a quantitative measure, the alignment between the modes
\begin{equation} A(\boldsymbol{\hat {q}}_1, \boldsymbol{\hat {q}}_2) = \frac {\boldsymbol{\hat {q}}_1^{\mathsf{H}} \unicode{x1D652}_{\!A} \boldsymbol{\hat {q}}_2}{\sqrt {\left (\boldsymbol{\hat {q}}_1^{\mathsf{H}} \unicode{x1D652}_{\!A} \boldsymbol{\hat {q}}_1\right )\left (\boldsymbol{\hat {q}}_2^{\mathsf{H}} \unicode{x1D652}_{\!A} \boldsymbol{\hat {q}}_2\right )}}\!, \end{equation}
based on an inner product defined by
$ \unicode{x1D652}_{\!A}$
. This measure is typically employed to assess the similarity between modes (Gudmundsson & Colonius Reference Gudmundsson and Colonius2011; Cavalieri et al. Reference Cavalieri, Rodríguez, Jordan, Colonius and Gervais2013; Pickering et al. Reference Pickering, Rigas, Schmidt, Sipp and Colonius2021) and is sometimes named correlation coefficient or normalised inner product. The alignment is bounded between 0 and 1, with 1 meaning ‘parallel’ and 0 meaning ‘orthogonal’ modes.
4. Data-driven analysis of dominant flow structures
In order to tackle objective (i) of this study, to identify the role of coherent structures in the broadband turbulent dynamics of the flow, we apply SPOD to identify coherent structures in the experimental flow data. In this section, we only consider the SPOD spectra. The respective mode shapes are shown in § 5, where they are compared with the RA results.
We first apply SPOD to the instantaneous pressure measurements from the downstream region of the bump. The block length of the SPOD was chosen as
$2\times 10^4$
snapshots, which yields a frequency resolution of
$\Delta f = 5\,\mathrm{Hz}$
(
$\Delta {\textit{St}}_h=0.023$
for the attached and
$\Delta {\textit{St}}_h=0.011$
for the separated case). The minimum resolvable convection velocity for streamwise-travelling waves is set by the spatial Nyquist criterion,
$u_{\textit{c,crit}}=2 {\textit{St}}\,\Delta x$
. For the dominant shedding frequencies of
$St=St_h/h=3.29$
and
$St=2.12$
in the attached and separated cases, respectively, and the normalised sensor spacing of
$\Delta x= 0.063$
, this corresponds to minimum resolvable convection velocities of
$u_{\textit{c,crit}} = 0.41$
and
$u_{\textit{c,crit}} = 0.27$
. In TSB flows, the expected convection velocities of coherent structures are typically in the range
$0.4 \lt u_{\textit{c}} \lt 0.6$
(Kiya & Sasaki Reference Kiya and Sasaki1983; Cherry et al. Reference Cherry, Hillier and Latour1984; Mohammed-Taifour & Weiss Reference Mohammed-Taifour and Weiss2016). Therefore, the expected coherent structures up to the dominant shedding frequency can be resolved using the pressure sensor data, especially in the separated case. The SPOD spectrum obtained from this dataset is shown in figure 6. In the upper panels, the sum over all eigenvalues at the respective frequency,
$\sum \lambda _i$
, is shown as a red line. Note that
$\sum \lambda _i$
is equal to the sum of the PSDs of the individual sensor signals. It can be seen that the leading eigenvalue approaches the red line in the highlighted frequency regimes, identified in the PSDs of single sensors (see figure 2). This is highlighted through the lower panels in figure 6, where the eigenvalues are divided by
$\sum \lambda _i$
. This representation shows the PSD share associated with the respective mode. Evidently, the prominent regimes identified so far also correspond to frequency ranges of the low-rank dynamics. In the medium frequency regime, approximately 30 % to 75 % of the signal’s PSD is represented by the leading SPOD mode in both cases. In the low-frequency regime, the leading mode accounts for approximately 40 % of the PSD in the attached – and 40 % to 65 % in the separated case. This observation shows that, while the low-frequency dynamics does not appear particularly prominent in the absence of flow separation (
$\textit{Re}=10^6$
), the present dynamics nevertheless appears to be low rank, indicating the presence of coherent structures.
Spectrum of the surface-pressure SPOD for the attached (a) and separated (b) cases, with the leading mode shown in black and subsequent modes in progressively lighter shades of grey. In the upper panels, the sum of all eigenvalues at each frequency is shown as a red line. In the lower panels, the eigenvalues are normalised by this value to indicate the PSD share of each mode. The blue-shaded region marks the medium-frequency regime, the purple-shaded region marks the low-frequency regime and tick marks at the top denote the frequency resolution.
The SPOD spectra of the PIV data from the different measurement windows for the attached (a) and separated (b) cases, with the leading mode shown in black and subsequent modes in progressively lighter shades of grey. Eigenvalues are normalised by
$\sum \lambda _i$
at each frequency to indicate the PSD share of each mode. The blue-shaded region marks the medium-frequency regime, the purple-shaded region marks the low-frequency regime and tick marks at the top denote the frequency resolution.
The SPOD is also performed on the PIV velocity time series in two streamwise and one spanwise measurement windows (PIV domains are shown in figure 1(b)). Here, the block length was chosen as
$200$
snapshots which yields a frequency resolution of
$\Delta f=1\,\mathrm{Hz}$
(
$\Delta {\textit{St}}_h=0.005$
for the attached and
$\Delta {\textit{St}}_h=0.002$
for the separated case). The Nyquist frequency is
$100\,\mathrm{Hz}$
, that is
${\textit{St}}_h=0.46$
for the attached and
${\textit{St}}_h=0.23$
for the separated case, and as such, the medium-frequency regime cannot be fully resolved for the separated case. Instead, we focus on the low-frequency regime with this analysis. Figure 7 shows the resulting normalised eigenvalue spectra. For the separated case, these reveal a very significant low-rank dynamics in the low-frequency regime. In all measurement windows, the leading mode approaches
$\approx 75\,\%$
of the PSD at the lowest resolvable frequency. Notably, the sub-leading mode in the spanwise window also shows clear separation from the subsequent modes. In the attached case, however, no low-rank dynamics is observed in the spanwise window. In the streamwise windows, the separation between the eigenvalues is also less pronounced than in the separated case. Nevertheless, the leading mode accounts for more than
$25\,\%$
of the PSD in the low-frequency regime, showing the presence of a low-rank dynamics. This observation is consistent with the previous observation that pressure fluctuations at low frequencies, in the absence of flow separation, are still low rank. The second dominant region with high gain separation in the medium-frequency regime (blue background), identified in the surface-pressure SPOD spectra in figure 6, is absent in the PIV SPOD spectra (figure 7). We attribute this to the aforementioned Nyquist rate limitation of the PIV measurements.
Regarding objective (i), we draw the conclusion that the medium-frequency dynamics in both cases is characterised by a low-rank dynamics driving coherent structures. The same holds for the low-frequency regime, especially in the separated case, and to a lesser degree also in the attached case.
The LSA eigenvalue spectra over spanwise wavenumber for the attached (a) and separated cases (b), coloured by
${\textit{St}}_h$
. The neutral-stability line
$\mathrm{Im}(\omega )=0$
is shown as a black line. The most unstable (i.e. maximum growth rate) mode is marked with a black cross.
5. Physics-based model for coherent structures
In order to address objective (ii) of this study, to compare the driving mechanisms of the coherent structures between the cases, we analyse the coherent dynamics of the flow in a physics-based way by employing linear mean field methods, specifically LSA and RA. We first investigate the global stability of the mean flow using LSA and subsequently employ RA to analyse non-modal mechanisms.
5.1. Linear stability analysis
Linear stability analysis is performed for each spanwise wavenumber
$n={\beta }/({2\pi })$
separately in the range of
$0\leqslant n \lt 16$
. Figure 8 shows the resulting eigenvalue spectra. Eigenvalues with
$\mathrm{Im}(\omega )\lt 0$
are stable, whereas those with
$\mathrm{Im}(\omega )\gt 0$
are unstable (recall the harmonic ansatz in (3.5)). Evidently, the attached case does not exhibit any global instabilities: all eigenvalues have negative growth rates and are located at approximately the same distance from the neutral-stability line. We note that the branch that appears to be growing for high
$n$
corresponds to spurious modes. In contrast, the separated case exhibits a distinct eigenvalue branch with numerically zero frequency (
${\textit{St}}_h\lt 10^{-10}$
) and elevated growth rates over
$2.5\leqslant n \lt 13$
. The positive growth rate is unexpected for a stationary mean-flow analysis and likely arises from model limitations, namely the spanwise-uniform mean-flow assumption and the eddy-viscosity modelling of the coherent Reynolds stresses. Crucially, the isolation of this branch from the rest of the spectrum indicates that, for this spanwise wavenumber range, the low-frequency dynamics is governed by a global mode.
The corresponding mode shape for the largest growth rate (marked by the black cross) is shown in figure 9. Note that the mode is real valued due to the frequency being zero. In the following sections, we show that this globally unstable zero-frequency mode provides a good explanation for the low-frequency dynamics observed in the separated case.
The LSA eigenmode with maximum growth rate for the separated case. Here,
$n=6.5$
(marked in figure 8 with a black cross).
5.2. Resolvent analysis
The previous analysis with LSA did not reveal any distinct eigenmodes in the attached case, although a low-rank dynamics has been identified in the data-driven analysis. We therefore employ RA in order to model coherent structures that originate from non-modal mechanisms.
The RA is performed for both cases. In the separated case, however, the modal instability identified in the previous section would lead to singularities in the RA gain spectrum, which would obfuscate, for example, the preferential spanwise wavenumber of the low-frequency dynamics. One approach to obtain an interpretable gain spectrum would be to apply a scaling factor to the eddy viscosity in order to stabilise the system, similarly as in Fuchs et al. (Reference Fuchs, Steinfurth, von, Jakob, Weiss and Oberleithner2026). However, we find that, for this case, the eddy-viscosity scaling deteriorates the linear models’ performance in explaining the data. We therefore take a different approach and apply discounted RA by adding an imaginary offset (discount factor) of
$\gamma = \mathrm{Im}(\omega ) = 0.5$
to the RA frequencies in the separated case. The chosen value for
$\gamma$
is slightly higher than the maximum growth rate of the unstable eigenvalues (see figure 8), effectively stabilising the system. It is important to note that, for each value of
$\gamma$
, a different RA spectrum is obtained: increasing the discount factor lowers the overall RA gains and flattens the spectrum (Rolandi et al. Reference Rolandi, Ribeiro, Yeh and Taira2024, Reference Rolandi, Smith, Amitay, Theofilis and Taira2025). In the present case, the corresponding RA mode shapes (in the region of interest) remain essentially unchanged for small discount factors, and only very large values lead to noticeable modifications. For the chosen
$\gamma =0.5$
, we demonstrate the negligible effect on the mode shapes in Appendix A. The discount factor is often interpreted as an effective time horizon
$t_{\textit{eff}}={1}/{\gamma }$
over which the structures grow (Jovanovic Reference Jovanovic2004; Rolandi et al. Reference Rolandi, Ribeiro, Yeh and Taira2024, Reference Rolandi, Smith, Amitay, Theofilis and Taira2025). For
$\gamma =0.5$
in the present case, this effective time horizon corresponds to
$t_{\textit{eff}}=2$
convective time scales based on
$L$
. Alternatively, the discounted RA spectrum can be interpreted as a slice through the pseudospectrum of the resolvent operator along a line in the complex-frequency plane specified by
$\gamma$
. We emphasise that discounting is applied solely to yield a finite, well-behaved spectrum, and the precise value chosen does not affect the qualitative interpretation of the results. In the attached case, all eigenvalues have a negative growth rate, so no discount factor is required. For a discussion of the relationship between the resolvent and the eigenvalues of the linear operator via the pseudospectrum, as well as the effect of discounting RA, we refer to Appendix A.
Resolvent gain as a function of Strouhal number and spanwise wavenumber for the attached (a) and separated cases (b), shown as a heatmap with a logarithmically scaled colour bar. Black markers indicate Strouhal number–spanwise wavenumber pairs discussed in the following sections.
A grid search over
$n$
and
$St_h$
is performed and contours of the resulting resolvent gain are plotted in figure 10. Note the logarithmic colour scaling, which spans multiple orders of magnitude and is different for the two cases. At high spanwise wavenumbers (
$n \gt 14$
), the resolvent modes display increasing levels of numerical noise and the gain separation between the leading and sub-leading resolvent modes drops off (not shown). Furthermore, a different mode family is picked up, which is located on the upstream face of the bump, where no time-resolved data for validation are available. These modes are excluded from the analysis, as indicated by the hatched area in figure 10. The black markers in the figure mark parameter combinations of
${\textit{St}}_h$
and
$n$
that will be further discussed in the following sections.
In both cases, a medium-frequency regime is identified. The highest gain in this regime is observed for non-zero
$n$
, but the gain remains elevated as
$n\to 0$
. The dominant frequency of this regime is
${\textit{St}}_h\approx 0.28$
for the attached case and slightly lower, at
${\textit{St}}_h\approx 0.18$
for the separated case (see square markers in figure 10), which is in line with the observations in the surface-pressure spectra presented in figure 2 (or the surface-pressure SPOD spectra in figure 6). The most notable difference between the two cases is seen at low frequencies,
${\textit{St}}_h\ll 0.1$
, down to
${\textit{St}}_h=0$
. Here, we see the influence of the distinct eigenmode in the separated case, which manifests as a broad regime of very high RA gain around
$n\approx 6.5$
. In the attached case, the corresponding eigenvalues at
${\textit{St}}_h=0$
are much more strongly damped, and the resolvent is therefore evaluated farther from any eigenvalue, so no comparable feature appears in the RA gain spectrum. However, there is still significant RA gain for
${\textit{St}}_h\lt 0.1$
and down to
${\textit{St}}_h =0$
in this case, especially for elevated
$n$
.
In this study, the resolvent gain is based on the TKE norm of the response mode and is thus indicative of the magnitude of the velocity fluctuations associated with the mode. In many experimental studies, however, the flow dynamics is investigated based on time-resolved wall-pressure measurements. For a discussion regarding the coherent surface-pressure fluctuations associated with the RA modes and their comparison with measurements from the experiment, we refer to Appendix B.
5.3. Validation of the linear analyses
In order to validate the linear analyses, we compute the alignment between the leading RA mode and the leading SPOD mode at the respective frequency. The alignment is shown in figure 11 as a heatmap over Strouhal number and spanwise wavenumber. It should be noted that the SPOD computed from streamwise PIV measurements is not associated with a specific spanwise wavenumber
$n$
, but instead reflects contributions from a range of
$n$
. Consequently, strong alignment is not expected for all individual
$n$
values. We consider the linear model valid if, at each frequency in the relevant regimes, at least one resolvent mode (over all
$n$
) has high alignment with the SPOD mode at this frequency. The maximum alignment at each frequency is expected at the dominant spanwise wavenumber of the dynamics in the flow. Because the PIV in FOV 1 and FOV 2 was not measured simultaneously, SPOD is performed for both areas separately. We therefore obtain two alignment values (one for FOV 1 and one for FOV 2) for each Strouhal number and spanwise wavenumber.
Alignment between RA modes and streamwise SPOD modes at respective frequency as a function of Strouhal number and spanwise wavenumber for the attached (a) and separated (b) cases. Upper and lower panels show the alignment with the SPOD from FOV 1 and FOV 2, respectively. Black markers indicate Strouhal number–spanwise wavenumber pairs discussed in the following sections.
In the attached case, high values of alignment are observed for
${\textit{St}}_h \lessapprox 0.1$
over a wide range of
$n$
. As the frequency exceeds
${\textit{St}}_h= 0.1$
, high alignment becomes limited to higher
$n$
. In the separated case, very high alignment is observed for
${\textit{St}}_h \lt 0.01$
with a preferential
$n$
around
$5$
. In the medium-frequency regime, around
${\textit{St}}_h=0.15$
, there is a rather confined regime of high alignment for
$n \approx 3$
. Generally, there is a trend that the alignment is high where the RA gain is high, confirming that the dominant coherent structures in the flow originate from a linear dynamics. However, for the low-frequency modes, there is a broad region of high alignment around the dominant spanwise wavenumbers, which reflects the fact that the corresponding mode shapes in the plane at
$z=0$
are relatively insensitive to
$n$
. For the medium-frequency modes, the alignment is confined to a narrower band of wavenumbers, especially for the separated case. We attribute this mostly to the higher sensitivity of the planar mode shapes to
$n$
in this regime, combined with the observation that the underlying mechanism appears to preferentially amplify these wavenumbers (see figure 10(b)).
For the separated case, we further compute the alignment between the RA modes and the
${\textit{St}}_h=0$
LSA mode with maximum growth rate at the respective spanwise wavenumber. This alignment is shown as a heatmap in figure 12. Evidently, as
${\textit{St}}_h \to 0$
, the RA and LSA modes converge in the range of spanwise wavenumbers where the eigenmode has a distinct growth rate. Figure 12 thus demonstrates that, in the separated case, RA reflects the modal instability at
${\textit{St}}_h=0$
throughout the low-frequency regime, with little contribution from non-modal effects.
For the attached case, comparing all LSA modes with
${\textit{St}}_h\approx 0$
against the RA mode at
${\textit{St}}_h=0$
and the respective
$n$
confirms that there is no significant alignment. The reader is referred to Appendix C for details.
Alignment between the zero-frequency LSA eigenmode and RA modes. Both modes share the same
$n$
but the RA mode frequency is indicated on the x-axis whereas the LSA mode is always at
${\textit{St}}_h=0$
. Contour lines show
$A=0.99$
(black) and
$A=0.999$
(white).
For representative frequencies and spanwise wavenumbers of the low- and medium-frequency regimes, the RA modes are shown in figure 13 and compared with the respective SPOD modes of the streamwise PIV windows in the downstream region of the bump. In this section, we only focus on SPOD modes of the streamwise PIV, which include the streamwise and vertical velocity components (
$u, v$
). A comparison of the RA mode with the SPOD modes from the spanwise SPIV data requires consideration of the three-dimensional structure of the modes, which will be addressed in § 6. The selected frequencies and spanwise wavenumbers are indicated in the RA gain spectrum (figure 10) and the RA–SPOD alignment heatmaps (figure 11) through diamond markers. Note that the SPOD was performed for both windows separately, as they were not measured simultaneously. For presentation purposes, the SPOD modes from both windows are still shown on the same axes, and visually framed with coloured lines for clarity. The low-frequency modes (
${\textit{St}}_h=0.01$
and
${\textit{St}}_h =0.002$
, for the attached and separated cases, respectively), shown in figures 13(a) and 13(b), resemble large-scale, streamwise-elongated structures. We use the term ‘streaky structures’ in the remainder of this manuscript in order to avoid confusion with the term ‘streaks’, which is often used specifically to refer to structures associated with the lift-up mechanism. Interestingly, we find qualitatively similar low-frequency modes in the attached and separated cases. In the attached case, the structure follows the attached shear layer along the wall, whereas the structure in the separated case is located within the free shear layer further away from the wall. This is notable because the low-frequency dynamics is often attributed to the presence of a TSB, which is not the case in the mean flow of the attached case. The medium-frequency modes (
${\textit{St}}_h=0.27$
and
${\textit{St}}_h =0.16$
, for the attached and separated cases, respectively), shown in figures 13(c) and 13(d), are likewise similar between the attached and separated cases. These modes resemble vortices with relatively small streamwise extent. In this case, the similarity is expected because this medium-frequency dynamics is typically associated with vortex shedding in the shear layer, which is present across both cases.
Mode shape comparison between RA and SPOD at selected frequencies and spanwise wavenumbers for attached (a,c) and separated (b,d) cases. The SPOD modes from both fields of view (FOVs) (independent measurements) are shown on the same axes, and visually framed with coloured lines for clarity. Three-dimensional isosurfaces of the RA forcing (cyan, magenta) and response (blue, red) modes are included for reference, with the black plane indicating the region used for the SPOD–RA comparison. Frequencies and spanwise wavenumbers are indicated in figures 10, 11 and 19 through diamond markers.
Based on these qualitative and quantitative comparisons between the LSA, RA and SPOD modes, we conclude that the linear modelling results qualitatively capture the relevant trends and mode shapes of the flows’ coherent dynamics with satisfactory accuracy. In Appendix B, we further show that RA accurately captures the coherent surface-pressure fluctuations on the downstream side of the bump.
5.4. Physical mechanisms
As the model is validated, we now turn back to objective (ii) and compare the physical mechanisms leading to the formation of coherent structures in both cases. Here, we focus on the low-frequency dynamics.
For the separated case, the low-frequency dynamics is associated with a distinct ‘steady’ (i.e.
${\textit{St}}_h=0$
), three-dimensional (i.e.
$n\neq 0$
) eigenmode. We therefore speak of a modal dynamics. In the study of laminar separation bubbles, such a mode was attributed to a centrifugal instability (Gallaire, Marquillie & Ehrenstein Reference Gallaire, Marquillie and Ehrenstein2007; Rodríguez et al. Reference Rodríguez, Gennaro and Juniper2013; Savarino et al. Reference Savarino, Poulain, Sipp and Rigas2024). In TSB flows, this type of mode was also previously identified (Cura et al. Reference Cura, Hanifi, Cavalieri and Weiss2024; Sarras et al. Reference Sarras, Tayeh, Mons and Marquet2024; Fuchs et al. Reference Fuchs, Steinfurth, von, Jakob, Weiss and Oberleithner2026) and it was suggested that a similar centrifugal mechanism is responsible. In the attached case, on the other hand, we did not identify any distinct eigenmode corresponding to a physically relevant feature. Yet, RA correctly captures the dynamics, with
$A\gt 0.75$
throughout the low-frequency regime (see figure 11). There are two possible interpretations for the dynamics in the attached case.
-
(i) The dynamics is modal but associated with a highly stable eigenvalue, which the LSA fails to isolate. In this case, the driving mechanism might be a similar centrifugal instability as in the separated case, acting on the geometry-induced curvature of the mean-flow streamlines.
-
(ii) The dynamics is non-modal; that is, it is not associated with a specific eigenmode but rather emerges from the non-orthogonality of several eigenmodes. A likely explanation for a non-modal dynamics would be the lift-up mechanism that generates streamwise-elongated low-frequency streaks through the transport of momentum across a shear layer, which is facilitated by the action of streamwise vorticity (Landahl Reference Landahl1980; Brandt Reference Brandt2014).
Note that the coherent structures generated through both of these mechanisms are very similar and are subsequently both described as low-frequency streaky structures, as introduced above. From the aforementioned results, we cannot draw a definite conclusion on whether a weakened version of the modal instability gives rise to the streaky structures in the attached case, or whether these originate through a qualitatively different mechanism like the lift-up effect. However, the absence of a distinct branch in the eigenvalue spectrum favours the latter explanation.
Here, we compare the cases at
$\textit{Re}=10^6$
(attached) and
$\textit{Re} =2\times 10^6$
(separated). We find that a further increase of the Reynolds number to
$4\times 10^6$
does not qualitatively change the dynamics of the separated flow. This analysis is available in Appendix D.
6. Finite-span effects on the low-frequency dynamics
Phase angle of the
$\hat {u}$
-component of the leading (a,b) and sub-leading (c,d) low-frequency SPOD modes in the spanwise plane for the attached (a,c) and separated cases (b,d). Transparency is set according to the magnitude of the mode with full transparency at zero magnitude and no transparency at half the maximum value.
In order to address objective (iii) of this study, to assess the role of the finite span and tunnel sidewalls on the dominant flow structures, we analyse the three-dimensional structure of the low-frequency dynamics. To this end, we consider SPOD of SPIV measurements in the spanwise plane as shown in figure 1.
The phase angle of the
$\hat {u}$
-component of the leading and sub-leading SPOD modes at
$f=1\,\mathrm{Hz}$
(the lowest non-zero frequency bin) are shown in figure 14. The phase fields are masked by the mode’s magnitude, with the transparency set inversely proportional to the absolute value of the mode. This depiction reveals the wave structure of the mode. In the case of travelling waves, the phase evolves continuously in space. Standing waves, on the other hand, are characterised by nodes, where the oscillation amplitude remains zero, and antinodes, where it reaches a maximum. Between nodes, the phase remains constant, whereas discontinuities occur at the nodes. In the attached case (a, c), the SPOD modes are very noisy and the mode cannot be as clearly extracted as from the side view data. In the separated case (b, d), on the other hand, the modes show clear evidence of a standing-wave structure. Here, the leading mode (b) corresponds to a standing-wave pattern with a node on the centreline, whereas the sub-leading mode (d) has an antinode at this position.
The observation of spanwise-standing waves constitutes a challenge for modelling these structures with RA because, by periodically expanding the RA mode in spanwise direction according to the ansatz in (3.5), a spanwise-travelling wave is obtained. In order to model the experimentally observed dynamics, we therefore apply a spanwise-standing-wave assumption (Fuchs et al. Reference Fuchs, Steinfurth, von, Jakob, Weiss and Oberleithner2026), where the RA mode is expanded as
\begin{equation} \begin{aligned} \boldsymbol{\hat {u}}^{\mathrm{SW}}_{\beta ,\omega }(x,y,z,t) &= \boldsymbol{\hat {u}}_{\beta ,\omega }(x,y)\exp (\mathrm{i}\beta [z-0.5L_{\textit{eff}}^{\text{SW}}] - \mathrm{i}\omega t) \\ &+ \boldsymbol{\hat {u}}_{-\beta ,\omega }(x,y)\exp (-\mathrm{i}\beta [z-0.5L_{\textit{eff}}^{\text{SW}}] - \mathrm{i}\omega t). \end{aligned} \end{equation}
Here,
$\boldsymbol{\hat {u}}_{\beta , \omega }$
denotes the velocity field associated with the leading resolvent response mode, and
$\boldsymbol{\hat {u}}_{-\beta ,\omega }$
is its spanwise-symmetric counterpart travelling in the opposite direction. The two modes are identical in all components except the spanwise velocity, which has the opposite sign due to the symmetry condition. An effective spanwise length scale for the standing-wave dynamics
$L_{\textit{eff}}^{\text{SW}}$
is introduced, and the term
$-0.5L_{\textit{eff}}^{\text{SW}}$
shifts the spanwise coordinate, which is convenient for imposing boundary conditions. We impose slip-wall boundary conditions (
$\hat {w}=0$
) at
$z=\pm 0.5L_{\textit{eff}}^{\text{SW}}$
(Fuchs et al. Reference Fuchs, Steinfurth, von, Jakob, Weiss and Oberleithner2026), which yields permissible wavenumbers
$\beta = {2\pi n^{\text{SW}}}/{L_{\textit{eff}}^{\text{SW}}}$
with
\begin{equation} n^{\text{SW}} = \{0.5, 1, 1.5, 2, 2.5, \ldots \}. \end{equation}
This model assumes that the hydrodynamic waves are reflected at the slip-wall boundaries without phase lag and resonate in the channel. We note that this model cannot satisfy no-slip boundary conditions, due to the phase relationships between the different components (Fuchs et al. Reference Fuchs, Steinfurth, von, Jakob, Weiss and Oberleithner2026). Fuchs et al. (Reference Fuchs, Steinfurth, von, Jakob, Weiss and Oberleithner2026), the wind tunnel width is the natural choice for the effective spanwise length scale, as their backward-facing ramp geometry is spanwise homogeneous and bounded by the wind tunnel sidewalls. In the present case, however, the spanwise profile of the bump and the highly three-dimensional flow enable multiple plausible choices for
$L_{\textit{eff}}^{\text{SW}}$
. For the results presented here, we chose
$L_{{eff}}^{\text{SW}}=L$
as in Fuchs et al. (Reference Fuchs, Steinfurth, von, Jakob, Weiss and Oberleithner2026), as the wind tunnel width represents the least ambiguous option. This leads to the slip-wall boundary conditions being placed at the sidewalls. However, we do not intend to imply with this choice that the wind tunnel sidewalls are necessarily the dynamically relevant features governing the formation of the standing-wave pattern. Alternative candidates include the bump shoulders or specific features of the mean flow, as discussed below. We note that for
$L_{\textit{eff}}^{\text{SW}}=L$
the wavenumber convention used in the previous sections
$n={\beta }/({2\pi })$
is recovered.
Alignment between the RA model and the low-frequency SPOD mode in the spanwise plane for the attached (a) and separated (b) cases. The markers show the alignment between standing-wave (SW) model and SPOD mode, and the dashed lines show the alignment of the travelling wave (TW) mode for comparison. Black markers and curves show the alignment with the leading, and red markers and curves with the sub-leading SPOD mode.
The alignment of both the travelling and standing-wave models with the leading and sub-leading SPOD modes in the low-frequency regime is shown in figure 15 for different spanwise wavenumbers. Half-integer wavenumbers
$n=\{0.5, 1.5, 2.5, \ldots \}$
correspond to standing waves with a node for
$\hat {u}$
on the centreline, whereas integer wavenumbers
$n=\{1, 2, 3,\ldots \}$
have an antinode at this position. At
$n=0$
, both the travelling wave and the standing wave reduce to a spanwise-constant mode.
In the attached case (see figure 15(a)), the results are inconclusive and both the travelling wave modes and the standing-wave model display low alignments with the spanwise SPOD modes. For modelling the leading SPOD mode, the standing-wave model with wavenumbers
$n=\{2, 3\}$
provides some improvement over the travelling wave modes. Yet, the respective alignments are still low, which is an expected result given the unclear structure of the SPOD modes (see figures 14(a) and 14(c)) and the low separation between the leading and sub-leading modes in the respective SPOD spectrum (see figure 7). In the separated case, however, the standing-wave model achieves far higher alignment with the SPOD mode compared with the travelling-wave modes, provided that the wavenumber corresponds to the same node position (i.e. either integer or half-integer
$n$
). Both the standing-wave pattern with a node at the centreline, observed in the leading SPOD mode, as well as the pattern with an antinode at the centreline, observed in the sub-leading SPOD mode, are well captured through the RA standing-wave model.
Comparison of the standing-wave RA model (contour lines) and SPOD mode (image behind contours) in the spanwise view of the bump. The leading SPOD mode is shown in (a) and the sub-leading SPOD mode is shown in (b). The colour scale for the SPOD modes is clipped at
$0.25$
times the respective maximum absolute value of
$\hat {u}$
to highlight the structures of the other components, which are much smaller in comparison. Grey lines visualise how the standing-wave systems fulfils slip-wall boundary conditions at the sidewalls:
$\hat {u}$
and
$\hat {v}$
have an antinode at sidewalls and
$\hat {w}$
has a node at the sidewalls. Here,
${\textit{St}}=0.002$
,
$\textit{Re}=2\times 10^6$
(separated). The reader is referred to the animated version of the figure.
The respective mode shapes, at the wavenumbers corresponding to the highest alignments, are displayed and compared with the respective SPOD modes in figure 16. The reader is referred to the animated version of this figure in the supplementary materials. This visualisation demonstrates the capability of the model to capture the main characteristics of the mode shape. This is particularly the case for the
$\hat {u}$
component that contains most of the modes’ energy. The dashed lines show how a standing-wave structure with the selected spanwise wavenumber fulfils the slip-wall condition at the sidewalls. However, there are some differences between the RA and SPOD mode shapes: most notably, the
$\hat {u}$
component of the SPOD mode displays a roughly triangular shape with vertices
$(z,y)$
around
$(\pm 0.125, 0.025)$
and
$(0, 0.1)$
, which is antisymmetric with respect to the node at the centreline but not symmetric with respect to the antinodes. The RA modes, on the other hand, are both antisymmetric with respect to the nodes and symmetric with respect to the antinodes, which is a consequence of the spanwise periodic ansatz for the RA modes (see (3.5)). Yet, alignments of
$A=0.87$
and
$0.80$
, for the leading and sub-leading modes, respectively, are obtained with this model, demonstrating its applicability to the present flow.
The relevant spanwise length scale for the standing-wave dynamics remains unclear, as the available data are insufficient to directly determine this value. If a different effective spanwise length scale is considered, the modes are simply represented by different wavenumbers based on that scale. Nevertheless, the modal structure provides a speculative argument suggesting that, if the standing waves result from a single effective spanwise length scale, this value would have to be relatively large. From the mode shapes in figures 14(b) and 14(d), we know that the leading mode requires a half-integer wavenumber, whereas the sub-leading mode requires an integer wavenumber in order to satisfy the node and antinode at
$z=0$
, respectively, under the assumed standing-wave ansatz. We can also deduce that the ratio of spanwise wavelengths (leading to sub-leading) is roughly
$1.5$
. If the effective spanwise length scale was close to the spanwise extend of the separation bubble, approximately
$0.18L$
, the leading mode would best be represented by
$n=0.5$
and the sub-leading mode by
$n=1$
(the wavelength of the sub-leading mode approximately coincides with the bubble width). This would give a wavelength ratio of
$1/0.5=2$
, which is in conflict with the data. Hence, the modes more likely correspond to higher wavenumbers on a larger length scale. For example,
$n=1.5$
and
$n=2$
give a wavelength ratio of
$1.33$
, which is in much better agreement with the data. This observation suggests that the effective spanwise length scale is at least two times the spanwise wavelength of the sub-leading mode, that is
$L_{\textit{eff}}^{\text{SW}}\geqslant 0.36$
.
The prominence of spanwise-standing-wave patterns in the low-frequency regime, similar to the findings by Fuchs et al. (Reference Fuchs, Steinfurth, von, Jakob, Weiss and Oberleithner2026), has broader implications for numerical simulations of separated flows. Many computational studies, for reasons of efficiency, assume periodic spanwise boundary conditions and simulate only a thin spanwise section of the flow (e.g. Balin & Jansen (Reference Balin and Jansen2021), Uzun & Malik (Reference Uzun and Malik2025)). These modelling choices can inadvertently filter out low-wavenumber modes, particularly those associated with large-scale breathing behaviour. Moreover, periodic boundary conditions restrict admissible spanwise wavenumbers to integer multiples of the fundamental wavenumber, excluding half-integer standing-wave patterns, which may arise from reflection of hydrodynamic waves from the channel sidewalls in configurations with finite span, and in the present case, account for up to 75 % of the PSD at low frequencies as measured in the spanwise SPIV (see figure 7). As a result, the low-frequency dynamics of TSBs may be underrepresented or misrepresented in such simulations, potentially affecting both qualitative understanding and quantitative model validation.
7. Conclusions
We have investigated the dynamics of attached and separated turbulent flows over the Boeing Gaussian Bump using both data-driven and physics-based approaches. Both flows exhibit a coherent low-frequency dynamics characterised by streamwise-elongated streaky structures. Notably, these low-frequency structures are not a unique signature of fully separated flow, challenging the conventional view that they arise solely from separation bubble breathing. Instead, they are already present in an attached flow state, so they may be a precursor or indicator of incipient separation.
In the separated flow, the dynamics is clearly modal (in the global sense), associated with a three-dimensional zero-frequency eigenmode that likely corresponds to a centrifugal instability. The instability produces a prominent spanwise-standing-wave pattern, which dominates the very-low-frequency spectral content. In contrast, the attached flow’s dynamics shows no evidence of a modal origin, and a non-modal amplification mechanism underlying the coherent streaky structures seems more likely. This delineation connects the Gaussian bump to recent work on APG-induced TSB breathing mechanisms. However, our results do not justify a definitive classification of the attached flow’s low-frequency dynamics as either modal or non-modal. Compared with the separated flow, this dynamics is weaker, and no standing-wave structure emerges.
The strong modal dynamics in the separated case helps explain persistent challenges in simulating this flow, even with high-fidelity CFD. The most prominent low-frequency dynamics is associated with spanwise wavelengths corresponding to roughly
$20\,\%$
to
$30\,\%$
of the wind tunnel width, exceeding the periodic span of many simulations. Additionally, the use of spanwise-periodic boundary conditions further restricts the resolved spanwise wavenumbers to integer multiples of the fundamental wavenumber. Small spanwise extents and periodic boundary conditions can therefore exclude the dominant low-frequency modes, especially half-integer standing-wave patterns resulting from sidewall reflections. These findings offer an explanation for discrepancies between simulations on spanwise-periodic domains and experimental measurements, contributing to a better understanding of the SBS and TSB dynamics. This study also highlights the potential of combining data assimilation with linear mean field methods to investigate the flow dynamics using mean-flow data from limited regions, an approach that can be applied to other flows. Open questions for future research include the mechanisms driving low-frequency streaky structures in the attached flow, the influence of spanwise confinement, and the full three-dimensional structure of the dynamics.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2026.11629.
Acknowledgements
The authors gratefully acknowledge P. Gray for providing the PIV snapshot data from the Smooth Body Separation Experiment, which was essential for validating the present analysis. The authors also thank P.S. Iyer for providing the wall skin friction data from their three-dimensional LES.
Funding
Funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - 506170981, 504349109. This work has been supported by the German Academic Exchange Service (DAAD).
Declaration of interests
The authors report no conflicts of interest.
Author contributions
This work builds upon R.K’.s Master’s thesis. R.K.: methodology, investigation, analysis, software, visualisation, writing–original draft, writing–review and editing. L.F.: methodology, analysis, writing–review and editing. G.R.: supervision, writing–review and editing. K.O.: supervision, funding acquisition, writing–review and editing. J.v.S.: conceptualisation, methodology, supervision, funding acquisition, writing–original draft, writing–review and editing. All authors approved the final manuscript.
Data availability statement
The data from the experiment at the University of Notre Dame are openly available in the NASA Turbulence Modeling Resource at https://tmbwg.github.io/turbmodels/OtherexpData/speedbumpsepexp.html. The data-assimilated mean flows are available at https://git.tu-berlin.de/pida-group. The code used for the linearised analyses is available at https://gitlab.com/felics-group/FELiCS.
Use of artificial intelligence (AI) tools
Artificial intelligence tools (DeepL and ChatGPT) were used solely for language refinement.
Appendix A. Effects of discounting resolvent analysis for the separated case
To highlight the relationship between RA gain and the eigenvalues of the linearised operator, and to demonstrate how discounting modifies the resulting spectrum, we compute a pseudospectrum slice (at a fixed frequency) of the resolvent operator for the separated case. In particular, we examine the slice at
${\textit{St}}_h=0$
, where the unstable eigenvalues are located. This pseudospectrum slice is shown in figure 17(a). The default RA gain spectrum corresponds to evaluating the resolvent along the real axis (i.e. at
$\mathrm{Im}(\omega )=0$
) for each
${\textit{St}}_h$
. This spectrum is shown as a heatmap in figure 17(b), and clearly shows the imprint of singularities at zero frequency and at spanwise wavenumbers
$n\approx 2.5$
and
$\approx 12.5$
, where the resolvent is evaluated very close to eigenvalues of the linear operator. These peaks obscure the underlying dependence on spanwise wavenumber at low frequencies. In contrast, the discounted RA gain spectrum shown in figure 10(b) is obtained by evaluating the resolvent at
$\mathrm{Im}(\omega )=0.5$
, which shifts the resolvent away from the eigenvalues and removes these singularities. The resulting gain distribution reveals a trend over
$n$
that closely matches the LSA growth rate.
Relationship between resolvent gain and eigenvalues of the linear operator for
$\textit{Re}=2\times 10^6$
(separated). Slice through the pseudospectrum at
$St_h=0$
(a) and gain heatmap at
$\mathrm{Im}(\omega )=0$
(b).
Alignment between the modes from discounted (
$\gamma =0.5$
) and default (
$\gamma =0$
) RA (a) and comparison of the mode shapes with the lowest alignment (b). Here,
$\textit{Re}=2\times 10^6$
(separated).
Although discounting dramatically affects the RA gain spectrum, the mode shapes remain very similar. This is demonstrated through figure 18(a), where the alignment between the modes from discounted (
$\gamma =0.5$
) and default (
$\gamma =0$
) RA is shown. The alignment was evaluated in the region
$0\leqslant x\leqslant 0.5$
and
$0\leqslant y \leqslant 0.2$
, which corresponds to the primary region of interest in this study. It should be noted, that different values would be obtained, should the alignment be computed over a different region. For most
${\textit{St}}_h$
and
$n$
,
$A\gt 0.98$
, but there are some modes with lower values, down to
$A\approx 0.85$
. Notably, these modes are not associated with a particularly high gain, and are therefore not of special interest. The mode pair with the lowest alignment (black marker in figure 18(a)) is shown in figure 18(b) for comparison. Although some differences between the mode shapes are clearly identified, these would likely not change the interpretation of the mode, and we therefore conclude that the effect of discounting on the RA mode shapes is negligible in the present case.
Appendix B. Surface-pressure signature of resolvent modes
In the experimental study of TSB flows, instantaneous surface-pressure measurements are often applied to gain insight into the flow dynamics (e.g. Weiss et al. (Reference Weiss, Mohammed-Taifour and Schwaab2015), Mohammed-Taifour & Weiss (Reference Mohammed-Taifour and Weiss2016), Gray (Reference Gray2023)). We therefore assess to what extent coherent structures in the flow translate into surface-pressure fluctuations. This is not immediately obvious: for example, we observe an additional regime of elevated gain in the separated case, around
$n=4$
and
${\textit{St}}_h=0.08$
, which is slightly below half the characteristic frequency of the medium-frequency regime in this case (see figure 10(b)). We refer to this as the intermediate regime. In this regime, the surface-pressure SPOD yields comparatively low separation between the leading and sub-leading modes (see figure 6); therefore, the observation of high RA gain is surprising.
Pressure signature of the resolvent modes at the bump surface as a function of Strouhal number and spanwise wavenumber for the attached (a) and separated cases (b), shown as a heatmap with a logarithmically scaled colour bar. Black markers indicate Strouhal number–spanwise wavenumber pairs discussed in the previous sections.
In order to facilitate a direct comparison of the RA spectrum with the surface-pressure SPOD spectra (figure 6), the pressure signature of the RA modes is computed by extracting the squared magnitude of the pressure component at the
$N$
sensor locations from the RA mode, and averaging over the number of sensors
\begin{equation} \tilde {p} = \frac {\sum _{i=1}^N|\hat {p}|^2(\boldsymbol{x}_i)}{N}, \end{equation}
where
$|\hat {p}|$
is the magnitude of the pressure component of the resolvent mode. Note that the TKE norm of the resolvent mode corresponds to the gain value,
$\sigma$
, at the respective
${\textit{St}}_h$
and
$n$
. The resulting pressure signature is shown in figure 19 as contours over
${\textit{St}}_h$
and
$n$
. In the attached case, the pressure signature is most significant in the medium-frequency regime. Interestingly, whereas the RA gain increases with
$n$
at this frequency, the response at lower
$n$
induces the strongest pressure signature at the sensor locations. In the separated case, the pressure signature displays pronounced low- and medium-frequency regimes. In the low-frequency regime, the pressure signature displays a similar dependence on the spanwise wavenumber as the resolvent gain, with a preferential
$n\approx 6.5$
, and vanishing for
$n\to 0$
. In the medium-frequency regime, the highest pressure signature is observed for
$n=0$
, a similar observation as in the attached case. Remarkably, the intermediate regime identified in the resolvent gain spectrum does not produce a distinct pressure imprint at the sensor locations. This highlights a critical point: the modes most amplified (in terms of the resolvent gain, based on TKE) are not necessarily those most observable in the pressure near the wall, and thus may not always be captured through surface measurements alone.
As a further step, we compare the shape of the modelled pressure fluctuation and the experimentally observed structures. As shown in figure 19, the largest pressure response is in the medium-frequency regime and associated with low spanwise wavenumbers. We therefore focus on the two-dimensional limit
$n=0$
for this analysis. Representative frequencies for the comparison were chosen according to the highest RA gain in the medium-frequency regime at
$n=0$
for the respective case. These are indicated in the RA gain and surface-pressure spectra (figures 10, 19) through square markers. Figure 20 shows the comparison of the RA pressure response mode, extracted at the bump surface, with the leading mode from the surface-pressure SPOD. For presentation purposes, the amplitude is scaled by its maximum value and the phase angle is unwrapped and set to zero at the position of the first sensor. The figure demonstrates the qualitative similarity of the mode shapes. A quantitative measure for the similarity of the modes is provided through the mode alignment. The pressure modes shown in figure 20 have alignment values of
$A=0.98$
for the attached and
$A=0.90$
for the separated case. These high alignments demonstrate the capability of RA to model the dominant coherent structure at this frequency.
Comparison of coherent surface-pressure fluctuations from measurements extracted using SPOD and from RA modelling for the attached (a) and separated (b) cases. Black plots and the left axis show the magnitude, normalised by its respective maximum value, and blue plots and the right axis show the unrolled phase angle. The bump geometry and positions of the pressure sensors (purple crosses) are included for reference. RA results correspond to a spanwise wavenumber of
$n=0$
.
Appendix C. Alignment between zero-frequency resolvent modes and eigenmodes
Alignment between all LSA modes at
${\textit{St}}_h\approx 0$
and the RA mode at
${\textit{St}}_h=0$
for the attached (a) and separated cases (b), shown as a heatmap over the spanwise wavenumber and LSA mode number.
In the LSA spectrum of the attached case, we did not identify any eigenvalue branch with distinctly elevated growth rates. This suggests that the dynamics identified with RA for this case is non-modal. To further support this hypothesis, we assess whether the RA modes are associated with a specific eigenmode.
The LSA solver computes 10 eigenmodes with close-to-zero frequencies for each
$n$
. In order to assess their similarity with the RA mode computed at
${\textit{St}}_h=0$
, we compare the respective modes based on their alignment. Specifically, at each
$n$
, we compute the alignment between all LSA modes
$({\textit{St}}_h\approx 0,n)$
with the RA mode
$({\textit{St}}_h=0, n)$
. The result is shown in figure 21. In the attached case, the alignment for all pairs is very low. We therefore conclude that the
${\textit{St}}_h=0$
RA mode identified in this case does not correspond to an eigenmode of the system. In the separated case, on the other hand, much higher alignment values are observed. In the range of approximately
$2\lt n\lt 15$
, a single LSA mode is very highly aligned with the RA mode at the respective
$n$
. This range of spanwise wavenumbers corresponds to the same range where the corresponding eigenvalue branch has a distinctly higher growth rate compared with the remaining eigenvalues (see figure 8). For
$n\lt 2$
or
$n\gt 15$
, there is no single eigenvalue with very high alignment. Instead, all eigenmodes are moderately aligned with the RA mode. We note that for these
$n$
, the growth rate of the respective eigenvalue is no longer elevated compared with the mode cloud (see figure 8).
Appendix D. Effect of Reynolds number on the separated flow
In their initial study, Williams et al. (Reference Williams, Samuell, Sarwas, Robbins and Ferrante2020) reported a regime of approximate Reynolds number invariance of the separated mean flow. Later experimental (Gray Reference Gray2023) and computational (Iyer & Malik Reference Iyer and Malik2023a
; Uzun & Malik Reference Uzun and Malik2025) studies mostly confirmed this observation. However, whereas the separation point and separated flow region are largely independent of the Reynolds number, differences in the boundary layer development have been reported (Uzun & Malik Reference Uzun and Malik2025). In order to assess the influence of the Reynolds number on the separated flow dynamics, we consider an additional flow case at
$\textit{Re}=4\times 10^6$
. The data-assimilated mean flow for this additional case was obtained through the same procedure as for the other cases. However, it was not included in our previous publication (Klopsch et al. Reference Klopsch, Fuchs, Rigas, Oberleithner, Von and Jakob2025) due to scope limitations.
Effect of increasing Reynolds number on the separated flow. (a) Comparison of mean-flow streamlines and eddy-viscosity contours between
$\textit{Re}=2\times 10^6$
and
$4\times 10^6$
, (b) LSA eigenvalue spectrum for
$\textit{Re}=4\times 10^6$
, (c) resolvent gain heatmap for
$\mathrm{Re}=4\times 10^6$
(d) and alignment between RA modes and SPOD modes at
$\textit{Re}=4\times 10^6$
.
Mean-flow streamlines and eddy-viscosity contours of this additional case and the separated case at
$\textit{Re}=2\times 10^6$
are compared in figure 22(a). Note that the same contour levels for the eddy viscosity are chosen for both cases. Evidently, the mean flows are very similar, but there are some differences: in the additional case, the flow in the separated region is slightly less tilted towards the wall. Moreover, the extent of the region with high eddy viscosity is larger in this case. Nevertheless, the mean flows are very similar, supporting the approximate Reynolds number independence of the mean flow (outside the boundary layer) in this Reynolds number regime. We proceed by performing LSA and RA on this additional mean flow in order to investigate whether these small differences change the dynamics of the flow predicted by the linear models. The LSA eigenvalue spectrum is shown in figure 22(b). Here we observe the same distinct branch of
${\textit{St}}_h=0$
eigenvalues as in the separated case. However, the maximum growth rate occurs at slightly lower spanwise wavenumbers. The RA gain spectrum is shown in figure 22(c). The spectrum is very similar to the separated case, but the highest gain in the low-frequency regime is shifted to slightly lower spanwise wavenumbers, reflecting the trend observed in the LSA. Nevertheless, the results of the linear analyses are remarkably similar between this additional case and the separated case. To validate the RA modes for this additional case, we perform SPOD on streamwise PIV snapshot data for this case as well, and compute the alignment between the RA and respective SPOD modes. The resulting alignment is shown in figure 22(d). For conciseness, the mean over both FOVs is shown. Very high alignments are observed in the low-frequency regime, demonstrating that this approximate Reynolds number independence is not an artefact of the linear models. Note that the Nyquist frequency of the SPOD is
${\textit{St}}_h=0.11$
in this case, which excludes most of the medium-frequency regime. The validation is therefore limited to the low-frequency regime. Based on this analysis we conclude that the dominant dynamics of the separated flow in the range
$2\times 10^6\leqslant Re\leqslant 4\times 10^6$
do not significantly depend on the Reynolds number.
Open access
