4.1 Characteristics of flow field

The interfacial mass transfer simulation was started using a fully developed turbulent flow field. This initial turbulent flow field was generated by running the simulation without calculating mass transfer. After a period of approximately $50L/U$ , the turbulence in the DNS domain (upper box in figure 1) was found to be fully developed, as verified by the statistically steady behaviour of the turbulent energy spectrum shown in figure 2. The mean energy spectrum, which was obtained by averaging the instantaneous energy spectrum at $z=2.75L$ over the interval $50L/U\leqslant t\leqslant 100L/U$ , clearly showed the existence of an inertial subrange with a $\unicode[STIX]{x1D705}^{-5/3}$ behaviour, similar to Flores, Riley & Horner-Devine (Reference Flores, Riley and Horner-Devine2017), as well as a broad dissipative range where $\unicode[STIX]{x1D705}>\unicode[STIX]{x1D702}$ . The same features were found in all spectra extracted at various other $z$ locations, indicating that the turbulent flow is very well resolved. Note that in the above $\unicode[STIX]{x1D705}$ is the wavenumber and $\unicode[STIX]{x1D702}$ is the Kolmogorov length scale. Subsequently, at $t=100L/U$ the full simulation was started by activating mass transfer.

Figure 2. The mean energy spectrum normalised by the Kolmogorov velocity $u_{\unicode[STIX]{x1D702}}$ and length $\unicode[STIX]{x1D702}$ scales. The dashed line indicates the $-5/3$ slope.

Figure 3. Flow statistics, ensemble averaged over $103L/U$ : (a) turbulent integral length scale $L_{11}$ ; (b) root-mean-square (r.m.s.) of horizontal $u$ and vertical $w$ velocity fluctuations; and (c) turbulent Reynolds number. The dotted lines in (a) and (c) correspond to $\overline{L_{11}}(z/L)\pm \unicode[STIX]{x1D70E}_{L}$ and $R_{T}(z/L)\pm \unicode[STIX]{x1D70E}_{R}$ , where $\unicode[STIX]{x1D70E}_{L}$ and $\unicode[STIX]{x1D70E}_{R}$ are the standard deviations of $L_{11}$ and $R_{T}$ , respectively.

The flow statistics shown in figure 3 were obtained by ensemble averaging the flow properties over the interval $50L/U\leqslant t\leqslant 153L/U$ . By combining time averaging with averaging in the large homogeneous horizontal directions, the size of this relatively short interval was deemed to be sufficient to obtain good-quality statistics. Figure 3(a) shows the variation in $z$ of the time-averaged longitudinal integral length scale $\overline{L_{11}}$ (where $\overline{\,\cdot \,}$  denotes averaging over a period in time), with

(4.1) $$\begin{eqnarray}L_{11}(z)=\int _{0}^{L_{x}/2}R_{11}(r,z)\,\text{d}r,\end{eqnarray}$$

where the longitudinal two-point correlation $R_{11}$ of the horizontal velocity $u$ is defined by

(4.2) $$\begin{eqnarray}R_{11}(r,z)=\frac{\displaystyle \int _{x=0}^{L_{x}/2}\int _{y=0}^{L_{y}}u(x,y,z)u(x+r,y,z)\,\text{d}y\,\text{d}x}{\displaystyle \int _{x=0}^{L_{x}/2}\int _{y=0}^{L_{y}}u^{2}(x,y,z)\,\text{d}y\,\text{d}x},\end{eqnarray}$$

in which $L_{x}\times L_{y}$ is the size of the horizontal plane. As in the grid-stirred experiments of, for example, Hopfinger & Toly (Reference Hopfinger and Toly1976) and Herlina & Jirka (Reference Jirka2008), the integral length scale in most of the lower part of the computational domain was found to increase linearly with increasing distance from the turbulent source. Furthermore, it can be seen that the characteristic integral length scale $L_{\infty }$ , defined by

(4.3) $$\begin{eqnarray}L_{\infty }=\overline{L_{11}}(z_{\infty }),\end{eqnarray}$$

