When bubbles bounce or stick

Domain: biology
Published: 2026-03-20
Editorial outcome: major_revision

Abstract
Bubble collisions with solid surfaces are critical in both natural and industrial contexts, ranging from gas absorption in bioreactors to marine aerosol production. Despite their ubiquity, the physical criteria underlying bounce dynamics remain unresolved. Here, through experiments and simulations, we map the phase diagram of rising bubbles impacting a wall in Galilei (Ga)–Bond (Bo) space, revealing four dynamic regimes: fully bouncing, underdamped non-bouncing, overdamped non-bouncing, and breakup. We find that bouncing is governed jointly by Ga and Bo, while underdamped dynamics depend solely on Ga. The initial rise distance modulates regime transitions only when shorter than five bubble radii, whereas longer rise distances in high Ga and Bo regions promote bubble breakup and suppress bouncing. We develop a unifying double-mass-spring-damper model, quantitatively matching the complete rebound and damped adhesion regime, and explain bouncing suppression in microgravity and low-viscosity fluids via energy dissipation analysis. Our work provides a unified framework that clarifies the governing mechanisms of bubble-impact dynamics, offering design principles for applications in chemical engineering, biomedicine, and environmental flows.

Introduction
Bubble dynamics in Newtonian fluids represent a fundamental challenge in interfacial science, with critical implications across natural phenomena and industrial applications1–3. Near solid boundaries, bubbles display rich behaviors, including bouncing, coalescence, fragmentation, and adhesion, arising from competing inertial, viscous, capillary, and thin-film forces4,5. Understanding and predicting these interactions is increasingly urgent across applications: thermal management, where gas-liquid interfaces constrain heat transfer in next-generation immersion cooling systems6; energy conversion, where bubble adhesion degrades electrolyzer efficiency7; and medicine, where the stability of ultrasound contrast agents depends on microscale bubble dynamics8. Despite this importance, predictive models for bubble–boundary interactions remain limited, especially in confined geometries relevant to emerging micro- and nanodevices.

The bouncing dynamics of bubbles near walls are largely governed by the conversion of surface energy, stored during interfacial deformation and bubble compression, into fluid kinetic energy. Surface tension, viscosity, wall properties, and surface contamination collectively modulate these interactions, as demonstrated experimentally9–12. As bubbles approach a boundary (smooth horizontal walls13, compound interfaces14,15, or more complex geometries such as inclined surfaces16 and curved boundaries17–19), they decelerate and may undergo multiple bounces before energy dissipation ultimately leads to adhesion to the wall. This evolution involves drainage of the intervening thin liquid film, the establishment of contact with the wall, and the formation of a stable three-phase contact line5. Here, we focus on the canonical impact of a single bubble on a rigid horizontal surface, which isolates the fundamental normal-incidence dynamics needed to derive universal scaling laws. Because impacts on complex geometries can often be decomposed into normal and tangential components, establishing the physics of bubble bouncing in this idealized setting is a prerequisite for predicting bubble behavior in realistic systems20,21.

High fluid viscosity and large interfacial deformations due to low surface tension strongly damp the bubble rise velocity and suppress bouncing22. While the Stokes number has been identified as a governing parameter for bubble and droplet collisions with submerged plates23–25, its a posteriori nature, which depends on the final collision state, limits its predictive utility under varying initial conditions. Recent work13 extended this framework to model oil droplet rebound by coupling a force-balance description with lubrication theory26,27. For gas-liquid systems, the Galilei (Ga) and Bond (Bo) numbers (alternatively, Ohnesorge and Bond number) effectively characterize the bubble rising dynamics and flow regimes, governing the transition from stable axisymmetric rising to asymmetric rising and breakup regimes28. However, a single bubble impacting a rigid horizontal surface poses distinct challenges because its inertia is far smaller than that of a droplet, and the dominant dissipation pathways differ. Existing single-mass-spring models23, though useful for specific zero-deformation contact scenarios, fail to generalize across parameter spaces. These limitations point to two critical gaps: (i) mechanistic models formulated in terms of fundamental dimensionless groups rather than phenomenological descriptors, and (ii) the overlooked role of short rise distances, Lrise. Most studies focus on bubbles reaching terminal velocity9,23,27,29, whereas micro/nanofluidic applications and secondary rebounds often occur under sub-terminal conditions. A systematic investigation of Lrise-dependent dynamics is therefore essential to improve predictive capabilities in confined flows.

Gas bubbles frequently encounter viscosity variations spanning orders of magnitude—from EPS(extracellular polymeric substances)-rich biological waters, where high viscosity strongly affects bubble rise speed, deformation, and residence time (see Fig. 1a), to industrial fluids where gas injection is essential for ecological safety, product quality, and even sustaining biological processes30–32. Motivated by these observations, we employ glycerol solutions to simulate viscous fluid conditions and investigate how fluid properties govern rising bubble–wall interactions (see Fig. 1a). Combining high-speed camera with direct numerical simulations, we construct a comprehensive phase diagram of bubble–wall collisions in Galilei–Bond space. We demonstrate that the transition between bounce and adhesion is controlled by the interplay of Ga and Bo, while overdamped behavior depends solely on Ga. Our key advance is a predictive double-mass-spring-damper model that quantitatively captures the transitions from complete rebound to damped adhesion. These results provide both fundamental principles for interfacial hydrodynamics and practical design rules for applications requiring controlled bubble adhesion (e.g., electrolysis cells) or precise impact management (e.g., biomedical contrast agents, aquaculture).

