Coalescence-induced jumping of unequal-sized droplets on superhydrophobic surfaces

With the rapid development of the electronic industry and continuous increase in integration density, the concern of device reliability due to thermal issues is increasing.^1–3 This trend has made the heat dissipation in small areas particularly prominent, thus requiring advanced heat dissipation techniques to be urgently addressed. Effective thermal management not only helps to ensure the performance and lifespan of electronic devices but also plays a positive role in energy efficiency and environmental protection.^4,5 In order to efficiently handle the excess heat generated, phase-change heat transfer, especially boiling and condensation, is considered to be an effective approach. Compared to traditional single-phase heat transfer, phase-change heat transfer embraces a higher heat transfer capacity, which is crucial in coping with the increasing heat load.^6,7 To address the issues arising from the increasing heat load, it becomes essential for cooling devices to employ efficient condensation methods as a viable solution. Therefore, it is necessary to have a deep understanding of the underlying mechanisms of the condensation process because the condensation process is considered to be a crucial component in phase-change heat transfer,⁸ including nucleation, growth, and droplet coalescence, which significantly affect the heat transfer performance. During condensation, the liquid coolant undergoes a phase change from gas to liquid and forms droplets on the surface. The key issues of nucleation, growth, and droplet coalescence determine the heat transfer performance.⁹ It is worth noting that droplet condensation is often regarded as the preferred way for high heat flux dissipation as it usually has a higher heat transfer performance than film condensation.¹⁰ However, in order to achieve sustainable droplet condensation, favorable nucleation sites need to be maintained, and the formed droplets need to be quickly removed. Wetting control plays a crucial role in this aspect. The studies have shown that by adjusting the surface wetting properties, nucleation and growth of droplets can be promoted while ensuring the rapid removal of droplets, thereby improving the heat transfer performance.¹¹ Superhydrophobic surfaces, capable of repelling water and preventing wetting, have gained significant attraction when utilized in the fields of self-cleaning,^12–15 anti-freezing,^16–19 microfluidic manipulation technology,^20–23 and heat transfer enhancement.^8,24–26 One of the fascinating phenomena observed on a superhydrophobic surface is the jumping behavior after coalescence, where two droplets coalesce and propel themselves off the surface in a jumping motion.²⁷ This jumping behavior has been studied extensively in recent years through experimental and numerical approaches.^28–32 

Previous investigations primarily highlighted the jumping behavior of equal-sized droplets on superhydrophobic surfaces,^33–35 with less attention given to unequal-sized droplets, which are frequently encountered in practical situations. Despite the extensive exploration of jump-off with a critical droplet size ratio,^36,37 the underlying physical mechanisms governing the specific dynamic behaviors await further elucidation. In particular, the association between droplet size, contact angle, and the jumping behavior of unequal-sized droplets has not yet been extensively studied.³⁸ Regarding the effect of contact angle and droplet size, an investigation has been attempted to probe the jump-off motion of a squeezed droplet within grooves.³⁹ However, the jump-off motion is attributed to the effect of solid–liquid interfacial tension between the surface and the droplet. In contrast, the jump-off motion of the coalesced droplet is not triggered by the compressive effect at the solid–liquid interface, whereas it depends on the fluid flow induced by the interfacial tension at the gas–liquid interface. The mechanisms governing these two situations in response to variations in the contact angle and droplet size are different. In addition, although the energy conversion mechanism of the coalesced droplet has been investigated previously,^40,41 a comprehensive understanding is still lacking, particularly for unequal-sized droplets, which is essential for improving the energy conversion efficiency for droplet coalescence in potential applications, such as dropwise condensation heat transfer.¹¹ 

