3.1. Flow topology

Figure 2(a) and 2(b) show snapshots of fully developed temperature fields obtained at ${\textit{Ma}}_h=0$ and $8.81\times 10^5$ , respectively. The figure illustrates the effect of Marangoni forces on the heat exchange mechanism between the cooled water surface and the warmer bottom boundary of the computational domain. Temperature isosurfaces at $T^*_1=a_{th}\langle T^*_s \rangle _{x,y}$ and $T^*_2= 2\langle T^*(z=0.5h) \rangle _{x,y} - T^*_1$ (where $a_{th}$ is given in the caption) were used to identify, e.g. the thermal plumes that are characteristic for buoyant instabilities. Irrespective of the presence of Marangoni forces, relatively few large buoyant rising plumes of similar size were found to emerge quasi-periodically from the bottom of the computational domain to the surface. While deeply penetrating plumes of similar size and frequency can also be observed to arise from the top in all simulations, their shape and the mechanism by which they form changes with ${\textit{Ma}}_h$ . In a fully developed RBM flow, heat fluxes at the surface and bottom of the computational domain need to be in equilibrium. Hence, at ${\textit{Ma}}_h=0$ the shape of these large plumes (that are purely buoyancy-generated) mirrors the situation at the bottom of the domain. Contrastingly, in the presence of Marangoni forces the local near-surface flow topology is significantly affected. For instance, at ${\textit{Ma}}_h=8.81\times 10^5$ areas containing a multitude of small-scale structures of relatively cold water emerge (cf. figure 2 b). Due to the entrainment of these Marangoni generated small-scale structures, with increasing ${\textit{Ma}}_h$ the shape of the large, deeply penetrating plumes becomes more and more irregular. The size of these small-scale structures tends to reduce as they become more numerous with increasing ${\textit{Ma}}_h$ . Because of their limited size, the buoyancy force acting on these structures is very small, so that the (horizontally) elongated plume-like structures were found to linger for a significant time at the surface forming a fine, non-uniform mesh.

Figure 3. Near-surface cross-sectional temperature distribution with fluctuating velocity vectors to identify counter-rotating vortices at one of the small plumes. The snapshot is from simulation with ${\textit{Ma}}_h=13.21\times 10^5$ in the plane $x/h=1$ at $t^*=0.00295$ .

As illustrated by the detailed cross-section in figure 3, underneath the joint edges of (typically small) neighbouring convection cells, pairs of counter-rotating vortices can be found lying side by side. The Marangoni-induced convection cells form due to surface-temperature gradients, which induce gradients in the surface tension. These gradients generate Marangoni forces that move warm (low surface tension) water along the surface to locations with relatively cool (high surface tension) water. While water moves along the surface it is gradually cooled due to evaporation. At the edges of convection cells, opposing streams of cooled water collide, generating small plumes of cold water that are actively pushed down into the upper bulk by the surface Marangoni forces. A more quantitative discussion of the Marangoni plumes will be presented in § 3.5.

The Marangoni effect on the near-surface flow topology is elucidated in figure 4, showing sequences of surface divergence ( $\beta ^*=-\partial w^*/ \partial z^*$ ) contour plots for cases M0, M1, M5, M10 and M15 (cf. table 1). A non-uniform distribution of large and small convection cells can already be observed at ${\textit{Ma}}_h=0.88\times 10^5$ (case M1). The average proportion of the surface occupied by small cells is observed to increase significantly with ${\textit{Ma}}_h$ (which is further quantitatively evidenced by the histograms in figure 5). As a result, on average the number of large cells becomes negligible (even though the area they occupy remains non-negligible) compared to the number of small convection cells. Obviously, as the number of cells increases, the characteristic length scale of the small cells tends to gradually reduce with ${\textit{Ma}}_h$ , which is in agreement with the dimensional analysis (2.13), as was previously reported by, e.g. Schwarzenberger et al. (Reference Schwarzenberger, Köllner, Linde, Boeck, Odenbach and Eckert2014). Note that to compute the area of individual convection cells, the cells were identified by isolated patches at the surface with a surface divergence $\beta ^* \geqslant -0.5\beta ^{*-}_{\textit{rms}}$ , where $\beta ^{*-}_{\textit{rms}}$ is the corresponding root mean square of all $\beta ^* \leqslant 0$ .

Figure 4. Sequences of scaled surface divergence $\beta ^*/\beta ^*_{\textit{rms}}$ contour plots extracted from cases M0, M1, M5, M10 and M15 with $[{\textit{Ma}}_h \,;\, \beta ^*_{\textit{rms}}]=[0\,;\, 0.44\times 10^4]$ , $[0.88\times 10^5\,;\, 0.80\times 10^4]$ , $[4.40\times 10^5\,;\, 2.71\times 10^4]$ , $[8.81\times 10^5\,;\, 6.51\times 10^4]$ , $[13.21\times 10^5\,;\, 9.4\times 10^4]$ , respectively. The first image in each row represents a snapshot near the end of a large-rising-plume event, characterised by an overall reduction in the number of cells and often by large regions completely devoid of small-scale Marangoni structures. The two subsequent images show the emergence of small-scale Marangoni plumes in the absence of large-plume events. In each row, the time interval between the first, second and third snapshots is $\Delta t^*=0.000115$ . The last snapshot is about 0.5 bulk time units apart from the first snapshot and illustrates the state immediately before the next large-plume event (see § 3.2 for the increase in $ \beta ^*_{\textit{rms}}$ with ${\textit{Ma}}_h$ ).

Figure 5. Distribution of convection cell footprint sizes from cases M5, M10 and M15 with ${\textit{Ma}}_h=4.40\times 10^5$ , $8.81\times 10^5$ and $13.21\times 10^5$ , respectively. Here $A_c/A_s$ is the proportion of surface area occupied by a convection cell.

