The marginal ice zone represents the periphery of the sea ice cover. In this region, the macroscale behaviour of the sea ice results from collisions and enduring contact between ice floes. This configuration closely resembles that of dense granular flows, which have been modelled successfully with the $\mu (I)$ rheology. Here, we present a continuum model based on the $\mu (I)$ rheology that treats sea ice as a compressible fluid, with the local sea ice concentration given by a dilatancy function $\varPhi (I)$. We infer expressions for $\mu (I)$ and $\varPhi (I)$ by nonlinear regression using data produced with a discrete element method (DEM) that considers polygon-shaped ice floes. We do this by driving the sea ice with a one-dimensional shearing ocean current. The resulting continuum model is a nonlinear system of equations with the sea ice velocity, local concentration and pressure as unknowns. The rheology is given by the sum of a plastic term and a viscous term. In the context of a periodic patch of ocean, which is effectively a one-dimensional problem, and under steady conditions, we prove this system to be well-posed, present a numerical algorithm for solving it, and compare its solutions to those of the DEM. These comparisons demonstrate the continuum model's ability to capture most of the DEM results accurately. The continuum model is particularly accurate for ocean currents faster than 0.25 m s$^{-1}$; however, for low concentrations and slow ocean currents, the continuum model is less effective in capturing the DEM results. In the latter case, the lack of accuracy of the continuum model is found to be accompanied by the breakdown of a balance between the average shear stress and the integrated ocean drag extracted from the DEM. Since this balance is expected to hold independently of our choice of rheology, this finding indicates that continuum models might not be able to describe sea ice dynamics for low concentrations and slow ocean currents.

Nonlinear electrokinetic phenomena, where electrically driven fluid flows depend nonlinearly on the applied voltage, are commonly encountered in aqueous suspensions of colloidal particles. A prime example is the induced-charge electro-osmosis, driven by an electric field acting on diffuse charge induced near a polarizable surface. Nonlinear electrohydrodynamic flows also occur in non-polar fluids, driven by the electric field acting on space charge induced by conductivity gradients. Here, we analyse the flows about a charge-neutral spherical solid particle in an applied uniform electric field that arise from conductivity dependence on local field intensity. The flow pattern varies with particle conductivity: while the flow about a conducting particle has a quadrupolar pattern similar to induced-charge electro-osmosis, albeit with opposite direction, the flow about an insulating particle has a more complex structure. We find that this flow induces a force on a particle near an electrode that varies non-trivially with particle conductivity: while it is repulsive for perfectly insulating particles and particles more conductive than the suspending medium, there exists a range of particle conductivities where the force is attractive. The force decays as the inverse square of the distance to the electrode and thus can dominate the dielectrophoretic attraction due to the image dipole, which falls off with the fourth power with the distance. This electrohydrodynamic lift opens new possibilities for colloidal manipulation and driven assembly by electric fields.

Linear and nonlinear contributions to the localization and dynamics of internal gravity waves in a stably stratified turbulent channel flow are investigated using data from direct numerical simulations (DNS). The classification into linear and nonlinear mechanisms is based on the resolvent formulation of the Navier–Stokes equations, which interprets velocity and temperature fluctuations (flow response) as the result of a linear operator (resolvent) acting on the nonlinear advection terms (forcing). Spatial and spatio-temporal power spectral densities computed from DNS data demonstrate that the stratified flow response is localized in spectral space and in the channel core, while the nonlinear forcing is broadband and spans up to the entire channel height. The localization of the velocity and temperature fluctuations in wavenumber and frequency is captured by the leading singular value of the resolvent operator. The wall-normal localization on the other hand results from a combination of linear dynamics and nonlinear forcing, and the latter is further examined using the cross-spectral density (CSD) tensor. Wall-normal subsets of the forcing CSD lead to flow responses that reveal a three-layer structure. The middle one hosts the critical layer of the gravity wave, and is termed the outer layer since it is flanked by an inner layer at the wall and the core region at the channel centre. Forcing within this outer layer generates the majority of the flow response in the channel core. Furthermore, a decomposition of the forcing CSD into velocity and temperature demonstrates that each imprints distinct phase relations on their associated responses, which lead to destructive interference and localization of the gravity waves in the channel core.