in which $z_{\infty }$ is the location where $\overline{L_{11}}$ is maximum, is approximately $5.67L$ . The characteristic integral length scale $L_{\infty }$ identifies the size of the largest eddies in the turbulent flow, which can only exist undisturbed at distances larger than $1L_{\infty }$ from the interface. Closer to the interface, inside the so-called surface-influenced layer, the turbulence starts to lose its isotropy (cf. e.g. Hunt & Graham Reference Hunt and Graham1978; Perot & Moin Reference Perot and Moin1995). Even though the depth of our DNS domain is somewhat smaller than $5.67L$ , a reasonable isotropy ( $|1-w_{rms}/u_{rms}|<0.1$ ) in most of the lower region ( $z/L<2.5$ ) is obtained (cf. figure 3 b). In agreement with Perot & Moin (Reference Perot and Moin1995) and Walker, Leighton & Garza-Rios (Reference Walker, Leighton and Garza-Rios1996), for example, figure 3(b) also shows an initial gradual decay in both horizontal and vertical fluctuations as the turbulence diffuses upwards. Close to the surface, a significant increase in decay rate of the vertical fluctuations $w_{rms}$ is observed, which in turn induces an increase in the horizontal fluctuations $u_{rms}$ due to the redistribution of the turbulent kinetic energy. More detailed discussions on this for free-slip and/or no-slip boundary conditions can be found elsewhere (e.g. Perot & Moin Reference Perot and Moin1995; Walker et al. Reference Walker, Leighton and Garza-Rios1996; Magnaudet Reference Magnaudet2003; Bodart, Cazalbou & Joly Reference Bodart, Cazalbou and Joly2010; McCorquodale & Munro Reference McCorquodale and Munro2017).

Typically, the turbulent Reynolds number for grid-stirred turbulence is defined by

(4.4) $$\begin{eqnarray}R_{T}=\frac{u_{\infty }\unicode[STIX]{x1D6EC}}{\unicode[STIX]{x1D708}},\end{eqnarray}$$

where $\unicode[STIX]{x1D6EC}=2L_{\infty }$ (e.g. Hopfinger & Toly Reference Hopfinger and Toly1976; Brumley & Jirka Reference Brumley and Jirka1987). This definition is also adopted here, where we used $u_{\infty }=u_{rms}(z_{\infty })$ for the characteristic velocity scale. The time average of the corresponding local turbulent Reynolds number

(4.5) $$\begin{eqnarray}R_{T}(z)=u_{rms}2\overline{L_{11}}/\unicode[STIX]{x1D708},\end{eqnarray}$$

shown in figure 3(c), was found to be approximately constant ( $R_{T}(z)\approx 1856$ ) in most of the lower DNS domain. In general, the turbulent flow statistics were similar to the far-field statistics in the experiments of Brumley & Jirka (Reference Brumley and Jirka1987) and Herlina & Jirka (Reference Jirka2008) and also to the statistics obtained at low to moderate $R_{T}$ in HW14, where a more extended comparison with literature was presented.

4.2 Statistics of turbulent scalar transport

As mentioned above, the scalar mass transfer calculations for $Sc=20$ and $500$ were activated at $t=100L/U$ using the initial conditions given in (3.2). Owing to the high $R_{T}$ , the transfer of saturated fluid from the surface to the bulk was relatively fast. This, together with the large horizontal extent of the computational domain, resulted in a speedy convergence of the normalised horizontally averaged turbulent mass flux statistics. At $t=126L/U$ time units, the total mass flux near the interface ( $0\leqslant \unicode[STIX]{x1D701}\leqslant 10\unicode[STIX]{x1D6FF}$ ) was found to be independent of $z$ . Scalar transport statistics were subsequently gathered in this quasi-steady regime from $t=126L/U$ until $153L/U$ .

As illustrated by the dotted lines in figure 3(a), the integral length scale of the turbulence $L_{11}$ varied significantly in time. Hence, it is expected that $L_{\infty }$ (and consequently $R_{T}$ ) – when calculated using a moving average over a relatively small time window – will also become time-dependent. Using this we obtained a significant range of turbulent Reynolds numbers for the verification of the small-eddy model (cf. § 4.4). Note that for the smaller time windows $L_{\infty }$ is the local maximum of $L_{11}$ nearest to the location where $L_{\infty }$ was obtained for the entire interval. The dependence of $L_{\infty }$ and $R_{T}$ on the size of the time interval can also be seen in table 1, where, when averaging over the entire interval $50L/U{-}153L/U$ , one has $L_{\infty }=5.67L$ and $R_{T}=1856$ . On the other hand, by averaging over the smaller interval, values for $L_{\infty }$ and $R_{T}$ of $4.54L$ and $1440$ , respectively, were obtained.

