Turbulent channel flow controlled by spanwise wall oscillations is studied using direct numerical simulations to improve how spanwise forcing reduces skin-friction drag. Harmonic wall oscillations generate a periodic transverse Stokes layer whose thickness $\delta$ is determined by the forcing period $T$. Although an optimal $T$ that maximises drag reduction is known to exist, its physical significance remains unclear. To elucidate it, we extend the spanwise Stokes layer by augmenting wall oscillation with an additional spanwise body force. In this formulation, $\delta$ and $T$ become decoupled and can be varied independently. The oscillating wall thus appears as a special and suboptimal case of spanwise forcing. Optimal performance is obtained for substantially smaller $T$ and larger $\delta$ than those of the classical Stokes layer. For the conditions examined, with Reynolds number and forcing amplitude held fixed, the maximum drag reduction increases by approximately one third, while the maximum net energy saving improves markedly from $-35\,\%$ to $+16\,\%$. These findings suggest that drag-reduction strategies based on spanwise forcing deserve renewed scrutiny: wall oscillation represents only one possible actuation method, and not necessarily the most effective one.

We consider free-surface flows driven by turbulence beneath the surface, particularly the strong free-surface turbulence (FST) regime, characterised by large Froude number ${\textit {Fr}}^2_T=\varepsilon /u_{\textit{rms}} {g}\,\gtrsim 0.1$. We study the surface layer, where air and water are highly mixed and turbulence modelling is challenging. We develop a definition of the surface-layer thickness $\delta _s$ based on the vertical derivative of intermittency $\gamma$ at the mean free surface $\bar {\eta }$, which, unlike previous definitions, is independent of the tail behaviour of $\gamma$. From direct numerical simulation (DNS) of statistically stationary, horizontally homogeneous strong FST, we show that scaling by $z^* = (z-\bar {\eta })/\delta _s$ collapses $\gamma$ across a wide range of ${\textit {Fr}}^2_T\in [0.03,0.3]$. The distribution more closely follows logistic rather than Gaussian behaviour. From the near-surface turbulence obtained from DNS, we make two general observations. First, we show that for strong FST there is minimal direct effect of the free surface on the isotropy, turbulence kinetic energy $\tilde {k}$ or dissipation rate $\varepsilon$ beneath the surface layer ($z^*\lt -0.5$). Instead, turbulence is only indirectly affected through the flux of kinetic energy into the surface layer. Second, we show that many relevant metrics within the surface layer ($z^*\in [-0.5, 0.5]$) collapse when appropriately scaled by $u_{\textit{rms}}^2=2\tilde {k}/3$ and $\varepsilon$ measured at $z^*=-0.5$. These observations suggest the possibility of a turbulence closure model which avoids direct modelling of $\tilde {k}$ and $\varepsilon$ in the surface layer. Towards this, we show that, across a wide range of ${\textit {Fr}}_T^2$, surface-layer thickness can be predicted by $\delta _s \approx 11.1 \,u_{\textit{rms}}^2 {g}^{-1}$ and energy flux into the surface layer by $W \approx 0.41 \,\varepsilon \delta _s$.

In this study, we develop a super-resolution (SR) model for homogeneous isotropic turbulence (HIT) inspired by the recently proposed low-inference-cost ResShift diffusion model. The training data are obtained from direct numerical simulation of two three-dimensional HIT cases with varying grid resolutions and Reynolds numbers ($ \textit{Re}_\lambda = 94$ and 173) to increase the model’s generalisability. The model is trained on two-dimensional snapshots rather than full three-dimensional fields, as training and inference on three-dimensional data would increase the computational cost significantly. Both the data from the whole domain and the data from a quarter of the domain are considered in the dataset to increase the diversity and quantity of training samples. This strategy also helps the model learn more localised flow structures and reduces dependence on global domain-specific patterns. The model is trained using single snapshots of velocity components for three upsampling factors of 4, 8 and 16. To assess the generalisability of the trained model, it is tested for flows under conditions different from those of the training data. Additionally, the high-resolution reconstruction of flow fields from low-resolution turbulent boundary layer data is performed to evaluate the model’s performance in anisotropic turbulence. The results show that the diffusion model presented in this study performs well in predicting the velocity field even for high upsampling factors, and unlike bicubic interpolation, convolutional neural network (CNN)- and U-Net-based models, it does not generate a visually blurry flow field when applied to high upsampling factors. It also outperforms bicubic interpolation, CNN- and U-Net-based models, as well as the traditional conditional denoising diffusion probabilistic model designed for SR, in predicting flow statistics. The model effectively extracts flow features, generates flow structures of varying sizes and shows strong performance in predicting vorticity. It also reproduces the energy spectrum at high wavenumbers with reasonable accuracy, indicating the recovery of small-scale structures often lost in coarse data. This capability is valuable for subgrid-scale stress estimation and helps improve the physical fidelity of large eddy simulation frameworks.

