2.1. Interfacial curvature

Figure 3(a) shows the probability density functions (PDFs) of the interfacial curvature $\kappa$ normalized by the bubble radius $D/2$ for bubbles that eventually break in Exp A. Here, $\kappa >0$ and $\kappa <0$ represent convex and concave interfaces, respectively; $t=0$ indicates the instant when the breakup occurs, and $t<0$ represents the moment before the bubble breakup; $\tau _\eta$ is the Kolmogorov time scale. The PDFs for the entire duration, from 10$\tau _\eta$ prior to the breakup to the breakup moment, follow a similar distribution, although a small difference over time can be observed. The peaks of all the PDFs are located at $\kappa D/2 \approx 0.9$, which is close to the undeformed spherical geometry with the mean curvature exactly at $\kappa D/2=1$. In addition, all the PDFs are positively skewed with a higher probability of finding interfaces with local positive curvatures, i.e. extruded outwards from the gas phase to the ambient liquid. This positive skewness is due potentially to the incompressibility of the inner gas. As part of the bubble interface is compressed by surrounding turbulence (which will be discussed in the following), the rest of the interface tends to extrude outwards. This extrusion might be not uniform but significantly sharp along certain directions with minimum resistance, resulting in $\kappa D/2\gg 1$.

Figure 3.

(a) The PDFs of the normalized interfacial curvature $\kappa D/2$ for all breaking bubbles in Exp A. Blue represents different times before the breakup. The red line represents the prediction based on the model (2.5). The solid part of the red line marks the range of $\kappa$ where the model is valid. (b) The same plots for Exp B.

Note that the negative $\kappa$ represents the local depression on the bubble interface, which is associated with the increase of local dynamic pressure. If this depression is considered as a result of the collision between the bubble interface and sub-bubble-scale eddies (Luo & Svendsen Reference Luo and Svendsen1996), then the left tail of the PDF can be understood and modelled by accounting for such interactions.

Let us picture a simple scenario in which a bubble with diameter $D$ encounters an energetic, sub-bubble-scale eddy with size $D_e< D$, as shown in the schematic of figure 3(a). To leading order, the local curvature can be assumed to scale with $1/D_e$, and the other part of the interface remains unchanged given the short interaction time. In order for the eddy to depress the local interface, the inertia of the eddy $\rho u_e^2$ must be comparable with the surface tension induced by the local curvature, i.e. $\sigma /D_e$, where $u_e$ is the velocity of the eddy. This relation leads to $\rho u_e^2>C_1\sigma /D_e$, where $C_1$ is a fitting parameter that will be discussed in detail later. Rearranging this equation leads to the minimum eddy velocity required to deform the interface:

(2.1)\begin{equation} u_{e,d}=\sqrt{C_1\sigma/(\rho D_e)}. \end{equation}

The conditional PDF of instantaneous eddy velocity $u_e$ in turbulence for a given eddy size $D_e$ can be expressed as

(2.2)\begin{equation} P(u_e\,|\,D_e)=3\sqrt{2}\,\epsilon_e^{2/3} D_e^{-1/3}\,P(\epsilon_e)/2, \end{equation}

where $\epsilon _e$ is the local energy dissipation rate at the eddy length scale (Qi et al. Reference Qi, Tan, Corbitt, Urbanik, Salibindla and Ni2022). The distribution of $\epsilon _e$ can be approximated by a log-normal function $P(\epsilon _e)=1/(\epsilon _e\sqrt {2{\rm \pi} \sigma _{\ln \epsilon }^2})\exp {[-(\ln {(\epsilon _e/\langle \epsilon \rangle )}+\sigma _{\ln \epsilon }^2/2)^2/{2\sigma _{\ln \epsilon }^2}]}$, given by the multi-fractal model (Kolmogorov Reference Kolmogorov1962; Meneveau & Sreenivasan Reference Meneveau and Sreenivasan1991). Here, $\sigma _{\ln \epsilon }^2=A+\mu \ln (L/D_e)$ is the variance; $A$ represents a large-scale variability, which is set at $A=0$ for convenience; $\mu \approx 0.25$ is the intermittency exponent; and $L$ is the integral length scale of turbulence.

Given the distribution of the instantaneous eddy velocity (2.2), the PDF of the size of the eddies that are sufficiently strong to depress the bubble interface can be expressed by

