Surface-temperature-induced Marangoni effects on developing buoyancy-driven flow

Abstract To investigate the initial development of the Rayleigh–Bénard–Marangoni (RBM) instability in a relatively deep domain, direct numerical simulations for a large range of Marangoni and Rayleigh numbers were performed. In the simulations, the surface was assumed to be flat and surface cooling was modelled by a constant heat flux. The small-scale dynamics of the flow and temperature fields near the surface was fully resolved by using a non-uniform vertical grid distribution. A detailed investigation of the differences in physical mechanisms that drive the Rayleigh- and Marangoni-dominated instabilities is presented. To this end, various properties such as the maturation rate of convection cells, the fluctuating kinetic energy and the surface characteristic length scale were studied. It was confirmed that buoyancy forces and surface-temperature-gradient-driven Marangoni forces enhance one another in promoting the development of the RBM instability. When using a relevant measure of the effective thermal boundary layer thickness as length scale, both the critical Marangoni and Rayleigh numbers, obtained for the purely Marangoni- and purely Rayleigh-driven instabilities, were found to be in good agreement with the literature.


Introduction
The present study is motivated by the problem of heat and mass transfer across the air-water surfaces of natural water bodies.An example of this is the transfer of heat and atmospheric gases (such as oxygen, carbon dioxide, methane) across the water surface in lakes and reservoirs.An accurate prediction of the transfer rate is important in order to produce a reliable global heat and greenhouse gas budget.Typically, wind shear is seen as the main source of turbulence generation that promotes interfacial heat and mass flux (Wanninkhof et al. 2009;Garbe et al. 2014).Other turbulence sources, especially buoyancy, are thereby usually neglected.Recently, buoyancy-driven heat and gas transfer induced by surface cooling, which is particularly important in lakes and ponds at low wind speeds, has attracted increasing attention (Rutgersson & Smedman 2010;Podgrajsek, Sahlee & Rutgersson 2014;MacIntyre et al. 2018).In this process, due to surface cooling (usually by evaporation), plumes and sheets of relatively heavy, cold (gas-saturated) surface water plunge down and are replaced by warm (unsaturated) bulk water, thereby promoting heat and gas transfer across the surface.
It is known that surface cooling not only results in unstable density gradients, but also may induce variation in surface tension due to local changes in surface temperature.Such variations in surface tension generate so-called Marangoni forces that induce flows from low-surface-tension (high-temperature) regions to high-surface-tension (low-temperature) regions.As a result of the buoyancy-induced and/or Marangoni-induced convective instabilities, convection cells are generated at the surface (e.g.Spangenberg & Rowland 1961;Maroto, Pérez-Muñuzuri & Romero-Cano 2007).The typical footprint of these convection cells indicates the presence of one or more high-temperature regions in their interior (Wissink & Herlina 2016), and shows that individual cells are separated by narrow regions of low temperature.Various theoretical and experimental studies have been carried out to investigate the onset of the underlying instability and the resulting pattern formation (see e.g.Pearson 1958;Chandrasekhar 1961;Bodenschatz, Pesch & Ahlers 2000;Nepomnyashchy, Legros & Simanovskii 2006).
Most of these studies were motivated by chemical engineering applications, and employed a relatively thin layer of fluid bounded from below by a solid wall (no-slip boundary).These studies showed that in a thin layer of fluid such horizontal temperature-gradient-induced Marangoni forces act together with buoyancy forces to promote mixing of cool surface water with warmer water from the (upper) bulk (e.g.Nield 1964;Toussaint et al. 2008).
Studies of this convective instability in non-shallow water bodies with a free surface, such as the work carried out by Spangenberg & Rowland (1961), Davenport (1972), Tan & Thorpe (1996) and Flack, Saylor & Smith (2001), are relatively few.Tan & Thorpe (1996) reported that surface-tension-driven convection usually only dominates in relatively thin layers of fluid (up to about 4 mm in thickness), while for layers thicker than 10 mm (as usually encountered in environmental engineering problems) buoyancy-driven convection dominates.The former would imply that, especially in a developing convective instability, Marangoni forces cannot be ignored as the maximum thickness of the thermal boundary layer (which is the most relevant thermal length scale) tends to be much smaller than 4 mm, irrespective of the water depth.
In previous three-dimensional unsteady numerical simulations of buoyancy-driven heat transfer across the air-water surface in relatively deep bodies of water, the physical effect of surface cooling was modelled by prescribing either a constant temperature (Wissink & Herlina 2016) or a constant heat flux (Handler et al. 1999;Fredriksson et al. 2016).In contrast to simulations with a constant surface temperature, in constant-heat-flux simulations variations in surface temperature typically occur which generate Marangoni forces.Despite their potentially significant contribution to the heat and mass transfer coefficients, these forces were usually ignored.
It should be noted that Marangoni forces can also be generated by local concentration variations of a component in a multicomponent fluid.An example of this can be found in evaporating binary droplets that contain water mixed with another fluid of different volatility and density.Here, evaporation of the more volatile component induces density differences as well as surface tension gradients.As a result, a competition between buoyant and Marangoni convection is obtained causing a variety of convection patterns (see Diddens, Li & Lohse (2021) and references therein).Another example are the surfactant-concentration-induced Marangoni forces studied by, for example, Khakpour, Shen & Yue (2011) and Wissink et al. (2017).In the latter study, a wide range of contamination levels was considered.It was shown that the interfacial mass transfer reduces monotonically with surfactant concentration, and a model was derived that predicts the interfacial mass transfer from a snapshot of the surfactant concentration distribution and turbulent flow characteristics of the water.
While surfactant-concentration-induced Marangoni forces were found to significantly reduce near-surface turbulence (and hence interfacial heat and mass transfer), the effect of the surface-temperature-gradient-induced Marangoni forces (studied here) would be to enhance the near-surface convection, thereby working in accord with the buoyant instability.To carry out a detailed investigation of the importance of Marangoni forces and their interaction with buoyancy forces in the initial development of the Rayleigh-Bénard-Marangoni (RBM) instability, fully resolved, three-dimensional, time-accurate numerical simulations were performed.In the simulations, evaporative cooling was modelled by a constant heat flux.Care was taken to ensure that the depth of the computational domain was sufficient to accurately represent the early stages of the development of the RBM instability in a deep body of water.The size of the computational domain as well as the interfacial heat flux were chosen in accordance with the typical set-up of laboratory experiments, such as those of Foster (1965).Direct numerical simulations (DNS) were performed for a wide range of Rayleigh (Ra L ) and Marangoni (Ma L ) numbers, including cases with Ra L = 0 and Ma L = 0, representing purely surface-tension-driven convection and purely buoyancy-driven convection, respectively.