From the sequences in figure 4, it can also be seen that, for ${\textit{Ma}}_h\geqslant 4.4\times 10^5$ , the number of small cells (and, hence, the surface area that they occupy) significantly fluctuates in time between two limiting states, one with relatively few cells and one with a large number of cells. The time interval between the two limiting states is determined by the characteristic period at which large plumes are generated by buoyancy forces at the bottom of the computational domain. Starting from the state with a relatively large number of (small) convection cells (roll cells), a (significant) reduction in the number of these cells is achieved through the surface interaction between warm buoyant plumes and Marangoni structures. On their way to the surface (depending on the strength of the plumes) the displaced water generated by these plumes either (i) sweeps near-surface structures immediately above to the side, resulting in rapidly spreading large convection cells, similar to the eruptions reported in Köllner et al. (Reference Köllner, Schwarzenberger, Eckert and Boeck2016) that are completely devoid of any small-scale Marangoni structures (cf. left most panes in figure 4), or (ii) gets entrained into the Marangoni cells, thereby (significantly) increasing their size (small eruptions). Both events result in an overall reduction of the number of cells.

After such reduction events, there is usually a period with no or very minor interactions between warm buoyant plumes and near-surface structures. During such periods the number of small cells quite rapidly increases due to a sequence of Marangoni-induced cell divisions. A possible explanation of the cell division process is provided below. It starts with pairs of counter-rotating vorticities lying side by side along opposite edges of the convection cell. Each vortex has a footprint consisting of a narrow strip of negative surface divergence at the joint edge of neighbouring (often) elongated convection cells, which are typically sandwiched between narrow areas of high positive surface divergence (cf. figure 4). If the cell is sufficiently large, these vortices will induce a pair of counter-rotating mirror vorticities (transporting fluid from the surface into the upper bulk). As these mirror vorticities become stronger, this process will eventually lead to the break up of the convection cell into smaller cells.

Below, the aforementioned observations will be further analysed and linked to the interfacial heat and mass transfer.

3.2. Flow statistics

Figures 6(a) and 6(b) show profiles of the normalised velocity fluctuations ( $\sqrt {\langle w^\prime w^\prime \rangle }/U_\kappa$ , $\sqrt {\langle u^\prime u^\prime + v^\prime v^\prime \rangle }/(2U_\kappa )$ ) as a function of $z/h$ , where $\langle \boldsymbol{\cdot }\rangle$ denotes averaging both in time and the (homogeneous) horizontal directions. The profiles shown are characteristic also for the other simulations, and illustrate that, irrespective of the Marangoni number, the vertical fluctuations in the bulk are roughly $1.5$ to $2.5$ times larger than the horizontal fluctuations. This indicates that the momentum transport in the bulk is predominantly in the vertical direction, which is due to the rising and sinking buoyant plumes (see also figure 2). The typical non-dimensional time scale for the large-scale fluctuations in the bulk – driven mostly by buoyancy – as estimated by

(3.1) \begin{equation} \tau _b^*=\tau _b U_\kappa /h=\frac {\kappa }{\int ^z_0 \sqrt {\langle w^\prime w^\prime \rangle } \mathrm{d}z} \end{equation}

was found to be about $0.002093$ in all simulations and, therefore, appears to not depend on ${\textit{Ma}}_h$ .

Figure 6. Flow statistics: (a) vertical velocity fluctuations and (b) horizontal velocity fluctuations as a function of $z/h$ extracted from simulations M0, M5, M10, M15, with ${\textit{Ma}}_h=0$ , $4.40\times 10^5$ , $8.81\times 10^5$ and $13.21\times 10^5$ , respectively. (c) Variation of surface divergence $\beta ^*_{\textit{rms}}$ and (d) surface kinetic energy $K_s$ as a function of ${\textit{Ma}}_h$ .

Figures 6(a) and 6(b) also show the significant enhancement in both the horizontal and vertical velocity fluctuations immediately below the surface with increasing ${\textit{Ma}}_h$ . These enhanced fluctuations are only able to penetrate thin layers adjacent to the surface (before the profiles intersect with the ${\textit{Ma}}_h=0$ profile). The thickness of these layers was found to be roughly $0.1h$ and $0.2h$ for the horizontal and vertical fluctuations, respectively. Further down, the flow is dominated by buoyancy.

The significant increase in the velocity fluctuations just below the surface is a direct consequence of the boundary condition (2.9): at the surface, the gradients in the horizontal fluctuations depend linearly on ${\textit{Ma}}_h$ . Hence, for Marangoni-dominated flows, with increasing ${\textit{Ma}}_h$ , the surface divergence ( $\beta ^*=-\partial w^*/ \partial z^*=\partial u^*/\partial x^*+\partial v^*/\partial y^*$ ) also increases linearly with ${\textit{Ma}}_h$ , which is evidenced in figure 6(c). As a result, the magnitude of the downward velocity of the Marangoni plumes at the surface increases, which is in agreement with (2.14). Note that while $\beta$ increases with ${\textit{Ma}}_h$ for surface-temperature-gradient-induced Marangoni forces, it decreases with increasing Marangoni number for surfactant-induced Marangoni forces. A detailed explanation of the physical mechanisms involved in the latter can be found in, e.g. Shen, Yue & Triantafyllou (Reference Shen, Yue and Triantafyllou2004), McKenna & McGillis (Reference McKenna and McGillis2004), Wissink et al. (Reference Wissink, Herlina, Akar and Uhlmann2017). The effects of these different mechanisms on interfacial mass transfer will be discussed in § 3.6.

In Wissink & Herlina (Reference Wissink and Herlina2023) it was shown that for the developing purely Marangoni driven instability, the total surface kinetic energy $\langle K_s \rangle$ linearly depends on the Marangoni number. When adding buoyancy (RBM instability), this linear dependency only persisted at sufficiently large Marangoni numbers. Similarly, also for fully developed RBM instabilities (with constant ${\textit{Ra}}_h$ ), a clear linear dependency of $\langle K_s \rangle$ on ${\textit{Ma}}_h$ was only found for larger ${\textit{Ma}}_h$ (cf. figure 6 d). For smaller ${\textit{Ma}}_h$ , buoyancy becomes increasingly important resulting in a gradual loss of linearity.

3.3. Mean temperature profile and heat transfer scaling