The mechanical response of elastic porous media confined within rigid geometries is central to a wide range of industrial, geological and biomedical systems. However, current models for these problems typically overlook the role of wall friction and particularly its interaction with confinement. Here, we develop a theoretical framework to describe the interplay between the mechanics of the medium and Coulomb friction at the confining walls for slow, quasistatic deformations in response to two canonical uniaxial forcings: piston-driven loading (i.e. an imposed effective stress at the top boundary) and fluid-driven loading (i.e. an imposed fluid pressure at the top boundary) followed by unloading. We find that, during compression, the stress field evolves according to a quasistatic advection–diffusion equation, extending classical poroelasticity results. The magnitude of friction is controlled by a single dimensionless number ($\mathcal{F}$) proportional to the friction coefficient and the aspect ratio of the confining geometry. During decompression, a portion of the solid matrix remains stuck due to friction, leading to hysteresis and to the propagation of a slip front. In piston-driven loading the frictional stress is directly coupled to the solid effective stress, leading to exponential damping of the loading and striking changes to the displacement field. However, this coupling limits the energy dissipated by friction. In fluid-driven loading the pressure gradient locally adds energy, decoupling elastic energy storage and frictional energy dissipation. The displacement remains qualitatively unchanged but is quantitatively reduced due to large energy dissipation. In both cases, friction can have a substantial impact on the apparent mechanical properties of the medium.

We investigate thermal boundary layer (BL) asymmetry in turbulent Rayleigh–Bénard convection (RBC) under both spherical and annular geometries using different BL theories. Unlike planar RBC, the spherical and annular configurations exhibit asymmetric thermal BLs near the inner and outer boundaries due to boundary curvature and non-uniform radial gravity. We generalise three BL frameworks – the Prandtl–Blasius BL model, the steady free-convective model and the fluctuating BL model – and apply them to both geometries. Direct numerical simulations (DNSs), based on the Oberbeck–Boussinesq equations, are performed in three-dimensional spherical RBC and three-dimensional annular RBC for various radius ratios ($\eta$), gravity profiles and also Prandtl numbers ($ \textit{Pr}$), to compare with the predictions of the extended BL models. We find that the BL asymmetries predicted by both the extended steady free-convective BL and the fluctuating BL agree well with DNS results, with the fluctuating BL model providing the best agreement for the mean temperature profiles. A force-balance analysis further shows that this better performance is consistent with the DNS observation that, in the wall-normal direction within the thermal BL, buoyancy is balanced primarily by the pressure-gradient force. This is consistent with the assumption underlying the steady free-convective and fluctuating BL models. Moreover, the fluctuating BL model explicitly accounts for the contribution of turbulent fluctuations to the heat flux, which further improves its agreement with the DNS mean temperature profiles. We derive analytical expressions for the bulk temperature and the thermal BL thickness ratio as functions of the radius ratio and gravity profile across different Prandtl-number regimes. These expressions are obtained by integrating the similarity thermal equation for both the inner and outer BLs using an approximate similarity streamfunction, and by closing the solutions through a heat-flux matching condition. The resulting leading-order expressions obtained from both the steady free-convective and fluctuating BL models are shown to be the same, and they agree well with DNS data. This analytical result provides a robust and practical tool for quantifying BL asymmetry in curved RBC systems.

We investigate the incompressible flow inside a two-dimensional square cavity, driven by the sliding motion of its four lids, all at the same speed and with facing lids moving in opposite directions. The problem has three symmetries: two mirror symmetries with respect to the diagonals and a $\pi$ rotation invariance about the centre of the cavity. The base flow, a steady state that has all three symmetries, is the unique solution at sufficiently low values of the Reynolds number ($ \textit{Re}$) and acts as a global attractor. At higher $ \textit{Re}$, it has become unstable and shares the phase space with a globally attracting space–time symmetric periodic orbit that, in addition to the rotational invariance, is also invariant under evolution over half a period followed by reflection about either of the diagonals. In between, a wealth of solution branches and intervening bifurcations mediate the transition process. In particular, a pair of steady states that break the mirror symmetries but are mirror-symmetry images of each other regulate the appearance and disappearance of a second space–time symmetric periodic orbit and a pair of asymmetric periodic orbits that are also mirror images of each other. The catalogue of instabilities includes both local (two pitchfork, two Hopf, a saddle-node and a cyclic fold) and global (two heteroclinic and one homoclinic) bifurcations. The sequence of transitions is explained in terms of a one-dimensional path through the parameter space of a codimension-four bifurcation: the double zero bifurcation with Z$_2$ symmetry and degeneracy of the third order terms.

