{'fraction_non_alphanumeric': 0.04315251062603559, 'fraction_numerical': 0.02326921691520784, 'mean_word_length': 4.4393939393939394, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Direct numerical simulation of turbulent mass transfer at the surface of an open channel flow

We present direct numerical simulation results of turbulent open channel flow at bulk Reynolds numbers up to 12000, coupled with (passive) scalar transport at Schmidt numbers up to 200. Care is taken to capture the very large scale motions which appear already for relatively modest Reynolds numbers. The transfer velocity at the flat, free surface is found to scale with the Schmidt number to the power"-1/2", in accordance with previous studies and theoretical predictions for uncontaminated surfaces. The scaling of the transfer velocity with Reynolds number is found to vary, depending on the Reynolds number definition used. To compare the present results with those obtained in other systems, we define a turbulent Reynolds number at the edge of the surface-influenced layer. This allows us to probe the two-regime model of Theofanous [Turbulent mass transfer at free, gas-liquid interfaces, with applications to open-channel, bubble and jet flows. Int. J. Heat Mass Transfer 19, 613--624, 1976], which is found to correctly predict that small-scale vortices significantly affect the mass transfer for turbulent Reynolds numbers larger than 500. It is further established that the root-mean-square of the surface divergence is, on average, proportional to the mean transfer velocity. However, the spatial correlation between instantaneous surface divergence and transfer velocity tends to decrease with increasing Schmidt number and increase with increasing Reynolds number. The latter is shown to be caused by an enhancement of the correlation in high-speed regions, which in turn is linked to the spatial distribution of surface-parallel vortices.


Introduction
Transfer of gases across a gas-liquid interface is of fundamental importance in many research fields, such as civil engineering and environmental science. The transfer rate is governed by a variety of processes, the most important being the hydrodynamic flow conditions at the liquid side, which are significantly affected by the surface condition. The latter includes surface contamination effects (e.g. Tsai and Liu, 2003;McKenna and McGillis, 2004;Wissink et al., 2017), surface roughness effects (e.g. Turney, 2016) and whitecap effects (e.g. Jähne, 2020;Brumer et al., 2017). In nature, turbulent flow typically generated by wind-shear, buoyancy and/or bottom-shear enhances the gas transfer rate. Wind-shear is usually the most effective driving mechanism and has been extensively studied (e.g. Plate and Friedrich, 1984;Komori et al., 1993;Wanninkhof et al., 2009;Garbe et al., 2014). At low wind-speeds, buoyancy and bottom-shear induced turbulence become more important. Particularly in deep water bodies, buoyancy effects due to surface cooling (e.g. by evaporation) often contribute significantly to the interfacial transfer of heat and greenhouse gases (e.g. Podgrajsek et al., 2014). In rivers or streams, bottom-shear induced turbulence is generally the 1 ℎ is the thermal diffusivity and is the molecular diffusivity). The main problem here is the very low diffusivity (high ) of many atmospheric gases in water, resulting in filaments with very large concentration gradients that need to be fully resolved. For this reason, DNS were usually performed at low to moderate Schmidt and/or Reynolds numbers. For example, Lakehal et al. (2003) and Tsai et al. (2013) performed DNS of interfacial mass transfer driven by surface-shear/wind for up to 10 and 7, respectively, while Schwertfirm and Manhart (2007) performed DNS of buoyancy-driven mass transfer for up to 49. Recently, DNS of high-realistic Schmidt numbers (e.g. ≈ 500 for O 2 and ≈ 600 for CO 2 in water) have been performed in buoyancy-driven flow by Wissink and Herlina (2016); Fredriksson et al. (2016) and in isotropic-turbulence driven flow by Wissink (2014, 2019). Note that the latter is often used to produce a near-surface turbulent flow field that, despite the lack of streamwise anisotropy, is used to mimic the flow field generated by bottom-shear, such as in open channel flow. Handler et al. (1999) were among the first to perform DNS of interfacial mass transfer in open channel flow, using a Prandtl number of = 2 and a friction Reynolds number of = / = 180, where is the friction velocity and is the channel height. It was observed that hairpin vortices were the dominant structures contributing to transfer of passive scalars in open channel flow. Nagaosa and Handler (2003), who performed DNS for = 150, 300 and = 1, showed that these hairpin vortices, which were generated at the bottom wall, evolved into ring-like vortices as they approach the free-surface. These ring-like vortices were observed to be present immediately below high mass transfer regions. Kermani et al. (2011) performed DNS at 300 and 0.71 ≤ ≤ 8. They focused on the quantification of the surface age (i.e. the time between two surface renewal events) and concluded that the random surface renewal model of Danckwerts (1951) is not applicable at small (young) surface age, where interfacial mass transfer is actually at its largest. More recently, Nagaosa and Handler (2012) performed DNS in open channel flow for 150 ≤ ≤ 600 with = 1. They evaluated the suitability of two characteristic time scales to approximate the renewal time in Danckwert's model. The first was based on the characteristic length and velocity scales at the free-surface, while the other was the reciprocal of the root-mean-square of the surface divergence. Both time scales were found to perform well. Magnaudet and Calmet (2006) used the computationally less demanding large eddy simulation (which only resolves the larger scales of motion) to be able to study interfacial mass transfer in open channel flow at a high friction Reynolds number of = 1280 and Schmidt numbers ranging from = 1 to 200. They confirmed the scaling of with −0.5 and found that the mass transfer was mostly driven by motions with a large streamwise extent.
In summary, because of the huge computational demands to fully resolve all scales of motion, until recently DNS of mass transfer in open channel flow was limited to ≤ 8 and ≤ 600, and relatively small computational domains. The latter making it impossible to capture the full extent of very large streamwise motions (VLSM). In this paper we present a detailed study of interfacial mass transfer for 180 ≤ ≤ 630 in combination with up to 200 using computational domains sufficiently large to capture the VLSM. The main objectives are (i) to study the effect of anisotropy on interfacial mass transfer, (ii) to confirm the scaling of the transfer velocity with also for larger Schmidt numbers, (iii) to show that turbulence dissipation is an important factor in the dynamical behaviour of surface divergence, (iv) to evaluate the validity of Theofanous' model for open channel flow, and finally (v) to produce detailed data on near surface scalar statistics.