Computational aspects
As mentioned above, this study is motivated by the problem of heat transfer across a water surface.In the range of temperatures considered, the water density varies approximately linearly with temperature.Because of the very small changes in density, the Boussinesq approximation is used to describe the flow induced by the unstable stratification near the surface.The non-dimensional Navier-Stokes equations are solved to determine the fluid motion.Following Balachandar (1992) the non-dimensionalisation uses a reference length scale L and velocity scale U = κ/L, where κ is the thermal diffusivity.The resulting continuity equation reads and the momentum equations are given by where x * 1 , x * 2 are the horizontal (x, y) directions, x * 3 is the vertical (z) direction, u * 1 , u * pressure, t * denotes time and Pr = ν/κ is the Prandtl number corresponding to the ratio of the momentum and thermal diffusivities.Note that the superscript ' * ' denotes non-dimensionalisation using L and U.The last term on the right-hand side of (2.2) represents the buoyancy force in the z direction, where δ i3 is the Kronecker delta, is the non-dimensional temperature, in which T B,0 is the initial temperature of the bulk and q = (∂T/∂z)| S is the constant temperature gradient at the water surface, and is the modified macroscale Rayleigh number in which g = −9.81m s −2 is the gravitational acceleration and α is the thermal expansion factor.The non-dimensional transport equation for the temperature T * is given by (2.5) Based on the model presented in Shen, Yue & Triantafyllou (2004), Wissink et al. (2017) produced a model for the Marangoni forces induced by gradients in the surfactant distribution using the ratio of the Marangoni and capillary numbers.Here, a similar model, which is equivalent to the model used by Pearson (1958), is employed for the Marangoni forces induced by gradients in the surface temperature, where the (modified) Marangoni number replaces the ratio of the Marangoni and capillary numbers, using σ and μ to denote the surface tension and dynamic viscosity, respectively.Note that for water at about 293.15 K, dσ/dT = −0.000151N m −1 K −1 , κ = 0.143 × 10 −6 m 2 s −1 , α = 0.000207 K −1 , Pr = 7.With the exception of dσ/dT (which was varied between −0.000151 and 0 N m −1 K −1 ), the macro Rayleigh and Marangoni numbers were based on the aforementioned values, together with the (arbitrarily chosen) fixed macro length scale L = 0.01 m.The range of parameters used can be seen in Appendix A. In all cases the size (in particular the depth) of the computational domain was sufficiently large so as to not affect the development of the instability.Hence, a priori a fixed macro length scale L was chosen to ensure that both Ra L and Ma L also become depth-independent.This is crucial when comparing different simulations with possibly different domain depths.The true characteristic length scale of a developing instability is, of course, time-dependent and related to the thermal boundary layer thickness, which is determined a posteriori.

Numerical method
To solve the governing equations presented above, the KCFlo solver described in Kubrak et al. (2013) was employed.The solver was used previously to study the influence of low-intensity (Herlina & Wissink 2014) and high-intensity (Herlina & Wissink 2016) isotropic turbulence on interfacial mass transfer as well as the influence of surfactants (Herlina & Wissink 2016;Wissink et al. 2017) and of a developing buoyant instability, where the interfacial cooling was modelled using a fixed temperature at the surface (Wissink & Herlina 2016).When required, the solver allows the usage of a dual, refined mesh to ensure that the steep gradients occurring in low-diffusivity scalar transport are fully resolved.For the discretisation of scalar convection, the fifth-order-accurate WENO scheme of (Liu, Osher & Chan 1994) is used, while for scalar diffusion a fourth-order-accurate central method is employed.The time integration of the scalar equations is performed using a three-stage Runge-Kutta scheme.
The discretisation of the three-dimensional incompressible Navier-Stokes equation is carried out using a fourth-order central kinetic energy conserving discretisation (Wissink 2004) of the convective terms, combined with a fourth-order-accurate discretisation of the diffusion.The second-order Adams-Bashforth method was employed for the time integration of the velocity field.For the spatial discretisation, the variables were arranged on a staggered Cartesian non-uniform mesh, where the scalars (pressure, temperature and concentration) are all defined in the middle of the grid cells while the velocities are defined in the middle of the cell faces.
After substituting the momentum equations into the continuity equation, a Poisson equation for the pressure was obtained which is solved using a conjugate gradient solver with simple diagonal preconditioning.The code was parallelised by splitting the computational domain into blocks of equal size, each allocated to its own processing core.Communication between the resulting processes was done by employing the standard message passing interface protocol.

Overview of simulations
A schematic of the computational domain is shown in figure 1.In all simulations, the domain was discretised using 200 × 200 × 252 grid points in the x, y and z directions, respectively.Further simulation parameters can be found in table 1.The mesh was stretched in the z direction in order to obtain a finer resolution near the interface with a node distribution that reads where n z is the number of nodes in the z direction and z(n z ) − z(0) = L z is the vertical extent of the computational domain.The stretching is controlled by the parameter ψ, which is set to ψ = 3 in all simulations.To illustrate the quality of the mesh, a grid refinement study is presented in Appendix B.
The boundary conditions used in the simulations are illustrated in figure 1.The much larger horizontal extent, which is typically encountered in natural water bodies, was modelled by employing periodic boundary conditions in the horizontal directions for Table 1.Overview of simulation parameters.The macroscale Ra L and Ma L are defined in (2.4) and (2.7), respectively.Note that the Ra L = 0 cases were obtained by setting α = 0.In total, 43 simulations were performed (cf.Appendix A).
all variables.To avoid fluid leaving the computational domain, at the top and bottom free-slip boundary conditions were used for the velocities.The non-dimensional heat flux at the top was fixed at ∂T * /∂z * = −1, while at the bottom the temperature was set at T * = 0. Initially, the velocity field was set to zero and the non-dimensional temperature in (2.3) was initialised by the analytical solution of the unsteady heat transfer in a semi-infinite domain with a constant heat flux given by (Carslaw & Jaeger 1959) The simulations started at t * = 1.4 × 10 −2 , corresponding to t = t 0 = 10 s, at which instance random disturbances that were uniformly distributed between 0 and 0.0257 were added to T * in (2.3) to trigger the RBM instability.Simulations were performed for Rayleigh numbers varying between Ra L = 0 and 212 000 and Marangoni numbers ranging from Ma L = 0 to 8150 (see table 1).
The bulk temperature T B is evaluated at z = 0.2L z , which was deemed to be sufficiently far away from the surface so as to not be influenced by the developing instability.Because snapshots were stored at intervals of 0.25 s, for presentation purposes it was decided to display the time in seconds and not in non-dimensional form.Blue and white colours correspond to relatively cold and warm fluid, respectively.Forces F g and F σ are gravity-and surface-tension-induced forces, respectively.Note that (a) and (b) are scaled differently, t 1 corresponds to the time at which the thermal boundary layer thickness is maximum and t 2 to the time at which the surface kinetic energy reaches its first peak.