Figure 7(a) shows time and horizontally averaged temperature profiles. All profiles are characterised by relatively thin thermal boundary layers at the surface and bottom, and have an extended constant (fully mixed) temperature region in the bulk. When increasing the Marangoni force, the surface temperature $\langle T^*_s \rangle$ can be seen to become warmer due to a more intense vertical mixing near the surface. Simultaneously, because of the increased transfer rate of relatively cold surface water into the bulk, the bulk temperature $\langle T^*_b \rangle$ becomes colder, and approaches an asymptotic value of about $-0.025$ (cf. figure 7 b). As can be seen in figure 7(c), this increased vertical mixing at the surface results in a thinning of the mean surface thermal boundary layer thickness $\langle \delta _T \rangle$ defined in (2.11).

Figure 7. (a) Normalised temperature profiles at various Marangoni numbers. (b) Variation of normalised surface temperature and bulk temperature as a function of ${\textit{Ma}}_h$ . (c) Non-dimensional mean thermal boundary layer thickness $\langle \delta _T/h \rangle$ .

Figure 8(a) shows the corresponding non-dimensional heat transfer velocity (Nusselt number), defined as

(3.2) \begin{equation} \textit{Nu}=\frac {q h}{\langle T_s-T_b\rangle } = \frac {h}{\langle \delta _T\rangle }, \end{equation}

as a function of $ ({\textit{Ma}}_h+\varepsilon {\textit{Ra}}_h )$ rather than ${\textit{Ma}}_h$ to reflect the fact that the Marangoni and buoyant instabilities reinforce one another (cf. e.g. Nield Reference Nield1964). At first sight, the data suggest a linear relation between ${\textit{Nu}}$ and $ ({\textit{Ma}}_h+\varepsilon {\textit{Ra}}_h )$ . However, for $ q\leqslant 0$ , it is to be expected that $\delta _T \rightarrow h$ due to thermal diffusion when ${\textit{Ma}}_h, {\textit{Ra}}_h \rightarrow 0$ , so that ${\textit{Nu}}\rightarrow {\textit{Nu}}_\kappa =1$ , where ${\textit{Nu}}_\kappa$ is the Nusselt number corresponding to the purely diffusive case. The latter is not predicted by a linear interpolation so that a power law relation,

(3.3) \begin{equation} ({\textit{Nu}}-{\textit{Nu}}_\kappa ) = a_N \left ({\textit{Ma}}_h+\varepsilon {\textit{Ra}}_h\right )^r, \end{equation}

is more applicable. Note that $\varepsilon =0.0016$ and $a_N$ (cf. figure 8 a) were determined empirically. In fact, the figure suggests the existence of two regimes. For the buoyancy-dominated regime, ${\textit{Ma}}_h$ can be seen to modulate the Rayleigh number, effectively replacing $ ({\textit{Ma}}_h+\varepsilon {\textit{Ra}}_h )$ by ${\textit{Ra}}_{\textit{eff}}$ . The obtained value of $r=0.25$ is consistent with the a priori estimate $\delta _R\propto \delta _T\propto {\textit{Ra}}_{\textit{eff}}^{-0.25}$ for purely buoyancy-driven flows (cf. § 2.2). Similarly, in the Marangoni-dominated regime, ${\textit{Ra}}_h$ slightly alters ${\textit{Ma}}_h$ , so that an effective Marangoni number ${\textit{Ma}}_{\textit{eff}}= ({\textit{Ma}}_h+\varepsilon {\textit{Ra}}_h )$ is obtained. Here, the value of $r=0.5$ was found to provide a very good interpolation of the data and is consistent with the a priori scaling $\delta _\sigma \propto \delta _T\propto {\textit{Ma}}_{\textit{eff}}^{-0.5}$ for purely Marangoni driven flows (cf. § 2.2). The existence of such a Marangoni-dominated regime is further evidenced in figure 8(b), where for ${\textit{Ma}}_h\geqslant 4.4\times 10^5$ , ${\textit{Ma}}_{\delta _T}={\textit{Ma}}_h (\delta _T/h)^2$ becomes approximately constant converging to a value of about $53$ . Note that the existence of such a dual regime (Rayleigh-dominated vs Marangoni-dominated) was also found in the developing RBM instability case (Wissink & Herlina Reference Wissink and Herlina2023). Furthermore, to facilitate comparison with experiments and other simulations, it could be convenient to recast the ${\textit{Ma}}_h+\varepsilon {\textit{Ra}}_h$ into $Bo_h^{-1}+\varepsilon$ , where $Bo_h={\textit{Ra}}_h/{\textit{Ma}}_h$ is the Bond number.

Figure 8. (a) Variation of non-dimensional heat transfer velocity ${\textit{Nu}}-{\textit{Nu}}_\kappa$ as a function of ${\textit{Ma}}_h+\varepsilon {\textit{Ra}}_h$ , with $\varepsilon =0.0016$ . (b) Variation of ${\textit{Ma}}_{\delta _T}$ with ${\textit{Ma}}_h$ .

3.4. Temperature fluctuations

As observed in figure 9(a), due to the action of Marangoni forces at the surface, significant differences in the $T^*_{\textit{rms}}=\sqrt {\langle T^{*\prime } T^{*\prime }\rangle }$ profiles can be seen in the upper bulk. These differences gradually decrease with distance to the surface, though their presence is still noticeable near $z=0$ . The temperature fluctuation at the surface, $T^*_{\textit{rms,s}}$ , and the peak in $T^*_{\textit{rms}}$ near the surface both become smaller with increasing ${\textit{Ma}}_h$ . This is associated with the variation of the temperature gradient $|{\langle \partial T^{*}/\partial z^{*} \rangle }|$ with $z/h$ shown in figure 9(b), where the locations of the $T^*_{\textit{rms}}$ peaks are indicated by the red crosses. It can be seen that for ${\textit{Ma}}_h \leqslant 1.76 \times 10^5$ , the location of these peaks is fully embedded in a region of relatively large $|{\langle \partial T^*/\partial z^* \rangle }|$ and, therefore, associated with a relatively large production of $T^*_{\textit{rms}}$ . In contrast, for ${\textit{Ma}}_h \geqslant 8.81\times 10^5$ , the peak in $T^*_{\textit{rms}}$ is located in a region of relatively small $ |{\langle\partial T^*/\partial z^* \rangle }|$ leading to a much reduced production of $T^*_{\textit{rms}}$ .