Numerical method
The mathematical formulation describing the process of mass transfer (e.g. dissolved gas transfer) in turbulent open channel flow comprises the incompressible Navier-Stokes equations for the flow and advection-diffusion equations for the scalar concentration field, which in non-dimensional form read * * = 0, (1) * * + * * * = − * * + 1 where * denotes the non dimensional form of , ( 1 , 2 , 3 ) = ( , , ) are the velocity components in the ( 1 , 2 , 3 ) = ( , , ) direction, respectively, is time, is pressure, is the Kronecker delta, is the dynamically adjusted forcing term added to the momentum equation to ensure a constant flow rate, is the scalar concentration and is the Schmidt number. The bulk Reynolds number is defined as = / , where is the height of the channel and is the bulk velocity. The latter is defined by = 1 ∫ 0 , where · and · denotes averaging over time and the homogeneous horizontal ( , ) directions, respectively. The scalar concentration is normalised by where ,0 is the initial concentration in the bulk and is the saturation concentration at the surface. The full set of governing equations was solved using the in-house KCFlo code (Kubrak et al., 2013). The momentum equations were discretised using a fourth-order central scheme for diffusion combined with a fourth-order kinetic energy conserving scheme for convection. The Poisson equation for the pressure was solved using a conjugate gradient solver, with simple diagonal preconditioning. The scalar advection-diffusion equations were discretised using a fifth-order weighted essentially non-oscillatory (WENO) scheme (Liu et al., 1994) for the convection and a fourth-order accurate central scheme for the diffusion. As the scalar diffusivities are much smaller than the momentum diffusivity, a dual-mesh approach (using a finer mesh for the scalar computation than for the flow computation) was used in order to accurately resolve the Batchelor scales in a computationally efficient manner.
The solutions of both flow and scalar fields were advanced in time using the second-order accurate explicit Adams-Bashforth scheme. Three scalar advection-diffusion equations with different Schmidt numbers were solved simultaneously, enabling a direct comparison of scalar transport processes driven by exactly the same background turbulence.
The constitutive parts of our numerical code were validated in Wissink (2004) for the flow solver and in Kubrak et al. (2013) for the scalar advection-diffusion. In the latter reference, the numerical treatment of the scalar transport process was compared to analytical solutions for pure advection as well as pure diffusion, and the convergence rate of the method was rigorously established. Furthermore, the effectiveness and convergence of our dual-meshing approach was demonstrated in e.g. Kubrak et al. (2013); Herlina and Wissink (2014); Wissink and Herlina (2016).