This study is devoted to the jumping behavior of unequal-sized droplets upon coalescence on a superhydrophobic surface, and it is numerically explored to provide a deep understanding of the association between contact angle, droplet size, and the jumping behavior and energy conversion mechanism. A continuous surface force model is used to describe interfacial tension effects, and a volume of fluid (VOF) coupled level-set method is utilized to track the interface boundaries, thereby capturing the evolving shapes of the coalesced droplet at different time instants. Variations in parameters during the droplet coalescence process, such as changes in the energy components involved, are extracted to better understand the energy conversion. The numerical model is first validated by comparing the numerical results with experimental data, including the droplet shapes at various moments. Further benefiting from the advantages of numerical methods, the evolution of the interface area over time is precisely captured to estimate the temporal variations in the surface energy. The evolution of other energies involved during the droplet coalescence process, such as oscillatory and translational kinetic energy, is also estimated to gain insights into the specific motion patterns. Oscillatory motion is quantitatively described by the axial length of the deformed droplet in three directions. For translational motion, the centroidal velocity and height are estimated to evaluate the jump-off performance of the coalesced droplet. Building upon such a comprehensive and executable evaluation approach, the present study explores the mechanism of altering the contact angle and droplet size on the jump-off performance of the coalesced droplet. The jump-off performance refers to the capability of the coalesced droplet to detach from a surface, which is closely related to the jump-off criteria. The jump-off suggests that the droplet detaches from the surface in an ascending manner and no longer touches the surface. In summary, the findings of the present research can provide guidance for the design and development of surfaces that exhibit enhanced jump-off performance of the coalesced droplet and offer insights into the establishment of rigorous criteria to judge the jump-off event.

A dimensionless parameter is helpful when investigating the effect of droplet size on physical quantities because it enables comparability among the droplet with different sizes. The energy is normalized by σ_lgr_t², and the interfacial area is normalized by 4πr_t². Here, the characteristic radius r_t is chosen as the radius of the liquid droplet corresponding to the average volume of two coalesced droplets, and r_t is equivalent to any radius of two coalesced droplets when two coalesced droplets are equal in size.

The centroidal velocity of two coalesced droplets is defined as $Ui=1Vdrop(∫VdropuidVcell)$⁠. The velocity and time scale are non-dimensionalized with the characteristic velocity and the characteristic time. According to the capillary-inertial scale, the characteristic time scale is defined as $tci=ρlrt3σ$⁠, and the corresponding characteristic velocity is defined as $uci=σρlrt$⁠.

In the computational domain, the coalescence of droplets on a wall in the three-dimensional framework is simulated, and the computational domain is shown in Fig. 1. The Navier-slip model is employed to explain the behavior of the contact line on the wall boundary, which permits movement of the contact line that is proportional to shear stress at the point of contact. The velocity of the contact line is calculated as $uTPCL=ε×∂ui∂nwall$⁠, with u_TPCL being the speed of the three-phase contact line (TPCL) and $∂ui∂nwall$ being the shear strain rate on the wall. The proportionality constant ε is equal to the size of a grid cell. By implementing a slip condition on the wall, the issue of infinite shear stress near the contact line is effectively mitigated.

The symmetry of the computational domain is used to reduce the number of grids in the numerical calculation, as shown in Fig. 1. In the computational domain, aside from the symmetric and wall boundaries, all other boundary conditions are set as pressure outlets, whose temperature and pressure are 293 K and 1 atm, respectively. The densities of water droplet and air are 998.2 and 1.2 kg/m³, with dynamic viscosities of 1.0 × 10⁻³ and 1.8 × 10⁻⁵ Pa s and thermal conductivities of 0.6 and 2.6 × 10⁻² W/(m K), respectively. At the initial moment, the gas phase and the liquid phase are stationary, and the droplet movement occurs due to the fluid flow triggered by the surface tension at the tangential interface of two droplets.

The discretization method of governing equations is based on finite volume, and the second order upwind scheme is used for spatial discretization. The pressure implicit with splitting of operator (PISO) scheme is used for pressure–velocity coupling. The numerical calculation for each step is considered to be convergent when the residuals of all equations are below 1.0 × 10⁻⁶.

In the computational domain, the motion of the coalesced droplet is the focus of attention; meanwhile, the velocity gradients at the interface are very large. Therefore, it is necessary to densify the grids for local refinement around the droplets, as shown in Fig. 1. In order to verify the grid independence, three cases with the number of grids of 1083k, 6125k, and 7043.75k are numerically investigated, where two droplets, with radii of 200 and 323 μm, are initially placed on the surface with a contact angle of 169° at the bottom of the computational domain. During the numerical calculation, the vertical kinetic energy of the coalesced droplet over time is monitored to check the grid independence, as shown in Fig. 2, because it is an important parameter to evaluate the jump-off performance of the coalesced droplet. In order to avoid the influence of the grid number on the estimation of the jump-off time and the vertical kinetic energy at the jump-off moment, the cases with different grid numbers are continued until the jump-off occurs, and then the jump-off time and the vertical kinetic energy at the jump-off moment are compared. According to the numerical results, the case with 6125k grids and the case with 1083k grids show differences of around 4% in dimensionless jump-off time and around 8% in dimensionless vertical kinetic energy, while the case with 7043.75k grids shows less than 1% difference in jump-off time and dimensionless vertical kinetic energy compared with those of 6125k grids. Therefore, the grid number of 6125k (66 grids on the diameter of the small droplet) is adopted in the numerical calculation to maintain enough accuracy without too much computational time and too many resources.