Figure 9. Profiles at various Marangoni numbers of (a) normalised temperature fluctuation, (b) normalised temperature gradient $|{\partial \langle T^* \rangle/\partial z^*}|$ , (c) $T^*_{\textit{rms}}/(T^*_b-T^*_s)$ , (d) $|{\partial ^2 \langle T^* \rangle/\partial z^{*2}}|$ . Variation of (e) surface-temperature fluctuations and ( f) boundary layer thickness with ${\textit{Ma}}_h$ .

To explain the latter, in figure 9(c), $T^*_{\textit{rms}}$ is scaled by $\langle T^*_b-T^*_s \rangle$ to remove any bias due to the ${\textit{Ma}}_h$ dependency of $\langle T^*_b-T^*_s \rangle$ . It can be seen that at the surface the scaled temperature fluctuation $\theta ^\prime = T^*_{\textit{rms,s}}/(T^*_b-T^*_s)$ only slightly reduces with increasing ${\textit{Ma}}_h$ for ${\textit{Ma}}_h \leqslant 2.65\times 10^5$ . Contrastingly, $\theta ^\prime$ values were found to increase with ${\textit{Ma}}_h$ for larger Marangoni numbers (cf. also figure 9 e). The latter indicates a major change in the surface flow dynamics as the RBM instability becomes dominated by Marangoni forces.

Furthermore, it can be observed that below the surface both $\theta ^\prime$ and the relative size of the peaks $\theta ^\prime _{p}-\theta ^\prime _{s}$ (where the subscripts $p, s$ correspond to the peak and surface, respectively) tend to increase with Marangoni number (cf. figure 9 c). This increase is closely linked to the decay rate of vertical temperature gradients expressed by $|{\langle \partial{^2} T^*/\partial {z^{*2}} \rangle }|$ shown in figure 9(d). At ${\textit{Ma}}_h=0$ , this decay is slowest and results in a relatively small peak in $\theta ^\prime$ located very close to the surface. The increase in the decay of the gradient $\left | \langle \partial T^*/\partial z^* \rangle \right |$ with increasing ${\textit{Ma}}_h$ results in an enhanced increase in $\theta ^\prime$ with distance to the surface, leading to an increase in $\theta ^\prime _{p}-\theta ^\prime _{s}$ and an approximately linear increase in the thickness $\delta ^*_{\textit{TT}}$ (which is the distance between the location of $\theta ^\prime _{p}$ – identified by the markers in figure 9 d – and the surface), at least for ${\textit{Ma}}_h\leqslant 4.4\times 10^5$ .

For the larger Marangoni numbers ${\textit{Ma}}_h \geqslant 8.81\times 10^5$ , the very rapid decrease in the gradient $|{\langle \partial T^*/\partial z^* \rangle }|$ in the near-surface region still leads to a further increase in $\theta ^\prime _{p}-\theta ^\prime _{s}$ . Simultaneously the location where the gradient becomes too weak to support any further growth in $\theta ^\prime$ approaches the surface, eventually resulting in a reduction in $\delta ^*_{\textit{TT}}$ . Note that the latter mechanism already affects the results at ${\textit{Ma}}_h=4.4\times 10^5$ , which can also be seen in figure 9( f), where $\delta ^*_{\textit{TT}}$ exhibits linear growth only for the smaller ${\textit{Ma}}_h \leqslant 2.64\times 10^5$ , while for larger ${\textit{Ma}}_h$ , this trend reverses. For comparison, also included in figure 9( f) are $\delta ^*_T$ , $\delta ^*_{\textit{wT}}$ (the distance between the peak in the $\langle w^{*\prime } T^{*\prime } \rangle$ profile and the surface) and finally $\delta ^*_{T10}$ , which is defined as the distance between the surface and the location where $ |\partial \langle T^* \rangle/\partial z^* |=0.1 |\partial T^*/\partial z^* |_s$ . It can be seen that for ${\textit{Ma}}_h \leqslant 2.64\times 10^5$ , $\delta ^*_{\textit{TT}}$ is significantly smaller than $\delta ^*_{T10}$ , indicating a non-negligible vertical gradient in $T$ for $0\leqslant \zeta \leqslant \delta ^*_{T10}$ , while for ${\textit{Ma}}_h \geqslant 8.81\times 10^5$ , $\delta ^*_{\textit{TT}}$ was found to be slightly larger than $\delta ^*_{T10}$ .

The above statistics can be linked to the flow dynamics discussed earlier in § 3.1 (cf. figures 2 to 4). Because of the very few (large) plumes emerging at ${\textit{Ma}}_h=0$ , the location of the peak in the horizontal temperature fluctuations is relatively close to the surface. For the smaller ${\textit{Ma}}_h\leqslant 4.40\times 10^5$ , the number of plumes somewhat increases with ${\textit{Ma}}_h$ causing the $T^*_{\textit{rms}}$ peak to move further away from the surface. This trend is reversed for the larger Marangoni numbers, where the average convection cell size significantly reduces and, hence, the size of the plumes formed at the edges of these cells becomes very small. Consequently, these small plumes lose buoyancy and tend to linger immediately under the surface for a significant time until either the larger of these plumes become sufficiently buoyant and moves into the upper bulk, or a significant buoyant event from below displaces these small convection cells, dragging them into the bulk. As a result, the $T^*_{\textit{rms}}$ peak approaches the surface with increasing ${\textit{Ma}}_h\geqslant 8.81\times 10^5$ .

3.5. Concentration profiles and Marangoni layer

As mentioned in § 2.4, after the fully developed flow and temperature fields were established, mass transfer calculations for four different Schmidt numbers ( ${\textit{Sc}}=16, 50 , 100, 200$ ) were activated starting with the initial concentration profiles given in (2.20). In all simulations, a statistically quasi-steady state for the mass concentrations was achieved within one bulk time unit.