Computational setup
As mentioned above, the present study focuses on low diffusivity mass transfer in a fully developed turbulent open channel flow. The size of the computational domain was × × in the streamwise ( ), vertical ( ) and spanwise ( ) directions, respectively (cf. figure 1). Three different domain sizes were used, 3 × × 3 , 12 × × 3 and 24 × × 6 . The total number of grid points employed in the simulations ranged from 4.72 · 10 6 to 1.22 · 10 10 . The Schmidt and bulk Reynolds numbers were varied from = 4 to 200 and = 2875 to 12000, respectively. A list of the simulations performed is provided in Table 1.
In all simulations, the surface was assumed to be flat, thereby neglecting any possible influence of waves on mass transfer. For both flow and scalar fields, periodic boundary conditions were used in the streamwise and spanwise directions. In the vertical direction, a free-slip boundary condition at the free-surface ( / = 1) and a no-slip boundary condition at the wall ( / = 0) of the channel were applied for the velocity field. For the scalars, at the surface a Dirichlet boundary condition with = , modeling the presence of the atmosphere, was employed, while at the wall a zero flux boundary condition ( / = 0) was used. To speed-up the formation of a turbulent concentration  are the number of grid points of the base mesh, while , , denote the refinement factors used for the finer scalar mesh (of the cases marked with superscript ) and Δ / is the time-window used for the flow statistics, where = / is the bulk time unit. Note that the time-window of the scalar statistics differ from Δ / , as will be explained in section 4.2.
boundary layer, the concentration field was initialised by * (1 − / ) = erfc 1 − / √︃ where * is the non-dimensional diffusion time, which was set to 10 for the smallest domain considered and 12 for the others. The spatial discretisation was performed on a non-uniform and staggered Cartesian mesh. The mesh was uniform in the homogeneous ( , ) directions and stretched in the vertical ( ) direction with refinement near the free-surface, located at ( ) = , and the wall, located at (0) = 0. The node distribution is given by where is the number of nodes in the -direction. The stretching is controlled by the parameter = + 1 − , where ( , ) was set to (2, 2) in the smallest domain, (3, 2) in the mid-sized domain and (2.3, 3) in the largest domain. The adequacy of the base mesh resolution for all simulations was confirmed by studying one-dimensional energy spectra of the velocity, as discussed in section 3.1.
As mentioned above, a dual-mesh strategy was employed to accurately resolve the evolution of the high scalar fields. The refined grid resolution was chosen such that (i) the vertical grid size Δ at the surface was less than 0.15 , where = / √ is the Batchelor scale and is the Kolmogorov scale, and (ii) the geometric mean of the grid cells Δ = 3 √︁ Δ Δ Δ was Δ < in the upper part of the domain ( / ≤ 0.55). These two conditions fulfil the criterion proposed by Grötzbach (1983), which ensures that the scalar mass transfer in all simulations is adequately resolved. One of the main findings reported in the investigations mentioned above is the appearance of very large coherent structures for large . To assess the suitability of the domain sizes used in the present simulations for capturing such large coherent structures, typical snapshots of the streamwise velocity fluctuations at / = 0.6 for G03, G06, G09 are shown in figures 2a, 2c and 2e, respectively. As can be seen in table 1, these cases have the highest for each of the three domain sizes considered. In all figures, high and low velocity elongated streaky structures, that are characteristic of open-channel flow, can be seen. The snapshots and corresponding two-point correlations  at / = 0.6 (figures 2b, 2d and 2f) indicate that the streamwise extent of these coherent structures was captured quite well in the largest domain and marginally in the mid-sized domain, but not in the smallest domain.

Domain size
For both = 365 and 630, the proper de-correlation of in the largest domain (shown for = 630 and / = 0.6 in figure 2f) was achieved in the spanwise direction for all / and in the streamwise direction for / ≤ 0.7. For / > 0.7 the minimum value for the streamwise was always smaller than 0.05, indicating a more marginal de-correlation. In contrast, for the lowest = 200, a de-correlation of in both horizontal directions was obtained for the entire depth. In the mid-sized simulations, the streamwise de-correlation was marginal for all cases, while in the smallest domain no de-correlation was achieved. As expected, a domain size of 3 × × 3 is too small to fully capture turbulent open-channel flow, even for as low as 190. Nevertheless, it will be shown in section 4.1 that for = 240 and 290 the interfacial mass transfer was found to be largely independent of the domain size.
Furthermore, in figure 2e a coherent structure of length ≈ 20 can be observed. In the literature, structures of length 10 are usually referred to as very large scale motions (VLSM), while the term large scale motions (LSM) is typically used for structures with a length of ≈ 1 − 3 . In open channel flow, the onset of the appearance of VLSM is still uncertain, i.e. to date such motions were confirmed to appear for ≥ 400 (Bauer et al., 2020), ≥ 550 (Wang and Richter, 2019) and experimentally at ≥ 700 (Peruzzi et al., 2020). The very large coherent structure seen in figure 2e is an example of such a very large scale motion that is found at = 630. Later it will be shown that in the present simulations VLSM could already be detected at = 365. Further confirmation of the appearance of VLSM based on e.g. pre-multiplied spectra and velocity fluctuations analyses, is presented in section 3.1.