In a vertical channel driven by an imposed horizontal temperature gradient, numerical simulations (Gao et al., Phys. Rev. E, vol. 88, 2013, 023010; Phys. Rev. E, vol. 91, 2015, 013006; Phys. Rev. E, vol. 97, 2018, 053107) have previously shown steady, time-periodic and chaotic dynamics. We explore the observed dynamics by constructing invariant solutions of the three-dimensional Oberbeck–Boussinesq equations, characterizing the stability of these equilibria and periodic orbits, and following the bifurcation structure of the solution branches under parametric continuation in Rayleigh number. We find that in a narrow vertically periodic domain of aspect ratio 10, the flow is dominated by the competition between three and four co-rotating rolls. We demonstrate that branches of three- and four-roll equilibria are connected and can be understood in terms of their discrete symmetries. Specifically, the $D_4$ symmetry of the four-roll branch dictates the existence of qualitatively different intermediate branches that themselves connect to the three-roll branch in a transcritical bifurcation due to $D_3$ symmetry. The physical appearance, disappearance, merging and splitting of rolls along the connecting branch provide a physical and phenomenological illustration of the equivariant theory of $D_3$–$D_4$ mode interaction. We observe other manifestations of the competition between three and four rolls, in which the symmetry in time or in the transverse direction is broken, leading to limit cycles or wavy rolls, respectively. Our work highlights the interest of combining numerical simulations, bifurcation theory and group theory, in order to understand the transitions between and origin of flow patterns.

First predicted by Richtmyer in 1960 and experimentally confirmed by Meshkov in 1969, the Richtmyer–Meshkov instability (RMI) is crucial in fields such as physics, astrophysics, inertial confinement fusion and high-energy-density physics. These disciplines often deal with strong shocks moving through condensed materials or high-pressure plasmas that exhibit non-ideal equations of state (EoS), thus requiring theoretical models with realistic fluid EoS for accurate RMI simulations. Approximate formulae for asymptotic growth rates, like those proposed by Richtmyer, are helpful but rely on heuristic prescriptions for compressible materials. These prescriptions can sometimes approximate the RMI growth rate well, but their accuracy remains uncertain without exact solutions, as the fully compressible RMI growth rate is influenced by both vorticity deposited during shock refraction and multiple sonic wave refractions. This study advances previous work by presenting an analytic, fully compressible theory of RMI for reflected shocks with arbitrary EoS. It compares theoretical predictions with heuristic prescriptions using ideal gas, van der Waals gas and three-term constitutive equations for simple metals, the latter being analysed with detailed and simplified ideal-gas-like EoS. We additionally offer an alternative explicit approximate formula for the asymptotic growth rate. The comprehensive model also incorporates the effects of constant-amplitude acoustic waves at the interface, associated with the D'yakov–Kontorovich instability in shocks.

Turbulent circular pipe flows subjected to axial system rotation are studied using direct numerical simulations (DNS) for a wide range of rotation numbers of $Ro_b = 0\unicode{x2013}20$ at a fixed Reynolds number. To ensure that energetic turbulent eddy motions are captured at high rotation numbers, long pipes up to $L_z = 180{\rm \pi} R$ are used in DNS. Two types of energy-containing flow structures have been observed. The first type is hairpin structures that are characteristic of the turbulent boundary layer developing over the pipe wall for both non-rotating and axially rotating flows. The second type is Taylor columns forming at moderate and high rotation numbers. Based on the study of two-point autocorrelation coefficients, it is observed that Taylor columns exhibit quasi-periods in both axial and azimuthal directions. According to the premultiplied spectra, Taylor columns feature one single characteristic axial length scale at the moderate rotation numbers but two at high rotation numbers. It is discovered that the axial system rotation suppresses the sweep events systematically and impedes the formation of hairpin structures. As the rotation number is increased, the turbulence kinetic energy held by Taylor columns enhances rapidly associated with significant increases in their axial length scales.

Aerodynamic breakup of vaporizing drops is commonly seen in many spray applications. While it is well known that vaporization can modulate interfacial instabilities, the impact of vaporization on drop aerobreakup is poorly understood. Detailed interface-resolved simulations were performed to systematically study the effect of vaporization, characterized by the Stefan number, on the drop breakup and acceleration for different Weber numbers and density ratios. It is observed that the resulting asymmetric vaporization rates and strengths of Stefan flow on the windward and leeward sides of the drop hinder bag development and prevent drop breakup. The critical Weber number thus generally increases with the Stefan number. The modulation of the boundary layer also contributes to a significant increase of drag coefficient. Numerical experiments were performed to affirm that the drop volume reduction plays a negligible role and the Stefan flow is the dominant reason for the breakup suppression and drag enhancement observed.