Figure 4. Scalar statistics: (a) normalised mean concentration profiles, (b) normalised r.m.s. profiles, and (c) normalised diffusive and turbulent scalar fluxes.

Table 1. Mean flow parameters and the interval $\unicode[STIX]{x0394}t$ over which they are determined. Scalar statistics were only gathered over the smaller of the two intervals.

Figure 5. Influence of $Sc$ on the horizontally averaged boundary layer thickness at $t=140L/U$ .

The resulting scalar statistics can be seen in figure 4. The vertical axes show the distance to the surface $\unicode[STIX]{x1D701}$ normalised by the mean concentration boundary layer thickness $\unicode[STIX]{x1D6FF}$ . The latter is defined by the distance to the surface where

(4.6) $$\begin{eqnarray}c_{rms}=\sqrt{\overline{\langle c^{2}\rangle -\langle c\rangle ^{2}}}\end{eqnarray}$$

is maximum. Note that $\langle \,\cdot \,\rangle$ denotes averaging in the homogeneous (horizontal) directions. Similarly to the lower- $R_{T}$ cases in HW14, the horizontally averaged boundary layer thickness $\unicode[STIX]{x1D6FF}$ was found to scale with $Sc^{-1/2}$ as illustrated for $t=140L/U$ in figure 5 (cf. further discussion in § 4.4).

Figure 4(a) shows the mean concentration profiles near the surface, normalised by calculating $(\overline{\langle c\rangle }-\overline{\langle c_{b}\rangle })/(c_{s}-\overline{\langle c_{b}\rangle })$ . Immediately underneath the surface all profiles are approximately linear, which is a clear indication that molecular diffusion dominates gas transfer. It can also be seen that in that region the normalised profiles all collapse. For similar $Sc$ , deeper down in the bulk – due to the more effective vertical mixing – a slightly increased concentration at higher $R_{T}$ is observed. This is in agreement with the trend found for lower $R_{T}$ in HW14 and the experiments of Herlina & Jirka (Reference Jirka2008).

The normalised concentration fluctuations $c_{rms}/(c_{s}-\overline{\langle c_{b}\rangle })$ , shown in figure 4(b), can be seen to grow rapidly from zero at the surface to a local maximum, where approximately the horizontal turbulent exchange of saturated and unsaturated fluid is most intensive, before declining again to zero in the well-mixed region of the bulk. The peaks observed in the normalised $c_{rms}$ (figure 4 b) were found to be approximately 0.3 in all cases, which is typical for numerical simulations with free-slip surface boundary conditions (e.g. Magnaudet & Calmet Reference Magnaudet and Calmet2006) and is also found in some experiments (e.g. Atmane & George Reference Atmane and George2002). Some other experiments, however, show $c_{rms}$ peaks as low as 0.1–0.2 (e.g. Chu & Jirka Reference Chu and Jirka1992; Herlina & Jirka Reference Jirka2008). These relatively low peaks could be a consequence of a significant surface contamination as shown in the DNS studies of, for example, Khakpour et al. (Reference Khakpour, Shen and Yue2011) and Wissink et al. (Reference Wissink, Herlina, Akar and Uhlmann2017).

Figure 4(c) shows the mean diffusive $\overline{\langle -D\unicode[STIX]{x2202}c/\unicode[STIX]{x2202}z\rangle }$ and turbulent fluxes $\overline{\langle c^{\prime }w^{\prime }\rangle }$ normalised by the mean diffusive mass flux $\overline{\langle \,j\rangle }$ at the interface, where

(4.7) $$\begin{eqnarray}j=\left.-D\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}z}\right|_{i}.\end{eqnarray}$$

The present results for $R_{T}\approx 1440$ with $Sc=20$ and $500$ are shown together with the GS500 results from HW14 with $R_{T}=507$ and $Sc=32$ . In agreement with the theory, the diffusive fluxes are maximum at the surface and rapidly decrease to zero in the bulk, while simultaneously the turbulent fluxes increase from zero at the surface such that the sum of the mean diffusive and turbulent fluxes remains constant at all $z$ locations. The latter indicates that the simulation has run sufficiently long to achieve a quasi-steady state. For all cases, the depth at which the normalised fluxes are in equilibrium is approximately the same ( $\unicode[STIX]{x1D701}\approx 0.65\unicode[STIX]{x1D6FF}$ ).