We perform direct numerical simulations of continuously growing broadband surface waves forced by a turbulent atmospheric boundary layer coupled with a developing underwater current. We resolve and analyse the multiscale space–time evolution of the waves by considering the wave spectrum in frequency and wavenumber space and describe the kinematics of nonlinear gravity–capillary waves under a current initially described by a viscous boundary layer and transitioning to turbulence at later times under the wind-wave forcing. The wave speed experiences a scale-dependent Doppler shift, with shorter waves shifted by currents closer to the surface, in agreement with the framework from Stewart & Joy (1974 Deep Sea Res. Oceanogr. Abstracts 21(12), 1039–1049). At low wave slopes, the wave energy concentrates along the linear dispersion relation. When the wave slope is high enough, we observe wave energy located in multiple branches associated with nonlinear bound harmonics travelling at the speed of a carrier mode. These nonlinear branches are well described by a generalized nonlinear dispersion relation that links each harmonic to the effective velocity of the carrier mode to which they are bound, and are found to be Doppler shifted with the carrier mode. The generalized Doppler-shifted nonlinear dispersion relation remains valid as the underwater current becomes turbulent, and the depth-varying mean current profile can be systematically reconstructed from the measured phase velocities from waves at different scales.

We show that spatio-temporal non-Markovianity of a Gaussian random synthetic velocity field is an essential property for modelling turbulent mixing. We demonstrate this using synthetically generated Gaussian incompressible velocity fields for passive scalar mixing. Including a separate velocity decorrelation time scale for each spatial scale (random sweeping) yields an essentially non-Markovian velocity field with a finite time memory decaying as $\tau ^{-6}$ (for a decaying spectrum) instead of an exponential decay (Markovian), which is obtained by including a constant time scale for all spatial scales, irrespective of the filtering function. We characterise the Lagrangian mixing statistics of both the Markovian and the non-Markovian synthetic fields and compare them against a corresponding incompressible direct numerical simulation (DNS). We show that the average pair dispersion is well captured by the non-Markovian fields across the ballistic, inertial and diffusive regimes. We also study diffusive passive scalar mixing in the Schmidt number range $\textit{Sc}\leqslant 1$ using the DNS and the synthetic fields. Both the synthetic fields recover the $-17/3$ scalar spectrum for low Schmidt numbers and inertial subrange in kinetic energy spectra. However, the mean fluctuation gradient magnitudes are severely under predicted by the Markovian synthetic fields compared with the non-Markovian synthetic fields. Additionally, the fluctuation gradients parallel to the mean gradient exhibit smaller skewness when stirred by the Markovian synthetic field compared with the non-Markovian fields. Finally, we show that the non-Markovian synthetic fields perform better in decaying scalar gradient simulations initialised by a concentrated sphere with high passive scalar concentration. Throughout, we compare our results with companion three-dimensional DNS to show the necessity of non-Markovianity in synthetic fields to capture mixing dynamics.

This work employs structured input–output analysis (SIOA) augmented by an eddy viscosity model (SIOA-e) to investigate turbulent flows over rigid and compliant walls. The SIOA-e framework demonstrates the capability in identifying both streamwise and spanwise dominant characteristic wavelengths for rigid wall turbulence. For compliant walls, the SIOA-e method predicts optimal compliant wall parameters associated with positive damping coefficients when minimizing input–output gain for near-wall cycle and very large-scale motions, respectively. The reduction of input–output gain due to the compliant wall is achieved by wall displacement resembling blowing and suction opposite to the wall-normal velocity of dominant streamwise vortices. However, optimized compliant wall parameters based on specific wavenumber–frequency combinations may amplify flow structures for other wavenumber–frequency pairs, potentially leading to an overall drag increase. For example, compliant wall parameters tuned for suppressing large-scale structures can affect both large- and small-scale structures. We also employ input–output analysis to predict convective velocity of wall displacement and pressure for turbulent flow over the compliant wall, and the predicted convective velocity of wall displacement is 0.53 times centreline velocity, which aligns well with recent experimental measurements.