Evolution of temperature field
At the surface, evaporative cooling (negative temperature gradient) results in the development of a thermal boundary layer, consisting of a cold, relatively heavy fluid on top of a relatively warm bulk fluid.Initially, the instability is triggered by adding random disturbances to the temperature field (cf.§ 2.2), resulting in small local pockets of relatively cold fluid.
In the case of purely buoyancy-driven instabilities, gravity causes these relatively heavy pockets of cold fluids to slightly sink, thereby replacing warmer fluid from the lower boundary layer.This process results in a convergent flow of water along the surface that is subjected to evaporative cooling.A schematic, showing the dynamics of this buoyant instability in a vertical cross-section through two neighbouring convection cells, can be seen in figure 2(a).The continuous accumulation of cooled surface water results in the development of a large pocket of cold, heavy water near the surface.At the time of this snapshot (t = t 1 -when the thermal boundary layer is maximum) the amount of cold water is sufficient to overcome the stabilising effect of diffusion as the buoyancy (gravity) forces (F g ) have become strong enough to pull the cold fluid away from the boundary layer down into the bulk.The snapshot in the right-hand panel of figure 2(a) shows the situation at the later time of t = t 2 -when the surface kinetic energy reaches its first maximum.Here, the cold fluid has been accelerated downwards and transformed into a cold water plume deeply penetrating into the bulk.
A schematic, showing the dynamics of the purely Marangoni-driven instability, can be seen in figure 2(b).The pockets of cold and warm fluid near the surface (generated by the addition of random disturbances to the temperature field) have relatively high and low surface tensions, respectively.This difference in surface tension results in Marangoni forces F σ (that act parallel to the surface) pushing surface water from the warm pockets towards the cold pockets.While travelling along the surface the relatively warm water is subjected to evaporative cooling.As a result, convection cells are formed due to the continuous supply of warm surface water from the interior that is cooled as it is transported towards the cell edges.At the edges, where two or more convection cells meet, opposing streams of cooled water result in a downward momentum, whereby cold fluid is pushed down into the upper bulk.In time, as the convection cells develop, the temperature gradients at the cell edges become stronger and stronger.At time t = t 1 the local Marangoni forces have grown sufficiently strong to overcome thermal diffusion.As a result, a significant amount of cold water is accumulated close to the surface at the edges of convection cells.In between t = t 1 and t = t 2 the convection cells rapidly mature, as the Marangoni forces further increase in strength, and become very effective in pushing cooled surface water into the upper bulk.
Note that the schematics shown in figure 2 illustrate the evolution of the instabilities driven by buoyancy and surface tension (Marangoni) at two important time instances t = t 1 and t = t 2 (see also § 4).
When the surface temperature is non-uniform, the two instabilities described above will typically act simultaneously.The effect of a changing Marangoni number Ma L on the evolution of the instantaneous surface temperature at a fixed Rayleigh number of Ra L = 11 000 is illustrated in figure 3.In time, at each Ma L it can be seen that convection cells rapidly mature as they become more clearly defined by an increase in the diameter of the areas with enhanced temperature resulting in a narrowing of the low-temperature regions separating the cells.
Below, it is shown that with increasing Marangoni number the rate at which these convection cells develop enhances.At t = t 0 = 10 s, in all simulations identical random disturbances were introduced to the temperature field.As a result, at t = 10.25 s the surface temperature fields were found to be very similar (cf.figure 3a-c).
When comparing temperature fields at t = 30 s, a clear difference can be seen in the state of development of the convection cells (cf.figure 3d-f ).Even though the vertical surface heat flux in all three simulations is the same, an earlier development (maturation) of well-defined convection cells at the larger Marangoni numbers was observed.This is especially clear for Ma L = 2700, where the convection cells are quite well defined when compared with the more diffusive temperature patterns obtained for the smaller Ma L .The observed enhanced maturation rate is explained by the promotion of RBM convection by the increasing Marangoni forces that originate from horizontal gradients of the surface temperature.
The remainder of figure 3 shows the further development of the convection cells at the two time instances t = t 1 (figure 3g-i) and t = t 2 (figure 3j-l); see also figure 2. It can be seen that the time needed to achieve both states reduces with increasing Ma L .This reflects the enhanced maturation rate induced by the increasing Marangoni forces (as observed at t = 30 s above).The effect of Marangoni forces is non-trivial, as evidenced by the significant reduction of the time at which the convection cells near the surface become well developed (mature), from t 2 = 143.75s at Ma L = 0 to t 2 = 53 s at Ma L = 2700.Furthermore, the size of the cells was found to become significantly smaller with increasing Ma L .A more quantitative discussion of the above observations is presented in the remainder of the paper.x/L where z S is the z coordinate at the surface and T B is the horizontally averaged temperature in the bulk.The onset of the RBM instability is identified by the time t = t 1 at which δ is maximum, where • denotes horizontal averaging.This time instance is an approximation of the time at which δ begins to deviate from the analytical solution for the purely diffusive case, which happens shortly before thermal plumes start plunging down.Note that the onset of the instability comes about much earlier than t = t 1 but is initially hidden by thermal diffusion.The development of δ for cases with Ra L = 11 000 and varying Ma L is shown in figure 4.This figure illustrates the importance of Marangoni forces in promoting the RBM instability as it shows that at a constant Rayleigh number, t 1 decreases significantly with increasing Ma L .
4.2.Surface temperature and fluctuating kinetic energy Figure 5(a,b) shows the temporal evolution of temperature fluctuations (T * rms = (T * − T * ) 2 ) and the fluctuating kinetic energy (K = ( u u + v v + w w )/2) evaluated at the surface for various Ma L at Ra L = 11 000.Note that the various values for Ma L were obtained by varying q dσ/dT.As mentioned above, the instability was triggered by adding random disturbances to the temperature field, which subsequently induced disturbances in the velocity field through buoyancy and/or Marangoni forces.In general, three phases were observed: (i) The transient phase, which is characterised by the diffusive damping of T * rms,S (note that the subscript 'S' indicates values at the surface).For Ma L > 0, there is an immediate connection between temperature and velocity gradients at the surface as prescribed by the model given in (2.6).This is reflected by the initial (instantaneous) jump in K S at the start of the simulation, and the subsequent diffusive damping of T * rms,S , leading to the damping of K S .With increasing Ma L , the initial damping of T * rms,S (and K S ) was found to significantly reduce.In contrast, in the absence of Marangoni forces (Ma L = 0), there is only an indirect connection (through buoyancy) between temperature gradients and velocity gradients at the surface.Consequently, despite the initial reduction in T * rms,S , the surface kinetic energy fluctuation was found to increase, with a growth rate that gradually reduced in time.This reduction in the growth rate comes to an end at approximately the same time as T * rms,S reaches its minimum.Note that for Ma L > 0, this buoyancy-induced gradual increase in K S is masked due to the aforementioned instantaneous jump in K S .(ii) Growth phase.After the transient phase, a simultaneous increase of T * rms,S and K S can be observed, indicating an increased growth in the buoyant and/or Marangoni instability.The growth rate of T * rms,S (and K S ) can be seen to strongly depend on the magnitude of Ma L .For instance, for Ma L = 0 the maximum T * rms,S and K S is achieved only after t ≈ 102 s, while for Ma L = 15 700, it is already achieved at t ≈ 18 s.(iii) The end phase begins after both T * rms,S and K S almost simultaneously reached their first (local) peak.At this stage, a quasi-steady state is not yet necessarily achieved, as for instance the number of convection cells and their topology may still change in time.
Note that in all simulations, the maximum thermal boundary layer thickness 4) is reached somewhat before the beginning of the end phase.Hereafter, in each simulation, the first local peak in K S is denoted by K S2 , while K S1 denotes K S at t = t 1 .length scale of the surface temperature (Λ S , defined in Appendix C).By comparing figures 3 and 6, it can be seen that Λ S (for at least up to t = t 1 ) is closely correlated to the average distance between two neighbouring convection cell centres:

Surface characteristic length scales
where N = N(t) is the number of cells in the computational domain at time t.Individual cells were identified by isolated patches at the surface with a surface divergence b ≥ b + rms , where b + rms is the corresponding root mean square of all b ≥ 0. Figure 6 shows that the surface characteristic length scale Λ S initially increases at least until t = t 1 , which corresponds to the time when the horizontally averaged thermal boundary layer thickness δ reaches its maximum (cf.figure 4).This growth in Λ S , and thus the reduction in N, indicates the formation of larger convection cells, for instance, due to the merging of two or more neighbouring cells, as was observed in figure 3.In the range t 1 < t < t 2 , however, the correlation between Λ S and λ S depends on which of the two forces (buoyancy or Marangoni) dominates.
For the smaller Ma L simulations, Λ S can be seen to somewhat reduce until at least t = t 2 .This small reduction is not necessarily only associated with the emergence of new convection cells, but possibly also with the enhanced maturation rate of initially very weak cells between t 1 and t 2 .In this phase, the dominant buoyancy forces in these simulations result in a significant acceleration of the falling plumes.This results in an enhanced updraft of warm fluid that leads to both a thinning of the boundary layer and an increase in footprint of well-defined convection cells at the surface, separated by relatively narrow sheets of falling cold fluid (cf.figure 3).
In contrast to the smaller Ma L simulations, for the larger Ma L simulations, where the Marangoni effect begins to dominate, the correlation between Λ S and λ S remains strong between t 1 and t 2 .In this period, Λ S significantly increases, which agrees with the observed reduction in the number of convection cells and an apparent increase in cell size (cf.figure 3i,l).
Note that as shown in figure 6, the growth rate of Λ S initially decreases with increasing Ma L , while simultaneously the maturation rate of the convection cells significantly enhances (cf.figure 3).As a result, significantly more (and smaller) convection cells are formed as Ma L increases.The associated decrease in maturation time is reflected in reductions in both t 1 as well as t 2 − t 1 .This can be explained as follows.In the buoyancy-dominated cases, significant time is needed to allow for a sufficient amount of cold water to accumulate near the surface so that cold plumes can be released.On the other hand, while Rayleigh forces always act in the vertical direction, Marangoni forces act parallel to the surface.Therefore, in the Marangoni-dominated cases, there is no need to achieve a significant thermal boundary layer thickness in order to obtain a well-developed instability.Given horizontal gradients in the surface temperature, Marangoni forces continuously push surface water to the edges of the convection cells, thereby sucking relatively warm, near-surface water upwards to the centre of these cells.The latter results in an increased temperature in the interior of the convection cells, which enhances their maturation rate by further strengthening the horizontal surface temperature gradients.