A novel fast-running model is developed to predict the three-dimensional (3-D) distribution of turbulent kinetic energy (TKE) in axisymmetric wake flows. This is achieved by mathematically solving the partial differential equation of the TKE transport using the Green's function method. The developed solution reduces to a double integral that can be computed numerically for a wake prescribed by any arbitrary velocity profile. It is shown that the solution can be further simplified to a single integral for wakes with Gaussian-like velocity-deficit profiles. Wind tunnel experiments were performed to compare model results against detailed 3-D laser Doppler anemometry data measured within the wake flow of a porous disk subject to a uniform free-stream flow. Furthermore, the new model is used to estimate the TKE distribution at the hub-height level of the rotating non-axisymmetric wake of a model wind turbine immersed in a rough-wall boundary layer. Our results show the important impact of operating conditions on TKE generation in wake flows, an effect not fully captured by existing empirical models. The wind-tunnel data also provide insights into the evolution of important turbulent flow quantities such as turbulent viscosity, mixing length and the TKE dissipation rate in wake flows. Both mixing length and turbulent viscosity are found to increase with the streamwise distance. The turbulent viscosity, however, reaches a plateau in the far-wake region. Consistent with the non-equilibrium theory, it is also observed that the normalised energy dissipation rate is not constant, and it increases with the streamwise distance.

Vertical thermal convection is a non-equilibrium system in which both buoyancy and shear forces play a role in driving the convective flow. Beyond the onset of convection, the driven dissipative system exhibits chaotic dynamics and turbulence. In a three-dimensional domain extended in both the vertical and the transverse dimensions, Gao et al. (Phys. Rev. E, vol. 97, 2018, 053107) have observed a variety of convection patterns which are not described by linear stability analysis. We investigate the fully nonlinear dynamics of vertical convection using a dynamical-systems approach based on the Oberbeck–Boussinesq equations. We compute the invariant solutions of these equations and the bifurcations that are responsible for the creation and termination of various branches. We map out a sequence of local bifurcations from the laminar base state, including simultaneous bifurcations involving patterned steady states with different symmetries. This atypical phenomenon of multiple branches simultaneously bifurcating from a single parent branch is explained by the role of $D_4$ symmetry. In addition, two global bifurcations are identified: first, a homoclinic cycle from modulated transverse rolls and second, a heteroclinic cycle linking two symmetry-related diamond-roll patterns. These are confirmed by phase space projections as well as the functional form of the divergence of the period close to the bifurcation points. The heteroclinic orbit is shown to be robust and to result from a 1:2 mode interaction. The intricacy of this bifurcation diagram highlights the essential role played by dynamical systems theory and computation in hydrodynamic configurations.

The friction drag of the axial flow along the outer surface of a cylinder varies with the cylinder radius and flow conditions. This study included direct numerical simulations of the axial turbulent flow along a circular cylinder under different conditions for obtaining the turbulence statistics and wall friction coefficient. Then the characteristics of velocity streaks were observed from a geometrical perspective of turbulence structures around the circular cylinder, and compared with the characteristics of the turbulence structures in a boundary layer on a flat plate. The results showed that the velocity streak spacing and the distance between the velocity streak and the cylinder surface in the viscous length scale do not vary substantively with the radius of the cylinder, and are the same as those of the turbulent flow along a flat plate. Therefore, they can be considered geometrical characteristics of the turbulence structure independent of the cylinder radius. Moreover, the friction coefficient per pair of high- and low-speed velocity streaks is the same as that of flat-plate turbulent flow, independent of the cylinder radius, and can be regarded as a dynamical characteristic for a pair of velocity streaks. Two equations were derived based on the characteristics of wall turbulence. The characteristics of the turbulence predicted by the two formulae were consistent with the simulation results. Consequently, we showed that the wall friction coefficient and number of the velocity streak pairs, which are statistical and structural characteristics of wall turbulence, can be predicted appropriately by specifying the radius Reynolds number.