(2.3)\begin{equation} P(D_e)\sim D_e^2\omega_c \int_{u_{e,d}}^\infty P(u_e\,|\,D_e)\,{\rm d} u_e. \end{equation}

This equation accounts for the interfacial area depressed by the bubble, which scales with ${\sim }D_e^2$ (as illustrated in figure 3a), and the frequency of the collision between the bubble and an eddy of size $D_e$ (Luo & Svendsen Reference Luo and Svendsen1996; Qi et al. Reference Qi, Tan, Corbitt, Urbanik, Salibindla and Ni2022), which can be estimated using

(2.4)\begin{equation} \omega_c\sim \langle\epsilon\rangle^{1/3}D^2D_e^{-11/3}. \end{equation}

Since the local curvature of the depression can be approximated by $\kappa \approx -2/D_e$, the PDF of the local negative curvature on the bubble interface can therefore be expressed as

(2.5)\begin{equation} P(\kappa)\sim\kappa^{-2}\,P(D_e). \end{equation}

The predicted PDF for $\kappa <0$ based on (2.5) is shown in figure 3(a) as the red solid curve, with only one fitting parameter $C_1=0.36$ to set the minimum eddy velocity $u_{e,d}$. It is evident that the model prediction agrees with the experimental results but only for a range of $\kappa$ because (2.5) works only for eddies within the inertial range. As a result, the predicted PDF extends only up to $\kappa D/2\approx -10$, corresponding to the eddy size of $D_e\approx 1$ mm, which is close to the lower limit of the inertial range (see Appendix C). In addition, the model can describe only local deformation so it cannot be used to predict the distribution of small curvature that corresponds to large deformation of the bubble size, i.e. $\kappa D/2\approx -1$. Nevertheless, the overall agreement for the scales considered suggests that the bubble local deformation is likely driven by the collision with sub-bubble-scale eddies at least statistically. In figure 3(b), the analysis above is also conducted for Exp B for bubbles with size 3–6 mm, which is around half of the bubble sizes in Exp A. The PDFs follow distributions similar to those in figure 3(a), and a similar good agreement between the model and experiments for the range of scales considered can still be established by setting the fitting parameter $C_1=0.14$ in the model.

Prior to the breakup, as shown in figure 2 by the green shaded area, bubbles experience strong deformation, including frequent stretching and depression, so the variation of curvature could be an indicator of the breakup dynamics. Figure 4(a) shows the time evolution of the averaged standard deviation of the curvature $\langle \sigma _{\kappa D/2}\rangle$ over all the breaking-bubble trajectories before the breakup occurs for Exp A. It is seen that $\langle \sigma _{\kappa D/2}\rangle$ remains approximately constant until approximately $t=-9\tau _\eta$, when $\langle \sigma _{\kappa D/2}\rangle$ begins to increase gradually, indicating a substantial variation of the interfacial curvature. The associated time scale $\tau _e\approx 9\tau _\eta$, as indicated by the blue shaded area, is then considered as the characteristic time scale of such a turbulence-driven inertial deformation process (Ruth et al. Reference Ruth, Aiyer, Rivière, Perrard and Deike2022) (blue shaded area in figure 2).

Figure 4.

(a) The time evolution of the averaged standard deviation of the curvature $\langle \sigma _{\kappa D/2}\rangle$ for Exp A. The blue shaded area marks the time scale $\tau _e$. (b) The same plots for Exp B. (c–e) The time evolution of the mean aspect ratio $\langle \alpha \rangle$ (solid lines) of breaking bubbles for Exp A, B and C, respectively. The blue shaded area marks the time scale $\tau _e$.