Scaling laws
In § 4, the effect of Marangoni forces on the maximum thermal boundary layer thickness (δ 1 ) and the surface kinetic energy fluctuation peaks (K S2 ) at a fixed Rayleigh number of Ra L = 11 000 was discussed.While δ 1 significantly decreases with increasing Ma L , K S2 was found to increase.Similar results were also observed for other fixed Ra L as shown in figure 7. It can be seen that with increasing Ma L all curves tend to approach the Ra L = 0 curves, indicating that the Marangoni forces begin to dominate the evolution of the instability.These Ra L = 0 curves show scalings of δ 1 ∝ Ma −0.50 L (figure 7a) and K S2 ∝ Ma L (figure 7b).Note that the Ra L = 0 cases were obtained by setting the thermal expansion factor α to zero, corresponding to water at a temperature of about 4 • C.
Similarly, for a number of fixed Ma L , the effect of Ra L on δ 1 and K S2 is summarised in figures 8(a) and 8(b), respectively.Here, it can be seen that for large Ra L , all curves tend to approach the Ma L = 0 curves, from which we can conclude that the Marangoni forces become negligible compared with the buoyancy forces.Note that the two Ma L = 0 curves indicate scalings of δ 1 ∝ Ra −0.25 L and K S2 ∝ Ra 0.5 L , respectively.To explain the scaling of δ 1 with Ma L at Ra L = 0, observed in figure 7 seen in figure 9(a) that for the entire range of Ma L investigated, both Ma δ 1 and Ma Λ S1 only fluctuate across relatively small ranges of 200 < Ma δ 1 < 250 and 410 < Ma Λ S1 < 440, respectively.Next, the Nusselt number where is obtained.By assuming that Ma δ 1 is constant (as somewhat indicated in figure 9a), is obtained, which is in excellent agreement with the scaling of δ 1 /L as seen in figure 7(a).
In figure 9(b), it can be seen that Ra δ 1 and Ra Λ S1 are bounded between 450 < Ra δ 1 < 550 and 10 300 < Ra Λ S1 < 15 700, respectively.Using a similar derivation as above, i.e.Nu δ 1 = L/δ 1 = [Ra L /Ra δ 1 ] 0.25 and assuming that Ra δ 1 is constant, it can be shown that δ 1 /L ∝ Ra −0.25 L , which is in close agreement with the slope of the Ma L = 0 curve in figure 8(a).Note that as Ra δ 1 is based on the maximum thermal boundary layer thickness that is reached slightly before the release of cold water plumes, it is thereby strongly related to the critical Rayleigh number below which the flow would be stable.
To discuss the scaling of the surface kinetic energy fluctuation peak K S2 with Ma L for Ra L = 0, the surface Reynolds number is introduced, which is based on the square root of the surface kinetic energy and the surface characteristic length scale evaluated at t = t 1 .In figure 10 Re S remains relatively close to 1.Because both Re S and ν are (approximately) constant, it follows that (5.5) Assuming that also Ma Λ S1 is constant (cf.figure 9a), from Nu Λ S1 = L/Λ S1 = (Ma L /Ma Λ S1 ) 0.5 , it immediately follows that √ Ma L ∝ 1/Λ S1 and, hence, K S1 ∝ Ma L .Analogously, for Ma L = 0, Re S was also found to be approximately constant (cf.figure 10b), such that (5.5) remains valid.Furthermore, by assuming that Ra Λ S1 is constant (as indicated in figure 9b), it follows that K S1 ∝ Ra 0.5 L .Note that K S1 was found to be approximately proportional to the surface kinetic energy peaks K S2 (not shown here).Hence, it can be concluded that K S2 ∝ Ma L and K S2 ∝ Ra 0.5 L , which are in agreement with the observations made in figures 7(b) and 8(b), respectively.

Relative importance of buoyancy and Marangoni forces
The Bond number expresses the relative importance of buoyancy and Marangoni forces in promoting interfacial flow, which results in a vertical heat exchange between the surface and the bulk.For sufficiently large Ma L , it was shown that Ma Λ S1 is approximately constant (cf.figure 9a) and Richardson number is proportional to Bo Λ S1 .This is evidenced in figure 11, which illustrates the domination of Marangoni forces approximately up to Ri T = 300 and Bo Λ S1 = 5.For larger Ri T , buoyancy effects become increasingly important as shown by the continuously increasing growth rate of Bo Λ S1 , indicating the presence of an asymptote at Ri T ≈ 1000.The presence of this asymptote is a consequence of the earlier observation ( § 5) that Ra Λ S1 is bounded by a relatively small interval about 1.3 × 10 4 for buoyancy-dominated flows, which combined with (6.2) implies that Ri T ≈ const. (6.3) Whether the RBM instability is Marangoni-or buoyancy-dominated can be detected, for example, by studying the evolution of the surface characteristic length scale (cf.figure 12).The dominating instability can be determined by the ratio r Λ = Λ S1 /Λ S2 .The smaller the value of r Λ < 1, the more the Marangoni forces dominate, while buoyancy forces progressively dominate with increasing r Λ > 1.Furthermore, r Λ ≈ 1 corresponds to the situation where buoyancy and Marangoni forces are in equilibrium, which is achieved for 300 < Ri T < 600.The lower bound of this range identifies the smallest Ri T at which Bo Λ S1 is no longer proportional to Ri T (cf.figure 11).Differences in the statistical properties between developing Rayleigh-Bénard (RB)and Bénard-Marangoni (BM)-dominated instabilities can also be seen in the surface temperature fluctuations (that are a measure of the convection cell strength) at t = t 2 (cf.figure 13a).For both purely Marangoni-and purely buoyancy-driven instabilities, a positive linear correlation T * rms,S = c T (T * B − T * S ) was obtained, with a constant of proportionality c T that was significantly larger in the purely buoyancy case.Furthermore, with increasing Ra L it was found that both (T * B − T * S ) and T * rms,S tend to become smaller.A similar observation was made for increasing Ma L at fixed Ra L , which is especially clearly visible in the figure at Ra L = 21 200 and at Ra L = 33 000.To further clarify these trends, in Appendix D detailed examples are presented for the case with varying Ra L at a fixed Ma L = 550, and for the case with varying Ma L at a fixed Ra L = 21 200.It can be calculated from figure 13(a) that T * rms,S /(T * B − T * S ) is approximately 0.23 for the Marangoni-dominated (low-Bo Λ S1 ) cases, while for the buoyancy-dominated cases it would be about 0.37.This is confirmed in figure 13 It can be seen that all results approximately collapse on one single curve, and that for Bo Λ S1 up to ≈1, the instability is Marangoni-dominated, while for Bo Λ S1 > 20 it is buoyancy-dominated.
A possible explanation of the significant reduction in T * rms,S for the Marangoni-dominated instabilities is provided in figure 14(a).This figure compares the horizontally averaged temperature profiles for zero buoyancy and zero surface tension obtained at t = t 2 with the corresponding case with realistic values for Ra L and Ma L for water at 293.15 K.Both profiles with Ma L = 8150 show an almost uniform local vertical temperature distribution with T * ≈ T * = −0.11(which is only slightly larger than T * S ≈ −0.13) evidencing extensive mixing.This mixing close to the surface is due to the horizontal acceleration of fluid at the surface by Marangoni forces.Where two opposing streams meet, (cold) surface water is pushed down into the boundary layer (even in the absence of buoyancy), while simultaneously warmer water from the boundary layer moves to the surface.The continuation of the above process then results in the accumulation of cold water plumes that remain near the surface due to the lack of buoyancy, explaining the aforementioned extensive mixing immediately underneath the surface.In contrast, due to buoyancy forces, in the purely buoyancy-driven instabilities the bulk of the mixing takes place significantly farther away from the surface (cf.figure 14a).As a result, both the temperature differences between the surface and the region of enhanced mixing, as well as the normalised T * rms,S are significantly larger (cf.figure 14b).Note that figure 14(a) illustrates that in realistic cases the Marangoni forces completely dominate the developing instability and, as a result, the time for the instability to develop effective vertical mixing (at t = t 2 ) is significantly smaller than without Marangoni forces.Because of this much smaller time t 2 , the near-surface water had been subjected to cooling for a significantly shorter period which explains the warmer temperature in the upper bulk compared with the buoyancy-dominated case.
The Marangoni forces, which are generated by gradients in the surface temperature, act to horizontally accelerate fluid, thereby increasing the turbulent kinetic energy K.As a result, in the Marangoni-dominated cases, K is maximum at the surface (cf.figure 15).In the absence of any further forcing below the surface, K very quickly reduces to zero in the downward direction.In contrast, in the pure buoyancy case the movement of the fluid at the surface is entirely generated by gravitational sources pulling heavy cold water downwards.Here, the local K will increase as soon as plumes of cold water start falling and are accelerated downwards, such that at t = t 2 the maximum K is not obtained at the surface but in the upper bulk.
To study the dynamics of the development of the purely buoyancy-and the purely Marangoni-driven instabilities, figure 16 shows instantaneous snapshots of the fluctuating temperatures (T * = T * − T * ) at t = t 1 and t = t 2 for Ra L = 70 800, Ma L = 0 and for Ra L = 0, Ma L = 2700.Note that these cases were selected because of their similar t 1 (time at which the maximum boundary layer is reached) and t 2 (time at which the first peak in the fluctuating kinetic energy is obtained).While Marangoni forces are generated by temperature gradients at the water surface, pushing cold water down at locations where convection cells meet, buoyancy forces work below the surface by pulling cold water into the bulk.The result of this can be seen in the pure buoyancy case, where at t = t 1 cold water has accumulated into very slow-moving primordial plumes that, in time, further develop and accelerate giving rise to the falling plumes observed at t = t 2 .In contrast, in the pure Marangoni case there is an increase in the variation in size of the relatively small convection cells, as can be seen at t = t 1 , while at t = t 2 the number of convection cells reduces.At both times, the cold water pushed down by Marangoni forces is visible in the y-z plane as sprout-like structures that remain very close to the surface.The penetration depth of these structures approximately scales with the average size of the convection cell footprints at the surface and, hence, with the characteristic length scale Λ S .Compared with the Marangoni-driven instability, at similar t = t 1 , the buoyancy-driven instability typically generates significantly larger convection cells, and hence a larger Λ S1 , while also the larger temperature range indicates an increased value of T * .The latter is in agreement with results shown in figure 13(a), while the former is in agreement with the results shown in figure 17.Note that the smaller structures underneath the surface of the pure Marangoni instability typically have a short turnaround time, resulting in a more intense mixing.