The coalescence behavior with a small droplet radius R_s of 200 µm, large droplet radius R_b of 323 µm, and contact angle of 169° is first investigated to validate the numerical model, and all parameters are the same as those in the experiments conducted by Yan et al.⁴² A slippery condition is employed to explain the three-phase contact line (TPCL) behavior. The VOF coupled level-set method is used for simulating the fluid flow and tracking the interface in this multiphase system because it can accurately capture the shape and movement of complex interfaces.⁴³ 

The numerical and experimental results are compared in Fig. 3, taking the case with a contact angle of 169°, R_s = 200 µm, and R_s/R_b = 0.62 as an example, where it can be observed that the numerical droplet shapes at different time instants agree well with the experimental findings (video S1).

Figure 3 illustrates a transient jumping behavior, and Fig. 4 depicts the temporal variations in the relevant physical quantities during the evolution of the phenomenon. In addition to velocity fields being displayed on the droplet interface, the variations in interfacial areas (liquid–solid and liquid–gas interfaces) and energy conversion with time during the coalescence are investigated to better interpret the jumping behavior [see Figs. 4(a) and 4(b)]. The total kinetic energy, including translational kinetic energy (TKE) and oscillatory kinetic energy (OKE), is shown in Fig. 4(c), and the dimensionless centroidal velocity is shown in Fig. 4(d).

Here, the hydrodynamic behavior of the coalesced droplet on a superhydrophobic surface is presented. Before the dimensionless time instant t* = 0.28, the droplets undergo coalescence to shape into a dumbbell geometry due to the fluid motion caused by surface tension. During this period, the liquid bridge undergoes rapid expansion, and the fluid flow decreases the liquid–gas interfacial area [see Fig. 4(a)]. As two droplets move toward the middle liquid bridge, the fluid flow also results in a decrease in the liquid–gas interfacial area [see Fig. 4(a)], and thereby the surface energy (SE) is released. About 45.1% of the released SE is converted into total kinetic energy, and the remaining portion is dissipated due to the viscous effect [see Fig. 4(b)]. At this moment, the liquid bridge is in a fast expansion stage, so there is a large local velocity gradient of the liquid bridge, as shown by the velocity vector in Fig. 3 at t* = 0.28, which causes a rapid increase in viscous dissipation [see Fig. 4(b)]. Since the droplet size is very small, the variation in gravitational potential energy is tiny [see Fig. 4(b)]. When t* is between 0.58 and 0.86, there is a pseudo-jump-off due to the fluid flow caused by the expanding liquid bridge, which makes the coalesced droplet temporarily leave the surface. At t* = 0.86, the expanding liquid bridge touches the underlying surface, and the further downward expansion of the liquid bridge is impeded, causing the bottom of the droplet to expand outward along the TPCL (see Fig. 3). Consequently, the liquid–solid interfacial area gradually increases [see Fig. 4(a)]. The two coalescing droplets continue to converge toward the middle, resulting in a continuous reduction in the liquid–gas interfacial area [see Fig. 4(a)], and the SE continues to be released and converted to total kinetic energy [see Fig. 4(b)]. At t* = 1.16, unlike two coalescing equal-sized droplets, the width of the liquid bridge between unequal-sized droplets gradually decreases, as shown in Fig. 3. At this moment, the conversion rate of the SE to kinetic energy is about 58.7%, and the large local velocity gradient at the liquid bridge causes a rapid increase in viscous dissipation. Due to the expansion of the TPCL, at t* = 1.44, the liquid–solid interfacial area increases and gradually approaches the maximum, as shown in Figs. 3 and 4(a). At t* = 1.72, the small droplet is completely engulfed where the velocity gradually turns outward, as shown in Fig. 3, and the liquid–gas interfacial area starts to increase, as shown in Fig. 4(a). At t* = 2.02, the droplet becomes nearly spheroidal (see Fig. 3), and the liquid–solid interfacial area experiences a local trough [see Fig. 4(a)]. A large amount of SE is released, resulting in an almost second peak of the total kinetic energy [see Fig. 4(b)]. At t* = 2.30, the total kinetic energy converts to SE, and the recovery of SE reaches the first peak while the total kinetic energy reaches the trough, as shown in Fig. 4(b). At t* = 2.60, the SE converts to total kinetic energy, on the contrary, and the total kinetic energy peaks for the third time at t* = 2.88. Subsequently, the second pseudo-jump-off occurs while the droplet moves upward at t* = 3.74, where the droplet undergoes a shape transformation from oblate to spheroid. Afterward, when the droplet changes from a spheroidal shape to a prolate shape, it touches the surface again at t* = 3.82 [see Fig. 4(a)], and its contact area with the surface reaches its maximum at t* = 4.10 [see the inset of Fig. 4(a)]. It should be noted that the real jump-off (t* = 4.60) means that the coalesced droplet no longer touches the surface during its upward motion [see Fig. 4(a) at t* = 10.00 for an example]. After the jump-off moment (t* > 4.60), the oscillation of the capillary wave shows the periodic characteristics with an oscillatory period of the SE and total kinetic energy of about 1.60 [see Fig. 4(b)], which is very close to the oscillatory period (π/2) of the equal-sized coalesced droplet jumping in the air.⁴⁴ Moreover, before the jump-off (t* < 4.60), there is also oscillation of the droplet, as can be seen in Fig. 4. At t* = 1.90, the coalesced droplet shapes into a spheroid for the first time [see the inset of Fig. 4(a)], which suggests that the SE reaches the trough for the first time [see Fig. 4(b)]. The second spheroidal shape [see the inset of Fig. 4(b)] and the second SE trough occur at t* = 2.90, as shown in Fig. 4(b). The third spheroidal shape [see the inset of Fig. 4(c)] occurs at the jump-off moment (t* = 4.60) when the SE is at the third trough [see Fig. 4(b)].