We consider the response of a flexibly mounted square prism placed in inertial-viscoelastic fluid flow with one degree-of-freedom in the cross-flow direction. Under these flow conditions, both inertia and elastic effects are significant. We model the system numerically using a two-way coupling scheme to simulate the interaction between the fluid and the spring–mass system at a Reynolds number of $\textit{Re}=200$ for two mass ratios of $m^* = 2$ and 20, and a Weissenberg number of $\textit{Wi}=2$, across a range of reduced velocities. We demonstrate that introducing fluid elasticity suppresses vortex-induced vibrations of square prisms, consistent with prior findings for circular bluff bodies. However, we find that fluid elasticity amplifies the galloping response in comparison with the response in a Newtonian fluid, leading to larger oscillation amplitudes and the onset of galloping at lower reduced velocities. The predicted enhancement in galloping is significant, particularly at low mass ratios, where no galloping is observed over the wide reduced velocity range tested for Newtonian fluids. We show that this enhancement of galloping is likely the result of the observation that the addition of viscoelasticity increases the magnitude of the rate of change of the transverse flow-induced force on the prism with increasing angle of attack of the incoming flow.

We investigate the influence of side-wall wetting on the linear stability of falling liquid films confined in the spanwise direction. A biglobal stability framework is developed, capturing inertia, viscosity, gravity, capillarity and geometric confinement. The base flow exhibits a curved meniscus and a streamwise velocity overshoot near the side walls. Linear stability analysis based on the Navier–Stokes equations is performed in two limiting regimes. In confined channels, where spanwise confinement stabilises moderate-wavenumber perturbations via side-wall boundary layers, wetting weakens this stabilisation; as the contact angle decreases, the neutral curves shift towards the unconfined one-dimensional limit, thus wetting acts as a relative destabilising mechanism. In contrast, in weakly confined channels where side-wall boundary layers do not provide confinement-induced stabilisation, wetting produces a net long-wave stabilisation ($k \rightarrow 0$), significantly increasing the critical Reynolds number. This effect strengthens as the contact angle decreases, indicating a competition between destabilising inertia and stabilising wetting-induced capillary forces. The predicted long-wave stabilisation effect is compared quantitatively with available experimental measurements, showing consistent trends and comparable magnitudes within the accessible parameter range. Perturbation eigenmode structures show that, in confined channels, the relative destabilisation is associated with near-wall vortical structures induced by the meniscus elevation and velocity overshoot, which reduce effective viscous damping. In contrast, in weakly confined channels, stabilisation is consistent with interface tensioning through strong anchoring of the perturbations at the side walls.

The interaction between acoustic waves and turbulent grazing flow over an acoustic liner is investigated using lattice-Boltzmann very-large-eddy simulations. A single-degree-of-freedom liner with 11 streamwise-aligned cavities is studied in a grazing flow impedance tube. The conditions replicate reference experiments from the Federal University of Santa Catarina. The influence of grazing flow (with a centreline Mach number of 0.32), acoustic wave amplitude, frequency and propagation direction relative to the mean flow is analysed. Impedance is computed using both direct (i.e. the in situ method) and model-fitting inference (i.e. the mode-matching) methods. The former reveals strong spatial variations; however, averaged values throughout the sample show minimal differences between upstream- and downstream-propagating waves, in contrast to what is obtained with the latter method. Flow analyses reveal that the orifices displace the flow away from the face sheet, with this effect amplified by acoustic waves and dependent on the wave propagation direction. Consequently, the boundary layer displacement thickness ($\delta ^*$) increases along the streamwise direction compared with a smooth wall and exhibits localised humps downstream of each orifice. The growth of $\delta ^*$ alters the flow dynamics within the orifices by weakening the shear layer at downstream positions. This influences the acoustic-induced mass flow rate through the orifices at equal sound pressure level, suggesting that acoustic energy is dissipated differently along the liner. The asymmetry of the flow field experienced by the acoustic wave, depending on its propagation direction, highlights the need to consider a spatially evolving turbulent flow when studying the acoustic–flow interaction and measuring impedance.