Note that, even though the Schmidt numbers in the simulation differ significantly, all normalised scalar statistics shown in figure 4(b,c) nearly collapse. When applying the present normalisation to the lower- $R_{T}$ cases in HW14, a similar data collapse is observed. This is illustrated by including the $R_{T}=507$ results (which are typical for $84\leqslant R_{T}\leqslant 507$ ) in figure 4. Based on this, we can conclude that the normalised results are almost invariant with respect to both $Sc$ and  $R_{T}$ .

4.3 Role of large and small scales

4.3.1 Surface divergence

Figure 6 depicts a comparison between the surface divergence ( $\unicode[STIX]{x1D6FD}=-\unicode[STIX]{x2202}w/\unicode[STIX]{x2202}z$ ) and the instantaneous mass transfer velocity for $Sc=20$ at an arbitrarily chosen time $t=140L/U$ . The $\unicode[STIX]{x1D6FD}$ contours (figure 6 a) clearly show the presence of convection cells corresponding to large structures of size ${\sim}L_{\infty }$ with positive surface divergence (upwelling of unsaturated fluid) separated by narrow regions of negative surface divergence (downwelling of saturated fluid). Compared to the surface divergence plot, the footprints of fine-scale turbulent structures can be more clearly identified in the corresponding instantaneous $K_{L}$ contour plot (for a more detailed comparison see figure 11). An important reason for this discrepancy is the much higher diffusivity of the fluid (affecting $\unicode[STIX]{x1D6FD}$ ) compared to the scalar (affecting $K_{L}$ ). An increase in $R_{T}$ leads to an increase in fine-scale turbulent structures at the surface. Because of the aforementioned difference in diffusivities, these structures will remain clearly visible for much longer in the $K_{L}$ contours than in the $\unicode[STIX]{x1D6FD}$ contours. This has a negative effect on the resulting time-averaged spatial correlation $\unicode[STIX]{x1D70C}(K_{L},\unicode[STIX]{x1D6FD})$ at high $R_{T}$ , as can be seen in figure 7, where in the present high- $R_{T}$ simulation the value of $\unicode[STIX]{x1D70C}(K_{L},\unicode[STIX]{x1D6FD})=0.78$ was lower than the values $\unicode[STIX]{x1D70C}(K_{L},\unicode[STIX]{x1D6FD})=0.89$ , 0.81 and 0.81 obtained from HW14 for $R_{T}=84$ , 195 and 507, respectively. This reduction of $\unicode[STIX]{x1D70C}(K_{L},\unicode[STIX]{x1D6FD})$ with $R_{T}$ also impacts the applicability of the surface divergence model at high $R_{T}$ . This is in agreement with the findings of Turney & Banerjee (Reference Turney and Banerjee2013), who observed a breakdown of the surface divergence model in the presence of small time scales. This will be further investigated below.

Figure 6. Snapshots at $t=140L/U$ of (a) surface divergence $\unicode[STIX]{x1D6FD}$ and (b)  $K_{L}$ for $Sc=20$ .

Figure 7. Time-averaged spatial correlation of the transfer velocity $K_{L}$ with the surface divergence $\unicode[STIX]{x1D6FD}$ for various $R_{T}$ .

4.3.2 Vortical structures

Figure 8. Snapshots of vortical structures identified using $\unicode[STIX]{x1D706}_{2}=-200$ (upper panels) and $\unicode[STIX]{x1D706}_{2}=-5$ (lower panels), normalised by $L_{\infty }$ and $u_{\infty }$ : (a,c) present DNS at $t=146.6L/U$ and (b,d) case GS500 from HW14 ( $R_{T}=507$ ).