Comparison with previous studies
In this section, the present DNS results -investigating the three-dimensional development of the RBM instability -are compared with those of previous studies, which mostly focused on either the RB or BM instability in thin layers of fluid.In such studies, the fluid depth d is commonly used as a fixed characteristic length scale.However, in non-shallow simulations, when studying the developing instability, d should be replaced by some measure of the thermal boundary layer thickness.When using the latter length scale, it is shown that the present results (obtained in a deep domain) agree well with those of previous studies.
Note that in all simulations, initially a fixed macro length scale of L = 0.01 m was used in combination with varying q, α and/or dσ/dT, in order to obtain a range of Ra L and Ma L values.The corresponding Crispation (Cr = ρνκ/σ δ 1 ) and Galileo (Ga = gH 3 /νκ) numbers ranged from 1.22 × 10 −7 to 1.95 × 10 −6 and from 7 × 10 4 to 2.85 × 10 8 , respectively.As in the present simulations Cr 1 and Ga 1, the flat surface assumption employed here is justified.depends non-monotonically on the Biot number (Bi = hd/κ c , where h is the heat transfer coefficient in W m −2 K −1 and κ c is the thermal conductivity) and weakly on the Prandtl number.For Bi = 1, they obtained Ma c ≈ 240, which is within the range 200 < Ma δ 1 < 250 obtained here at Ra L = 0 (cf.figure 9a).A similar result was obtained in the stability analysis of Smith & Davis (1983), who found that Ma c ≈ 250 at Pr = 7 for a fluid layer subjected to a constant horizontal temperature gradient (with return flow conditions).
It is also interesting to compare the characteristic size of the convection cells at the critical time t = t 1 .Using (4.4) and the results shown for Ra L = 0 in figure 17, it can be derived that λ S /δ 1 , which is associated with the non-dimensional critical wavelength, is about 2.8.This is in agreement with Smith & Davis (1983), who found that λ S /d ≈ 2.5, as well as with the experimental results 2.1 ≤ λ S /d ≤ 2.5 obtained by Bénard (cf.Pearson 1958). Furthermore, recently, Toussaint et al. (2008) found a value of λ S /d ≈ 2.6 in their experiments studying instabilities induced by evaporation during the drying of a thin layer of polymer solution.