Results
Dynamics of Bubble-Wall Impacts.
We delineate the bouncing dynamics of bubble–wall impacts by first examining air bubbles of radius of R = 2mm rising in glycerol–water mixtures of varied concentration (see full details in Methods). Each bubble is released and impacts on a rigid horizontal wall surface after rising a distance Lrise = 10R = 20mm, during which the trajectory and shape are shown to be axisymmetric. To better describe the bubble’s shape and position, we define Ht and Hb as the distances from the top and bottom vertices of the bubble to the wall, respectively. The bubble’s width and height are denoted by DH (horizontal diameter) and DV (vertical diameter). These quantities are obtained through image analysis (see supplementary S1). We conduct three cases with different concentrations of =0.1,0.5,0.8 and quantify the bubble dynamics via the temporal evolution of Ht and Hb (see Fig. 1b). In low-concentration solution (=0.1), the bubble’s speed at t =0ms (defined when the bubble bottom is about 4.25mm from the surface) is 0.27m/s with an aspect ratioDH/DV2.73. As the bubble approaches the wall, lubrication pressure in the draining film increases and the bubble is compressed, with a maximum vertical compression and horizontal spread at t =15.6ms (3.49). We denote the current Hb value as Hb,min. The bubble then recoils under surface tension, moves downward, and fully detaches by t = 35.2ms (denoted as full bounce, FB), with the Ht value at this instant as the rebound height Ht,max. Subsequent impact processes will no longer cause a bounce. Increasing the glycerol concentration to moderate ( = 0.5) suppresses the rebound. The initial speed and deformation are reduced (0.18m/s, and 1.34). Subsequently, an underdamped-like response ensues with repeated oscillations, and the apex stays attached to the wall throughout (denoted as underdamped non-bounce, UD-NB). At high concentration =0.85, the initial speed and deformation decrease further (0.07m/s and 1.18), and the motion becomes monotonic: from t = 0ms to 50ms, the bubble rises until the apex nearly touches the wall, after which the bottom point slowly rises and the bubble creeps to equilibrium without any recoil (denoted as overdamped non-bounce, OD-NB).

FB exhibits more than one positive peak for Ht, whereas Ht = 0 for UD-NB and OD-NB after impact. Accordingly, UD-NB exhibits an underdamped response: Hb first falls below the distance from the equilibrium bottom vertices to the wall Hb,eq = 2.3mm, then overshoots and undergoes damped oscillations about Hb,eq until settling, whereas OD-NB decreases monotonically toward Hb,eq. We identify different bubble phases using the maximal rebound height Ht,max and the minimal bottom height Hb,min after the first impact (see Fig. 1c). For bubbles with a radius of 2mm, the transition from FB to UD-NB near c,1 = 0.4750.025, evidenced by Ht,max  0 for FB. And the subsequent transition from UD-NB to OD-NB follows near c,2 = 0.7750.025, marked by Hb,min  Hb,eq for UD-NB. Crucially, the transitions from FB, UD-NB to OD-NB also persist at R = 0.7, 1mm, but the transitions are not uniquely set by concentration. Rather, they arise from the coupled action of gravity, viscosity, and capillarity—a coupling we next interrogate using the Galilei number and Bond number, which quantify gravitational-to-viscous and gravitational-to-capillary force ratios, respectively.

Phase diagram
The experiments, however, access only a limited region of parameter space, spanning 0.1 < Bo < 1 and 1 < Ga < 230 for bubbles of radius R = 0.7-2 mm. To interrogate the gravity–viscosity–capillarity coupling more broadly, we complement the measurements with direct numerical simulations of the two-phase flow with surface tension, sweeping 10-4 < Bo < 500 and 1 < Ga < 1000 to explore the bouncing dynamic of the rising bubble28. This combined approach enables comprehensive analysis across parameter regimes inaccessible to experiments alone. We systematically analyse bubble impact dynamics after a fixed rise distance of Lrise = 10R, focusing on three behavioural regimes introduced above (see Fig. 2a). To determine the validity range of axisymmetric rise dynamics at Lrise = 10R, we validate the axisymmetric–asymmetric rise transition boundary reported in Ref. 28 using three-dimensional simulations at Bo=0.1, Ga=100 and 1000. In these cases, the horizontal deviation of the bubble’s center of mass (rs) remains negligible (rs < 0.03R, see Fig. 2g and supplementary movie 2). For bubbles in low viscosity liquid (Ga > 40), the axisymmetric model remains valid for Bo < 0.1, beyond which asymmetric rising bubbles develop and are simulated by 3D model (upper-right of the dashed line in Ga-Bo space, see Fig. 2a).

The FB regime occurs at high Ga and low Bo. Figure. 2b plots the bubble’s maximum bounce height Ht,max. It is zero at the regime boundary and increases with distance from it, peaking at 8.03 when Bo = 0.001, Ga = 1000. For Bo < 0.01, Ht,max is essentially independent of Bo (see Fig. 2b). Transitioning to lower Ga or higher Bo yields UD-NB region, where bubbles exhibit damped oscillations about equilibrium position Hb,eq (see Fig. 3c). OD-NB happens in low Ga region. In addition, sufficiently large Ga and Bo lead to a breakup regime, in which the cap shape of the bubble formed during rising cannot be maintained upon impact with the wall. Under the action of drainage forces, the center of the cap gradually thins and eventually ruptures, forming a toroidal shape along with several small satellite bubbles (see Fig. 2e, f). Remarkably, the UD-NB to OD-NB transition depends exclusively on Ga (≈5.5). The observed threshold shift toward higher Ga (≈9) at extreme Bo values may reflect that low surface tension effects promoting dimple formation (high Bo), or that computational resolution limits detecting minimal oscillations (low Bo).

In the experimental study, we primarily vary Galilei numbers via changing solution concentration at fixed bubble radius. Here we examine how the Bond number controls the transitions at fixed Galilei number (Ga = 100). Figure 2c–e presents three representative cases with Bo = 0.1, 1.6, and 10 (see also the supplementary movie 2). At Bo = 0.1 (Fig. 2c), the bubble exhibits a full bounce after impact. The smaller Bo (i.e., stronger capillary dominance) limits deformation and thus promotes bouncing compared with Fig. 1b. Increasing to Bo = 1.6 suppresses bouncing (Fig. 2d): the bubble deforms more strongly, develops a pronounced dimple, and undergoes an underdamped non-bounce (UD-NB). A further increase to Bo = 10 leads to breakup upon impact. In this case, the bubble motion remains nearly axisymmetric: the film ruptures only at the center, while the rim retains a coherent toroidal shape (Fig. 2e and supplementary movie 2). Conversely, when significant deviation from axisymmetry occurs (rs > R, see Fig. 2g), as in the case Ga=250, Bo=2.5, the toroidal structure ruptures on one side, yielding a C-shaped ring (see Fig. 2f and supplementary movie 2). Importantly, we focus on cases in which the bubble remains intact during ascent and ruptures only upon impact with the substrate, in contrast to the peripheral and central breakup modes observed during free rise by Tripathi et al.28. Our observation suggests that the solid wall promotes rupture by destabilizing the cap during impact. Finally, these results show that increasing Bo alone at fixed Ga is insufficient to trigger the overdamped non-bounce (OD-NB) regime; achieving OD-NB requires stronger viscous damping rather than merely weaker capillarity.