Spanwise vortex instability and the growth of secondary hairpin-like vortical structures in the wake of an oscillating foil are investigated numerically at Reynolds number 8000 in a range of chord-based Strouhal number ($0.32 \le St_c \le 0.56$). The phase-offset ($\phi$) between the heaving and pitching motion is $\phi = 90^\circ$. The wake at the lowest $St_c$ (0.32) is characterized by a single system of streamwise hairpin-like structures that evolve from the core vorticity outflux of the secondary leading edge vortex (LEV) over the foil boundary. The primary LEV features spanwise dislocations, but it does not reveal substantial changes advecting downstream. Increasing $St_c$ beyond 0.32 reveals that the transition in spanwise instability characterizes the deformation of primary LEV cores, which subsequently transforms to hairpin-like secondary structures. At higher $St_c$, stronger trailing edge vortices (TEVs) grow in close proximity to the primary LEVs, which contributes to an enhanced localized vortex compression and tilting near dislocations. This phenomenon amplifies the undulation amplitude of primary LEVs, eventually leading to vortex tearing. The larger circulation of TEVs with increasing $St_c$ provides an additional explanation for an accelerated vortex compression that coincides with a faster transition of spanwise LEV instability to secondary hairpin-like structures in the wake.

Microscopic irregularity (roughness) of bounding surfaces affects macroscopic dynamics of fluid flows. Its effect on bulk flow is usually quantified empirically by means of a roughness coefficient. A new approach, which treats rough surfaces (e.g. parallel plates) as random fields whose statistical properties can be inferred from measurements, is presented. The mapping of a random flow domain onto its deterministic counterpart, and the subsequent stochastic averaging of the transformed Stokes equations, yield expressions for the effective viscosity and roughness coefficient in terms of the statistical characteristics of the irregular geometry of the boundaries. The analytical nature of the solutions allows one to handle surface roughness characterized by short correlation lengths, a challenging feature for numerical stochastic simulations.

In spray cleaning, a multitude of small drops, violently accelerated by a high-speed gas stream, strike a dirty surface. This process is extremely effective: very little dirt can resist it. This is true even for dirt particles whose characteristic size is less than 100 nm. Spray cleaning is classically modelled by balancing particle adhesion with either inertial stress or viscous shear near the surface, the latter being calculated using droplet size and velocity as the characteristic length and velocity. This results in dimensionless numbers that are often well below one, suggesting that the mechanical stress exerted on the surface by the drop impact that detaches the particle is not well captured. Using quantitative nanoscale measurements, we show that the remarkable efficiency of spray cleaning results from the forced spreading of each droplet on the surface, which generates an unsteady and inhomogeneous shear confined to a boundary layer entrained in the wake of the liquid–solid contact line. In the very first moments of impact, the boundary layer is extremely thin, yielding a gigantic stress: the contact line of the spreading droplets sweeps all the surface particles away. We propose a quantitative model of spray cleaning based on this unsteady surface stress, which agrees well with (i) experimental data obtained with spray droplets of $34\ \mathrm {\mu }$m mean radius impacting the surface to be cleaned at a mean velocity ranging between 30 and 70 m s$^{-1}$ and contamination by nanoparticles of varying nature and shape and (ii) data in the literature on spray cleaning.

Simulations of elastic turbulence, the chaotic flow of highly elastic and inertialess polymer solutions, are plagued by numerical difficulties: the chaotically advected polymer conformation tensor develops extremely large gradients and can lose its positive-definiteness, which triggers numerical instabilities. While efforts to tackle these issues have produced a plethora of specialized techniques – tensor decompositions, artificial diffusion, and shock-capturing advection schemes – we still lack an unambiguous route to accurate and efficient simulations. In this work, we show that even when a simulation is numerically stable, maintaining positive-definiteness and displaying the expected chaotic fluctuations, it can still suffer from errors significant enough to distort the large-scale dynamics and flow structures. We focus on two-dimensional simulations of the Oldroyd-B and FENE-P equations, driven by a large-scale cellular body forcing. We first compare two positivity-preserving decompositions of the conformation tensor: symmetric square root (SSR) and Cholesky with a logarithmic transformation (Cholesky-log). While both simulations yield chaotic flows, only the latter preserves the pattern of the forcing, i.e. its fluctuating vortical cells remain ordered in a lattice. In contrast, the SSR simulation exhibits distorted vortical cells that shrink, expand and reorient constantly. To identify the accurate simulation, we appeal to a hitherto overlooked mathematical bound on the determinant of the conformation tensor, which unequivocally rejects the SSR simulation. Importantly, the accuracy of the Cholesky-log simulation is shown to arise from the logarithmic transformation. We also consider local artificial diffusion, a potential low-cost alternative to high-order advection schemes. Unfortunately, the artificially enhanced diffusive smearing of polymer stress in regions of intense stretching substantially modifies the global dynamics. We then show how the spurious large-scale motions, identified here, contaminate predictions of scalar mixing. Finally, we discuss the effects of spatial resolution, which controls the steepness of gradients in a non-diffusive simulation.