The extreme heat fluxes characteristic of hypersonic flows significantly limit the flight envelope of hypersonic vehicles. The role of hydrodynamic instability and the onset of laminar-to-turbulent boundary layer transition is of notable importance. The effect of streaks on the suppression of planar (second Mack mode) instabilities has been previously investigated, but a potentially passive and non-intrusive control method has not been established yet. Recent work shows that streaks can be generated through a spanwise variation in surface temperature. This method exploits the aerothermodynamic characteristics of the flow, and therefore promises to be robust. This work uses direct numerical simulations to determine and quantify the effectiveness of this novel control method in the suppression of second Mack mode instability for a hypersonic boundary layer over a flat plate. The computational analyses cover a range of Mach numbers, 4.8–6, and wall temperature ratios representative of both wind tunnel testing and flight scenarios. Among the range of configurations investigated, the energy of the second Mack mode is reduced by up to approximately 60 % by the steady streaks. The streak wavelength parameter plays a significant role in the stabilisation benefits. For a Mach 6 configuration, for the most linearly amplified second Mack mode disturbance frequency, nearly optimum performance is achieved for a spanwise wavelength of approximately 8–10 times the local boundary layer thickness. These findings open new avenues for controlling hypersonic boundary layers and offer valuable guidance for future experimental campaigns aimed at validating this novel control strategy.

Attached cavitation is likely the most common form of developed hydrodynamic cavitation, yet the reason for its dominance remains unclear. From the experimental side, a natural approach is to seed controllable nuclei and observe their evolution. We propose a laser-based on-demand nucleation method that generates micro- and nanobubbles as nuclei in Venturi flows, enabling unprecedented spatio-temporal control of hydrodynamic cavitation inception. For single-bubble cases, we find that attached cavitation occurs when the bubble surface enters the boundary layer of the channel where the pressure is below the vapour pressure. Based on it, we construct a phase diagram of cavitation regimes as a function of cavitation number and non-dimensional wall distance. Extending to multiple bubbles, assuming a random spatial distribution of nuclei within the laser-illuminated region, we develop a simple model to estimate the probability of attached cavitation. Results show that, at typical cavitation numbers, only a few bubbles suffice for attached cavitation to occur with nearly 100 % probability. Our finding provides new insights into why nuclei in hydrodynamic processes tend to develop into attached cavitation.

Bubble pairs are effective modulators of liquid jets. We investigate the jetting of an air bubble driven by a laser-induced cavitation bubble using high-speed imaging, compressible volume-of-fluid (VoF) simulations and theoretical analysis. Three distinct jet types emerge, depending on the stand-off distance $\gamma$ and size ratio $\eta$ between the bubbles. Jet formation proceeds through two stages: an initial shock-induced acceleration followed by flow focusing on the concave liquid–air interface. We derive scaling relations, $V_0=1.1 p_0R_0/(\rho cR_l)((\gamma (1+\eta )-1)/\eta )^{-1.6}$ for the shock-driven stage and $V_m={}(1+(0.8-0.5\gamma )\eta ^{0.75})V_0$ for the flow focusing stage in the strong jet regime, both of which agree closely with experimental and numerical measurements. Here, $V_0$ and $V_m$ denote the velocity increments associated with shock-wave-induced acceleration and flow focusing stages, respectively. The variables $p_0$, $R_0$, $\rho$, $c$ and $R_l$ represent the initial pressure and radius of the cavitation bubble, the fluid density, the speed of sound in the liquid and the maximum volume-equivalent radius of the cavitation bubble, respectively. A $(\eta ,\gamma )$ phase diagram delineates the weak, strong and explosive jets, with regime boundaries accurately captured by the theoretically derived transitions.

Boiling and bubble injection are effective strategies for enhancing heat transfer between solid boundaries and a working fluid in numerous industrial applications, including nuclear reactors and molten metal processing. Motivated by this, we conduct direct numerical simulations of a vertical, turbulent, differentially heated, bubble-laden channel flow. The Prandtl number $\textit{Pr}$, kept identical in both phases, is varied across three representative values – $0.07$ (liquid metals), $0.7$ (vapour) and $7$ (water) – to span thermal transport regimes across three orders of magnitude. The simulations are conducted at a friction Reynolds number $\textit{Re}_\tau =150$, void fraction $\alpha =5.4\,\%$ and a density ratio $\rho _r=0.1$ (defined as the bubble-to-carrier density). The bubbles substantially alter the hydrodynamic structure of the flow, amplifying turbulent fluctuations and mixing. Their interaction with the thermal boundary layers disrupts the characteristic streaky structures near the heated walls, fragmenting them into smaller and more chaotic patterns. To elucidate this mechanism, we examine the bubble-induced modifications to the temperature field and show that temperature becomes decorrelated from velocity. Consequently, the heat-transfer enhancement arises primarily from an increase in convective heat flux driven by intensified wall-normal velocity fluctuations. The thermal boundary layer is markedly thinned, and the Nusselt number nearly doubles across all examined cases.