In addition to the curvature, another quantity that reflects the deformation of bubbles is the aspect ratio $\alpha$. Figure 4(c) shows the time evolution of the mean aspect ratio $\langle \alpha \rangle$ of all the breaking bubbles for the same dataset prior to the breakup moment. Here, the aspect ratio $\alpha =2r_{maj}/D$ is defined as the ratio between the semi-major axis $r_{maj}$ and the bubble spherical-equivalent radius $D/2$. In figure 4(c), it is evident that $\langle \alpha \rangle$ far from the breakup moment remains almost constant with $\langle \alpha \rangle >1$, suggesting that bubbles experience strong temporal oscillation (as indicated by the green shaded area in figure 2). However, since those fluctuations do not follow the same frequency or phase, averaging the signals over many bubbles results in a nearly constant $\langle \alpha \rangle$. Close to the breakup moment, $\langle \alpha \rangle$ exhibits clear growth, implying that $\alpha$ for all bubbles is likely to increase during this time period. A method is applied to extract this time scale by calculating the intersection between two lines: the dotted line that captures the early plateau, and the dashed line that is fitted over the range of $\langle \alpha \rangle$ from the half-height to the peak, as illustrated in figure 4(c). A time scale approximately $7.8\tau _\eta$ can be then extracted as indicated by the blue shaded area. This time scale is similar to $\tau _e\approx 9\tau _\eta$ obtained from figure 4(a), and is thus considered as the same time scale associated with the bubble inertial deformation driven by turbulence.

The same analysis of $\langle \sigma _{\kappa D/2}\rangle$ and $\langle \kappa \rangle$ is repeated for Exp B with 3–6 mm bubbles, as shown in figures 4(b,c), respectively. Most of the previous discussion regarding Exp A remains the same here, and a shorter inertial deformation time scale $\tau _e\approx 6\tau _\eta$ can be obtained consistently from the time evolution of both $\langle \sigma _{\kappa D/2}\rangle$ and $\langle \kappa \rangle$. For Exp C, since the bubble is too small so that the reliable reconstruction of the curvature is not possible, $\tau _e$ is determined only based on the time evolution of $\langle \alpha \rangle$ (figure 4e).

2.2. The inertial deformation time scale

Following a procedure similar to that discussed above, the time scale of the inertial deformation, i.e. $\tau _e$, is extracted and shown in figure 5 for all the datasets. In the case when $\tau _e$ obtained based on the time evolution of $\langle \sigma _{\kappa D/2}\rangle$ and $\langle \kappa \rangle$ is slightly inconsistent (e.g. figures 4a,c), the averaged $\tau _e$ is used.

Figure 5.

The inertial deformation time scale $\tau _e$ as a function of the bubble size. Squares represent the experimental data from all the datasets. The solid lines are the predictions based on the model (2.8). The purple dashed lines are the turn-over times of the bubble-sized eddies. The black dashed line represents the second-mode natural oscillation time scale of the bubble.

Given the fact that $\tau _e$ is driven by surrounding turbulence, two simple approaches are available to model $\tau _e$. The first approach is the Kolmogorov–Hinze (KH) framework (Kolmogorov Reference Kolmogorov1949; Hinze Reference Hinze1955), based on which the time scale of the inertial deformation process $\tau _e$ can be estimated using the turn-over time scale of the bubble-sized eddy, i.e. $\tau _D= D/(\sqrt {C_2}(\langle \epsilon \rangle D)^{1/3})$, which is shown in figure 5 as purple dashed lines for various $\langle \epsilon \rangle$. The other approach is to associate $\tau _e$ with the resonance oscillation of the bubble, which could also lead to breakup (Risso & Fabre Reference Risso and Fabre1998). The time scale of this resonance oscillation is given by $\tau _2=2{\rm \pi} \sqrt {\rho D^3/(96\sigma )}$ (Lamb Reference Lamb1879), which is shown in figure 5 as black dashed lines. In figure 5, it is evident that both approaches overestimate $\tau _e$ compared to the experimental data.

The alternative way to model $\tau _e$ is by following an approach similar to that in the previous work by Qi et al. (Reference Qi, Tan, Corbitt, Urbanik, Salibindla and Ni2022), in which the sub-bubble-scale eddy contribution to the breakup is incorporated, and the breakup criterion is set by two relationships based on the inertia and time scale of the sub-bubble-scale eddy. Following this work, the minimum requirement of the eddy velocity $u_{e,b}$ (or equivalently, the eddy kinetic energy) to break the bubble can be written as

(2.6)\begin{equation} u_{e,b}(D_e,D)=\max\left(\sqrt{\sigma/(\rho D_e)},\sqrt{96\sigma/(4{\rm \pi} \rho D^3 D_e^{-2})}\right)\!. \end{equation}

Given $u_{e,b}$, considering the distribution of eddy velocity $P(u_e\,|\,D_e)$ of an eddy (2.2), the probability $P'_b$ for this eddy to break the bubble is expressed following the integration