Figure 10(a) shows (passive) scalar concentration profiles at ${\textit{Sc}}=200$ . The shape of the mean profiles can be seen to significantly change with increasing Marangoni number. As typical for a buoyancy-dominated flow, at ${\textit{Ma}}_h=0$ the profile shows a gradual decrease in concentration when moving away from the surface. Contrastingly, for ${\textit{Ma}}_h\geqslant 4.40\times 10^5$ , the profiles are characterised by the appearance of a local minimum in the concentration very close to the surface. This minimum is followed by a local maximum, indicating the approximate location of the edge of the relatively thin Marangoni layer.

Figure 10. Near-surface concentration profiles (a) for various ${\textit{Ma}}_h$ at ${\textit{Sc}}=200$ , (b) for various ${\textit{Sc}}$ at ${\textit{Ma}}_h=13.2\times 10^5$ (case M15).

Figure 11. Near-surface cross-sectional concentration distribution for ${\textit{Sc}}=200$ and ${\textit{Ma}}_h=13.21\times 10^5$ in the plane $x/h=1$ at (a) $t^*=0.00224$ , (b) $t^*=0.00295$ . (c) Detailed view of the Marangoni plume identified by the box in (b), together with fluctuating velocity vectors.

An explanation of these features can be seen in figure 11, showing cross-sectional snapshots at $x/h=1$ of the concentration at ${\textit{Sc}}=200$ focusing on the near-surface region of simulation M15 ( ${\textit{Ma}}_h=13.21\times 10^5$ ). The snapshots shown were extracted at $t^*=0.00224$ (figure 11 a) and $t^*=0.00295$ (figure 11 b). The corresponding convection cell footprints at the surface can be found in figure 4. Both snapshots show a multitude of plumes of varying sizes. Most of the plumes are small with a limited penetration depth and are associated with the small convection cells observed in figure 4 that are typical for Marangoni-dominated flows. Figure 11(c) shows a detailed cross-section through the Marangoni plume identified by the rectangle in figure 11(b), together with fluctuating velocity vectors. As illustrated in the figure, because of the specific topology of this characteristic Marangoni plume, which is connected to the surface by its stalk, a clearly identifiable minimum in concentration is obtained between the surface and the lobes of the plume. The Marangoni plume typically starts its life as a downward moving jet of cold, saturated flow that induces two counter-rotating surface-parallel vortices on either side of its symmetry plane (cf. also figure 3). As it approaches its maximum penetration depth, the two vortices move warm, unsaturated ambient water upwards around the lobes and finally inwards immediately below the surface. This process is responsible for the appearance of a local minimum in the mean near-surface concentration profile.

Figure 12. Profiles of the normalised concentration fluctuations $c^*_{\textit{rms}}/(c^*_s-c^*_b)$ for ${\textit{Sc}}=16-200$ as obtained in simulations M0, M1, M5, M10, M15.

The presence of the Marangoni sublayers in the mean concentration profiles at ${\textit{Sc}}=200$ and ${\textit{Ma}}_h\geqslant 4.40\times 10^5$ are visible owing to the combination of very low scalar diffusivity and a sufficiently high Marangoni number, which ensures that scalar diffusivity and buoyancy effects, respectively, do not mask the characteristic features of the Marangoni-dominated region in the concentration profiles. In contrast, as depicted in figure 10(b) (showing the mean concentration profiles obtained at ${\textit{Ma}}_h=13.21\times 10^5$ for ${\textit{Sc}}=16, 50 ,100, 200$ ), the local minima in the profiles that are located very close to the surface gradually vanish as the diffusivity increases with reducing Schmidt number. Figures 10(a) and 10(b) also clearly show the (expected) significant reduction in the concentration boundary layer thickness (steepening of the concentration gradient ( $\partial c^*/\partial z^*$ ) at the surface) both with increasing ${\textit{Ma}}_h$ and ${\textit{Sc}}$ . The former due to stronger fluctuations introduced to the upper bulk by the Marangoni forces and the latter by reduced scalar diffusivity.

In figure 12 the effect of increasing ${\textit{Ma}}_h$ and ${\textit{Sc}}$ on the concentration fluctuation ( $c^*_{\textit{rms}}= \sqrt { \langle c^{*\prime } c^{*\prime } \rangle }$ ) profile is shown. A significant difference can be observed between the profiles of the two most extreme cases studied here (M0 with ${\textit{Ma}}_h=0, Sc=16$ versus M15 with ${\textit{Ma}}_h=13.21\times 10^5, Sc=200$ ). While in the former case, only one peak near the surface is observed, two peaks can be clearly identified in the latter case. The distinct characteristics of these profiles for the different regimes are discussed in more detail below.

For cases M0 and M1 ( ${\textit{Ma}}_h \leqslant 0.88\times 10^5$ ), the $c^*_{\textit{rms}}$ profiles show a well-defined maximum (identified by ‘x’). The distance between this (local) maximum and the surface corresponds to the concentration boundary layer thickness $\langle \delta _c \rangle /h$ , which is found to gradually decrease with increasing Schmidt number.

For cases M5–M15 ( ${\textit{Ma}}_h \geqslant 4.40\times 10^5$ ) and ${\textit{Sc}} \geqslant 100$ , two well-defined $c^*_{\textit{rms}}$ peaks were observed. The surface-nearest peak was found to coincide with the edge of the concentration boundary layer, of which the thickness

(3.4) \begin{equation} \frac {\langle \delta _c \rangle }{h} =\left |{\frac {c^*_s-c^*_b}{\left . \partial c^* / \partial z^* \right |_s}}\right | \end{equation}

can be seen to reduce with increasing ${\textit{Sc}}$ . The thickness, $\delta _M$ , of the Marangoni layer is defined by the distance between the surface and the penetration depth of the Marangoni plumes. It was found that $\delta _M$ approximately corresponds to the distance between the surface and the location of the lower peak (identified by the dashed grey line). The figure clearly shows that this location only depends on ${\textit{Ma}}_h$ and is ${\textit{Sc}}$ independent. Note that for the cases M5–M15, the $c^*_{\textit{rms}}$ profiles at ${\textit{Sc}}=50$ already suggest the presence of two peaks of which the lower one coincides with the dashed grey line.