Flow characteristics 3.1 Flow statistics
All simulations were started from fully developed turbulent flow fields. The shown statistics, if not stated otherwise, were obtained by averaging in the homogeneous ( , ) directions and in time (cf. table 1). Figure 3a shows that the mean streamwise velocity profiles in all simulations agree well with the law of the wall, including a logarithmic region + = 1 ( + ) + , with standard values for the von Karman constant ( = 0.40) and for the intercept ( = 5). The one-dimensional streamwise energy spectrum of the velocity components for the simulation with the largest Reynolds number (G09) is shown in figure 3b. As also observed for all other simulations, the existence of an inertial subrange, indicated by the −5/3 power law, can be clearly seen. Furthermore, no energy pile-up at high wavenumbers was observed, demonstrating that the smallest scales of motion were well resolved in all simulations. In addition, it was found that = 0.166 0.88 , which is in agreement with the results for closed channel flow as shown by Pope (2000) and Lee and Moser (2015). Figure 4 shows contours of the longitudinal velocity fluctuation spectra in the streamwise ( ) and the spanwise ( ) direction from the simulations performed in the 24 × × 6 domain (G07-G09). The energy spectra were premultiplied by the wavenumber and are shown as a function of the wavelength and distance from the wall. As expected, for increasing Reynolds number the amount of energy at smaller wavelengths was found to increase, especially near the wall of the channel. In the streamwise spectra (figures 4a, 4b, 4c), an energy peak at a wavelength of / ≈ 3 was identified in all three cases. This peak is associated with LSM. It should be noted that in the spectra, the location of the peaks with higher wavelengths becomes less accurate due to the limited size of the computational domain = 24 , = 6 . Even though the location may not be entirely correct, the energy peaks observed at 10 , which extend over virtually the whole channel height, indicate the presence of VLSM for ≥ 365 (G08, G09). When examining the spanwise spectra (figures 4d, 4e, 4f), energy peaks at / ≈ 1, which relate to LSM, were detected in all three cases. High energy values at / ≥ 2 (typical for VLSM) were observed for = 365 and = 630 when / ≥ 0.3 and 0.1, respectively. The root mean square (r.m.s.) of the velocity fluctuations for the simulations performed using the mid-and large sized domains are shown in figure 5. As expected, a decrease in the boundary layer thickness at the wall with increasing Reynolds number was observed. The increase in the peak of the streamwise component is in good agreement with the values reported in Kim et al. (1987); Pope (2000); Wang and Richter (2019). For / ≥ 0.3, the + profiles for < 365 were found to nearly collapse on one curve, while for ≥ 365 + tends to grow with . This increase in + tends to be associated to the presence of VLSM (Kim and Adrian, 1999;Del Álamo et al., 2004;Hoyas and Jiménez, 2006). Figures 5b and 5c show that for / ≥ 0.3, similar + profiles and similar + profiles were obtained for all Reynolds numbers. At the surface, vertical fluctuations are damped and the turbulent kinetic energy is redistributed in the horizontal directions, which explains the increase observed near the surface in the + and + profiles. The thickness of the surface influenced layer was determined by the distance between the surface and the location, ∞ , where ( ) = ( ). can be seen as a measure of the Reynolds stress anisotropy, where the lower limit = 1 corresponds to 2D flow, while isotropic flow would result in = 3/2. In open channel flow, the former is achieved at the free-surface. With increasing distance from the surface ( − > 0), the flow becomes less anisotropic and hence the magnitude of increases until it reaches a maximum at the edge of the surface influenced layer ( = ∞ ). Here, in all simulations, ∞ was found to be located somewhere between 0.6 and 0.75 (cf. table 2), such that 0.25 ≤ ≤ 0.4 . When comparing simulations of the same domain size, the surface influenced layer was found to decrease with increasing . Figure 6 shows the integral length scales for the velocity components in the homogeneous  Figure 7 shows contours of the streamwise-averaged -fluctuation ( / ), together with streamwise-averaged velocity vectors of the two components in the plane. Large areas with high and low speed streamwise flow can be observed, which extend almost from the wall to the free-surface of the channel. Downward moving flow is typically present in the high streamwise velocity areas, while in the low-speed areas the flow tends to move upwards. It was observed in figure 2a that these high and low speed areas extend over a significant streamwise portion of the channel, and are related to VLSM. Figure 8a shows contours of streamwise velocity fluctuations / in the plane at / = 0.5. Superimposed on this plot are small-scale vortical structures in the interval 0.5 ≤ / ≤ 1 that are visualised using the normalised -criterion,