(2.7)\begin{equation} P'_b=\int_{u_{e,b}}^\infty P(u_e\,|\,D_e)\,{\rm d} u_e. \end{equation}

In order to estimate the time scale, it is assumed that the time required for a bubble to be broken by an eddy with velocity $u_e$ is given by $D/(2u_e)$. The factor 2 in the denominator comes from the observation that typically, the formation of the neck involves retraction of the interface simultaneously from both sides. Considering the distribution of eddy velocity $P(u_e\,|\,D_e)$ in (2.2), the breakup probability $P'_b$ from (2.7) and the collision frequency $\omega _c$ (2.4), the expected inertial time scale, i.e. $\tau _e$, can be expressed by

(2.8)\begin{equation} \tau_e=\frac{\displaystyle\int^D_{10\eta} \int_{u_{e,b}}^\infty D/(2u_e)\,P(u_e\,|\,D_e)\,P_b' \omega_c \,{\rm d} u_e \,{\rm d} D_e}{\displaystyle\int^D_{10\eta}P_b' \omega_c \,{\rm d} D_e}, \end{equation}

where the contribution from all the sub-bubble-scale eddies in the inertial range, from approximately ${\sim }10\eta$ ($\eta$ being the Kolmogorov length scale) to the bubble size, are incorporated. Eddies in the dissipation range are not considered here as they have negligible energy to break the bubble, and (2.2) and (2.4) no longer hold. Note that in this model, no free parameter is involved.

It is worth emphasizing the difference between (2.8) and the bubble lifetime model reported in Qi et al. (Reference Qi, Tan, Corbitt, Urbanik, Salibindla and Ni2022). The bubble lifetime model in Qi et al. (Reference Qi, Tan, Corbitt, Urbanik, Salibindla and Ni2022) considers the averaged time required for the bubble to encounter an eddy that eventually leads to a successful breakup, whereas (2.8) is the weighted average of the time between a successful collision event and the final breakup.

The prediction based on (2.8) for various $\langle \epsilon \rangle$ is shown as blue solid lines in figure 5. Compared to the classical Kolmogorov–Hinze framework, it is evident that the predicted $\tau _e$ agrees better with the experimental data, further emphasizing the role of energetic, sub-bubble-scale eddies in the inertial deformation process. It is noted that the Kolmogorov–Hinze framework (purple lines) could be modified to achieve reasonable agreement with the data by using a different eddy size. Specifically, the eddy size can be assumed to be $\beta D$, where $\beta$ is a coefficient, and the resulting eddy turnover time becomes: $\tau_{\beta D}= \beta D/(\sqrt {C_2}(\langle \epsilon \rangle \beta D)^{1/3})$). To bring the purple lines down to the experimental data in figure 5, it is evident that $\beta$ has to be smaller than 1. This result highlights our point that the inertial deformation time scale could be dominated by the turn-over time scale of a sub-bubble-scale eddy, consistent with the underlying mechanism of (2.8) although only one sub-bubble scale is considered here. Unfortunately, given the range of $We$ covered in this work, it is difficult to conclude which model (modified Kolmogorov–Hinze framework or (2.8)) better predicts the inertial deformation time scale. It is also worth noting that (2.8) slightly over-predicts $\tau _e$ at $D=11$ mm, and it is likely driven by the buoyancy effect. In particular, the Eötvös number $Eo=\rho g D^2/\sigma$ for such bubble size is around $Eo=16.8\gg 1$. This effect likely destabilizes the bubbles and results in an accelerated breakup process, especially for large bubbles with relatively low turbulence intensity.

2.3. Interfacial velocity

In addition to $\tau _e$, a different time scale can be extracted from the interfacial velocity statistics. Figure 6(a) shows the PDFs of the interfacial velocity $u_f$ for 3–6 mm breaking bubbles over different times in Exp B. Here, $u_f>0$ and $u_f<0$ represent outward and inward interfacial velocity, respectively, from the perspective of bubbles. As a result, $u_f>0$ indicates the interface moving from the gas phase towards the outer liquid phase, and vice versa. In this figure, $u_f$ is normalized by its standard deviation $\sigma _{u_f}'$. It is evident that the PDFs for $t<-2\tau _\eta$ collapse well and follow a similar distribution. However, the PDF at $-2\tau _\eta < t<0$ is slightly off, with a higher left tail and a lower right tail. In addition, all the PDFs clearly deviate from a Gaussian distribution, which is shown in figure 6(a) by the black dashed line. The much higher probability of finding extreme interfacial velocity implies that the interfacial velocity is more intermittent compared to the single-phase fluctuation velocity as shown in figure 6 by the grey solid line. Moreover, the PDFs of $u_f$ are all negatively skewed, indicating that the inward interfacial velocity is substantially stronger compared to the outward. This strong inward interfacial velocity could be linked to the distribution of the interfacial curvature, which shows clear positive skewness (see figure 3). As the majority of the interface is convex with $\kappa D/2\gg 1$, the resulting surface tension pointing inwards from liquid to gas tends to significantly accelerate the contraction of the interface, leading to stronger negative interfacial velocity.