Purely buoyancy-driven case
In this section, the critical Rayleigh numbers Ra δ 1 obtained in the present simulations are compared with the critical Rayleigh numbers Ra c presented in the literature.Comparable to Ra δ 1 , Ra c is calculated using some measure of the effective thermal boundary layer thickness as length scale.Sparrow, Goldstein & Jonsson (1964) conducted a stability analysis for linear and nonlinear temperature profiles and various upper and lower boundary conditions.The largest Ra c was found for the upper rigid surface boundary condition with fixed temperatures at both upper and lower surfaces.It was indicated that Ra c reduces (i) when replacing the rigid by a free upper surface boundary condition, (ii) when replacing the fixed upper surface temperature by a constant heat flux and (iii) when using a nonlinear instead of a linear temperature profile.For nonlinear temperature profiles and rigid, fixed temperature upper and lower surfaces, Ra c was found to vary between 560 and 595.A further reduction in Ra c for a free upper surface with a constant heat flux is expected, and would be in line with the present results, where Ra δ 1 was found to vary between 450 and 550.
The aforementioned range is also in accordance with the result Ra c = 500 obtained by Davenport (1972) for their n-octanol, deep-pool, free-surface case, considering a nonlinear temperature profile combined with a constant heat flux (cf.figure VII-2 in Davenport 1972).To compare the present results with the experiments of Foster (1965) and Spangenberg & Rowland (1961), the following modified Rayleigh numbers are introduced: present results obtained at Ma L = 0. Note that for the present DNS, t c was approximated by t = t 1 and β by (7.4) where t d = 10 s is the time at which random disturbances were added to the temperature field (cf.§ 2.2).It can be seen that when using Ra γ , the present data agree well with those of Foster (1965) and Spangenberg & Rowland (1961).
The close agreement between the present DNS and the experimental data might be somewhat coincidental, as it is unlikely that the boundary and initial conditions are identical.For example, the chosen initial random disturbance level might or might not be a close match to the experimental condition.At the same time, even in well-controlled laboratory conditions, it is nearly impossible to maintain a perfect Ma L = 0 condition, which implies a surface without any contamination.As shown in Herlina & Wissink (2016) and Wissink et al. (2017), even a slightly contaminated surface, due to the presence of, for example, dust, tracer particles or surfactants, will lead to a significant damping of near-surface turbulence.Hence, in the laboratory, horizontal surface-temperature-gradient-induced Marangoni forces (promoting near-surface turbulence) are likely to be more or less negated by contamination-induced Marangoni forces.
It is worth mentioning that the Ra β L values in the experiments of Foster (1965) are contained within the range of Ra β L values used in the present DNS.Also, Ra γ was found to be similar to Ra δ1 , demonstrating that (7.3) provides a good estimation for δ 1 (see also figure 4).Nield (1964) proposed the existence of a tight but imperfect coupling between the buoyancy and Marangoni instabilities, which can be expressed by are identified as the minimum of Ra δ 1 at Ma L = 0 and Ma δ 1 at Ra L = 0, respectively (cf.figure 9a,b).It should be noted (as also mentioned previously) that the values of Ra δ 1 and Ma δ 1 are affected by the prescribed initial conditions used in the simulations, especially the initial level of random disturbances added to the temperature field combined with the initial thermal boundary layer thickness (which were kept constant in the simulations).

Coupled buoyancy-driven and Marangoni-driven instabilities
Note that the study of the specific effects of the initial conditions is beyond the scope of this paper.Nevertheless, as can be seen in figure 19, the tight coupling between Rayleigh and Marangoni instabilities, such as found by Nield (1964), also approximately holds in our transient simulations.Hence, also in the present (deep-domain) developing case, it was confirmed that both instabilities enhance one another.

Conclusions
This paper presents a fully resolved, three-dimensional numerical study of the initial development of the RBM instability, i.e. up to the time when the surface kinetic energy reaches its first maximum (at t = t 2 ).It is important to note that the mechanisms of the surface-tension-gradient-driven (BM) and the buoyancy-driven (RB) instabilities are different, though they do enhance one another.The latter is evidenced by the increased maturation rate of convection cells obtained at fixed Rayleigh number, Ra L , when increasing the Marangoni number, Ma L .Similarly, such an increased maturation rate was also observed at fixed Ma L when increasing Ra L .Three phases could be identified in the initial evolution of the RBM instability: (i) the transient phase, (ii) the growth phase and (iii) the end phase.The transient phase is identified by the diffusive damping of the disturbances added to the temperature field, as reflected by the initial reduction in T * rms,S .In this phase, the immediate connection between surface temperature and velocity gradients (through the prescribed boundary condition) was found to result in a significant jump in K S at the start of the simulation only in the presence of Marangoni forces.In the growth phase, both T * rms,S and K S simultaneously increase, with a growth rate that significantly enhances with increasing Ma L .The end phase begins approximately at t = t 2 when both T * rms,S and K S have reached their first 962 A23-23 Both the RB-and BM-dominated instabilities are characterised by an initial rapid reduction in the number of convection cell centres and an increase in surface characteristic length scale, Λ S .When using (4.4) for the average distance between adjacent convection cells (λ S ), a good correlation was found between λ S and the surface characteristic length scale (Λ S ) during the transient phase (up to t = t 1 ).For t > t 1 , this positive correlation only persists in the BM-dominated cases, where Λ S grows significantly.In contrast, in the RB-dominated cases this correlation is lost, as Λ S tends to decline in the growth phase (between t = t 1 and t = t 2 ).
The relatively wide range of Ra L and Ma L simulated allowed an examination of a number of scaling laws.For fixed Ra L , the following were found: (i) At Ra L = 0 and 550 ≤ Ma L ≤ 8150, Ma δ 1 , Ma Λ S1 and Re S were shown to remain approximately constant, and as a consequence, δ 1 was found to scale with Ma −0.5 L and K S2 /(κ 2 /L 2 ) with Ma L .(ii) For fixed Ra L > 0, with increasing Ma L , δ 1 and K S2 were observed to asymptotically approach their respective Ra L = 0 curve.
When varying Ra L at fixed Ma L , the following were shown: (i) At Ma L = 0 and for 5700 ≤ Ra L ≤ 212 000, Ra δ 1 and Re S were observed to be approximately constant.As a consequence, δ 1 was shown to scale with Ra −0.25 L and K S2 /(κ 2 /L 2 ) with Ma 0.5 L .(ii) For fixed Ma L > 0, with increasing Ra L , δ 1 and K S2 asymptotically approach their respective Ma L = 0 curve.
Note that the 'critical' Rayleigh and Marangoni numbers, 450 < Ra δ 1 < 550 and 200 < Ma δ 1 < 250, were obtained for the specific (fixed) initial conditions (e.g.initial temperature profile and initial disturbance level) employed in the present simulations.Any changes in these initial conditions may affect these parameters.
The RB and BM instabilities were found to be in approximate equilibrium when the ratio r Λ = Λ S1 /Λ S2 ≈ 1. Flows dominated by the BM instability are characterised by a small r Λ < 1, while increasing values of r Λ > 1 indicate that the flow instability becomes more and more RB-dominated.
Positive linear correlations T * rms,S ∝ (T * B − T * S ) were obtained for both purely Marangoni-and purely buoyancy-driven instabilities, where both (T * B − T * S ) and T * rms,S tend to become smaller with increasing Ma L and Ra L , respectively.In general, when both Ra L and Ma L are non-zero, instabilities were found to be Marangoni-dominated for Bond numbers Bo Λ S1 up to ≈1 and buoyancy-dominated for Bo Λ S1 > 20.The significant reduction in T * rms,S for the Marangoni-dominated simulations was attributed to the significant mixing induced immediately underneath the surface, in contrast to the buoyancy-dominated simulations, where the bulk of the mixing was induced much deeper down.The differences in mixing locations are explained by the fact that buoyancy tends to promote the formation of cold water plumes that are pulled down by gravitational forces into the deeper bulk.In contrast, Marangoni forces horizontally accelerate fluid inducing flows along the surface from low-to high-surface-tension areas, associated with relatively warm and cold surface temperatures, respectively.The collision of two or more opposing In one of the simulations, a heat flux of 46.5 W m −2 was combined with dσ/dT = 0.000151 N m −1 K −1 , corresponding to the surface tension gradient for water at ≈293.15 K.As a result, it was found that the Marangoni forces completely dominate the instability to such a degree that the mean temperature profiles up to t = t 2 were very nearly identical irrespective of the presence of buoyancy forces.Consequently, in any numerical study of buoyancy-driven interfacial mass transfer, in which surface cooling is modelled by a constant heat flux, Marangoni forces should be taken into account.(corresponding to the mesh used in the actual simulations) and R2. Figure 20(b) shows the horizontally averaged temperature profile from R1 against z/L at t = 30 s.It can be seen that there are at least nine grid points inside the thermal boundary layer, which is sufficient to fully resolve it.
To further illustrate the quality of the mesh employed, in figure 21 instantaneous profiles of the normalised temperature at t = 30 s are shown.The profiles were extracted at x/L = 2.5 with vertical locations of z/L = 4.6 and 4.975.Given that the perturbation was introduced at t = 12 s and was given sufficient time (until t = 30 s) for the instability to develop, the results obtained on the three meshes are generally in good agreement.This is especially true when comparing the temperature profiles of simulations R1 and R2, which almost perfectly overlap.Hence, the 200 × 200 × 252 mesh was deemed to be sufficiently fine to fully resolve all simulation results presented in this paper.Note that while in the first step the grid was refined by a factor of 2, in the second step a further refinement of about 1.3 in the z direction and 1.5 in the x,y directions was employed.The latter was chosen to save computing time, and was judged to be adequate due to the high order of accuracy of the WENO solver.
Additionally, while post-processing the actual simulation results, further evidence was obtained by showing that the ratio Figure 1.Schematic of computational domain.