Figure 13. Variation of the Marangoni layer thickness $\langle \delta _M \rangle /h$ with ${\textit{Ma}}_h$ . For comparison, also shown are the estimates $\langle \delta _M \rangle /h \approx a_1 \sqrt {A_c}/h$ and $\langle \delta _M \rangle /h \approx a_2 \delta _{\sigma }/h \propto {\textit{Ma}}_h^{-0.5}$ , with $a_1=0.1$ and $a_2=1.6$ . Here $\langle A_c \rangle$ denotes the approximate convection cell footprint area and $\delta _\sigma$ is the a priori estimate of the Marangoni layer (2.13).

Figure 13 shows that the thickness of the Marangoni layer $\langle \delta _M \rangle$ decreases with increasing ${\textit{Ma}}_h$ according to the power law $\langle \delta _M \rangle \propto {\textit{Ma}}_h^{-0.5}$ , which is in agreement with the a priori determined scaling of the Marangoni length scale $\delta _\sigma$ (2.13). Also shown is an estimate of $\langle \delta _M \rangle$ using $\langle \delta _M/h \rangle \approx 0.1 \sqrt {A_c}/h$ , where $A_c$ is the approximate convection cell footprint area. For the larger Marangoni numbers, this estimate provides a reasonably good prediction of the Marangoni layer thickness, confirming the close relation between the characteristic Marangoni convection cell footprint and the thickness of the Marangoni layer.

3.6. Scaling of mass transfer velocity

Previously it was shown that due to increased mixing near the surface, the non-dimensional heat transfer, expressed by the Nusselt number, significantly increases with Marangoni number. The same is expected for the mass transfer velocity

(3.5) \begin{equation} k_L = \frac {\left . - \!D \partial c / \partial z \right |_s} {c_s-c_b} \end{equation}

of the passive scalars, representing atmospheric gases. This expected increase in $K_{\!L}=\langle k_L \rangle$ is in stark contrast to what was observed at contaminated surfaces, where the surfactant-concentration-gradient-induced Marangoni forces tend to inhibit mass transfer. In Wissink et al. (Reference Wissink, Herlina, Akar and Uhlmann2017) it was shown that $n$ in $K_{\!L} \propto \textit{Sc}^{-n}$ gradually decreases from $n=0.667$ for severely contaminated surfaces to $n=0.5$ for surfactant-free surfaces. In the current simulations, the power $n$ (obtained using the results for ${\textit{Sc}}\geqslant 50$ ) was found to continuously reduce further with increasing Marangoni number from $n\approx 0.5$ at ${\textit{Ma}}_h=0$ to $n\approx 0.438$ at ${\textit{Ma}}_h=13.2\times 10^5$ (see figure 14 a). This reduction in $n$ confirms the promotion of interfacial mass transfer by surface-temperature-gradient-induced Marangoni forces to values clearly below $n=0.5$ . The results presented make it likely that even lower values of $n$ are possible for ${\textit{Ma}}_h \gt 13.2 \times 10^5$ . In addition, at each ${\textit{Ma}}_h$ the unique scaling of $K_{\!L}$ with ${\textit{Sc}}$ (shown in figure 14(a) for ${\textit{Sc}}\,\geqslant \,50$ ) can likely be extrapolated to ${\textit{Sc}}\gt 200$ , as supported by previous DNS results performed for ${\textit{Sc}}$ up to $500$ (e.g. Herlina & Wissink Reference Herlina and Wissink2014; Wissink & Herlina Reference Wissink and Herlina2016; Herlina & Wissink Reference Herlina and Wissink2019).

Figure 14. (a) Scaling of the normalised transfer velocity $K_{\!L}/U_\kappa \propto \textit{Sc}^{-n}$ for cases M0, M1, M5, M10, M15, with ${\textit{Ma}}_h=0$ , $0.88\!\!\times \!\!10^5$ , $4.40\!\!\times \!\!10^5$ , $8.81\!\!\times \!\!10^5$ and $13.21\!\!\times \!\!10^5$ , respectively. The obtained values for $n$ are based on the higher ${\textit{Sc}}\geqslant 50$ cases, indicated by the solid lines. (b) Variation of the normalised transfer velocity $K_{\!L}$ and $Sh$ as a function of $({\textit{Ma}}_h+\varepsilon {\textit{Ra}}_h)$ . The solid line represents $K_{\!L}\textit{Sc}^n/U_\kappa = 0.23 ({\textit{Ma}}_h+\varepsilon {\textit{Ra}}_h)^{0.5}$ and $({\textit{Sh}}-{\textit{Sh}}_D)\textit{Sc}^{n-1} = 0.033 ({\textit{Ma}}_h+\varepsilon {\textit{Ra}}_h)^{0.5}$ , while the dashed line represents $K_{\!L}\textit{Sc}^n/U_\kappa = 6.12 ({\textit{Ma}}_h+\varepsilon {\textit{Ra}}_h)^{0.25}$ and $({\textit{Sh}}-{\textit{Sh}}_D)\textit{Sc}^{n-1} = 0.88 ({\textit{Ma}}_h+\varepsilon {\textit{Ra}}_h)^{0.25}$ , where $\varepsilon =0.0016$ and $n$ can be found in the legend of (a).

This overall trend can be explained using the eddy diffusivity model of Ledwell (Reference Ledwell1984) that reads

(3.6) \begin{equation} K_{\!L}\approx \left ( \int _0^\infty [D+D_1]^{-1} d\zeta \right )^{-1}\!\!, \end{equation}

and is valid for $D$ sufficiently small. Here, the eddy diffusivity $D_1$ is modelled by the product of the mixing length ${L}(\zeta ) \approx \zeta$ (distance to the surface) and the mixing velocity $w_{\textit{rms}}(\zeta )$ (the root mean square of the vertical velocity) evaluated near the surface, which is responsible for the vertical transport of mass. By subsequently approximating the behaviour of $w_{\textit{rms}}(\zeta )$ immediately below the surface using $w_{\textit{rms}}(\zeta ) \approx a \zeta ^\gamma$ (where $a$ is a constant), the eddy diffusivity

(3.7) \begin{equation} D_1(\zeta )\approx a \zeta ^{1+\gamma } \end{equation}