Figures 5(a) and 5(b) illustrate the energy conversion and two forms of kinetic energy at the jump-off moment at different contact angles, where the case with the small droplet radius R_s of 200 µm, R_s/R_b = 0.62, and contact angle of 169° serves as the baseline. In addition, Figs. 5(c) and 5(d) illustrate the temporal centroidal velocity and centroidal height as $h=1Vdrop(∫VdropzdVcell)$ at different contact angles.

When the contact angle decreases from 169° to 160° (video S2), about 2.3% more SE is released; however, with an increase of 5.5% in viscous dissipation, the total kinetic energy decreases by 7.9%, as shown in Fig. 5(a), leading to an 8.7% decrease in OKE and a 3.8% decrease in TKE, as shown in Fig. 5(b). However, the 3.8% reduction in TKE is mainly manifested in the decrease in the upward velocity W* [see Fig. 5(c)], leading to a significant reduction in centroidal height [see Fig. 5(d)]. It should be noted in Fig. 5(d) that the centroidal height is nearly the maximum and that the liquid–solid interfacial area is the minimum at t* = 4.60. However, the liquid–solid interfacial area does not reach zero at this moment, so the coalesced droplet fails to jump off.

When the contact angle increases from 169° to 180° (video S3), the increased viscous dissipation (3.4%) is greater than the increased SE (2.2%), so the total kinetic energy decreases by about 1.6%, as shown in Fig. 5(a). Although the total kinetic energy only decreases by 1.6%, the OKE increases by 4.3%, and the TKE decreases significantly by 35.1% [see Fig. 5(b)]. The decrease in the TKE is within the XOY plane (the non-jump-off direction), particularly for the U* in the X-axis direction [see Fig. 5(c)]. When the jump-off occurs at t* = 4.60 [see the inset of Fig. 5(d)], the upward velocity W* in the Z-axis direction (the jump-off direction) shows a 14.2% decrease compared to that of 169°, as shown in Fig. 5(c).