A Double-Mass-Spring-Damper Model for Predicting Bubble Bouncing Dynamics
To elucidate bubble bouncing mechanisms, we develop a double-mass-spring-damper model, extending prior droplet models33 by incorporating viscous drag and gravitational effects absent in single-mass approximations24. The system comprises a spring (initial length d = 2R, stiffness k) representing bubble deformation, two effective masses (m1, m2) at the bubble's lower (x1) and upper (x2) interfaces (Fig. 3). This framework physically captures bubble-specific dynamics. Unlike droplets with concentrated mass, bubbles exhibit negligible internal mass, with motion dominated by fluid added mass (m* =(4/3)R3lCm) and viscous damping. The bubble-wall interaction is a two-phase interaction: the approach phase and the bouncing phase. During the approach phase, rising bubble deformation converts potential energy into surface/kinetic energy (spring-like compression). In the subsequent bouncing phase, stored surface energy drives fluid kinetic energy (spring rebound analogy). Key innovations include explicit viscous dissipation terms for thin-film hydrodynamics, gravity-dependent terms for buoyancy effects, and dual-mass formulation capturing axisymmetric interface motions.

We focus on the bouncing phase, during which the top point of the bubble remains in contact with the wall. This implies that both the position and velocity of the upper mass m2 are zero (x2 = 0, dx2/dt = 0). Consequently, the viscous drag acting on the entire bubble is solely attributed to the lower mass m1. The viscous force is expressed as F = Cv = 6Rdx1/dt, where  = 2/3 + R/(4h) reflects the combined contribution of bulk drag and the classical lubrication limit for a mobile bubble27. The elastic force generated by the gas–liquid interface is modelled as Fk = k(x1 − x2 − d), where k = (8/5) is the effective stiffness of the bubble, derived from the surface energy stored during small deformations24. The buoyancy force Fg = (4/3)R3lg is assumed to act equally on both masses, while the added mass, representing the inertia of the displaced fluid, is asymmetrically distributed between two masses (m1 + m2 = m*, m1 = m2 = m*).

Considering these forces, the Newtonian equations of motion for the two-mass system are governed by the combined effects of buoyancy, viscous drag, elastic (spring) force, and reaction forces at the wall, which can be expressed as equations (1) and (2). Equation (1) can be modified as equation (3), whose solution (x1) depends on the discriminant (C/m1)2 − 4k/m1 and initial conditions. When d2x2/dt2  0, Fwall  0, m2 moves away from the wall, which gives a sufficient condition for the bubble to bounce as equation (4).

m1 d2x1/dt2 = −Fg + k(x1 − d) − C dx1/dt (1)
m2 d2x2/dt2 = −Fg + k(x1 − d) + Fwall (2)
where Fwall is the reaction force the wall puts on m2 when it contacts the wall.
d2x1/dt2 + (C/m1) dx1/dt + (k/m1)x1 = −(Fg + kd)/m1 (3)
d2x2/dt2 = −Fg + k(x1 − d)  0 (4)

To study the bouncing process of the bubble with a compression of 2 (see Fig. 3b), we determine the initial conditions: x1(t = 0) = d − 2, dx1/dt(t = 0) = 0. This compression arises from the combined effects of inertial and buoyant forces, and Fg  k(2) holds. When (C/m1)2 − 4k/m1  0, x1 = x*1 and dx*1/dt  0, the bubble is balanced or slowly approach but not exceed the equilibrium position (x*1(t = ) = d − Fg/k), which provides the criterion from OD-NB to UD-NB (equation (5)). In contrast, when (C/m1)2 − 4k/m1  0, x1 = x**1 and x**1 has a maximum value of max(x**1) = x**1(t = /)  d − Fg/2k, meaning that bubbles oscillate. When the bubble bounce off, Fwall(t) = (2k − Fg)e/ − Fg  0 satisfies, which gives the criterion for bouncing bubbles (equation (6)). The details of the derivations are given in supplementary S6.

Ga = Gac (5)
Ga = Gac exp(ln(Bo1/2)/) (6)
where Gac = 5.5  0.5 and  = 3.0  0.1 are estimated as an approximate range and fitted to the simulation and experimental data in Fig. 2a. Two criteria for fully bouncing and underdamped non-bouncing exhibit good agreement with simulation and experimental data. When Ga → , Bo  4. However, in the region of Ga > 40, Bo ≈ 1, we observe that large rising distances (Lrise > 10R) induce non-axisymmetric ascent and impact trajectories. These instabilities suppress bubble bouncing, thereby establishing Bo = 1 as the upper boundary of the FB regime (see Fig. 2a).

Energy analysis
Beyond the equations of motion, we quantify bubble bouncing behavior through complete energy budget analysis. We examine the energy evolution during the entire impact process, from initial wall contact (t0) to detachment (t1), with the system's energy balance given by:
E0 = Ek(t) + Eg(t) + E(t) + E(t) (7)
where Ek(t) and E(t) are the kinetic and viscous dissipation energies for the whole fluid, with E(t0) = 0. Eg(t) and E(t) are the change of gravitational potential energy and surface energies, with Eg(Vcenter = 0) and E(t=0)(spherical gas-liquid surface) as the energy base to ensure that change of energies is always positive for the impact event. Each energy component is normalized by the initial total energy E0 at t = t0.

We first examine how liquid viscosity, characterized by the Galilei number (Ga = 1 ~ 100), suppresses bubble bouncing in a microgravity environment, where the Bond number (Bo = 0.001) ensures capillarity dominates over gravity. The initial time is defined as the moment when the bubble center is at 1.1R from the wall (t0 = t(zcenter = 1.1R)). Figure. 4a, b depict the energy balance evolution for Ga = 100 at Lrise = R and 10R, respectively. Initially, the system energy E0 consists of kinetic energy (0.84E0 for Lrise = R; 0.97E0 for Lrise = 10R) and gravitational potential energy (0.16E0 for Lrise = R; 0.03E0 for Lrise = 10R). As the bubble deforms, this energy transfers into surface energy until the bubble center’s velocity (Vcenter) reaches zero, marking maximal deformation. At this point, most energy is stored as surface energy (0.78E0 for Lrise = R; 0.84E0 for Lrise = 10R), while a small fraction remains as kinetic energy (0.033E0 and 0.045E0, respectively), indicating residual fluid motion. During retraction, surface energy reconverts to kinetic energy, recovering 0.63E0 (Lrise = R) or 0.76E0 (Lrise = 10R) at take-off. Most of the remaining energy is dissipated and little stored in the deformed bubble’s surface. Snapshots indicate that dissipation primarily occurs near the bubble and within the thin liquid film between the bubble and the wall. When Ga is reduced to 3 (see Fig. 4c, d), the bubble exhibits no bouncing for either Lrise = R or 10R. Here, viscous dissipation increases significantly, while the gravitational potential as well as kinetic energies decrease monotonically and gradually to zero, with minimal change in surface energy. Specifically, more than 85% of the energy for Lrise = 10R (97% for Lrise = R) is dissipated as the bubble moves from 1.1R to R. Subsequently, kinetic energy decays without significant conversion to surface or gravitational potential energy, indicating a complete absence of UD-NB phase.