In figure 8, vortical structures in part of the computational domains are identified using the $\unicode[STIX]{x1D706}_{2}$ criterion of Jeong & Hussain (Reference Jeong and Hussain1995). To allow a direct comparison of the structures obtained in the present DNS and GS500 (HW14), the $\unicode[STIX]{x1D706}_{2}$ eigenvalues were normalised using the relevant characteristic turbulent length ( $L_{\infty }$ ) and velocity ( $u_{\infty }$ ) scales. Figure 8(a,b) shows that, because of the larger $R_{T}$ , significantly more vortical structures are found in the present simulation than in the GS500 simulation, while the thickness of the structures – which depends on the Kolmogorov length scale – is smaller. The lack of vortical structures, as observed in large areas of the GS500 simulation, indicates that the instantaneous turbulence is not uniformly distributed. Consequently, it is to be expected that the instantaneous mass transfer velocity will also not be uniformly distributed but will be elevated in regions with increased turbulence.

Figure 9. (a) Variation of the mean Kolmogorov length scale with depth. (b) Detailed view of (a) immediately beneath the surface.

To further explore what actually happens close to the surface, in figure 8(c,d) the vortical structures in the region $0\leqslant \unicode[STIX]{x1D701}/L_{\infty }\leqslant 0.04$ are visualised using $\unicode[STIX]{x1D706}_{2}=-5$ . Again, the number and size of the vortical structures in the present simulation differ significantly from the ones in GS500. As explained below, because of the free-slip boundary condition at the top, the axes of the structures close to the surface are either parallel or orthogonal to the interface. Figure 9 shows the Kolmogorov length scale $\unicode[STIX]{x1D702}$ for the present simulation as a function of the distance to the surface. Between $\unicode[STIX]{x1D701}/L_{\infty }\approx 1$ (near the bottom) and $\unicode[STIX]{x1D701}/L_{\infty }\approx 0.05$ , $\unicode[STIX]{x1D702}$ was found to increase as the turbulence diffusing from below dissipates. As seen in figures 8(a,c) and 3(b), when further approaching the surface, the flow rapidly becomes more and more two-dimensional. In upwelling regions, the increase in horizontal velocities towards the surface stretches initially ‘randomly’ oriented vortical structures and aligns them horizontally. Because of the large size of these upwelling regions (compared to the downwelling regions), this horizontal alignment is very common and causes vortical structures (on average) to become thinner towards the surface. As the diameter of the vortical structures scales with $\unicode[STIX]{x1D702}$ (Jimenez et al. Reference Jimenez, Wray, Saffman and Rogallo1993), the thinning mentioned above explains the decrease in $\unicode[STIX]{x1D702}$ observed in figure 9 between $\unicode[STIX]{x1D701}/L_{\infty }\approx 0.05$ and $0.008$ . Closer to the surface, the horizontal velocities become constant due to the free-slip boundary condition. As a result, the horizontal stretching no longer increases, which is an explanation of the slight increase in $\unicode[STIX]{x1D702}$ observed between ${\approx}0.008$ and  $0$ . In the relatively small downwelling-dominated regions, the downward velocity increases with distance from the surface, thereby stretching and strengthening the local surface-attached (surface-normal) vortical structures (SAVS). Note that the SAVS observed in the downwelling regions are seeded by relatively weak surface-normal vorticity transported upwards in the upwelling regions, where they attach to the surface and are subsequently moved towards the downwelling regions.

While the SAVS promote mixing in the horizontal direction, the structures with axes parallel to the interface promote the exchange of fluid between the upper bulk and the surface. Their relation with mass transfer is studied in figure 10(a), by combining contours of $K_{L}$ with a cross-section of the vortical structures from figure 8(c) at $\unicode[STIX]{x1D701}/L_{\infty }=0.01$ . The horizontally aligned slender vortical structures can be found near the edges of areas with high $K_{L}$ , adjacent to locations where saturated fluid is transported downwards. The SAVS, on the other hand, are generally located in areas with relatively low $K_{L}$ . To distinguish the SAVS from the horizontally aligned vortical structures, in figure 10(b), the $\unicode[STIX]{x1D706}_{2}$ isolines are combined with contours of the surface-normal vorticity component $\unicode[STIX]{x1D714}_{z}$ . In contrast to the horizontally aligned vortices, SAVS are marked by elevated levels of $\unicode[STIX]{x1D714}_{z}$ in their entire interior.

Figure 10. Isolines of $\unicode[STIX]{x1D706}_{2}=-5$ at $\unicode[STIX]{x1D701}/L_{\infty }=0.01$ and $t=146.6L/U$ with (a) contours of $K_{L}$ for $Sc=20$ and (b) contours of $|\unicode[STIX]{x1D714}_{z}|$ . Note that $\unicode[STIX]{x1D706}_{2}$ , $K_{L}$ and $\unicode[STIX]{x1D714}_{z}$ were normalised using $L_{\infty }$ and $u_{\infty }$ .