Figure 6.

(a) PDFs of the interfacial velocity $u_f$ for 3–6 mm breaking bubbles in Exp B. Purple solid lines represent the different times before the breakup. The black dashed line marks the Gaussian distribution. The grey solid line shows the PDF of one component of the single-phase fluctuation velocity $u_x^\prime$ normalized by its standard deviation. (b) The time evolution of the averaged standard deviation of the interfacial velocity $\langle \sigma _{u_f}\rangle$. The purple shaded area marks the time scale $\tau _c$.

In addition to the PDFs, the time evolution of the standard deviation of the interfacial velocity, averaged over all the breaking-bubble trajectories $\langle \sigma _{u_f}\rangle$, is shown in figure 6(b). Here, $\langle \sigma _{u_f}\rangle$ remains almost constant at the beginning, indicating that the distribution of $u_f$ does not experience significant changes. At $t\approx -6\tau _\eta$, approximately one inertial deformation time scale earlier, before the breakup ($\tau _e\approx 6\tau _\eta$), $\langle \sigma _{u_f}\rangle$ begins to grow. Such growth suggests that the inertial deformation process of the bubble could affect the distribution of the interfacial velocity. Later, at approximately $t\approx -2\tau_\eta$ to $-1\tau _\eta$, another substantial increase of $\langle \sigma _{u_f}\rangle$ is observed. This sudden increase leads to a new time scale $\tau _c\approx \tau_\eta$ to $2\tau _\eta$, as indicated by the purple shaded area in figure 6(b). Note that this time scale $\tau _c$ is significantly shorter than the corresponding inertial deformation time scale ($\tau _e\approx 6\tau _\eta$), indicating that a distinct physical process occurs right before the breakup.

The mechanism that leads to the sudden increase of $\langle \sigma _{u_f}\rangle$ is likely linked to the necking process before the neck finally pinches off, as illustrated by the purple shaded area in figure 2. With a sufficient amount of deformation, the neck shrinks rapidly right before the breakup, which is accompanied by a fast and local inward interfacial velocity. This interfacial velocity eventually leads to the sudden increase of $\langle \sigma _{u_f}\rangle$ in figure 6(b), and the rise of the left tail of the PDF in figure 6(a) during $-\tau _c< t<0$.

Based on the discussion above, the time scale $\tau _c$ of the final pinch-off can be modelled using the capillary time scale that has been discussed before (Villermaux Reference Villermaux2020; Rivière et al. Reference Rivière, Ruth, Mostert, Deike and Perrard2022; Ruth et al. Reference Ruth, Aiyer, Rivière, Perrard and Deike2022). The model of the capillary time scale follows $\tau _c=\rho ^{1/2} \sigma ^{-1/2} \delta ^{3/2}/(2\sqrt {3})$, where $\delta$ is the typical width of the bubble neck that eventually becomes unstable and leads to the rupture of the interface. Here, $\delta$ can be estimated using the size of the smallest daughter bubble resulting from the breakup (Rivière et al. Reference Rivière, Ruth, Mostert, Deike and Perrard2022). In our experiments, the neck width and the daughter bubble size vary from case to case. The mean neck width is estimated of the order of $\delta \sim O(1)$ mm for bubbles with $D\gtrsim 3$ mm. Using $\delta =1$ mm as the neck width leads to an estimation $\tau _c\approx 1.5\tau _\eta$, which is reasonably close to the range of $\tau _c$ measured from the experiment (figure 6b).

2.4. Breakup time scales and separation