We numerically investigate the steady and unsteady wakes of three-dimensional permeable disks over Reynolds number ($\textit{Re}$) range 100–300 and Darcy number ($Da$) range $10^{-9}$–$10^{-3}$. For disks with low permeability ($Da\le 8\times 10^{-5}$), the dynamical transition route is the same as that of impervious disks, with the critical $\textit{Re}$ for all bifurcations increasing with decreasing permeability. In contrast, for disks with high permeability ($Da\ge 2\times 10^{-4}$), all unsteady bifurcations are suppressed, and the wake remains in a steady regime throughout the $\textit{Re}$ range considered. Interestingly, at moderate $Da$, permeability gives rise to two previously unreported flow regimes. The first is the ‘SVR breathing’ regime, occurring at $Da\approx 10^{-4}$ and $\textit{Re}\approx 200$, and is attributed to the subharmonic lock-in between two distinct unsteady dynamics: the shedding of hairpin vortices and the low-frequency unsteadiness of the near-wake recirculation regions. The second is the ‘intermittency’ regime, which occurs at $Da\approx 1.5\times 10^{-4}$, $\textit{Re}\approx 200$; the wake alternates irregularly between two periodic modes with orthogonal planes of symmetry. Future work might include verifying whether intermittency arises from the energy competition between two modes, as the vortices lack sufficient energy to sustain stable single-mode harmonic oscillations. These findings demonstrate that permeability can fundamentally alter wake dynamics and introduce new wake structures that do not occur on an impervious disk.

The drag wake of a towed inclined 6 : 1 prolate spheroid in unstratified and stratified ambients is investigated experimentally using stereoscopic particle image velocimetry. Tow speed, stratification strength and inclination angle are varied independently, resulting in a parameter space spanning Reynolds numbers $\textit{Re} = (1.25{-}20) \times 10^3$, Froude numbers $\textit{Fr}= U/\textit{ND} = 2{-}32$ and $\infty$, and inclination angles $\theta = 0^\circ$ and $20^\circ$. Measurements are repeated at each parameter combination to obtain converged wake statistics for $3 \leqslant x/D \leqslant 40$. Unstratified measurements provide a baseline experimental dataset for inclined spheroids that has not previously been reported. In the absence of stratification, inclination generates persistent wake asymmetries and a net vertical impulse that deflects the wake trajectory. Although inclined configurations exhibit larger initial wake heights than axisymmetric cases, the early wake evolution collapses when scaled by an effective body diameter, indicating that this increase is geometric in origin. Regular vertical velocity protrusions are observed in inclined wakes, with a characteristic spacing that depends on Reynolds number but shows no measurable dependence on Froude number. At sufficiently low Froude number, buoyancy influences the near-body flow, and modifies the wake trajectory and streamwise velocity profiles. For $\textit{Re} = 5000$, wake heights for both axisymmetric and inclined configurations collapse across stratification strengths when scaled by an effective diameter. In this regime, the wake trajectory exhibits oscillations with period $2\pi /N$, in agreement with previously reported stratified wake dynamics.

We perform causal analysis on the low-dimensional Galerkin model for shear flow developed by Moehlis et al. (New J. Phys., vol. 6, 2004, 56). Our method integrates both equation-based analysis and the proposed Galerkin-based Granger causality (GGC) to investigate the effect of the nonlinear terms on the dynamics. Two types of quadratic interactions are identified: a fully triadic interaction and a modulated two-mode coupling. The propagation of these interactions through the nonlinear dynamics leads to a directed cause-and-effect network. Furthermore, the relative importance of each mode amplitude on the dynamics of the target mode is quantified. This analysis provides a deeper understanding of the nonlinear dynamics and distills control opportunities. To demonstrate the applicability of the proposed GGC to realistic flows where Galerkin projection is impractical, a turbulent lid-driven cavity flow is further studied. We foresee applications of the proposed causal analysis framework as valuable tools for Galerkin modelling – guiding investigations of modal causality, prediction uncertainty, model-order reduction and control design.