Figure 11. Zoomed view of figure 6 combined with isolines of surface-normal vorticity magnitude at $|\unicode[STIX]{x1D714}_{z}|\,(L_{\infty }/u_{\infty })=14,~34,~51$ and $68$ . (a) Contours of $\unicode[STIX]{x1D6FD}$ and (b) contours of $K_{L}$ .

Figure 11 shows a detailed view of the $\unicode[STIX]{x1D6FD}$ and $K_{L}$ contours from figure 6 combined with isolines of relatively high-intensity surface-normal vorticity $|\unicode[STIX]{x1D714}_{z}|$ . As mentioned above, because of the stretching of SAVS, it is expected to see such isolines only in downwelling areas. When this downwelling is very strong, SAVS undergo significant stretching and quickly become exceedingly thin, leading to their dissipation within a relatively short time. SAVS in areas of very weak downwelling (as seen in figure 11 a) tend to survive for a long time in the absence of large-scale turbulent motion. Because of the extended exposure time at the surface, the fluid inside these structures has become highly saturated due to diffusive mass transfer, as indicated by the local low $K_{L}$ values in figure 11(b). While the correlation of these SAVS with low $K_{L}$ values is very good, the correlation with low $\unicode[STIX]{x1D6FD}$ values is rather poor. This indicates that the local correlation $\unicode[STIX]{x1D70C}(\unicode[STIX]{x1D6FD},K_{L})$ will be rather small. Note that the above observations are generally confirmed in all recorded snapshots.

The dynamics of vortical structures in open-channel flow and their influence on interfacial mass transfer has been discussed by, for example, Pan & Banerjee (Reference Pan and Banerjee1995), Nagaosa (Reference Nagaosa1999) and Nagaosa & Handler (Reference Nagaosa and Handler2003). Pan & Banerjee (Reference Pan and Banerjee1995) discovered the existence of relatively large ring-like vortices at the surface of open-channel flow that are associated with vertical mixing. In the interior of the ring, unsaturated fluid from the bulk is transported towards the surface (splat), while immediately outside of the ring, saturated fluid is transported down into the bulk (antisplat). Nagaosa & Handler (Reference Nagaosa and Handler2003) found that the ring-like vortical structures originate from near the bottom of the channel and start their life as hairpin vortices. In our case, in the absence of bottom shear and associated vortical structures like hairpin vortices, the vortical dynamics of the turbulent flow changes significantly. Individual upwellings are usually relatively small but tend to cluster together into larger areas separated by narrow regions with strong downwelling. The topology of the clusters is very complex and rapidly changes in time. As discussed above (see also figure 10), the vortices parallel to the surface, which are responsible for vertical mixing, tend to lie along the edges of upwelling regions and sometimes resemble parts of ring-like vortices, similar to the ones previously observed by, for example, Nagaosa & Handler (Reference Nagaosa and Handler2003). As mentioned earlier, the SAVS only promote horizontal mixing and hence do not directly contribute to the turbulent vertical mass transfer, which is in agreement with the findings of Nagaosa (Reference Nagaosa1999).

4.3.3 Spectra of turbulent mass flux

Figure 12. (a) Premultiplied cospectra of turbulent mass flux $c^{\prime }w^{\prime }$ at $\unicode[STIX]{x1D701}=5\unicode[STIX]{x1D6FF}$ . Cumulative cospectra of $c^{\prime }w^{\prime }$ at (b)  $\unicode[STIX]{x1D701}=5\unicode[STIX]{x1D6FF}$ and (c)  $\unicode[STIX]{x1D701}=\unicode[STIX]{x1D6FF}$ .

As illustrated above, the mass transfer is directly related to the presence of vortical structures immediately underneath the surface. With increasing $R_{T}$ , the quantity of vortical structures increases, while their diameter typically decreases. To further investigate the relative importance of large- and small-scale structures for vertical mass transfer, time-averaged spectra of the turbulent mass flux $c^{\prime }w^{\prime }$ are plotted as a function of the normalised wavenumber $\unicode[STIX]{x1D705}L_{\infty }$ (cf. figure 12). The present results, obtained at $R_{T}=1440$ , are compared to the results of GS500 at $R_{T}=507$ and GS80 at $R_{T}=84$ (see HW14).