Figure 2 illustrates two time scales associated with the breakup process, including the inertial deformation time $\tau _e$ and the capillary time scale $\tau _c$. In addition, the bubble lifetime $\tau _l$, which is defined as the time interval between two consecutive breakups, is also shown, as indicated by the red shaded area. During this bubble lifetime, bubbles experience complicated deformation and oscillation due to the continuous bombardment by eddies of various sizes (green shaded area), until the event that the bubble encounters intense eddies, leading to strong inertial deformation followed by the neck thinning, which leads ultimately to the breakup.

To further compare the bubble lifetime $\tau _l$ with the other time scales, the experimental data from previous works are included. Figure 7 shows the measured bubble lifetime $\tau _l$ by Vejražka et al. (Reference Vejražka, Zedníková and Stanovský2018) (red symbols) and Martínez-Bazán et al. (Reference Martínez-Bazán, Montanes and Lasheras1999) (blue symbols) as functions of the Weber number. Here, $\tau _l$ is calculated following $\tau _l=1/g_b$, where $g_b$ is the bubble breakup frequency that can be measured experimentally (Håkansson Reference Håkansson2020). Both datasets show distinct trends. Specifically, Vejražka et al. (Reference Vejražka, Zedníková and Stanovský2018) reported a decrease in bubble lifetime as $We$ increases, while Martínez-Bazán et al. (Reference Martínez-Bazán, Montanes and Lasheras1999) suggested a slightly increasing $\tau _l/\tau _\eta$. This difference in the Weber number dependence marks the decades-long debate on whether the lifetime of bubbles increases or decreases as a function of $We$.

Figure 7.

The comparison among different time scales associated with bubble breakup predicted by models, including the bubble lifetime $\tau _l$ (red solid line), the inertial deformation time scale $\tau _e$ (blue solid line), and the capillary time scale $\tau _c$ (grey solid line), as functions of $We$. The red and blue symbols represent the experimental data of $\tau _l$ by Vejražka et al. (Reference Vejražka, Zedníková and Stanovský2018) and Martínez-Bazán et al. (Reference Martínez-Bazán, Montanes and Lasheras1999), respectively. The red and blue shaded areas mark $We<10$ and $We>10$.

In addition to the experimental results, the solid red line in figure 7 also shows the prediction by a bubble lifetime model (Qi et al. Reference Qi, Tan, Corbitt, Urbanik, Salibindla and Ni2022) for 3 mm bubbles at various $\langle \epsilon \rangle$. This model already considers the contribution from sub-bubble-scale eddies, i.e. bubble breakups are accelerated by the collisions with sub-bubble-scale energetic eddies with the bubble lifetime dominated by the frequency of those extreme events. The model agrees well with the experimental data for $We\lesssim 10$; however, it decays continuously and underestimates the lifetime for $We\gtrsim 10$. Note that this model does not account for the fact that the lifetime $\tau _l$ cannot drop indefinitely as bubbles still need at least one inertial deformation time scale $\tau _e$ to complete the breakup process. Figure 7 shows the predicted $\tau _e$ by (2.8) for 3 mm bubbles as the blue solid line. The result by (2.8) is consistent with experimental data for $\tau _l$ at the large Weber number limit. Note that both experimental datasets can be explained by considering the maximum between the predicted bubble lifetime (Qi et al. Reference Qi, Tan, Corbitt, Urbanik, Salibindla and Ni2022) and the predicted inertial deformation time scale (2.8), which cross over each other at $We\approx 10$. The agreement suggests that different trends observed in the experiments were the result of the transition of different time scales at play.

This result also highlights the role of the Weber number in determining the time scale separation for breakup in turbulent two-phase flows. For small $We$, since the majority of collisions between bubbles and surrounding eddies lead only to deformation but not breakup, the lifetime of the bubble is dominated by the long oscillation stage, and the scale separation between the lifetime and the final inertial deformation time scale is significant. Such scale separation depends only on $We$. For large $We$, the two time scales become one, and no scale separation can be observed. This is analogous to the role played by the Reynolds number in controlling the scale separation in single-phase turbulence where the separation between the Kolmogorov scale and integral scale diminishes as $Re$ decreases. The capillary time scale $\tau _c$ (Rivière et al. Reference Rivière, Ruth, Mostert, Deike and Perrard2022) (grey solid line in figure 7), however, remains nearly a small constant regardless of $We$ as it is not affected by surrounding turbulence.