is obtained. For free-slip and no-slip boundary conditions, $\gamma =1$ and $\gamma =2$ , respectively, can be obtained analytically using a Taylor series expansion of $w^\prime$ near the surface. Employing these values of $\gamma$ in (3.6) leads to the scaling $K_{\!L}\propto D^{n}$ (equivalent to $K_{\!L}\propto \textit{Sc}^{-n}$ ), where $n=1/2$ for free-slip and $n=2/3$ for no-slip boundary conditions (as shown in Ledwell (Reference Ledwell1984)). As mentioned above, Wissink et al. (Reference Wissink, Herlina, Akar and Uhlmann2017) found a smooth transition between these two states with increasing surface pollution. Following this approach, this smooth transition can be retrieved using $1 \leqslant \gamma \leqslant 2$ in (3.7). Similarly, in the present simulations, eddy diffusivities $D_1$ obtained for positive $\gamma \leqslant 1$ can be employed to predict the observed gradual changes in $n$ (cf. figure 14 a). Overall, by combining all results, the present data suggest that $n$ increases from $\approx 0.438$ at ${\textit{Ma}}_h = 13.2 \times 10^5$ to $n=0.5$ at ${\textit{Ma}}_h=0$ , while a further increase to $n=0.667$ is obtained in the study of surfactant-induced Marangoni forces. Note that in the latter case the transport of surfactant-free (high surface tension) bulk water to the surface is counteracted by surfactant-gradient-induced Marangoni forces. This is opposite to what is happening in our present simulations, where the transport of relatively warm (low surface tension) bulk water is promoted by surface-temperature-gradient-induced Marangoni forces.

Figure 14(a) also shows that for constant ${\textit{Sc}}$ , the transfer velocity increases with ${\textit{Ma}}_h$ . This increase correlates with the steepening of the corresponding surface concentration gradients at ${\textit{Sc}}=200$ (and the thinning of the concentration boundary layer) with increasing ${\textit{Ma}}_h$ , shown in figure 10(a). This enhancement is also illustrated in figure 14(b), which shows the equivalent power law relations

(3.8) \begin{equation} K_{\!L} \textit{Sc}^{n}/U_\kappa = a_{K} \left ({\textit{Ma}}_h+\varepsilon {\textit{Ra}}_h\right )^{r} \end{equation}

and

(3.9) \begin{equation} ({\textit{Sh}}-{\textit{Sh}}_D) \textit{Sc}^{n-1} = a_{S} \left ({\textit{Ma}}_h+\varepsilon {\textit{Ra}}_h\right )^{r}, \end{equation}

where $Sh=K_{\!L}h/D$ is the Sherwood number, $Sh_D=1$ is the Sherwood number for the purely diffusive case, $n$ is displayed in figure 14(a), $a_K$ , $a_S$ are the coefficients of proportionality (cf. caption of figure 14 b) and, as in (3.3), $r=0.25$ in the buoyancy-dominated regime, $r=0.5$ in the Marangoni-dominated regime and $\varepsilon =0.0016$ . Using this scaled transfer velocity the data points for various Schmidt numbers can be seen to collapse and evidence is provided of the existence of the two regimes. Note that because of the extremely expensive computational cost, the (high) Schmidt number mass transfer statistics were gathered over a time interval of only two bulk time units and only for a selection of ${\textit{Ma}}_h$ (cf. table 1), while heat transfer statistics were gathered for all ${\textit{Ma}}_h$ over time intervals exceeding $10$ bulk time units. Despite the fewer data points, the evidence for a buoyancy-dominated regime in the mass transfer case is still quite clear.

3.7. Surface structures and quantities

As can be seen in figure 5, due to the strength of the Marangoni forces at ${\textit{Ma}}_h\geqslant 4.40\times 10^5$ , the vast majority of convection cells is very small compared to the surface area of the computational domain. Only a few large cells (induced by buoyancy forces) can be observed at each time instance. A comparison of the number of convection cells at the surface with the horizontally averaged instantaneous transfer velocity at ${\textit{Sc}}=200$ is shown in figure 15. As the vast majority of the cells are small, this figure effectively illustrates that the number of small cells is strongly correlated ( $\rho (\langle k_L/U_\kappa \rangle _{x,y},N_c) \geqslant 0.79$ ) with the surface-averaged transfer velocity. Correlations of comparable quality, as shown in the figure for ${\textit{Sc}}=200$ , were also obtained for ${\textit{Sc}} \geqslant 50$ . This indicates that, on average, the contribution of the small cells to the transfer velocity (relative to their surface area) is significantly larger than that of the large cells. This is in accordance with what was observed comparing small to larger, uniformly distributed convection cells for the developing RBM instability (Wissink & Herlina Reference Wissink and Herlina2023).

Figure 15. Variation in time of the number of convection cells $N_c$ and the scaled transfer velocity $\langle k_L / U_\kappa \rangle _{x,y}$ at ${\textit{Sc}}=200$ for cases M5, M10 and M15. Also shown are the corresponding correlation coefficients $\rho (\langle k_L \rangle _{x,y},N_c)$ .

Figure 16. Snapshots showing contours of the surface divergence ( $\beta ^*=\beta h/U_\kappa$ ) at the surface and $T^*$ , $c^*_{200}$ , at distances of $\zeta /h=0.001100$ and $\zeta /h=0.0005487$ to the surface of simulations M0 with ${\textit{Ma}}_h=0$ (upper panes) and M15 with ${\textit{Ma}}_h=13.21\times 10^5$ (lower panes), respectively. The contours are extracted at an arbitrary time $t^*=0.002286$ ( $t/\tau _b = 1.1$ ).

Figure 17. Correlation coefficients between (a) instantaneous surface divergence and temperature fields $\rho (\beta ^*,T^*)$ , and (b) instantaneous $\beta ^*$ and scalars $\rho (\beta ^*,c^*_{\textit{Sc}})$ as functions of Marangoni number ${\textit{Ma}}_h$ .