The energy budgets at take-off in the bubble bounce regime for Bo = 0.001 are summarized in Figs. 4e, f for Lrise = R and 10R, respectively. The influence of Ga on energy conversion is primarily reflected in viscous dissipation and kinetic energy. A discrepancy arises at low Ga between Lrise = R and 10R, as the delayed contact initiation excludes a significant portion of the energy dissipated during bubble impact. By adjusting the contact start time to an earlier point t0 = t(zcenter = 1.5R), this discrepancy is resolved (details can be found in supplementary S7), and the energy budget for Lrise = 10R (see Fig. 4f) aligns with that of Lrise = R. A comparison of the take-off energy budgets for Lrise = R and Lrise = 10R shows nearly identical energy distributions, regardless of the initial rising distance. However, as Ga decreases, a larger fraction of kinetic energy is lost to viscous dissipation, while surface and gravitational potential energies remain nearly constant. This shift in energy conversion is the key mechanism suppressing bubble bouncing.

We then examine how surface tension and gravity (quantified by the Bond number, Bo = 0.001–1) modulate bubble bouncing in low-viscosity fluids at fixed Galilei number (Ga = 40). Energy analysis commences when the bubble’s center is positioned at R from the wall (t0 = t(zcenter = R)). For (Ga, Bo, Lrise) = (40, 0.01, R) and (40, 0.01, 10R), viscous dissipation remains substantial despite high Ga, owing to pronounced bubble deformation (Fig. 5a, b). Increasing Bo to 1, beyond Boc, bouncing is fully suppressed (Fig. 5c, d). Instead, the bubble undergoes damped oscillations characterized by cyclic energy conversion: kinetic and gravitational potential energy transform into surface energy during compression, peaking at maximum deformation (Vcenter = 0) with residual fluid motion. The system asymptotically approaches hydrostatic equilibrium through successive energy transfers.

Energy budgets at take-off (Figs. 5e, f) reveal that increasing Bo enhances viscous dissipation and gravitational potential energy, while reducing kinetic energy. Notably, unlike the high-dissipation regime (Fig. 4), substantial energy persists as gravitational potential energy during bouncing attempts. This demonstrates that buoyancy (Bo) and viscous stress (Ga) act synergistically to inhibit bouncing, even when inertial forces dominate.

Effect of rising distance on bubble dynamics
Figure. 6a maps the regime transitions in Ga–Bo space at six characteristic rising distances (Lrise = 0.1R to 20R), revealing two key findings. Firstly, when the bubble starts far from the wall (Lrise > 5R), the two transition boundaries converge, demonstrating negligible influence of rise distance beyond 5R. Secondly, when Lrise is shorter than R, the transition boundaries shift toward higher Ga values. For the cases at Lrise > 5R, the parameters Gac = 5.5  0.5 and  = 3.0  0.1 in equations (5) and (6) are obtained by fitting simulation data, while Gac = 8.0  0.5 and  = 1.0  0.1 for the case at Lrise = 0.1R. A notable phenomenon is observed in the region of Ga > 40, Bo≈1 (see Fig. 6a, b and supplementary movie 5), where large rising distance induces complex non-axisymmetric ascent trajectories and alter the impact dynamics. Consequently, cases that fall within the FB regime at short rising distances (Lrise = R) are observed to transition into the UD-NB regime as the rising distance increases to 20R (at Ga=250, Bo=1). Notably, a mixed regime emerges at Bo < 10−3 and Ga < 10, where the rebound height becomes negligible (~0.001R) and bounce and overdamp non-bounce criteria coalesce. UD-NB manifests as micro-bounces with equilibrium position near 2R (Fig. 3c). We plot the maximum velocity (Vmax) distribution in Ga-Bo space for Lrise = 60R, 5R in supplementary Fig. S7, demonstrating a nonlinear relationship between maximum velocity and the three governing parameters (Ga, Bo, and Lrise), which can be effectively predicted using machine learning approaches34,35. The regime boundaries align with the maximum velocity contours, establishing flow dynamics as the primary determinant of regime shifts24. Specifically, the transition occurs when Vmax reaches ~85% of the terminal velocity, defined as the maximum velocity reached at Lrise = 60R (see supplementary S6).

We examine two representative cases (Ga = 20, Bo = 0.1 and Ga = 6, Bo = 0.1) to quantify how rising distance (Lrise = 0.1R–20R) governs regime transitions. Figure. 6c presents the temporal evolution of the center velocity and top point coordinates, with key snapshots for Lrise = 0.5R and 5R shown in Fig. 6e. The bubble reaches the terminal velocity at Lrise = 20R, while achieving approximately 85% of terminal velocity at Lrise = 5R. Notably, for Lrise = 5R, 10R and 20R, the velocity profiles during approach and bouncing stages remain nearly identical, suggesting that once the impact velocity reaches a critical threshold, bubbles exhibit consistent bouncing behavior regardless of rise height. For the case of Ga = 6, Bo = 0.1 (see Fig. 6d), Lrise not only affects rebound height but also the approaching stage. Insufficient rising distance prevents complete film drainage, causing the bubble to fail to cross its equilibrium position (OD-NB), where the oscillation of the bottom points occurs but remains confined over the equilibrium position. Further reducing Ga to 4 eliminates oscillations entirely, producing a strictly overdamp model without any oscillatory motion.