Flow structures
where is the second invariant of the velocity gradient tensor (Hunt et al., 1988) and is the root mean square of determined using temporal and spatial averaging in , . In figure 8a, it can be seen that the vast majority of the small-scale vortical structures is present inside large low-speed streaks that extend toward the surface (cf. also figure 7). The strong concentration of these small structures in low-speed streaks can be explained as follows. Low-speed streaks form when relatively slow moving, highly turbulent flow from the lower part of the boundary layer is ejected to the upper, less turbulent part of the boundary layer (e.g. Komori et al., 1982). This larger turbulence intensity is responsible for the concentration of small scale structures observed inside the low-speed streaks.
When the small vortices approach the surface they either align with or become orthogonal to the surface, which are referred to by surface-aligned and surface-attached vortices, respectively. Figure 8b shows an instantaneous snapshot of the upper part ( / ≥ 0.9) of the channel. Most of the surface-attached vortices were found above downwelling motions corresponding to high-speed streaks. As mentioned above, in low-speed streaks most of the vortical structures are present. These structures tend to align with the surface due to the shear generated underneath the divergent flow at the surface. In these low-speed regions, the surface aligned vortices are often ring-shaped, which according to Nagaosa and Handler (2003) started their life as hairpin vortices from near the wall of the channel.
The implications of the above on interfacial mass transfer will be discussed in section 4.4. Figure 9 shows typical contours of the concentration and qualitatively depicts the interaction between the turbulent open channel flow and the scalar transport at = 7 (a),(c) and = 100 (b),(d). The thickness of the concentration boundary layer , which will be defined in equation (12) below, depends on both the Reynolds and the Schmidt number. Comparing the vertical crosssections ( , planes) shown in figures 9a and 9b, it can be clearly seen that an increase in Schmidt number at constant results in much finer concentration filaments in the bulk and a significantly reduced boundary layer thickness. Note that for a shear-free surface, the thickness scales with −0.5 (cf. figure 10a) so that =7 ≈ 3.78 =100 . This scaling is not taken into account in figures 9c and 9d, which shows surface-parallel planes at a distance of 0.0003 to the surface, corresponding to 0.019 and 0.071 for = 7 and 100, respectively. Hence, the range of the scalar concentrations in figure 9d had to be increased by a factor of ≈ √︁ 100/7 = 3.78 to obtain similar contours. Nevertheless, small differences can still be observed locally due to differences in diffusion in the horizontal directions.

Instantaneous results
Note that the above reduction in at a fixed is due to the increase in interfacial mass transfer resistance with increasing Schmidt number. At a fixed , the increase in turbulence in the bulk associated with an increase in results in improved mixing with a reduction of , which in this case promotes mass transfer.