The seminal Bolgiano–Obukhov (BO) theory established the fundamental framework for turbulent mixing and energy transfer in stably stratified fluids. However, the presence of BO scalings remains debatable despite their being observed in stably stratified atmospheric layers and convective turbulence. In this study, we performed precise temperature measurements with 51 high-resolution loggers above the seafloor for 46 h on the continental shelf of the northern South China Sea. The temperature observation exhibits three layers with increasing distance from the seafloor: the bottom mixed layer (BML), the mixing zone and the internal wave zone. A BO-like scaling $\alpha =-1.34\pm 0.10$ is observed in the temperature spectrum when the BML is in a weakly stable stratified ($N\sim 0.0018$ rad s$^{-1}$) and strongly sheared ($Ri\sim 0.0027$) condition, whereas in the unstably stratified convective turbulence of the BML, the scaling $\alpha =-1.76\pm 0.10$ clearly deviated from the BO theory but approached the classical $-$5/3 scaling in isotropic turbulence. This suggests that the convective turbulence is not the promise of BO scaling. In the mixing zone, where internal waves alternately interact with the BML, the scaling follows the Kolmogorov scaling. In the internal wave zone, the scaling $\alpha =-2.12 \pm 0.15$ is observed in the turbulence range and possible mechanisms are provided.

This work reports an experimental study of the turbulent entrainment into the planar wake of a circular cylinder, exposed to various turbulent backgrounds, from the near- to the far-field. The background turbulence features independently varying turbulence intensity and integral length scale, thereby rendering different turbulent/turbulent interfaces (TTIs) between the background and the primary flow (wake). Combined, simultaneous particle image velocimetry and planar laser induced fluorescence measurements were conducted to quantify the entrainment characteristics across these various TTIs at an inlet Reynolds number of 3800. The primary focus was on understanding how turbulent entrainment evolves spatially in conjunction with the rapid development of the large-scale coherent vortices in the planar wake, and how such evolution is affected by the background turbulence. It is found that TTIs can establish two layers when the background turbulence is sufficiently intense, which distinguishes TTIs from the turbulent/non-turbulent interface (TNTI). The two layers are underpinned by different physical mechanisms but have the same thickness and appear to scale with the local Kolmogorov length scale after the wake spreading transition position (Chen & Buxton, J. Fluid Mech., vol. 969, 2023, A4). It is also found that the probability density functions of the entrainment velocity for both TTIs and a TNTI display power law tails, which are associated with extremely large entrainment velocities occurring more frequently than for a Gaussian process. These intermittent, extreme entrainment velocities make a remarkable contribution to the mean entrainment velocity, particularly in the near wake, which leads to a much higher mean entrainment velocity than farther downstream, for both a TNTI and the TTIs. Conditionally averaged analysis reveals that these extreme events of the entrainment velocity are directly associated with intense enstrophy structures close to the interface.