Methods
Impact experiments
The apparatus (supplementary Fig. S1a) enables controlled bubble generation and release from a position 10R below a quartz plate, with the dynamics recorded at 5000 fps (M220M Revealer high-speed camera). Bubbles (R = 0.7–2.0 mm) are generated through needles of six gauges (34G–10G; inner diameters range 0.06–2.7 mm) in glycerol solutions (0–80%), with gas injection precisely regulated by a syringe pump. The wettability of the quartz plate has been shown not to affect the first impact of bubbles (See the supplementary movie 6), because a liquid film prevents the formation of a three-phase contact line between the bubble and the plate36. However, superhydrophobic surfaces facilitate the removal of the liquid film, so a three-phase contact line forms on superhydrophobic quartz plates during the subsequent impact, which contrasts sharply with the multiple bounces of bubbles on superhydrophilic surfaces. Here, we use superhydrophilic quartz plates for all cases.

Image analysis
We quantify bubble morphology using ImageJ-processed binary masks (see supplementary Fig. S1a), extracting horizontal (DH) and vertical (DV) diameters to calculate the equivalent radius (R = 3DH2DV / 2), while tracking vertical positions of the top (Ht) and bottom (Hb) interfaces with sub-pixel resolution throughout impact events.

The aqueous solutions of glycerol
To vary the viscosity, we use standard aqueous glycerol solutions at different concentrations. The relevant fluid properties are determined from the glycerol concentration and temperature. Experiments are conducted at 29 ± 1°C due to the heating from the illumination. For each solution, the liquid density, viscosity, and surface tension are characterized (see supplementary Fig. S1c), allowing Ga and Bo to be varied systematically. The corresponding Ga and Bo for bubbles with radii of 0.7, 1.0, and 2.0 mm are shown in supplementary Fig. S1d.

Direct numerical simulations
We employ an arbitrary Lagrangian–Eulerian (ALE) formulation to perform direct numerical simulations of axisymmetric bubble-rise dynamics and impact. For non-axisymmetric cases, we use a three-dimensional volume of fluid (VoF) method in Basilisk following Tripathi et al.28

All physical quantities are non-dimensionalized using R, V0 = gR, l, l, R/g = as the characteristic length, velocity, viscosity, density, and time scales, respectively. After fixing the gas-liquid density ratio (a/l = 0.001) and the viscosity ratio (a/l = 0.01), the two remaining dimensionless control parameters are the Galilei number Ga=l gRR/l and the Bond number Bo=lgR2/, which fully specify the gas–liquid system28,35. Accordingly, dimensional quantities (e.g., the gravitational force g = (0, 0, −1), density l = 1, a = 0.001 and viscosity l = 1/Ga, a = 0.01/Ga) of the gas and liquid can be determined. In addition, the initial bubble position, defined as the distance from the bubble apex to the wall (Lrise), must be specified to obtain a unique solution.

  u = 0 (8)
(u/t + u   u) = −p + 2u +   ((u)T) + g (9)

The ALE approach is chosen for the axisymmetric cases because it provides high interface resolution at relatively low computational cost, enabling faithful reconstruction of the bubble evolution. The incompressible Navier–Stokes equations (equations (8) and (9)) are solved within the ALE framework to obtain the velocity (ua, ul) and pressure (pa, pl) fields in both the gas (bubble) and the surrounding liquid domains. The gas-liquid interface S partitions the system into two distinct domains: Ω₁ (air bubble) and Ω₂ (surrounding liquid). Within the ALE framework, the interface is treated in a Lagrangian manner and advected by the local flow. Specifically, the interface dynamics are governed by the normal mesh velocity umesh, which satisfies a mobile boundary condition (ua = ul) in liquid systems without lubricant27.

A balance of velocity, pressure and stress needs to be ensured at the boundary conditions. To track the interface position, the kinematic boundary condition is applied. Under the assumption of no mass transfer across S, the velocity is continuous across the interface37 (see equation (10)). The difference in pressure between the two sides of the interface is caused by surface tension and satisfies equation (11).
umesh  ni = ua  ni = ul  ni (10)
(pa − pl)ni + (l(ul + (ul)T) − a(ua + (ua)T))  ni = ni (11)
where, ni, the unit normal vector of i-th point on S, is positive from the gas towards liquid.   ni =  is the mean curvature, and the determination of  requires that no topological changes occur at the interface 37. The dimensionless expression for surface tension is =1/Bo, without considering the Marangoni effect (i.e., the gradient of tangential surface tension at the interface). The fluid motion changes the position of the interface, and then the Yeoh smoothing model (Stiffening factor set to 10) is used to determine the positions of the grid in the gas and fluid regions38.

The axisymmetric bubble rising dynamics are simulated using the arbitrary Lagrangian-Eulerian (ALE) framework implemented in the commercial finite element analysis software COMSOL Multiphysics®. The laminar flow module is employed to solve the velocity and pressure fields in the fluid domain, while the moving mesh technique dynamically tracks the interface position. To maintain mesh quality, automatic remeshing is implemented to prevent mesh distortion and avoid negative volume elements. Additional details regarding numerical simulation setup, postprocessing, mesh independence and validation are conducted in supplementary S2.

The VoF approach is employed to handle cases exhibiting significant latera displacement and complex topological changes. The Navier-Stokes equations (equations (8) and (9)) are solved using the open-source code Basilisk39. The solver utilizes a generic time loop with a CFL-limited time step. The velocity advection term is computed using the Bell–Colella–Glaz second-order unsplit upwind scheme, while the piecewise linear interface calculation (PLIC) method is employed for geometric interface reconstruction40,41. The interface is tracked by VoF marker function f, which evolves according to the advection equation:
f/t + u   f = 0 (12)
where f = 1 denotes the gas phase (bubble) and f = 0 for liquid phase. Under the one-fluid approximation, the local density  and viscosity  are defined as:
 = fa + (1 − f)l (13)
 = fa + (1 − f)l (14)