Figure 12(a) shows the cospectra premultiplied by the wavenumber $\unicode[STIX]{x1D705}$ . Compared to the present simulation at $R_{T}=1440$ , the energy-containing range at $R_{T}=507$ moves to smaller values of $\unicode[STIX]{x1D705}L_{\infty }$ . At $R_{T}=84$ , only a relatively narrow band of energy-containing scales was observed, which coincides with the largest energy-containing scales observed at $R_{T}=507$ . Hence, it can be concluded that at low $R_{T}$ large scales dominate the turbulent mass flux, while at high $R_{T}$ the contribution of the small scales becomes more important, supporting the two-regime mass transfer concept of Theofanous et al. (Reference Theofanous, Houze and Brumfield1976). Near the critical $R_{T}$ , both small and large scales were found to contribute significantly to the overall turbulent mass transfer.

To further illustrate the contribution of the different scales to the mass transfer, the cumulative cospectra of $c^{\prime }w^{\prime }$ are shown in figure 12(b,c) at $\unicode[STIX]{x1D701}=5\unicode[STIX]{x1D6FF}$ and $\unicode[STIX]{x1D6FF}$ , respectively. It can be seen that with increasing $R_{T}$ the smaller scales become more and more important for the vertical turbulent mass transport. For instance, at $\unicode[STIX]{x1D701}=5\unicode[STIX]{x1D6FF}$ , the contribution of relatively large scales with wavenumbers $\unicode[STIX]{x1D705}L_{\infty }\leqslant 10$ to the total energy of the turbulent flux decreases from $70\,\%$ for $R_{T}=84$ to $55\,\%$ for $R_{T}=507$ and finally to $34\,\%$ for $R_{T}=1440$ . Higher up in the boundary layer, at $\unicode[STIX]{x1D701}=\unicode[STIX]{x1D6FF}$ , the presence of the free-slip surface causes the flow to become increasingly two-dimensional. This is evidenced by the change in orientation of the vortical structures when approaching the surface as mentioned above (cf. figure 8 a,c). The inverse energy cascade, that is typical for two-dimensional turbulence, explains the increased importance of the larger scales to the vertical mass transfer when approaching the surface (cf. figure 12 c).

4.4 Scaling of transfer velocity with $Sc$ and $R_{T}$

Figure 13. Influence of $Sc$ on $K_{L}$ . Data shown are from the snapshot at $t=140L/U$ .

In agreement with the experiments of Chu & Jirka (Reference Chu and Jirka1992), Atmane & George (Reference Atmane and George2002) and Herlina & Jirka (Reference Jirka2008), for example, in the HW14 simulations it was found that the concentration boundary layer thickness $\unicode[STIX]{x1D6FF}$ depends on the level of isotropic turbulence diffusing from below. This turbulence promotes the interchange of unsaturated fluid from the bulk with saturated fluid from the surface. For constant $Sc$ , the boundary layer thickness (normalised by $L_{\infty }$ ) decreases with increasing turbulence level, resulting in a steeper concentration gradient at the surface and therefore an increased normalised mass transfer velocity. On the other hand, for constant $R_{T}$ , it was shown that $\unicode[STIX]{x1D6FF}$ scales with $Sc^{-1/2}$ (cf. figure 5). The same scaling was found for $K_{L}$ , as can be seen in figure 13. This implies that, for constant $R_{T}$ and free-slip surface boundary conditions, the mass transfer velocity $K_{L}$ is proportional to $\unicode[STIX]{x1D6FF}$ . Note that the same instantaneous scaling behaviour $K_{L}\propto Sc^{-1/2}$ was consistently found for all $t\geqslant 126L/U$ . Also, this scaling is in agreement with the theoretical scaling found for free-slip surface conditions (e.g. Ledwell Reference Ledwell1984; Jähne & Haussecker Reference Jähne and Haussecker1998).