This is attributed to the fact that the deformation of the droplet on the surface is less drastic at a contact angle of 180° before the jump-off. The comparison between the deformation of the coalesced droplet at contact angles of 169° and 180° is shown in Fig. 6. It can be observed that the contraction and expansion of the coalesced droplet shape occur mainly in the XOY plane because the spheroidal axial length changes dramatically (the difference between a peak and a trough is almost 1.0) in the X and Y directions. Conversely, the spheroidal axial length curve shows the smallest amplitude (the difference between a peak and a trough is less than 0.3) in the Z direction. The amplitudes of the solid curve for the results of a contact angle of 169° are generally larger in the X and Y directions, as shown in Fig. 6, showing that the deformation of the coalesced droplet is more drastic at a contact angle of 169°. It signifies a more significant interference of the symmetry of the oscillatory process at a contact angle of 169°, resulting in a decreased oscillatory kinetic energy and increased translational kinetic energy, as shown in Fig. 6. Moreover, after the liquid bridge touches the surface (t* = 0.86), the further expansion along the TPCL on the surface is less, as proved by the smaller liquid–solid interfacial area [see the inset of Fig. 5(d)]. Therefore, the minimum deformation of the droplet on the surface at a contact angle of 180° resulted in more OKE and less TKE, which implies that a higher contact angle might not be necessarily beneficial for the jump-off. It is also demonstrated that two equal-sized droplets at a contact angle of 165° showed a larger jumping velocity than that at 180° upon coalescence.⁴⁵ 

A higher contact angle of the surface makes the droplets have less pinning and spread, indicating a potentially better performance of droplet motion. For example, in the groove-induced droplet jumping and vibration-induced droplet shedding, the jump-off performance exhibits a positive correlation with the contact angle.^39,46 However, the jump-off performance of the coalesced droplet displays a different correlation with the contact angle. Specifically, the droplet jump-off performance is superior at a contact angle of 169° when compared to the case of 180°. This can be attributed to the increased complexity of dynamic phenomena triggered by the interfacial tension at the gas–liquid interface, which differs from the inducement of the compressive effect³⁹ and the vibration effect.⁴⁶ 

Droplet size is also an indispensable factor in droplet coalescence behavior. Figures 7(a) and 7(b) show the energy conversion and two forms of kinetic energy in different droplet sizes at the jump-off moment. The case with the small droplet radius R_s of 200 µm, R_s/R_b = 0.62, and a contact angle of 169° serves as the baseline. Moreover, Figs. 7(c) and 7(d) illustrate the variation in the centroidal velocity and centroidal height with time in different small droplet sizes, respectively.

When increasing the small droplet radius from 200 to 400 µm (video S4), about 3.2% more SE is released, which mainly results in a 4.2% increase in viscous dissipation, while the total kinetic energy is almost the same [see Fig. 7(a)]. However, there is a significant difference in the two forms of kinetic energy, with a 5.3% increase in the OKE and a 28.9% decrease in the TKE [see Fig. 7(b)]. The 28.9% reduction in the TKE is mainly in the jump-off direction, showing a significant decrease in the upward velocity W* [see Fig. 7(c)]. Therefore, the lowest centroidal height in all cases can be seen from Fig. 7(d), and the liquid–solid interfacial area remains positive [see the inset of Fig. 7(d)], which results in the non-jump-off. In the other cases with smaller droplet sizes, the centroidal height rises higher, resulting in a more pronounced reduction in the liquid–solid interfacial area, eventually reaching zero with the occurrence of the jump-off.

As the radius of the small droplet decreases from 200 to 100 µm (video S5), there is an insignificant change in energy conversion, with the energy variation being less than 1% [see Fig. 7(a)], but with a 5.3% decrease in the OKE and a 25% increase in the TKE [see Fig. 7(b)]. The increased TKE is mainly in the jump-off direction velocity W* [see Fig. 7(c)], resulting in a higher centroidal height [see Fig. 7(d)].

As the radius of the small droplet decreases further from 200 to 50 µm (video S6), the release of the SE increases negligently, with the viscous dissipation increased by 3.0%, causing an 8.5% decrease in the total kinetic energy [see Fig. 7(a)]. The OKE and the TKE variations are more pronounced than the total kinetic energy, with a 15.5% decrease and a 31% increase, respectively [see Fig. 7(b)]. The higher TKE results in a 45.8% larger velocity W* in the jump-off direction, which is the largest among all cases [see Fig. 7(c)]. Meanwhile, the centroidal height is also the highest among all cases [see Fig. 7(d)].