To ensure the temporal stability of interface propagation, a constant Courant–Friedrichs–Lewy (CFL) number of 0.5 is strictly maintained41. Further details regarding the computational setup and grid resolution are provided in supplementary S3.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The paper investigates, using multiphase DNS and experiments, bubble-wall impacts and subsequent bouncing dynamics of the bubble. The study elucidates how fluid properties influence the dynamics and also provide a mechanistic mass-spring-damper model to predict the dynamics.\nThe paper is fairly well written and the results are analyzed and explained thoroughly. The subject of the paper is interesting and worthwhile to investigate. I have, however, concerns that I would like the authors to carefully consider before I can recommend publication in Nature Communications. The main comments are:\n1. Although the findings may be useful for understanding and predicting bubble bouncing dynamics, it is not clear how the present results can be used to predict bubble bouncing dynamics in a general (and perhaps more realistic) bubbly flow system with more complex geometries (for example if a bubble impacts a wall that is not horizontal but at an angle) or at larger Ga and Bo numbers. Such limitations should be clearly stated. In relation to this, it should be discussed how the studied problem setup (a bubble rising a short distance and impacting a horizontal wall) is relevant for the applications mentioned in the introduction. This would make it easier for the readers to understand when and where the results are applicable.\n2. Related to lines 62-65 and Fig. 1: How and why is phytoplankton EPS relevant to the present study? There are many (organic and non-organic) water solutions with viscosities more than hundred times that of pure water. That the liquid viscosity influence the bubble dynamics is not surprising and seems like an odd motivation for the present work.\nThe rest of the study in mainly presented in a general manner, using non-dimensional numbers to represent the fluid properties. It is therefore somewhat confusing why such a specific liquid is mentioned both in the introduction and in Fig. 1.\n3. According to the Ga-Bo phase-space reported by Tripathi et al. (2015), the bubble dynamics are axisymmetrical in the region of about Bo < 10-20 and Ga < 30-40. However, on line 132-133 the authors claim that the axisymmertic model is valid for Bo < 1 and Ga > 40. Is this claim based on the experiments or substantiated by other means? This should be carefully explained.\nWithin the axisymmetrical region of the phase-space, it is of course reasonable to use an axisymmetric assumption when performing the DNS. On the other hand, outside of these ranges asymmetric bubble shapes and oscillatory motions may significantly influence the bouncing dynamics. For these reasons, it is difficult to understand the physical relevance of the DNS performed up to Ga=100 for any Bo. In lines 133-134, the authors mention that the results are indeed unphysical when Bo > 1 (although I believe this is more likely the case for Ga>40) but that they use the same axisymmetric framework to preserve regime boundary continuity. This is not a satisfactory argument for including unphysical results in the subsequent analysis and unfortunately makes the study less trustworthy. The discrepancies introduced by the axisymmetric assumption on the results should be carefully checked and discussed and, if needed, such results should be removed.\n4. I like that the effect of Lrise is thoroughly analysed. However, in relation to Fig.6, it is not entirely clear what “condition1” and “condition2” refer to in Fig. 6a. I believe this is connected to Eq. 5 and 6 although this does not seem to be explicitly stated. It is also not clear what the markers in Fig. 6a represent. Are these simulation cases at the transition border or something else?\nIn relation to Fig. 6c, I believe the influence of Lrise on Vcenter should also be reported for the highest Ga in your study that likely need the longest rising distance to reach terminal velocity.\nMinor comments:\nLines 40-42: According to the phase-space reported by Tripathi et al. 2015, the Ga and Bo governs the bubble dynamics at all relevant ranges, not just at low Ga and Bo. Also, what is meant here by “symmetric acceleration”?\nLines 67: What is meant by “monotonic behavior” and why does that depend only on Ga?\nLines 88-90: Seems to be a duplicate sentence here.\nLines 118-119: The symbol in the Bo-interval is not displayed correctly. Also what is the Ga-range of the experiments?\nLines 325-329: I believe also the Lrise parameter is needed to determine the simulation system.\nLines 337-338 and Eq. 10: Is the “normal mesh velocity” equal to the interface velocity or what is this physically representing? Also, the velocity is normally continuous across the interface as long as there is no mass transfer across the interface. A lubricant fluid in between to immiscible fluids would also satisfy a no-slip condition.\nRelated to section S3: The refinement level used in COMSOL is perhaps not the most relevant quantity to report grid resolutions with. I would recommend stating all resolutions in the number of grid points used per bubble diameter (the initial or spherical equivalent bubble diameter). This would apply also for Fig. S3 where the legends are not easy to interpret.\nReferences:\nTripathi, M. K., Sahu, K. C., & Govindarajan, R. (2015). Dynamics of an initially spherical bubble rising in quiescent liquid. Nature communications, 6(1), 6268.", "reviewer_2": "The paper investigates, using multiphase DNS and experiments, bubble-wall impacts and subsequent bouncing dynamics of the bubble. The study elucidates how fluid properties influence the dynamics and also provide a mechanistic mass-spring-damper model to predict the dynamics. The paper is fairly well written and the results are analyzed and explained thoroughly. The subject of the paper is interesting and worthwhile to investigate. I have, however, concerns that I would like the authors to carefully consider before I can recommend publication in Nature Communications. The main comments are:", "reviewer_3": "The paper provides interesting information concerning the rebound of bubbles from a horizontal rigid surface. The authors conducted experimental studies in the range of Bo (0.05–1) and 1 < Ga < 100. The results of simulations performed using the Arbitrary Lagrangian–Eulerian (ALE) method in COMSOL Multiphysics® were obtained for 10⁻⁴ < Bo < 20 and 1 < Ga < 100. Based on the numerical simulation results, an analysis of different forms of energy is also presented. Moreover, an interesting model of two masses connected by a linear spring and a linear damper is introduced.\nIn the introduction and in the description of the publication (line 121), the authors refer to asymmetric bubble rising. Asymmetric bubble rising occurs when bubbles deviate from a straight vertical path, often due to flow disturbances, leading to zigzag, spiral, or oscillating bubble motions and changes in shape. In such cases, bubble rotation may also occur.\nHowever, in all examples presented in the paper, the bubble remains axisymmetric. This fact contradicts the statement regarding asymmetric bubble rising.\nThe experimental and numerical parameter ranges differ significantly, and the simulation results do not show velocity fields near the moving bubble, which makes it difficult to assess whether the model correctly describes asymmetric bubble rising. I suspect that, during the impact of an asymmetric bubble with a rigid plate, both the spring and the damper in the two-mass model should not be linear.\nIn my opinion, this inconsistency between the experimental and numerical ranges makes it difficult to accept the general conclusion stated in the paper: “Our work resolves decades of ambiguity in bubble-impact dynamics.”\nI believe that the paper should be revised in such a way that experimental and numerical results are directly compared within the same ranges of the Bond and Galileo numbers. The remaining results may be presented as a hypothesis, but with a justification of the differences in bubble dynamics which, upon contacting a flat surface, are not axisymmetric.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": ""}, {"round": 2, "reviewer_comments": {"reviewer_1": "The reviewer thanks the authors for their detailed responses and revisions addressing my previous comments, which have helped to improve the clarity and quality of the manuscript. Still, in response to the rebuttal letter and the revised manuscript I have the following remarks:\nComment 3. The manuscript provides insufficient details regarding the new three-dimensional DNS methodology. Key information is missing, including the numerical code used, the specific multiphase DNS approach (e.g., VOF or alternative interface-capturing methods), and the spatial and temporal resolution.\nIt is also unclear which cases shown in Fig. 2 correspond to two-dimensional versus three-dimensional simulations. Are all cases above the dashed line obtained from fully resolved 3D DNS? This should be explicitly stated and clearly defined in the figure caption and/or main text. In particular, are all reported cases indeed performed using fully resolved 3D DNS?\nAt high Galileo numbers, accurately resolving the thin boundary layers at the bubble interface requires very fine grid resolutions. Bubble dynamics are known to be highly sensitive to under-resolved interfacial boundary layers, as the vorticity generated in these regions strongly governs the bubble trajectory. Therefore, a grid-independence (or grid-convergence) study for the highest-Ga cases is necessary to demonstrate that the adopted resolution is sufficient across the entire phase space considered.\nComment 3c. The authors note that for Ga > 40 and Bo ≈ 1, large rising distances lead to complex, non-axisymmetric ascent trajectories that modify the impact dynamics. I would expect similar behavior to occur for Bo > 1 as well. Is this not observed in your simulations? If not, please clarify the underlying reasons, or discuss whether this regime was examined.", "reviewer_2": "In the revised version, the authors adequately addressed my concerns. Therefore, I am now favorable to accept the paper for publication.", "reviewer_3": "In my opinion, the authors have provided comprehensive responses to both my and co-reviewer comments. The authors have adjusted the ranges of the presented experimental data to the simulation results (Fig. 2). They have also discussed the limitations of applicability of the model proposed in the paper, which describes the motion of two masses connected by a spring and a damper.\nAlthough the paper focuses on bubble collisions with solid surfaces, the simulation results also illustrate the process of bubble flow in a liquid, taking into account instabilities in their motion.\nI believe that, in its current form, the paper is suitable for publication in the Nature Communications.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "Response: We sincerely appreciate the reviewer’s encouraging evaluation and constructive comments. We have carefully considered each suggestion and revised the manuscript accordingly. Detailed responses are provided below.\nComment 3. The manuscript provides insufficient details regarding the new three-dimensional DNS methodology. Key information is missing, including the numerical code used, the specific multiphase DNS approach (e.g., VOF or alternative interface-capturing methods), and the spatial and temporal resolution.\nIt is also unclear which cases shown in Fig. 2 correspond to two-dimensional versus three-dimensional simulations. Are all cases above the dashed line obtained from fully resolved 3D DNS? This should be explicitly stated and clearly defined in the figure caption and/or main text. In particular, are all reported cases indeed performed using fully resolved 3D DNS?\nAt high Galileo numbers, accurately resolving the thin boundary layers at the bubble interface requires very fine grid resolutions. Bubble dynamics are known to be highly sensitive to under-resolved interfacial boundary layers, as the vorticity generated in these regions strongly governs the bubble trajectory. Therefore, a grid-independence (or grid-convergence) study for the highest-Ga cases is necessary to demonstrate that the adopted resolution is sufficient across the entire phase space considered.\nResponse: We thank the reviewer for highlighting the need for greater technical detail regarding our numerical implementation and apologize for the insufficient information in the original manuscript. We have substantially expanded the methodology section to include details of the numerical code, interface-capturing method, boundary conditions, and spatial resolution. In addition, we have clarified the distinction between the 2D and 3D simulations in Figure 2 and included a grid-independence study in the Supplementary Material.\n(a) On the 3D-DNS Methodology:\nThe three-dimensional direct numerical simulation (3D-DNS) framework employed in this study is implemented using the open-source solver Basilisk, following the numerical approach of Tripathi et al.¹. In response to the reviewer’s request, we have added a detailed description of the numerical implementation in the Methods section of the revised manuscript (Revised manuscript, from Page 12 Lines 397), as presented below:\n“The VoF approach is employed to handle cases exhibiting significant latera displacement and complex topological changes. The Navier-Stokes equations (equations Error! Reference source not found. and Error! Reference source not found.) are solved using the open-source code Basilisk2. The solver utilizes a generic time loop with a CFL-limited timestep and employs the Bell–Colella–Glaz advection scheme for interface reconstruction3,4. The interface is tracked by VoF marker function f , which evolves according to the advection equation:\n  f t + u   f = 0 (1)\nwhere f = 1 denotes the gas phase (bubble) and f = 0 for liquid phase. Under the one-fluid approximation, the local density and viscosity  are defined as:\nf a (1 f ) l    = + − (2)\n= f +(1− f) (3)\na l\nTo ensure the temporal stability of interface propagation, a constant Courant–Friedrichs–Lewy (CFL) number of 0.5 is strictly maintained4. Further details regarding the computational setup and grid resolution are provided in Supplementary S3.”\nand (Revised Supplementary Materials, from Page 8 Lines 155):\n“S3. 3D VoF simulations using Basilisk\nThe three-dimensional (3D) simulations are performed using the volume of fluid (VoF) method within the open-source solver Basilisk2. The computational setups, including domain size, initial bubble position and the boundary conditions for velocity and pressure, remain consistent with the two-dimensional model. We employ an adaptive mesh refinement (AMR) scheme3, which dynamically refines the regions exhibiting significant variations in interfacial curvature and velocity gradients, with the minimum grid size set to  m in = R 3 2 . Grid adaptation is based on minimizing the error estimated in the volume of fluid tracers ( f ), velocity field with tolerances of 10−3 and 10−2, respectively, using the wavelet algorithm3. To model an ideal wetting substrate, a volume fraction boundary condition of f = 0 is imposed at the wall, at which finer mesh (  m in 4 ) is employed. This constraint maintains a thin lubrication film with a thickness of  8 between the bubble interface and the solid surface to prevent min numerical contact line1. A constant Courant–Friedrichs–Lewy (CFL) number of 0.5 is strictly enforced to ensure the temporal stability of the interface reconstruction4.”\n(b) Clarification on Figure 2 (2D vs. 3D):\nWe thank the reviewer for this clarification. The dashed line in Figure 2 indeed denotes the onset of path instability (axisymmetric versus non-axisymmetric regimes). We confirm that all cases above this line are obtained from fully resolved 3D DNS, and Figures 2c–f are derived from 3D simulations. The Figure 2 caption has been revised to explicitly clarify this point (Revised manuscript, from Page 15 Lines 453):\n“Figure 2 | Phase diagram. (a) Regime map of impacting behavior in both experiments (empty symbols) and simulations (filled symbols) at L = 10R, which rise includes the fully bouncing (FB, blue), underdamped non-bouncing (UD-NB, red), overdamped non-bouncing (OD-NB, black) and breakup domains (magenta pink). Black stars represent the case shown in (c-e). The black dashed line indicates the boundary for an axisymmetric rising bubble. All cases above the dashed line are obtained from 3D simulation. (b) Normalized maximum rebound height H *t,m a x = H t,m a x R , where the inset shows its variation with Ga at Bo = 10-4, 10-3, and 10-2. (c-e) Representative 3D snapshots of (c) the FB case (Ga=100, Bo=0.1), (d) the UD-NB case (Ga = 100, Bo = 1.6), and (e) the breakup case (Ga=100, Bo=10) (see also the Supplementary movie SM2). (f) 3D snapshots of the breakup case at Ga=250, Bo=2.5, where the bubble exhibits C-shaped ring shape. (g) The radial distance ( r s R ) of the center of the bubble measured in the horizontal plane from the original location for the cases of (c-f).”\n(c) Grid Independence Study:\nWe agree that resolving the boundary layer at high Ga numbers is critical for capturing accurate vorticity dynamics and trajectories. Prior to our parametric analysis, we performed a rigorous grid-independence study at Ga=1000, Bo=0.1,1, and Lrise=10R. We tested grid resolutions of  m in = R 1 6 ， R 3 2 ， R 6 4 . Our results confirmed that  m in = R 3 2 is sufficient to achieve grid-independent results (the temporal evolution of the bubble's vertical velocity V ). The mesh size is consistent with the size used in z Tripathi et al.1. We have now included these validation results in the Supplementary Material to demonstrate the reliability of our adopted resolution across the phase space in (Revised Supplementary Materials, from Page 8 Lines 168):\n“We conduct a grid independence study using two representative parameter sets from the non-axisymmetric regime in Fig. S5. We compare the temporal evolution of the bubble’s vertical velocity (V ) using three minimum grid sizes of z  m in = R 1 6 ， R 3 2 ， R 6 4 with approximately 500,000, 3,500,000, and 15,000,000 cells. As shown in Fig. S5, the profiles are in excellent agreement, despite minor deviations during the rebound stage. Consequently, the  = R 32 case is deemed min sufficient to capture physics and adopted for the 3D results to optimize computational resources.”\nFigure S1 | Mesh independence for 3D rising bubbles. (a) (Ga, Bo, Lrise) = (1000, 1,10R), (b) (Ga, Bo, Lrise) = (1000, 0.1,10R). In all panels, the vertical velocity (V ) of z bubble center and time are normalized by V = gR and = R/ g , separately.\nComment 3c. The authors note that for Ga > 40 and Bo ≈ 1, large rising distances lead to complex, non-axisymmetric ascent trajectories that modify the impact dynamics. I would expect similar behavior to occur for Bo > 1 as well. Is this not observed in your simulations? If not, please clarify the underlying reasons, or discuss whether this regime was examined.\nResponse: We appreciate the reviewer’s insightful comment. You are correct that this phenomenon is prevalent for Bo > 1. However, since Bo ≈ 1 represents the boundary condition for the onset of this behaviour, our analysis in Figure 6 focused on variations near this critical regime. In our response to Comment 3c, we highlighted the influence of rising height on the bubble bouncing boundary and clarified the changes in the bouncing mechanism induced by increasing rising height.\nRegarding whether the same behaviour occurs for Bo > 1, in this regime bubbles do not bounce, regardless of the rising distance. Instead, they undergo breakup, which can be broadly classified into two types:\na) In-flight breakup: The bubble ruptures centrally during the ascent, forming the “central breakup region V” mentioned in Tripathi et al. 1\nb) Impact-induced breakup: The bubble remains intact during the ascent but ruptures centrally upon impact with the surface (as shown in Fig. 2e, f). This can be further categorized into symmetric (Fig. 2e) and asymmetric breakup (Fig. 2f).\nIn summary, the transition in impact modes induced by rising distance (specifically from bounce to non-bounce) observed in the regime of Ga > 40 and Bo ≈ 1 does not occur for Bo > 1. Undoubtedly, compared with short rising distances, longer rising distances significantly influence the ascent trajectory and bubble morphology, which in turn affect the impact behavior. However, in the section “Effect of rising distance on bubble dynamics,” our primary objective was to examine the influence of rising distance on the phase boundaries. Therefore, we did not further distinguish the specific rupture modes within the breakup regime in that section."}, {"round": 3, "reviewer_comments": {"reviewer_1": "The reviewer thanks the authors for their detailed responses and for the revisions addressing the previous comments, which have significantly improved the clarity and overall quality of the manuscript. I have only one minor remark concerning lines 399–401 of the revised manuscript: in Basilisk, interface reconstruction and advection are typically performed using geometric fluxes based on a PLIC method, rather than the Bell–Colella–Glaz advection scheme. Subject to this minor clarification, I recommend the manuscript for publication."}, "author_response": "Response: We sincerely thank the reviewer for the careful reading and for pointing out this important clarification regarding the advection and interface reconstruction schemes in Basilisk. We agree with the reviewer that interface reconstruction is performed using the geometric PLIC method rather than the Bell–Colella–Glaz (BCG) scheme. Our original sentence was ambiguous and may have caused misunderstanding.\nIn the revised manuscript, we have now corrected the description to clearly distinguish the numerical treatment of the velocity advection term from the interface reconstruction method. The revised text is provided (Revised manuscript, from Page 12 Lines 397),\n“The solver utilizes a generic time loop with a CFL-limited time step. The velocity advection term is computed using the Bell–Colella–Glaz second-order unsplit upwind scheme, while the piecewise linear interface calculation (PLIC) method is employed for geometric interface reconstruction40,41.”"}]

Decision letter