Understanding the solutal convection is a crucial step towards accurate prediction of CO$_2$ sequestration in deep saline aquifers. In this study, pore-scale resolved direct numerical simulations (DNS) are performed to analyse the scaling laws of the solutal convection in porous media. The porous media studied are composed of uniformly distributed square or circular elements. The Rayleigh numbers in the range $426 \le Ra \le 80\,000$, the Darcy numbers in the range $1.7\times 10^{-8} \le Da \le 8.8\times 10^{-6}$ and the Schmidt numbers in the range $250 \le Sc \le 10^4$ are considered in the DNS. The main time, length and velocity scales affecting the solutal convection are classified as the diffusive scales ($t_I$, $l_I$ and $u_I$), the convective scales ($t_{II}$, $l_{II}$ and $u_{II}$) and the shut-down scales ($t_{III}$, $l_{III}$ and $u_{III}$). These scales determine the pore-scale Rayleigh number $Ra_K$ and the Rayleigh number $Ra$. Based on the DNS results, the scaling laws for the transient dissolution flux are proposed in the different regimes of convection. It is found that the onset time of convection ($t_{oc}$) and the period of the flux-growth regime ($\Delta t_{fg}$) vary linearly with the convective time scale $t_{II}$. The merging of the original plumes and the re-initiation of the new plumes start in the same period, resulting in the merging re-initiation regime. Horizontal dispersion plays an important role in plume merging. The dissolution flux $F$ and the time since the onset of convection $t^{\ast }$ have an $F / u_{II} \sim (t^{\ast }/ t_{II})^{-0.2}$ scaling in the later stage of the merging re-initiation regime. This scaling is caused by the merging of the low-wavenumber plumes. It becomes more pronounced with decreasing porosity and leads to the nonlinear relationship between the Sherwood number $Sh$ and $Ra$ when the domain is not high enough for the plumes to fully develop. According to the DNS results, a regime is expected that is independent of both $Ra$ and $Ra_K$, while the dimensionless constant flux $F_{cf}/u_{II}$ in this regime decreases with decreasing porosity. When the mean solute concentration reaches approximately 30 %, convection enters the shut-down regime. For large $Ra$, the dimensionless flux in the shut-down regime follows the scaling law $F/u_{III}\sim (t/t_{III})^{-1.2}$ at large porosity ($\phi =0.56$) and $F/u_{III}\sim (t/t_{III})^{-1.5}$ at small porosity ($\phi =0.36$ or $0.32$). The study shows the significant pore-scale effect on the convection. The DNS cases in this study are mainly for simplified geometries and large $Ra_K$. This can lead to uncertainties of the obtained scaling coefficients. It is important to determine the scaling coefficients for real geological formations to predict a real CO$_2$ sequestration process.

Developing a model to describe the shock-accelerated cylindrical fluid layer with arbitrary Atwood numbers is essential for uncovering the effect of Atwood numbers on the perturbation growth. The recent model (J. Fluid Mech., vol. 969, 2023, p. A6) reveals several contributions to the instability evolution of a shock-accelerated cylindrical fluid layer but its applicability is limited to cases with an absolute value of Atwood numbers close to $1$, due to the employment of the thin-shell correction and interface coupling effect of the fluid layer in vacuum. By employing the linear stability analysis on a cylindrical fluid layer in which two interfaces separate three arbitrary-density fluids, the present work generalizes the thin-shell correction and interface coupling effect, and thus, extends the recent model to cases with arbitrary Atwood numbers. The accuracy of this extended model in describing the instability evolution of the shock-accelerated fluid layer before reshock is confirmed via direct numerical simulations. In the verification simulations, three fluid-layer configurations are considered, where the outer and intermediate fluids remain fixed and the density of the inner fluid is reduced. Moreover, the mechanisms underlying the effect of the Atwood number at the inner interface on the perturbation growth are mainly elucidated by employing the model to analyse each contribution. As the Atwood number decreases, the dominant contribution of the Richtmyer–Meshkov instability is enhanced due to the stronger waves reverberated inside the layer, leading to weakened perturbation growth at initial in-phase interfaces and enhanced perturbation growth at initial anti-phase interfaces.

The resolvent analysis reveals the worst-case disturbances and the most amplified response in a fluid flow that can develop around a stationary base state. The recent work by Padovan et al. (J. Fluid Mech., vol. 900, 2020, A14) extended the classical resolvent analysis to the harmonic resolvent analysis framework by incorporating the time-varying nature of the base flow. The harmonic resolvent analysis can capture the triadic interactions between perturbations at two different frequencies through a base flow at a particular frequency. The singular values of the harmonic resolvent operator act as a gain between the spatiotemporal forcing and the response provided by the singular vectors. In the current study, we formulate the harmonic resolvent analysis framework for compressible flows based on the linearized Navier–Stokes equation (i.e. operator-based formulation). We validate our approach by applying the technique to the low-Mach-number flow past an airfoil. We further illustrate the application of this method to compressible cavity flows at Mach numbers of 0.6 and 0.8 with a length-to-depth ratio of $2$. For the cavity flow at a Mach number of 0.6, the harmonic resolvent analysis reveals that the nonlinear cross-frequency interactions dominate the amplification of perturbations at frequencies that are harmonics of the leading Rossiter mode in the nonlinear flow. The findings demonstrate a physically consistent representation of an energy transfer from slow-evolving modes toward fast-evolving modes in the flow through cross-frequency interactions. For the cavity flow at a Mach number of 0.8, the analysis also sheds light on the nature of cross-frequency interaction in a cavity flow with two coexisting resonances.