It is not sufficient to merely study the centroidal velocity to compare the jumping performance, and more mechanism discussion is needed.⁴⁷ In addition, employing the characteristic number Oh for investigating the viscous dissipation and establishing its relationship with droplet size is also not sufficient.^36,37,48 Numerical results indicate that the dimensionless viscous dissipation is relatively insensitive to variations in droplet size. In correlation with this, in the dimensionless energy curve, negligible differences are observed in surface energy and kinetic energy following the coalescence of equal-sized droplets.⁴⁹ Therefore, the mechanism of size effect on droplet behavior needs more in-depth understanding, such as the translational and oscillatory kinetic energy. Compared with the case with the small droplet radius R_s of 200 µm, the case with the small droplet radius R_s of 50 µm has a 31.0% larger translational kinetic energy and a 15.5% smaller oscillatory kinetic energy at the jump-off moment [see Fig. 7(b)], while the surface energy increases by only 0.3% and the viscous dissipation increases by only 3.0% [see Fig. 7(a)].

The physical mechanism driving the behavior of droplets in air after coalescence is the conversion between surface energy and kinetic energy, which is accompanied by increasing viscous dissipation over time, indicating a decreasing magnitude of energy conversion. During the coalescence process on the surface, the downward extension of the coalesced droplet is hindered, leading to the disruption of its longitudinal symmetry. In the subsequent process of recoiling toward the center after extension, the fluid flow at the bottom of the droplet on the surface is restricted from flowing upward, resulting in an upward rise in the centroid of the droplet. If the rising height is significant enough to make the coalesced droplet bottom detach from the surface, the event of jump-off occurs, which is desired for droplet removal from the surface. In summary, the jump-off of the coalesced droplet is caused by the non-wetting surface breaking the symmetric oscillation of the coalesced droplet.⁴⁴ When the equal-sized droplets coalesce in the air, significant oscillations in physical quantities, such as axial length, occur with notable periodic characteristics.⁴⁴ It is worth noting that, even for two unequal-sized droplets, the oscillation of physical quantities during the coalescence shows the same periodic characteristics, regardless of the contact angle and droplet size [see Fig. 5(c), the inset of Figs. 5(d) and 7(c), and the inset of Fig. 7(d)]. During the oscillation after droplet coalescence, pseudo-jump-off may exist, but this is not the expected consequence since the droplet can still touch the surface afterward and undergo pinning instead of truly jumping off. To distinguish the pseudo-jump-off, a new criterion is proposed, which is based on the condition that the centroid of the droplet reaches a threshold height with an upward velocity, indicating a sufficient jump-off performance. Therefore, a jump-off criterion is suggested as the centroidal height of the coalesced droplet at any moment is larger than the radius of the sphere estimated from the summed volumes of the coalesced droplets and as the upward velocity W* is positive. Such a criterion is validated by the numerical results, as shown in Fig. 8.

In conclusion, this study provides a comprehensive and deep analysis of the jumping behavior of two unequal-sized droplets on a superhydrophobic surface in terms of the interfacial area, energy conversion, kinetic energy, and centroidal velocity. The findings indicate that simply increasing the contact angle does not necessarily benefit the jump-off because the minimum deformation of the droplet on the surface at a contact angle of 180° leads to a higher proportion of oscillatory kinetic energy and a lower proportion of translational kinetic energy. As the small droplet radius decreases from 200 to 50 µm, the kinetic energy decreases by 8.5%; however, the translational kinetic energy increases by 31%, which results in a 45.8% larger upward velocity. A criterion is proposed to judge whether the jump-off behavior occurs or not to help design superhydrophobic surfaces with controllable jumping behavior.

See the supplementary material for S1–S6, which are videos of the coalescence-induced behavior of unequal-sized droplets on superhydrophobic surfaces.

This research was supported by the National Natural Science Foundation of China, under Contract No. 51976117, and the National Key Research and Development Program of China, under Contract No. 2018YFA0702300.

The authors have no conflicts to disclose.

Ting-en Huang: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (equal); Software (lead); Validation (equal); Visualization (equal); Writing – original draft (lead). Peng Zhang: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (lead); Investigation (equal); Project administration (lead); Resources (lead); Supervision (lead); Writing – review & editing (lead).

The data that support the findings of this study are available from the corresponding author upon reasonable request.