Statistics of scalar transport
As mentioned in section 3.1 all simulations were started from fully developed turbulent flow fields with the scalars initialised by (5). After a transient period needed to ensure that the scalar statistics were quasi-steady, scalar averaging was carried out using a time-window of Δ / (see table 2).
The thickness of the diffusive concentration boundary layer is identified using As illustrated in figure 10a for simulation G07, in all simulations was found to scale with −0.5 which is in agreement with the theoretical prediction for a shear-free interface (e.g. Ledwell, 1984). Included in this plot are the thicknesses of the Kolmogorov sublayer and the Batchelor sublayer at the interface. Except for = 7, it was found that < < , which is in agreement with Herlina and Wissink (2014, hereafter HW14). Figure 10b shows the variation of √ / with the bulk Reynolds number . It can be seen that for = 4000 and 5000, the normalised boundary layer thickness is nearly independent of the computational domain size. The best fit through the data points was found to be √ / ∝ −0.67 . This will be discussed further in section 4.3, where the scaling will be linked to the transfer velocity. Figure 11 shows normalised mean vertical profiles of the concentrations, the r.m.s. of concentration fluctuations and the mass fluxes at various Reynolds numbers (G07, G08 and G09). As discussed above, the boundary layer thickness depends on both the molecular diffusivity ( ) and the Reynolds number ( ). Thus, it is expected that the vertical profiles of the normalised mean scalar quantities exhibit self-similarity when the vertical ( − ) direction is normalised by . (Note that the profiles for the different Schmidt numbers also collapse). In all simulations, the magnitude of / was found to be virtually identical to the distance between the surface and the location at which the r.m.s. of the concentrations reaches its maximum. Hence, the peak in figure 11b is located at ( − )/ = 1. The maximum /( − ) values were ≈ 0.3, which is in agreement with previous numerical (Magnaudet The total averaged vertical mass flux, = + , comprises a diffusive component and a turbulent component ( ) = .
(15) Figure 11c illustrates that dominates at the surface, as is damped due to the two-dimensionality imposed by the free-slip boundary condition. With increasing distance to the surface the contribution of to the total mass flux reduces, while at the same time the contribution of becomes increasingly important until it entirely dominates . It can also be seen that at ( − )/ 0.65 the diffusive and turbulent mass fluxes become equal. Furthermore, it was found that in the range 2 ( − )/ 10, the turbulent mass flux agrees reasonably well with the mass flux at the surface ( ).
From the surface ( = ) down to a depth of ( − ) ≈ 2 , the normalised mean vertical profiles of the scalar quantities shown in figure 11 were found to nearly collapse with the profiles reported in Herlina and Wissink (2019, hereafter HW19) (not included in the plot), where the interfacial mass transfer was driven by isotropic-turbulence diffusing upwards from the opposite boundary. For the highest Reynolds number case (G09), the agreement between the open channel profiles produced here and in the case driven by isotropic turbulence remains reasonably well up to a distance of at least 10 from the free-surface. This indicates that the strong anisotropy present in the velocity field of turbulent open channel flow (which increases with increasing Reynolds number) does not affect the shape of the normalised profiles at least up to the aforementioned depth.

Scaling of mass transfer velocity
The local, instantaneous transfer velocity is defined as For open channel flow with no wind-shear, the mean transfer velocity is commonly parameterised using the bulk Reynolds number or the friction Reynolds number . In figure 10 it was found that / scales with −0.5 and −0.67 . Hence, as the normalised transfer velocities / obtained here scale with −0.5 and −0.33 . The former scaling was found to be valid only for clean (surfactant-free) surfaces (Wissink et al., 2017). The power = −0.33 in the latter scaling is somewhat smaller than the one found in the DNS of Nagaosa and Handler (2012, hereafter NH12), who obtained ∝ −0.25 in small-box simulations at = 1. When using the friction Reynolds number rather than , / was found to scale with −0.23 , see figure 12a. This power of = −0.23 is similar to the powers found in open channel flow experiments of Moog and Jirka (1999); Lau (1975); Gulliver and Halverson (1989) where ranges between −0.19 and −0.29. Theofanous et al. (1976) proposed to use the turbulent Reynolds number , defined by as a measure of turbulence characteristics that is independent on the way that turbulence is generated, e.g. by wind-shear, bottom-shear, or buoyancy. In (19), ∞ is the turbulence intensity and ∞ is the streamwise integral length scale, taken here as and evaluated at the edge of the surface influenced layer = ∞ , see table 2. Note that could not be computed for the smallest domain size since no decorrelation of in the streamwise direction was achieved. As can be seen   (19). Note that all values were obtained by averaging over a time-window Δ / in which the scalar statistics are quasi-steady, see table 1. in table 2, our data fall near and above the critical turbulent Reynolds number = 500, that in the two-regime model of Theofanous et al. (1976) separates the large-eddy ( ≤ 500) from the small-eddy ( ≥ 500) regime. Figure 12b shows that in our open channel flow simulations 0.5 / ∞ scales with 0.35 −0.25 when > 500. The normalised / ∞ values and hence its scaling with are in close agreement with the results found in the absence of mean shear by HW19, which confirms that for > 500 the turbulent characteristics are dominated by small eddies (i.e.

Flow structures and mass transfer
Previous authors have shown that the surface divergence model of McCready et al. (1986), ∝ √ , where is the r.m.s of the surface divergence = − / | = , generally provides a good prediction of the mass transfer velocity, in open channel flow (Sanjou et al., 2017;Turney and Banerjee, 2013) as well as in turbulence without mean shear (HW19, McKenna and McGillis, 2004). The approximately linear variation of √︁ /( ) with √︁ / , shown in figure 13, confirms that the surface divergence model performs reasonably well also in the present simulations. For interfacial mass transfer driven by isotropic turbulence, it was found that despite the good correlation between and the mean transfer velocity , the correlation between the local, instantaneous transfer velocity ( , , ) and the instantaneous surface divergence tends to deteriorate with (Herlina and Wissink, 2019). As explained below, in our present simulations no such deterioration of the time-averaged correlation ( ( , , ), ( , , )) with an increase in turbulence level was observed. Hereafter, the aforementioned correlation will be referred to as ( , ).
Figure 14a depicts the variation of ( , ) with obtained in the present simulations, combined with results from HW19. It can be seen that initially, for isotropic turbulence ( = 20) the correlation ( , ) tends to reduce with increasing and appears to reach a plateau around ≈ 0.8. In contrast, in the large-box simulations of the present open channel flow study at = 16, the correlation ( , ) approaches a somewhat smaller value of approximately 0.72 from below. As can be seen in figure 14b, a similar improvement of ( , ) with was also observed in the present small-box simulations and, somewhat less clear, also in the mid-sized simulations. Note that this is in agreement with the trend observed in small-box simulations reported in NH12 at = 1. The good agreement achieved between the small and large box simulations is likely due to the incomplete de-correlation in the streamwise direction in the small box, which mimics the presence of large streamwise structures as can only be fully captured in much larger computational domain. Contrastingly, in the mid-sized box a marginal streamwise de-correlation was achieved causing the maximum size of the streamwise structures to become significantly smaller than in the large box simulations. Note that a more detailed study of the effect of the domain size on mass transfer is beyond the scope of this paper.
The Schmidt number effect on the correlation coefficient ( , ) is shown in figure 14c. The correlation was found to decrease with increasing Schmidt number. This can also be seen in figure 14b when comparing the present results with NH12. For = 1 momentum and scalar diffusivity are identical and the correlation is not influenced by differences in diffusive time scales. With increasing scalar diffusivity, the difference in diffusive time scales between momentum and scalar increases, resulting in a reduced correlation between and . To investigate the slight increase of ( , ) with the Reynolds number -observed in open channel flow simulations -in more detail, table 3 reports the correlation ( , ) separately for low-speed regions ( < − 0.5 ( )) and high-speed regions ( > + 0.5 ( )), where ( ) denotes the standard deviation of . It was found that ( , ) is approximately independent of in the low-speed regions, while it increases with in the high-speed regions. As a result, an overall increase in the correlation was observed with increasing turbulence Reynolds number. This increased correlation was linked to the distribution of surface parallel vortical structures (SPVS) close to the surface. These structures are responsible for the vertical mixing of fluid close to the surface and are relatively stable. In Run low-speed high-speed  Table 3: Correlation coefficient, ( , ), at = 16 for = 3200 (G07), = 6400 (G08) and = 12000 (G09), together with the corresponding fractions of the free-surface area covered by high ( > + 0.5 ( )) and low ( < − 0.5 ( )) speed regions.
an idealised case, they would look like a surface parallel ring vortex that continuously brings up unsaturated fluid in its centre to the surface leading to local maxima in and , while on its outer side the now partially saturated fluid is returned to the bulk leading to local minima in and as illustrated in figures 15a and c. Hence, these SPVS typically contribute to a good correlation ( , ). For low , these SPVS are mainly present in the low-speed regions, while for larger this distribution was observed to become more uniform (possibly due to increased mixing), see figures 15b and d. As a result, ( , ) was found to increase inside the high-speed regions, and hence overall. The above is consistent with the observation that ( , ) was found to be markedly higher in the low-speed regions than in the high-speed regions.
By comparing figures 15a and 15c it can also be seen that the typical size of the pattern (and hence the surface divergence pattern) reduces with increasing Reynolds number. This is consistent with the results shown in figure 16, where the integral length scale was found to decrease with increasing . Because of the reasonably good correlation between surface divergence and mass transfer, it is inferred that turbulence structures of similar size as the integral length scale of the surface divergence are of key importance to the overall mass transfer. Additionally, above it was found that / ∞ ∝ −0.25 , which indicates that mass transfer is significantly influenced by very small scale turbulent structures. To explore this further, first the importance on the vertical mass transfer of turbulence structures of various size is investigated in figure 17. Here, the normalised premultiplied spectrum of is shown as a function of the normalised wavenumber . Also indicated in the figure are the locations of the wavenumber associated with ∞ ( ∞ = 2 / ∞ ) and the wavenumber associated with ( = 2 / ). It can be seen that in all three cases, is associated with much smaller turbulence structures (higher wavenumbers) than ∞ . The peaks in the premultiplied energy spectra were located in between the wavenumbers ∞ and . Also, as expected, the peaks tend to move towards larger wavenumbers (smaller scales) at higher turbulence levels (from ≈ 10 for = 465 to ≈ 20 for = 1581 to ≈ 30 for = 2833). For the cases > 500, the peaks were located at wavenumbers significantly larger than ∞ . Second, in figure 18, for each simulation the approximate location of the aforementioned peaks in the streamwise 1D turbulent energy spectrum is shown. It can be seen that the wavenumber associated with each peak is located close to the Taylor microscale ( ), where dissipative effects are no longer negligible.
Third, also shown in figure 18 are the exact locations of . In each of the simulations, was found to be markedly larger than the wavenumber associated with the Taylor microscale , and hence small scale turbulence dissipation significantly affects surface divergence. This is in agreement with the theoretical considerations of McCready et al. (1986) and further evidences the association of the surface divergence integral length scale ( ) with small rather than large turbulence scales. The above confirms the importance of small (dissipative) scales on interfacial

Conclusions
Mass transfer across the clean, flat surface of a fully-developed turbulent open channel flow was studied in detail using large-scale direct numerical simulations at various Reynolds (180 ≤ ≤ 630) and Schmidt numbers (4 ≤ ≤ 200). Turbulent flow statistics were found to be in good agreement with the literature (e.g. Wang and Richter, 2019;Peruzzi et al., 2020). By analysing pre-multiplied spectra and velocity fluctuations, very large scale motions (VLSM), detected previously for ≥ 400, were identified already for ≥ 365 in the simulations with the largest computational domain (24 × × 6 ). In agreement with the literature, small scale vortical structures were mainly found to be located inside the large (long-lived) low-speed streaks. Towards the surface, these vortices were observed to become either aligned or orthogonal to the surface. By studying the near-surface region in detail, it was found that when approaching the surface these small scale vortices, which originate from the lower part of the turbulent boundary layer, become more uniformly distributed in the spanwise direction at high . Close to the surface, the vertical profiles of ( − )/( − ), /( − ), and / , / were found to collapse after normalising the -coordinate using ( − )/ . In addition, these profiles were also found to collapse with those reported in HW19, which  Table 1). The arrows indicate the locations of associated with the approximate peaks in * , while the lines indicate : , −5/3 slope; , location of = 2 / ; , location of = 2 / .
were obtained using isotropic (rather than anisotropic) turbulence diffusing from below. The normalised mean transfer velocity was found to scale with −0.5 for Schmidt numbers up to = 200, which is in agreement with the theoretical prediction for free-slip (clean) surfaces. The scaling of the normalised with Reynolds number was found to vary with −0.5 / = −0.33 , −0.5 / = −0.23 and −0.5 / ∞ = ∞ −0.25 . The two former scalings are relatively close to the data reported in literature. For > 500, the latter scaling (with ∞ = 0.35) is in very good agreement with HW19 and the small eddy model of Banerjee et al. (1968); Lamont and Scott (1970). Hence, the ansatz of Theofanous et al. (1976) that the scaling of −0.5 / ∞ with −0.25 is independent on the way the turbulence is generated was found to be supported for the case of bottom-shear induced turbulence which can be replaced by isotropic turbulence diffusing from below without affecting the scaling. Note that in the open channel flow analysis ∞ and ∞ are determined at the edge of the surface influenced layer ( = ∞ ), which was identified as the location where ( ) = ( + + ) /( + ) is maximum.
Furthermore, it was found that √︁ /( ) scales with 0.47 √︁ / , which supports the validity of the surface divergence model for all and considered. A detailed study of the correlation ( , ) of the local mass transfer velocity and the surface divergence showed that it improves with increasing Reynolds number, but deteriorates with increasing . Using conditional averaging, it was found that improvements in ( , ) with increasing Reynolds number were entirely due to increases in the correlation inside high-speed regions, while the correlation in the low-speed regions remained almost constant (and significantly higher than in the high-speed regions). The observed improvement in correlation was due to the surface parallel vortical structures (SPVS) -which are responsible for vertical mixing -becoming more uniformly distributed close to the surface with increasing Reynolds number.
Finally, the locations of the wavenumbers , and associated with the Taylor microscale, the surface divergence integral length and the approximate location of the peaks of the premultiplied spectral density, respectively, were determined for = 465, 1581, 2833. It was found that and almost coincide, while was significantly larger for all Reynolds numbers and well outside of the inertial range. Such that it can be concluded that the surface divergence, and hence the mean transfer velocity (which is well correlated with ), are notably affected by the turbulence dissipation and hence the small length scales. This again supports the validity of the small eddy model for the prediction of for 500 in open channel flow.