Free-slip boundary conditions are applied to represent a clean (surfactant-free) surface. Contaminated surfaces, on the other hand, can be modelled, for example, by using a no-slip boundary condition to represent some of the physics pertinent to severe contamination (Herlina & Wissink Reference Herlina and Wissink2016) or by explicitly solving the surfactant transport at the surface and modelling the influence of surfactants on the near-surface velocity field (Shen, Yue & Triantafyllou Reference Shen, Yue and Triantafyllou2004; Khakpour et al. Reference Khakpour, Shen and Yue2011; Wissink et al. Reference Wissink, Herlina, Akar and Uhlmann2017). As seen in Wissink et al. (Reference Wissink, Herlina, Akar and Uhlmann2017), the scaling behaviour of the mass transfer velocity with the Schmidt number smoothly changes from $Sc^{-1/2}$ for clean surfaces to $Sc^{-2/3}$ for severely contaminated surfaces.

Figure 14. Variation of the normalised $K_{L}$ with $R_{T}$ . For the free-slip boundary condition, the present DNS results are combined with results from HW14 obtained at low to moderate $R_{T}$ . In addition, results from HW16, using a no-slip boundary condition, are also shown. Note that the lines represent the fitted scaling laws.

As mentioned above, the present high- $R_{T}$ , high- $Sc$ DNS study builds on our previous DNS of interfacial mass transfer driven by isotropic turbulence diffusing from below (see Herlina & Wissink Reference Herlina and Wissink2014, Reference Herlina and Wissink2016; Wissink et al. Reference Wissink, Herlina, Akar and Uhlmann2017). Figure 14 shows the normalised mass transfer velocity as a function of the turbulent Reynolds number $R_{T}$ ,

(4.8) $$\begin{eqnarray}\frac{K_{L}}{u_{\infty }Sc^{-n}}=cR_{T}^{-q},\end{eqnarray}$$

where $n=1/2$ and $n=2/3$ for the free-slip and no-slip surface boundary conditions, respectively, while $q$ depends on the size of the dominating turbulent scales, as discussed below. Note that, after using $u_{\infty }$ and $\unicode[STIX]{x1D6EC}$ for non-dimensionalisation, it can be seen that for $n=q=1/2$ equation (4.8) is equivalent to (1.3), while for $n=1/2$ and $q=1/4$ it is equivalent to (1.4) (cf. Theofanous et al. Reference Theofanous, Houze and Brumfield1976). The data shown in figure 14 originate from our previous (HW14) and present DNS of interfacial mass transfer driven by isotropic turbulence with free-slip surface boundary conditions. While HW14 only considered low to moderate $R_{T}$ values, our present DNS for the first time combines high- $R_{T}$ turbulence ( $R_{T}\approx 1440$ ) with low-diffusivity mass transfer ( $Sc$ up to $500$ ), which is typical for dissolved oxygen in water. In agreement with Theofanous (Reference Theofanous, Brutsaert and Jirka1984), two regimes can be identified in figure 14. For turbulent Reynolds numbers smaller than $R_{T}\approx 500$ , the HW14 numerical results show a scaling of the normalised mass transfer velocity with $R_{T}^{-1/2}$ , which supports the large-eddy model of Fortescue & Pearson (Reference Fortescue and Pearson1967). For larger $R_{T}$ , our present DNS results – obtained by combining horizontal averaging with ensemble averaging over time periods of $30L/U$ – clearly show the presence of a different scaling law where $(K_{L}/u_{\infty }Sc^{-1/2})\propto R_{T}^{-1/4}$ , which is in agreement with the small-eddy model of Banerjee et al. (Reference Banerjee, Scott and Rhodes1968) and Lamont & Scott (Reference Lamont and Scott1970). The same $K_{L}$ dependence on $R_{T}^{-1/4}$ was reported in the experimental work of Herlina & Jirka (Reference Jirka2008). The data shown in figure 14 are complemented by results from our no-slip simulations (cf. HW16). As in the free-slip case, both the small- and large-eddy regimes can be clearly identified. It can be seen in Herlina & Wissink (Reference Herlina and Wissink2016) that the fitted line through the free-slip numerical results for high $R_{T}$ also agrees well with the upper bound of the experimental findings of McKenna & McGillis (Reference McKenna and McGillis2004) obtained using a clean surface. Note that, even though we expect a smooth transition between the two regimes, it was found that the line through the present free-slip DNS data also intersects the result obtained in HW14 at $R_{T}=507$ . Further investigations would be needed to determine a detailed picture of the transitional regimes around $R_{T}\approx 500$ for both the free-slip and the no-slip case.