After discussing temporal correlations, we now focus on the spatial correlations of surface divergence, temperature and concentrations. Figure 16 shows instantaneous contours of $\beta ^*$ at the surface and $T^*$ , $c^*_{200}$ adjacent to the surface for simulations M0 (upper panes) and M15 (lower panes). While at ${\textit{Ma}}_h=0$ (case M0) the correlation between $\beta ^*$ and $T^*$ is quite good, the quality of this correlation gradually reduces with increasing ${\textit{Ma}}_h$ (see figure 17 a), resulting in a rather poor correlation at ${\textit{Ma}}_h=13.2\times 10^5$ (case M15). This deterioration can be explained by the fact that surface tension forces directly depend on surface-temperature gradients. These forces become negligibly small in areas with virtually uniform high or low temperature. Horizontal surface-temperature gradients appear every time warm water moves towards the surface. Subsequently, whilst the generated Marangoni forces move this warm water towards the convection cell edges, it is cooled due to evaporation. At the edges, opposing Marangoni-induced currents from neighbouring convection cells collide and push this now relatively cold surface water downwards ( $\beta ^*\lt 0$ ). As a result, the locations of negative surface divergence remain well correlated with areas of relatively cold water. However, the correlation between areas of warm water and positive surface divergence becomes quite weak. This is because the local Marangoni forces generate strong longitudinal vortices along the edges of convection cells. These vortices not only push cold surface water down but also transport water from the upper bulk upwards. Due to intense mixing and relatively large heat diffusion, this water is relatively cold compared to the surface water in the middle of the convection cells, thereby adversely affecting $|{\rho (\beta ^*,T^*)}|$ .

In contrast to the weak correlation between $\beta ^*$ and $T^*$ for ${\textit{Ma}}_h \gt 0$ mentioned above, the correlation between areas of negative surface divergence (which acts to transport saturated surface water downwards) and areas of relatively high scalar concentration ( $c^*_{\textit{Sc}}$ ) is relatively good (cf. figure 16). For the buoyancy-dominated simulations at ${\textit{Ma}}_h=0$ and $0.88\times 10^5$ , the correlations were found to be independent of ${\textit{Sc}}$ . In contrast, the Marangoni-dominated simulations ( ${\textit{Ma}}_h\geqslant 4.40 \times 10^5$ ) show a clear reduction in the correlation coefficient with ${\textit{Sc}}$ . Initially, starting at ${\textit{Ma}}_h=0$ the correlation was found to worsen with increasing ${\textit{Ma}}_h$ before becoming virtually constant for ${\textit{Ma}}_h \geqslant 4.40 \times 10^5$ .

Figure 18. Normalised heat and mass transfer velocities as functions of the surface divergence $\beta ^*_{\textit{rms}}=\beta _{\textit{rms}} h^2/\kappa$ : (a) ${\textit{Sh}}-{\textit{Sh}}_D$ vs $\beta ^*_{\textit{rms}}$ with $a_\beta =1.14, 1.68, 2.43,$ for ${\textit{Sc}}=50, 100, 200$ , respectively, and (b) ${\textit{Nu}}-{\textit{Nu}}_\kappa$ vs $\beta ^*_{\textit{rms}}$ .

Above, the instantaneous spatial correlation coefficients between $-\beta ^*$ and $c^*_{200}$ as well as $\beta ^*$ and $T^*$ were shown to be $\geqslant 0.40$ for all ${\textit{Ma}}_h\geqslant 0$ . As the near-surface concentration and temperature are directly related to their respective transfer velocities, it is to be expected that the Sherwood number $Sh=K_{\!L} h / D$ and the Nusselt number ${\textit{Nu}}$ can both be linked to $\beta ^*_{\textit{rms}}=\beta _{\textit{rms}}h/U_\kappa$ . Based on McCready, Vassiliadou & Hanratty (Reference McCready, Vassiliadou and Hanratty1986), it immediately follows that $Sh \propto \sqrt {\beta ^*_{\textit{rms}}}$ . This is confirmed in figure 18(a), which shows results obtained from the mass transfer simulations at ${\textit{Sc}}=50, 100, 200$ as well as their interpolations

(3.10) \begin{equation} {\textit{Sh}}-{\textit{Sh}}_D = a_\beta \sqrt {\beta ^*_{\textit{rms}}}, \end{equation}

with $a_\beta =1.14, 1.68, 2.43$ for ${\textit{Sc}}=50, 100, 200$ , respectively. It can be seen that, for each ${\textit{Sc}}$ , the results obtained for the larger $\beta ^*_{\textit{rms}}$ (corresponding to the higher ${\textit{Ma}}_h$ ) are in very good agreement with the interpolating curve. Contrastingly, for the smaller $\beta ^*_{\textit{rms}}$ , the mass transfer results show a slight deviation, which can be seen to become more pronounced in the corresponding Nusselt number plot (cf. figure 18 b). In this figure, for high $\beta ^*_{\textit{rms}}$ (large ${\textit{Ma}}_h$ ), the results closely match the

(3.11) \begin{equation} {\textit{Nu}}-{\textit{Nu}}_\kappa = 0.51\,{\left (\beta ^*_{\textit{rms}}\right )}^{0.5} \end{equation}

curve, while the small $\beta ^*_{\textit{rms}}$ (small ${\textit{Ma}}_h$ ) results obey a different scaling law:

(3.12) \begin{equation} {\textit{Nu}}-{\textit{Nu}}_\kappa = 9.85\,{\left (\beta ^*_{\textit{rms}}\right )}^{0.2}. \end{equation}

These different scalings can be explained by the fact that the surface divergence $\beta ^*_{\textit{rms}}$ (which represents the surface velocity frequencies) scales with $U_\sigma / \delta _\sigma\! \propto\! {\textit{Ma}}_h$ for Marangoni-dominated flows and with $U_{\!R}/\delta _R\propto {\textit{Ra}}_h^{1.25}$ for fully developed buoyancy-dominated flows (cf. (2.13)–(2.16)). Hence, combined with the results presented in figure 8(a) and (3.3), it immediately follows that for a fully developed flow, the power in the scaling $({\textit{Nu}}-{\textit{Nu}}_\kappa ) \propto (\beta ^*_{\textit{rms}})^p$ yields $p=0.2$ in the buoyancy-dominated regime and $p=0.5$ in the Marangoni-dominated regime.