Figure 2 .
Figure 2. Schematic showing forces driving (a) purely buoyancy-and (b) purely Marangoni-induced flow patterns.Blue and white colours correspond to relatively cold and warm fluid, respectively.Forces F g and F σ are gravity-and surface-tension-induced forces, respectively.Note that (a) and (b) are scaled differently, t 1 corresponds to the time at which the thermal boundary layer thickness is maximum and t 2 to the time at which the surface kinetic energy reaches its first peak.

Figure 5 .
Figure 5. Evolution of (a) surface temperature fluctuations T * rms,S and (b) surface fluctuating kinetic energy K S for various Ma L at Ra L = 11 000.The cross symbol identifies t 1 , which is the time at which the first peak in δ is reached, and the line colours blue to cyan indicate low-to high-Ma L cases, respectively (cf.figure 4).

Figure 6 .
Figure6.Evolution of the surface characteristic length scale for various Ma L at Ra L = 11 000.The markers × and + correspond to t 1 and t 2 , respectively.

Figure 7 .
Figure 7. Variation in (a) maximum thermal boundary layer thickness δ 1 (and corresponding Nusselt number Nu δ 1 ; cf.(5.2)) and (b) maximum surface fluctuating kinetic energy K S2 with Ma L for various fixed Ra L .

Figure 8 .
Figure 8. Variation in (a) maximum thermal boundary layer thickness δ 1 (and corresponding Nusselt number Nu δ 1 ; cf.(5.2)) and (b) maximum surface fluctuating kinetic energy K S2 with Ra L for various fixed Ma L .

Figure 9 .Figure 10 .
Figure 9. Variation of (a) Ma δ 1 and Ma Λ S1 with Ma L at Ra L = 0 and (b) Ra δ 1 and Ra Λ S1 with Ra L at Ma L = 0.

Figure 12 .Figure 13
Figure 12.Temporal evolution of the surface characteristic length scale at various Ri T .The markers × and + correspond to Λ S1 and Λ S2 , respectively.

Figure 17 .
Figure 17.Variation of the surface characteristic length scale Λ S1 at t = t 1 with δ 1 .Virtually all results are within the envelope indicated by the dashed lines, defined by Λ S1 = 1.4δ 1 and Λ S1 = 2.2δ 1 .

Figure 20
Figure 20.(a) Ratio between grid size and Batchelor scale from R0, R1 and R2 as a function of z/L at t = 30 s.(b) Profile of the horizontally averaged normalised temperature T * from R1 at t = 30 s.

Figure 23 .
Figure 23.Variation of T * rms with (T * B − T * S ): (a) for a variety of Ra L (as shown) at Ma L = 550 and (b) for a variety of Ma L (as shown) at Ra L = 21 200.
Foster (1965)urely diffusive thermal boundary layer thickness, see (4.2), at the observed critical time t = t c .The variation of Ra γ with Ra β L , computed using the β and t c values reported inFoster (1965), is shown in figure18, where it is compared with the Figure 18.Critical Rayleigh numbers Ra γ (and Ra δ 1 ) versus Ra β L .
Even though the convection cells are now well developed, the flow at the surface is still changing and not yet in a quasi-steady state.Note that the growth phase is the main focus of this paper. https://doi.org/10.1017/jfm.2023.263

Table 3 .
Simulations performed for the grid refinement study.Ratio ΔR is defined in (B5).
table 3). Figure 22.Vertical profiles of the maximum ratio between grid size and Batchelor scale evaluated over the entire simulation time.The parameters of each simulation are listed in table 2.