Direct numerical simulation of gas transfer across the air–water interface driven by buoyant convection

A series of direct numerical simulations of mass transfer across the air–water interface driven by buoyancy-induced convection have been carried out to elucidate the physical mechanisms that play a role in the transfer of heat and atmospheric gases. The buoyant instability is caused by the presence of a thin layer of cold water situated on top of a body of warm water. In time, heat and atmospheric gases diffuse into the uppermost part of the thermal boundary layer and are subsequently transported down into the bulk by falling sheets and plumes of cold water. Using a specifically designed numerical code for the discretization of scalar convection and diffusion, it was possible to accurately resolve this buoyant-instability-induced transport of atmospheric gases into the bulk at a realistic Prandtl number ( $\mathit{Pr}=6$ ) and Schmidt numbers ranging from $\mathit{Sc}=20$ to $\mathit{Sc}=500$ . The simulations presented here provided a detailed insight into instantaneous gas transfer processes. The falling plumes with highly gas-saturated fluid in their core were found to penetrate deep inside the bulk. With an initial temperature difference between the water surface and the bulk of slightly above $2$ K, peaks in the instantaneous heat flux in excess of $1600~\text{W}~\text{m}^{-2}$ were observed, proving the potential effectiveness of buoyant-convective heat and gas transfer. Furthermore, the validity of the scaling law for the ratio of gas and heat transfer velocities $K_{L}/H_{L}\propto (\mathit{Pr}/\mathit{Sc})^{0.5}$ for the entire range of Schmidt numbers considered was confirmed. A good time-accurate approximation of $K_{L}$ was found using surface information such as velocity fluctuations and convection cell size or surface divergence. A reasonable time accuracy for the $K_{L}$ estimation was obtained using the horizontal integral length scale and the root mean square of the horizontal velocity fluctuations in the upper part of the bulk.

A series of direct numerical simulations of mass transfer across the air-water interface driven by buoyancy-induced convection have been carried out to elucidate the physical mechanisms that play a role in the transfer of heat and atmospheric gases. The buoyant instability is caused by the presence of a thin layer of cold water situated on top of a body of warm water. In time, heat and atmospheric gases diffuse into the uppermost part of the thermal boundary layer and are subsequently transported down into the bulk by falling sheets and plumes of cold water. Using a specifically designed numerical code for the discretization of scalar convection and diffusion, it was possible to accurately resolve this buoyant-instability-induced transport of atmospheric gases into the bulk at a realistic Prandtl number (Pr = 6) and Schmidt numbers ranging from Sc = 20 to Sc = 500. The simulations presented here provided a detailed insight into instantaneous gas transfer processes. The falling plumes with highly gas-saturated fluid in their core were found to penetrate deep inside the bulk. With an initial temperature difference between the water surface and the bulk of slightly above 2 K, peaks in the instantaneous heat flux in excess of 1600 W m −2 were observed, proving the potential effectiveness of buoyant-convective heat and gas transfer. Furthermore, the validity of the scaling law for the ratio of gas and heat transfer velocities K L /H L ∝ (Pr/Sc) 0.5 for the entire range of Schmidt numbers considered was confirmed. A good time-accurate approximation of K L was found using surface information such as velocity fluctuations and convection cell size or surface divergence. A reasonable time accuracy for the K L estimation was obtained using the horizontal integral length scale and the root mean square of the horizontal velocity fluctuations in the upper part of the bulk.

Introduction
During the past decade, gas transfer across the air-water interface has received increasing interest because of its importance to the global greenhouse gas budget. At the air-water interface, heat and gas fluxes are dominated by molecular diffusion, which can be enhanced significantly by convective mixing. In nature, three sources of near-surface convective mixing can be identified. The first and most studied one is mixing by wind-shear-induced turbulence at the water surface. The second is turbulence generated at the bottom of a stream or by ocean currents that subsequently diffuses upwards to the surface where it enhances turbulent mixing. The third and least studied mechanism is due to a buoyant-convective instability, which is responsible for mixing during low-wind-speed conditions. In the ocean, free convection due to buoyant instabilities is often caused by vertical temperature and/or salinity gradients. This free convection induces near-surface turbulence, which enhances the exchange of heat and greenhouse gases (in particular carbon dioxide) between the atmosphere and the oceans.
Wind-shear-induced mixing has attracted significant attention because turbulence introduced immediately at the water surface is a very effective mechanism for promoting gas transfer. Especially larger wind speeds can induce significant wave growth and wave breaking, which further aid the mixing process (MacIntyre et al. 1999). Hence, in most cases the predictive estimation of heat and greenhouse gas fluxes only takes wind shear into account, which often results in an overor under-prediction of the actual gas flux, especially in low-wind conditions (e.g. Nightingale et al. 2000;Crusius & Wanninkhof 2003;Smith et al. 2011).
Sheltered inland water bodies and equatorial ocean regions are examples of regions with low-wind-speed conditions that typically have large air-water CO 2 concentration gradients. These regions contribute significantly to the global CO 2 budget (Feely et al. 2002), both as a major source of climatically important gases (in lakes and reservoirs: Cole et al. 2007;MacIntyre et al. 2010) as well as a major sink for CO 2 (in the equatorial ocean: Ishii et al. 2014). MacIntyre et al. (2010) showed that for wind speeds lower than 5 m s −1 the gas transfer velocity K L is independent of wind speed and the buoyancy-driven flux dominates. In the ocean Soloviev & Klinger (2001) and Rutgersson, Smedman & Sahlee (2011) have reported that in low-wind-speed conditions buoyant convection due to surface cooling contributes significantly to the surface gas exchange.
Another important example of buoyant-convective mixing in nature is the absorption of oxygen into lakes or reservoirs driven by surface cooling. At night, the gassaturated surface layer is cooled and becomes slightly heavier than the water underneath. Small fluctuations may trigger an instability, resulting in cold oxygensaturated surface water plunging down and being replaced by warmer unsaturated fluid from below. This process results in the re-aeration of sheltered bodies of water, which is needed for biodegradation of organic wastes.
In order to improve current models used for the prediction of gas fluxes, there is a need for a detailed understanding of buoyant-convective mixing processes. To explore the effect of these processes requires an in-depth investigation of the interaction between the near-surface flow field and the interfacial mass flux, which -in contrast to the relatively abundant empirical quantification of the gas transfer velocity -is currently still lacking.
The gas flux across the air-water interface j z is defined by where K L is the transfer velocity, C s is the saturated gas concentration at the interface and C b is the gas concentration in the bulk. The transfer velocity K L is governed by the complex interaction of molecular diffusion and turbulent transport on either side of the surface and also depends on surface conditions as well as the chemical reactivity of the dissolved gas. At a constant temperature Henry's law states that at equilibrium the solubility of an atmospheric gas in water increases linearly with its partial pressure in air, p a , so that C s = p a /H, where H denotes Henry's constant. (Note that the evaluation of C s is important in order to accurately determine j z .) The partition coefficient of a gas between the atmosphere and water is defined by K = C g /C s , where C g is the gas concentration in the air. (For an ideal gas, K is related to H by K = H/RT, where R is the gas constant and T is the temperature.) As shown by Liss (1973), the total transfer velocity K L can be related to the individual gas transfer coefficients (k L and k g for the liquid and gas phase, respectively) by This 'resistance-in-series model' shows that the ratio k L /K k g determines whether the concentration boundary layer on the air or on the water side controls the gas transfer process (offers more resistance). With increasing solubility (which may be caused by a decreasing temperature) the control of the atmospheric gas transfer process will gradually shift to the air side. Environmentally important gases typically have a low (O 2 , CO, NO) to moderate (CO 2 ) solubility. Because of this and their significantly lower diffusivities in water than in air (hence k L k g ), for such gases k L /K k g is small and their transfer process at the air-water interface is concentrated in an extremely thin boundary layer (10-1000 µm) on the liquid side (see Liss 1973;Jähne & Haussecker 1998). It should be noted that CO 2 is a chemically reactive gas in water and can be hydrated at the water surface, which chemically enhances the transfer velocity (see Wanninkhof & Knox 1996).
The small size of the concentration boundary layer on the liquid side makes it extremely difficult to carry out direct measurements of the diffusive (−D ∂C/∂z) and turbulent (w C ) vertical mass fluxes, where C is the concentration of the solute, w is the velocity in the vertical (z) direction and D is the molecular diffusion coefficient (double primes denoting fluctuations with respect to statistical averages). Together, the diffusive and turbulent mass fluxes determine the total mass flux j z in the vertical direction. The mass flux averaged over horizontal planes is a function of depth (z) and time (t), as given by Because it is difficult to perform measurements inside the tiny layer of gas-saturated fluid at the water surface, previous measurements focused mainly on measuring the overall transfer velocity K L . Equations (1.1) and (1.3) imply that to gain a fundamental understanding of the overall transfer velocity K L it is necessary to resolve both the diffusive and turbulent mass transfer terms. Numerous studies proposed empirical models which relate the transfer velocity K L to global parameters such as the bulk velocity of a stream, the bottom slope and friction (e.g. Tsivoglou & Wallace 1972;Chu & Jirka 2003). As mentioned above, for the case where buoyancy is likely to dominate, K L has been mostly related to wind shear alone (e.g Liss & Slater 1974;Nightingale et al. 2000;Crusius & Wanninkhof 2003). This has led to considerable inaccuracies in the estimation of K L , especially for wind speeds below 5 m s −1 .
Furthermore, previous studies on buoyant-convective-induced flow were mainly concerned with the heat flux and flow patterns in a fluid layer confined by a heated plate at the bottom and a cooled plate at the top (Amati et al. 2005;Verzicco & Sreenivasan 2008;Schumacher 2009), where no-slip boundary conditions apply. In many applications in nature such conditions do not apply, and as the current work is motivated by the interfacial mass transfer at a gas-liquid interface, a flat surface assumption with a free-slip boundary condition is far more appropriate to model the water surface.
While Rayleigh-Bénard convection between two solid boundaries has been studied quite extensively, buoyant convection induced by a cooled water surface has received much less attention. Spangenberg & Rowland (1961) were among the first who visualized the surface and vertical temperature patterns during evaporative cooling using the schlieren imaging technique. They observed that the cool surface water forms a net-like pattern when it collects along lines as it plunges down, forming vertical sheets. Similar observations were made by Katsaros et al. (1977). The infrared images of Volino & Smith (1999), however, did not reveal any net-like pattern at the surface. Instead, they only observed regions with vortical structures, which could be explained by their relatively high Rayleigh number and/or the presence of surface contamination. Flack, Saylor & Smith (2001), for instance, observed that surface contamination can significantly affect the surface temperature pattern. A detailed investigation on how the surface pattern in buoyancy-driven flow is affected by surface contamination, Rayleigh number or Schmidt number is still lacking.
Only a limited number of quantitative investigations of this problem have been performed. Deadorff, Willis & Lilly (1969) measured the mean temperature profile of the bulk water and horizontal length scales, while Katsaros et al. (1977) conducted temperature measurements within the boundary layer. The latter found that the boundary layer profile was highly nonlinear due to unsteadiness. Bukhari & Siddiqui (2006) used particle image velocimetry (PIV) to measure the turbulent structure below the surface, while Volino & Smith (1999) performed simultaneous temperature and PIV measurements at several horizontal and vertical cross-sections of the fluid. They observed that the surface temperature did not always correlate well with the velocity field.
While all buoyant-convective flow investigations studied the evolution of flow structures, only a limited number related buoyant convection to gas exchange across the air-water interface. For instance, Eugster et al. (2003) and MacIntyre et al. (2010) measured the CO 2 exchange rate between the atmosphere and a lake. Their quantification of the global K L confirmed the importance of buoyancydriven convective mixing in relation to interfacial gas transfer. Schladow et al. (2002) performed gas transfer measurements in the laboratory using laser-induced fluorescence (LIF), where the buoyant-convective instability was triggered by surface cooling. The synoptic two-dimensional (2D) PIV-LIF measurements of Jirka, Herlina & Niepelt (2010) showed that the oxygen-rich fluid that is transported by cold plumes from the surface into the bulk does initially form mushroom-like structures, which eventually deform into an anchor-like shape. To further elucidate the actual dynamics of the falling plumes, three-dimensional (3D) visualizations are needed. Also, the accuracy of the laboratory measurements was insufficient to resolve the minute fluctuations in C and w , which is necessary to quantify directly the turbulent mass flux and to relate the near-surface turbulence to the mixing properties.
Direct numerical simulations (DNS) related to heat/mass transfer across the free surface of an open channel flow have been performed by various research groups (Nagaosa 1999;Yamamoto, Kunugi & Serizawa 2001). These simulations showed a correlation between the vortices ejected from the bottom region and the near-surface concentration field, and proved the usefulness of DNS in obtaining detailed statistics of the scalar fields. Other researchers, such as Kunugi & Satake (2002), simulated the gas transfer promoted by wind-shear-generated turbulence. They were able to establish a relationship between low-speed streaks and the transfer of CO 2 from the interface into the bulk region. Magnaudet & Calmet (2006) used large-eddy simulations (LES) to analyse the concentration and flow field in the near-surface region of an open channel flow at a high Reynolds number. They reported large-scale turbulent upwelling and downwelling motions. Because of the high Reynolds number and the relatively large size of the computational domain, the numerical resolution of the very thin diffusive sublayers was only marginal. Kermani & Shen (2009) used DNS to study characteristics of interfacial mass transfer driven by free-surface turbulence. They focused on the quantification of the surface age related to the surface renewal models pioneered by Higbie's penetration theory (Higbie 1935) and further elaborated by Danckwert's random surface renewal model Danckwerts (1951). These models assume that the turbulence in the bulk region of the fluid transports the low-gas-saturated fluid portion up to the surface where gas transfer takes place. Hasegawa & Kasagi (2009) performed a hybrid DNS/LES numerical simulation of a coupled air-water turbulent flow and associated mass transfer to investigate the relation between the interfacial mass transfer and the surface divergence. Khakpour, Shen & Yue (2011) performed DNS with turbulent shear flow to study the effect of a surfactant-contaminated surface on the gas transfer process. It should be noted that, due to limitations in computing capacity, the above DNS studies were mostly conducted at low Schmidt numbers (<10). Recently, Herlina & Wissink (2014) carried out a series of DNS to study the effect of bottom-shear-induced turbulence on interfacial gas transfer at realistic Schmidt numbers. To resolve the scalar flow field for Schmidt numbers up to Sc = 500, a dual-mesh strategy was implemented together with a fifth-order-accurate weighted essentially non-oscillatory (WENO) scheme for scalar convection (see Kubrak et al. 2013).
This paper reports on the results obtained in a series of DNS calculations of interfacial gas transfer driven by buoyant convection. The heat and mass transfer at a Prandtl number of Pr = 6 and Schmidt numbers up to Sc = 500, which are relevant for the transport of atmospheric gases in water, were calculated using a slightly adapted version of the code by Kubrak et al. (2013). The complex mass transfer process is elucidated by fully resolving the thin concentration boundary layer at the air-water interface as well as the thin filaments present in the deeper bulk region. This allowed us to carry out an in-depth investigation into the dynamics of the falling sheets and plumes that are responsible for the heat and mass transfer from the surface into the bulk. Also, the growth rate of the convection cells in time was investigated as well as the relation between the flow field, the temperature and the gas concentration. Finally, the dependence of the transfer velocity K L on the Schmidt number, Rayleigh number, typical surface information (e.g. the average size of the convection cells and the surface divergence) and typical bulk information (e.g. bulk integral length scale and root mean square (r.m.s.) of the velocity) was assessed.

Numerical method
The problem under consideration is the gas transfer across the air-water interface driven by buoyant convection caused by temperature differences between the water surface and the bulk. The initial temperature difference is assumed to be small; in our case T b,0 − T s was 2.2 K, where T b,0 and T s are the temperatures at the bulk and the surface, respectively. In this region the density is approximately a linear function of the temperature and the Boussinesq approximation is used to describe the effect of the small variations in density on the fluid flow.
The fluid motion is governed by the non-dimensional Navier-Stokes equations. As in Balachandar (1992) the flow is driven by an unstable temperature gradient and the non-dimensionalization is carried out using a length scale L and a velocity scale U = κ/L, where κ = 0.146 × 10 −6 m 2 s −1 is the thermal diffusivity of water at 298.15 K.
Note that this study only concerns the initial development of the Rayleigh instability, hence the length scale is chosen to be fixed at a value of L = 0.00925 m that is independent of the actual depth of the computational domain. The resulting continuity equation reads ∂u i ∂x i = 0 (2.1) and the momentum equations are given by where x 1 = x, x 2 = y are the horizontal directions and x 3 = z is the vertical direction, u 1 = u, u 2 = v and u 3 = w are the velocities in the x, y and z directions, respectively, p is the pressure and t denotes time. The last term on the right-hand side represents the buoyancy force in the z direction where δ i3 is the Kronecker delta, is the non-dimensional temperature, Pr = ν/κ = 6 is the Prandtl number corresponding to the ratio of the momentum and the thermal diffusivities of water at 298.15 K and is the macroscale Rayleigh number, in which α is the thermal expansion factor in K −1 , g = 9.81 m s −2 is the gravitational acceleration and T = (T b − T s ) is the temperature difference between the bulk and the surface. Note that because of the adiabatic boundary conditions at the bottom, T b depends on time. In (2.2), Ra L | t=0 is calculated using The non-dimensional transport equation of the temperature T * is given by Similarly, the passive scalar transport is governed by the dimensionless 3D convectiondiffusion equation of the non-dimensional scalar concentration C * , with C b,0 and C s defined analogously to T b,0 and T s , respectively. In the remainder of the paper, T and C rather than T * and C * will be used to refer to the non-dimensional temperature and scalar concentration, respectively. To improve clarity, time will be displayed in seconds. Unless mentioned otherwise, T b and C b will be evaluated at d/2, where d is the domain depth. For the simulations, the full set of governing equations mentioned above was solved using the code described in Kubrak et al. (2013). The solver was developed specifically with the aim to resolve the details of low-diffusivity scalar transport, which is governed by a very thin boundary layer at the air-water interface. For the scalar transport a fifth-order-accurate WENO scheme (Liu, Osher & Chan 1994) for the convection was combined with a fourth-order-accurate central method for the diffusion in order to accurately capture high concentration gradients. For the time integration of the scalar a three-stage Runge-Kutta scheme was employed. To solve the incompressible 3D Navier-Stokes equations for the fluid flow, a central finite-difference approach with a fourth-order-accurate discretization of the diffusion and a fourth-order-accurate kinetic-energy-conserving discretization of the convection (Wissink 2004) was used. The spatial discretization was performed on a non-uniform mesh using a staggered variable arrangement, where the scalar concentrations, temperature and pressure were all defined in the middle of the grid cells. The time integration was performed using the second-order-accurate Adams-Bashforth method. The Poisson equation for the pressure, which is obtained after substituting the discretized momentum equation into the continuity equation, was solved using a conjugate gradient solver with a simple diagonal preconditioning.
As the scalar diffusivity is much smaller than the momentum diffusivity, significantly more grid points are needed to accurately resolve the evolution of the scalar. To deal with this in a computationally efficient manner, a dual-meshing strategy was employed in which the scalar concentration field was solved on a finer mesh than the velocity field.

Overview of simulations
The set-up of the numerical simulations was inspired by earlier experiments performed at Karlsruhe Institute of Technology (KIT) (Jirka et al. 2010). In the experiments, the water temperature was held fixed at 293.15 K and oxygen was stripped out. The body of water was subsequently exposed to a cold gas at its surface, leading to a temperature difference between bulk and surface. As a result, a thin thermal boundary layer of cool water formed adjacent to the air-water interface, with an even thinner layer of gas-saturated water at the top. As the cool water is slightly heavier than the water in the bulk region below, small disturbances may trigger a Rayleigh instability by which the cold water (partially saturated with gas) from above will penetrate the bulk region below, thereby interchanging unsaturated water with saturated water. While the experiments were performed in a 50 × 50 × 65 cm 3 tank with a water depth of 42 cm, the computational domain in the simulations only covered part of the physical domain, as illustrated in figure 1. The domain and grid sizes used in the various simulations are given in table 1. The mesh is stretched in the z direction in order to obtain a finer resolution near the interface. The node distribution is given by  Run In all simulations the Prandtl number was set to Pr = 6 (typical for water at 298.15 K), f RS is the refinement factor for the scalar mesh, α T is the expansion factor and Ra L | t=0 is the initial macroscale Rayleigh number. Note that the conclusions presented in this paper are based on results mostly from BC3.
for k = 1, . . . , n z − 1, with where n z is the number of nodes in the z direction. The stretching is controlled by the parameter σ , which is set to σ = 3 in all simulations. As mentioned above, a dual-meshing strategy is employed in which the gas concentration field is solved on a finer mesh (dual grid) than the velocity and temperature fields. The boundary conditions used in the simulations are illustrated in For water at 298.15 K this corresponds to temperature differences between the water surface and the bulk of T = 1.48 K and 2.22 K, respectively. Initially at t = 0 the velocity field was set to zero. Below the top boundary, the concentration field was set to C = 0 and the temperature field to T = 1. After allowing the thermal and concentration boundary layers to develop for t = 9.6 s, small random disturbances (uniformly distributed between T/ T = 0 and 0.020) were added to the temperature field to trigger the instability.

Grid refinement study
According to the Grötzbach criterion (Grötzbach 1983), a scalar is sufficiently resolved if the geometric mean of the grid size, where L B is the Batchelor scale defined by for the temperature and mass transport, respectively. It can be seen in table 2 that in all simulations the geometric mean of the grid size is less than πL B , and near the air-water interface the vertical size of the grid cells is significantly smaller than the thermal and concentration boundary layer thicknesses, so that all employed meshes fulfill the Grötzbach criterion.
To further ensure an accurate resolution of all near-surface details of the buoyantconvective instability down to the Batchelor scale, a rigorous grid refinement study for the base mesh was performed. To facilitate this study, in a separate simulation, similar to VR2 (see table 1), random disturbances were added to the temperature field at t = 9.6 s. This simulation was subsequently run for a further 1.9 s after which the resulting perturbed temperature field was saved and used as the initial condition (starting the simulations at t = 11.5 s) for the temperature in VR1, VR2 and VR3. This procedure allows a direct comparison of instantaneous flow and temperature fields calculated on different meshes. VR1 has the coarsest mesh with 280 × 280 × 184 points, the typical mesh used in VR2 is approximately 1.4 times  Hence, it can be concluded that the mesh used in VR2 is sufficiently fine to resolve the buoyancy-induced instability.
As the scalar diffusivity is much smaller than the momentum and thermal diffusivities, a refined mesh is used to be able to accurately resolve the scalar concentration for Schmidt numbers up to Sc = 500. To test which refinement factor is needed, a series of simulations (SR1, SR2, SR3) with the same 400 × 400 × 256 base mesh and refinement factors of f RS = 1, 2, 3 were performed. In figure 3 a snapshot at t = 52.9 s of the scalar concentration (Sc = 500), obtained in SR1 and SR3 at x/L = 2.5, is presented. The snapshot shows the ability of the employed code to resolve regions where steep gradients occur without introducing spurious oscillations, including in the deeper bulk region where the mesh size becomes relatively coarse (see figure 1). In figure 4 the cross-sectional profiles at x/L = 2.5, z/L = 4.38 obtained in SR1, SR2 and SR3 are compared. The resolution is significantly improved when increasing the scalar mesh refinement factor from 1 to 2, while an increment from 2 to 3 yields only a slight further improvement. In the near-surface region, which is our actual area of interest, the results using f RS = 2 and f RS = 3 are nearly identical. Thus, a refinement factor f RS 2 is deemed to be sufficient to accurately resolve the scalar transport.

Restriction on simulation time due to domain size
The simulations aim to elucidate the physics of gas transfer across the air-water interface driven by buoyant convection in deep calm bodies of water. To fully describe this phenomenon, the size of the computational domain would have to be large enough to allow all cold water plumes to move downwards until they naturally lose their negative buoyancy. Unfortunately, the current computational resources are insufficient to fully simulate this problem. Hence, it was decided to concentrate on the initial development of the Rayleigh instability and the early stages of transition, including the formation of convection cells and their influence on interfacial gas transfer. The maximum simulation time, during which the physics does not depend on the size of the computational domain, was determined by performing simulations using computational domains with a variety of sizes.
To determine the maximum simulation time during which bottom effects are negligible, two simulations were performed with different depth sizes (BC1 and BC2) using exactly the same disturbance to the initial temperature field at t = 9.6 s. As illustrated by the snapshots in figure 5(a,b), the surface temperature contours of simulations BC1 and BC2 up to t = 57.7 s were found to be identical even after the first plumes reached the bottom of the domain of simulation BC1 at t = 48.1 s. When comparing figure 5(c,d), a very slight difference can be seen close to the bottom of BC1 in the snapshots at t = 57.7 s showing temperature contours in the plane x = 2.5L from simulations BC1 and BC2. Higher up in the domain, however, the contours were found to be in very good agreement. Figure 6(a) shows that also the w peak versus time profiles are identical up to t = 57.7 s, where w peak is the maximum of the vertical w 2 profile, with · denoting horizontal averaging. Hence, in the free-fall regime the vertical velocity in BC1 and BC2 does not depend on depth. This regime starts at t ≈ 29.3 s and ends when the bottom effects begin to influence the flow dynamics. As illustrated in figure 6(b), the effect of the bottom on the total heat flux is felt first at t = 63.5 s, which is approximately 15.4 s after the plumes reached the bottom of the shallower domain (BC1). The results also suggest that for The results from the simulations with a larger horizontal domain (BC3 and BC4) are also included in figure 6. Because the initial disturbances added to the temperature field were not identical (though of the same level), these results do not show the same degree of agreement as in BC1 and BC2. By extrapolating these results, it is to be expected that, after the plumes reached the bottom and the fall velocity becomes constant, the convection cell size would continue to increase in size provided there are no horizontal restrictions. As the focus in this paper is on gas transfer induced by a buoyant-convective instability in deep waters, the remainder of the paper will mainly focus on the near-surface fluid flow and the scalar transport process before bottom effects become noticeable. Unless mentioned differently, the results presented below are all taken from the largest simulation BC3.

Visualization of thermal structures
The sequence of snapshots in figure 7 shows isosurfaces of the temperature field from BC3. It illustrates the 3D transport process driven by a buoyant-convective instability that occurs due to surface cooling. At first, heat transfer is dominated entirely by molecular diffusion, leading to a thickening of the thermal (cooler) boundary layer (not shown here). In this unstably stratified configuration, small fluctuations added to the temperature field at t = 9.6 s trigger an instability that is characterized by the formation of thin falling 'sheets' of saturated cold water. To replenish the falling cold water, warmer water from the bulk starts to flow upwards, resulting in the formation of convection cells. In the interior of the cells, one or more sources of upwelling warm water exist. Once the warm water reaches the surface, it moves radially outwards approximately parallel to the surface. While flowing along the surface, the warm water is cooled and slowly becomes increasingly heavy until it starts falling down, usually in the form of vertical sheets (see also Spangenberg & Rowland 1961) along the edges of the cells at a relatively high velocity, as can be seen by the colouring of the isosurface. At intersections of three or more convection cells, the added quantity of cold water increases the local sink speed, which results in the formation of mushroom-like plumes that are clearly visible in the snapshots at t = 38.9 s and t = 42.8 s. As the plumes sink, they tend to rotate about their axis, thereby distorting the mushroom-like shape and changing it into something resembling a blade. The rotational movement appears to enhance the stretching and downward movement of the blade-shaped plumes. The latter is reflected in figure 6(a), where w peak (after a transient period of t ≈ 43.3 s) does appear to increase somewhat in time as long as the bottom is not reached. As the convection cells grow, separate plumes merge into larger plumes and separate sheets into larger sheets. Note that a sheet typically spans the joint edge of two convection cells, while a blade originates from a plume that forms at the intersection of three or more convection cells and typically has a much smaller width and larger penetration depth than a sheet.

Top structures: net-like pattern, size of convection cells
Snapshots of the non-dimensionalized temperature distribution just beneath the top boundary, at a distance of z s = 0.0029L from the surface, are shown in figure 8. The convection cells discussed above are responsible for the formation of this typical netlike pattern at the surface. The snapshots show that the convection cells tend to grow in time, as illustrated by the fact that the average cell size at t = 46.6 s is significantly smaller than at t = 87 s. The increased area of high surface temperature that can be seen inside many individual cells at t = 87 s indicates that the upwelling is stronger than at t = 46.6 s. This is indirectly confirmed by the isosurfaces shown in figure 7, where the downward-falling plumes are much more pronounced at t = 87 s than at t = 46.6 s. As a result, a stronger upwelling motion would be needed to replenish the water at the surface.
Snapshots at t = 46.6 s and 67.8 s of the fluctuating temperature T − T in various z planes are shown in figure 9. It can be seen that the same net-like pattern that is observed at the surface is also identifiable in the upper bulk to depths of z s ≈ 0.5L and z s ≈ 0.8L at t = 46.6 s and 67.8 s, respectively. For larger depths, the convection cells lose their coherence because of the limited penetration depth of the sheets of falling cold water (for the smaller structures) and the interaction between the larger structures as they fall deeper into the bulk (cf. figure 8). The figure also shows the significant change in the temperature fluctuation intensity with depth. The intensity changes dramatically within a distance of 0.1L from the surface. Further down, the intensity decreases again as the bulk becomes more mixed by the interacting falling plumes.
The existence of net-like patterns at the surface has been known for a long time -see e.g. the qualitative description in Spangenberg & Rowland (1961). The well-resolved time-accurate data produced in the present DNS allows a quantitative estimation of the average convection cell size at the surface (L C ). Here L C is estimated by twice the average integral length scale at the surface, where L x L y is the geometric mean of the integral length scales in the x and y directions (L x and L y ) based on the u and v velocities, respectively. Figure 10(a) depicts the average convection cell size L C from BC3 estimated using (4.1). Also shown in the plot is the alternative estimation of the cell size L CT , which is calculated similarly to L C , after replacing the interfacial u and v velocities by the temperature in the grid plane immediately below the surface. It should be noted that the structures found in the interfacial ∂w/∂z contours and in the temperature contours seen in  Based on this, it is expected that L C and L CT are of comparable size. This is confirmed by the fact that the estimation of both L C and L CT using (4.1) is found to correspond well with the mean size of the convection cells displayed in figure 8. Figure 10(a) shows that the horizontal integral length scales L C and L CT tend to grow in time, which is in agreement with observations made in figures 8 and 10, where sequences of clearly identifiable structures show a significant growth of the convection cell size in time. A comparison of the growth of the cell size obtained in BC3 and BC5 with Rayleigh numbers Ra L | t=0 = 34 000 and Ra L | t=0 = 23 000, respectively, is given in figure 10(b). In both cases, the cell size begins to increase rapidly approximately 7 s after the first plumes start to break away from the thermal boundary layer, which happens at t ≈ 29.3 s for BC3 and at t ≈ 36.1 s for BC5. Furthermore, figure 10(a) shows that the integral length scale in BC3 keeps increasing, indicating a continued growth of the convection cells until the end of the simulation. It is expected that this growth will end only when one convection cell spans the entire computational domain. Figure 11 contrasts the diffusivities of the temperature and the scalar C Sc=500 by combining contours of T with isolines of C Sc=500 . A good correlation in the middle of the downward-falling sheets is obtained between the cold water from the surface and high scalar concentrations, which is maintained even in the head of the mushroom. As the difference in diffusivities of T and C becomes less, areas of low T and high C will tend to coincide more and more. This is illustrated in figure 12 by the almost perfect negative correlation between T and C Sc=20 of ρ(T, C Sc=20 ) ≈ −0.92 compared to values of ρ(T, C Sc=500 ) between −0.455 and −0.495.

Correlation of T and C
It is interesting to study the evolution of the scalar distribution of C Sc=500 in more detail. In an idealized situation, at the location where two equally strong convection cells meet, we would expect the formation of a symmetric 2D sheet of falling saturated water to form in the middle of the planar jet (corresponding to the coldest water). Because of the free-slip boundary conditions, the near-surface flow will be almost uniform. Hence, it is likely that, at the locations where convection cells meet, the relative scaling of the thermal and concentration boundary layer thicknesses (see figure 11) would be at least partially conserved in the falling sheets of cold water. Indeed, a relative scaling of δ C /δ T ≈ (Sc/Pr) −0.5 was retrieved in the centre of the falling sheets, where δ T and δ C are the thicknesses of the thermal and concentration layers inside the sheets. This almost perfect conservation of the relative boundary layer thicknesses in the falling sheets implies that the almost uniform flow extends further down than just the near-surface region.

Thermal boundary layer, Nusselt and Rayleigh numbers
Before the instability sets in, the colder temperature at the surface drives molecular diffusion and the thermal boundary layer thickness increases with time according to δ = (πκt) 1/2 , where κ is the thermal diffusivity (see e.g. Turner 1979). Figure 13(a) shows the temporal variation of the thermal boundary layer thicknesses δ cf and δ q . The former is identified by the distance from the surface where the temperature fluctuation is maximum, while the latter is defined by where subscript i denotes the interface. The small disturbances added to the temperature field at t = 9.6 s were found to result in an almost instantaneous growth of the Rayleigh instability. This is proven by the fact that, immediately after the addition of the disturbances, organized fluid movements (although very small) at the top surface started to appear. Figure 13(a,b) illustrates that, at least up to t = 27 s, diffusion dominates as the growth of the boundary layer thickness adheres to δ = (πκt) 1/2 . At approximately t = 29.3 s the horizontally averaged boundary layer thickness reaches its maximum before starting to decrease abruptly. This instance matches the time when the cold fluid collected along the border of small cells first starts to plunge down in the form of thin vertical sheets and plumes. As in Howard (1964), the thermal boundary layer Rayleigh number Ra δ is defined by (4.3) Figure 13(b) shows a detailed plot of the thermal boundary layer evolution without and with disturbances for two different expansion factors, α T = 0.000572 and 0.000381, as used in BC3 and BC5, respectively. Based on the boundary layer thickness reached at the time when fluid elements first break away from the boundary layer in BC3 and BC5, critical Ra δ values of approximately 2000 and 1800 were found, respectively. The initial evolution of the thermal boundary layer thickness in all simulations with identical α T and initial disturbance level was found to be the same. It should be noted that in our case T is the temperature difference between the water surface and the bulk, while in Howard (1964) T is the temperature difference between the lower hot and upper cold surfaces. For the case of two rigid boundaries, Howard (1964) estimated a critical Ra δ of the order of 1000.
In figure 13(a,b) it can also be seen that for BC3 immediately after t = 29.3 s both the mean boundary layer thickness and Ra δ rapidly decrease and finally reach plateaus at approximately 0.083L and 18, respectively. Figure 13(c) shows the temporal variation of T = T b − T s , where · denotes horizontal averaging, and the macroscale Rayleigh number Ra L . The initial development of the instability and its affect on gas transfer across the air-water interface gives rise to time-dependent Rayleigh numbers and a time-dependent Nusselt number given by where H L = q/ T is the heat transfer coefficient and the dimensional heat flux q is defined by where κ c is the thermal conductivity. From the definitions above it follows that Nu = (Ra L /Ra δ ) 1/3 (4.6) (e.g. Howard 1964). Figure 14 shows the relation between Nu and Ra δ . The data points are calculated using instantaneous horizontally averaged data. As the boundary layer thickness is correlated directly to the interfacial heat flux, the relation Nu ∝ Ra −1/3 δ holds for all times, including the diffusion-dominated regime at the beginning of the simulation and the regime where the effects of the limited depth become important at the end of the simulation. In the free-fall regime, after t ≈ 41.3 s, Ra δ becomes virtually constant so that with a value of c Ra in the region between 0.35 and 0.4 (see figure 15a). Some time after the free-fall regime, additional advective motions begin to affect the surface heat flux, resulting in a reduction in c Ra . Because the process is unsteady, the constant of proportionality c Ra in (4.7) was found to be sensitive to the distance from the surface, z b , where T b was evaluated (figure 15a,b for the case BC3). In the simulations with a domain depth d = 10L (BC2 and BC3), for instance, c Ra slightly increased by ≈0.02 when z b was decreased from 5L to 2.5L. To allow a consistent comparison between simulations with different depths (figure 15c), for this case it was decided to fix z b to a value of 2.5L, resulting in a constant of proportionality of c Ra = 0.39. In the present simulations, temperature was only fixed at the top. To compare the present c Ra to the free-slip results of Herring (1963), who used fixed temperatures at top and bottom, the different definitions of T need to be taken into account. As Herring's T is twice as large, his c Ra needs to be pre-multiplied by 2 1/3 , giving a value of ≈0.39, which is identical to the value obtained here. Figure 16 shows horizontally and time-averaged near-surface profiles of (a) the temperature T and (b) the temperature fluctuation T for various 9.6 s time intervals between t = 41.3 and 87 s. Note that the purpose of the time averaging over these short intervals is to slightly improve the statistics. When plotted against z s /δ cf , both the normalized mean temperature T b − T / T b − T s and the normalized r.m.s. temperature T / T b − T s profiles can be seen to collapse reasonably well. This illustrates that, in the quasi-steady state part of the free-fall regime, which is reached after a transient period (further detailed in § 4.5), the instantaneous horizontally averaged scalar profiles near the surface are approximately time-invariant (quasi-steady state). Furthermore, figure 16(a) shows that the horizontally averaged temperature eventually becomes quasi-depth-independent in the upper part of the bulk region. Figure 17 shows the mean profiles of (a) the velocity fluctuations and (b) the scalar quantities T and C averaged over the entire free-fall period. The horizontal velocity fluctuations u and v reach their maximum of 82 % of w peak at the surface where they affect the gas transfer to the unsaturated water flowing up from the bulk. For z s > 0.7L the horizontal fluctuations become almost constant. The vertical fluctuations reach their peak value at z s ≈ 1.9L.

Vertical profiles
In figure 18(a) the normalized r.m.s. profiles of the temperature (Pr = 6) and concentrations (Sc = 20 and 500) can be seen to agree reasonably well with a peak value of approximately 0.25. The normalized diffusive and turbulent vertical mass fluxes shown in figure 18(b) were also found to collapse nicely. As expected, the turbulent vertical mass fluxes were observed to quickly increase with depth. Simultaneously, the diffusive mass flux reduced from its maximum value at the surface so that the sum of the two mass fluxes (in the region with quasi-steady horizontally averaged statistics) remained constant for all z s in the upper bulk. At z s ≈ 0.7δ cf the turbulent and diffusive mass fluxes were found to be equal, while for larger z s the turbulent vertical mass flux dominated the total mass flux. Finally, for z s 2δ cf the diffusive mass flux was found to be zero while the turbulent mass flux reached its maximum value of (D∂ C /∂z| i ) and fully dominated the total vertical mass flux. This result supports the usage of the eddy-correlation measurement technique (at the liquid side) in the bulk to determine the gas flux at the interface. 4.5. Estimation of H L and K L In figure 6(b) the evolution of the temporal interfacial (dimensional) heat flux from BC3 was shown. With an initial temperature difference between the water surface (a) Velocity fluctuations: -· -, u /w peak ; ---, v /w peak ; --, w /w peak . (b) Temperature and scalars: --, temperature; ---, C Sc=20 ; -· -, C Sc=500 . and the bulk (T b,0 − T s ) of 2.2 K, instantaneous peaks in the heat flux in excess of 1600 W m −2 were observed. However, in time, the heat flux was found to slowly decrease (while the convection cell size increased). In BC5, where T b,0 − T s was reduced to 1.48 K, the maximum instantaneous heat flux decreased to approximately 1000 W m −2 . The reduction in the initial temperature difference from BC5 resulted in a thicker thermal boundary layer and hence a lower heat flux. The basic physical mechanism, however, is the same. Mainly because of the large initial temperature difference between the bulk and the water surface, in this study the heat fluxes were larger (though of similar order of magnitude) than most reported values observed in the oceans or lakes (e.g. Shay & Gregg 1984MacIntyre, Romero & Kling 2002). It is clear that a one-to-one comparison between the heat fluxes obtained in the present simulations with the ones observed in nature is difficult. Apart from the additional uncertainties in nature, such as the presence of surface contamination that is known to reduce the interfacial mass transfer significantly, an accurate determination of the surface temperature in field and laboratory measurements is extremely difficult to achieve. In most experiments T s is assumed to be identical to the air temperature some distance above the water surface. In still air, however, the formation of a very thick thermal boundary layer on the air side will result in a significant reduction of the temperature difference between the water surface and the bulk in a matter of seconds. On the other hand, evaporation effects (not considered here, as the effects are confined to the surface, where we assume a constant temperature) may result in a surface temperature that is lower than the air temperature some distance above the interface. In the present simulations we were able to fully resolve the initial (transient) development of the Rayleigh instability with relatively large heat fluxes and steep gradients that occur in nature when the water surface is suddenly cooled. Figure 19 shows the horizontally and time-averaged gas transfer velocity and boundary layer thickness as a function of the Schmidt number. It can be seen that both K L and δ cf perfectly scale with Sc −0.5 . Time averaging was performed from t = 41.3 s to t = 87 s. For t 87 s bottom effects became non-negligible, though the scaling of both K L and δ cf with Sc −0.5 was found to persist. This is in agreement with the K L and δ scaling found in our previous DNS of grid-stirred turbulence-driven gas transfer (Herlina & Wissink 2014) and is typical for cases with free-slip boundary conditions at the interface (perfectly clean water surface). Note that for contaminated surfaces the free-slip condition does not hold at the surface (see e.g. Khakpour et al. 2011).
The Sherwood number, defined as An important aim of this study is to correlate surface information to the gas transfer velocity K L . As can be seen in figures 7 and 8, the near-surface flow in the present buoyancy-driven gas transfer simulations is fully dominated by large structures. Based on this observation, it can be concluded that the large-eddy model of Fortescue & Pearson (1967) should be able to provide an accurate estimation of K L . According to this model, the surface renewal rate r can be estimated using velocity (U r ) and length (L r ) scales that are typical for the larger structures in the flow, so that r = U r /L r . Buoyancy-driven flows exhibit a net-like structure of convection cells at the surface (see figure 8). The average diameter of these convection cells, L C (see figure 10), estimated by (4.1), appears to be the appropriate choice for L r . Near the interface the flow in the convection cells is almost horizontal so that U r can be based on the velocity fluctuations u surf at the surface, defined by u surf = √ u v (see figure 20a). Hence, a renewal rate of r = u surf /L C is obtained, and K L and H L can be predicted using where D and κ are the scalar and thermal diffusivities, respectively. In figure 20(b) the predictions given in (4.10) are plotted together with K L (for Sc = 20 and 500) and H L . For t 38.5 s a very good agreement between the aforementioned equations and both K L and H L is obtained using a value of 1.1 for both c K and c H . This good agreement is not only observed in the quasi-steady state part of the free-fall regime but also persists well after the plumes have reached the bottom of the domain and the unsteadiness introduced by the falling plumes starts to accumulate inside the computational domain. Only during the initial phase of the simulation, where diffusive mass transfer dominates the scaling of K L and H L , is (4.10) not applicable. Note that when replacing L C by L CT in (4.10) an approximation of K L is obtained that is of equal quality for t 50 s and slightly worse for t < 50 s.
As proposed by McCready, Vassiliadou & Hanratty (1986), when very detailed surface information is available so that the surface divergence can be calculated, K L and H L can be estimated using where β is a measure of the surface divergence strength (see Turney & Banerjee 2013) given by (4.12) In figure 20(c) a constant of proportionality of c β = 0.59 in (4.11) is shown to give very good approximations. Note that Turney & Banerjee (2013) recently found that surface divergence theory is not applicable when a portion of the near-surface motions holds small time scales. The time scales in the buoyancy-driven flow presented here, however, are sufficiently long for surface divergence theory to be applicable. Figure 20(a) shows a comparison of w peak and u surf against time. It can be seen that initially the qualitative behaviour of w peak is very similar to that of K L shown in figure 20(b,c). In the free-fall regime for t 48.1 s, w peak begins to grow while K L reached a plateau. The reason for the very good initial agreement is that w peak is located near the surface so that w peak ≈ w , which is strongly related to the r.m.s. of the surface divergence.
By comparing figure 20(b) and (c) it can be seen that the quality of the estimation of H L and K L using (4.10) or (4.11) is similar for t 87 s. For t > 87 s, when bottom effects start to become important, a slightly better estimation is obtained when using the surface divergence model (4.11). Figure 21(a,b) shows estimations of the transfer velocities in various simulations using (4.10) and (4.11), respectively. The plotted data are time-averaged over the quasisteady part of the free-fall regime. Also displayed in figure 21 Above, we have shown that a very good time-accurate approximation of K L can be obtained using surface information. When no such information is available, it might be possible to use subsurface information instead, such as the vertical turbulent mass flux (see figure 18) or the velocity fluctuations combined with the integral length scale as detailed below. Analogously to the grid-stirred case of Herlina & Wissink (2014), can be employed to estimate the transfer velocity, where the bulk turbulent Reynolds number is given by R T = 2L ∞ u /ν, (4.14) L ∞ is the integral length scale and u is the r.m.s. of the horizontal velocity fluctuations. Note that (4.13) is a variant of the large-eddy model of Fortescue & Pearson (1967). To carry out a proper calculation of K L using (4.13), bulk measurements need to be performed within the convective mixing layer. In the present simulations, a very low R T 50 was obtained. As illustrated in figure 22, a reasonable time accuracy was achieved using a coefficient of proportionality of c RT = 1.8, which is comparable to the value of 1.6 found in the grid-stirred case.  . Estimations of H L and K L using data that were time-averaged over the quasisteady part of the free-fall regime: (a) using (4.10), (b) using (4.11). The GS data points are from the numerical simulations by Herlina & Wissink (2014); the HJ and KG results are from experiments performed by Herlina & Jirka (2008) and McKenna & McGillis (2004), respectively.

Conclusions
Large-scale DNS of mass transfer across a flat clean surface driven by buoyant convection have been carried out. The simulations are motivated by previous experiments performed at KIT (Jirka et al. 2010). The aim of the DNS is to produce and analyse highly resolved data leading to a more in-depth and detailed understanding of the physical mechanisms that play a role in buoyancy-driven gas transfer across the air-water interface. DNS calculations allow a full control of all boundary conditions, including the temperature and possible contaminations of the water surface. In the experiments, however, such a detailed control or knowledge of surface boundary conditions is never achieved. For instance, in the experiments by Jirka et al. (2010), the exact surface temperature could not be measured. Furthermore, in still air the large difference in thermal diffusivities of air and water will result in the formation of a very thick thermal boundary layer on the air side of the interface. As a result, the initial temperature difference between the water surface and the bulk will significantly reduce in a matter of seconds. Evaporation, on the contrary, could result in an interfacial temperature that is lower than the air temperature above. Another problem faced by experimentalists is unwanted surface contamination that tends to inhibit the free movement of water particles along the water surface. The latter also reduces the vertical fluxes that drive the gas transfer across the air-water interface.
The simulations were performed using a code especially developed to perform highly accurate simulations of mass transfer at realistic (high) Schmidt and Prandtl numbers on a computationally feasible mesh while avoiding spurious oscillations of the scalar quantities (Kubrak et al. 2013). The purpose of the simulations was to investigate buoyancy-driven gas transfer across the air-water interface. The presented results are all transient and range from the onset of the Rayleigh instability until the limited depth of the domain begins to affect the gas flux at the surface. The Rayleigh instability was triggered by adding random disturbances to the temperature field at t = 9.6 s. The overall effect of the developing instability was found to be largely the same. In those cases where exactly the same disturbance field was employed (BC1 and BC2), it was possible to perform a one-to-one comparison of instantaneous flow, temperature and concentration fields. These fields were found to be identical until the plumes in the shallower simulation, BC1, reached the bottom of the computational domain at t ≈ 48.1 s, while it took a further time lag of t = 15.4 s before the surface heat flux was affected. Hence, to describe this evolving Rayleigh instability, depth-independent Rayleigh numbers were used, based on either the thermal boundary layer thickness or a fixed length scale L. For larger temperature differences between the surface and the bulk (corresponding to larger Rayleigh numbers), a faster-growing instability was obtained. Initially the thermal boundary layer thickness was fully determined by thermal diffusion. Eventually, the growing Rayleigh instability resulted in the vertical convective flux becoming strong enough to stop a further increase of the thermal boundary layer thickness. The critical boundary layer Rayleigh number, Ra δ , at maximum boundary layer thickness was found to be Ra δ ≈ 2000 for the simulations with an expansion factor of α T = 0.572 × 10 −3 and Ra δ ≈ 1800 for α T = 0.381 × 10 −3 .
The Rayleigh instability is seen to give rise to the formation of relatively small convection cells that tend to grow in time. At the water surface these cells are visible as a net-like structure. In the interior of the cells one or more upwellings of low-saturated warm fluid from the bulk can be found. This fluid subsequently flows radially outwards while it is cooled by the cold air and gets increasingly saturated with atmospheric gases. At the boundaries of the cells, the cooled, now heavier, fluid subsequently falls down, forming thin vertical sheets of cold fluid. In the centre plane of these sheets, where the fluid is coldest, the highest concentration of atmospheric gases was found. The scaling of the thermal boundary layer thickness δ T = √ Sc/Pr δ C with the concentration boundary layer thickness, δ C , was found to be conserved in the sheets. The mushroom-like structures that are often reported in the literature were found to be formed at the intersection of three or more convection cells, leading to a local increase in the amount of cold fluid plunging down deep in the bulk at a relatively high velocity. Momentarily heat fluxes were obtained in excess of 1600 W m −2 . The above illustrates the effectiveness of buoyant-convective gas and heat transfer. The ratio of the gas transfer and the heat transfer velocities, K L /H L , was found to be equal to (Pr/Sc) 0.5 . Hence, even when only information on H L is available, estimates of K L for Schmidt numbers between Sc = 20 up to (at least) Sc = 500 can be obtained.
To model the renewal rate r (Danckwerts 1951), it is possible to use the r.m.s. of the surface divergence β as initially proposed by McCready et al. (1986). The best fit for K L ∝ √ Dβ was obtained using a constant of proportionality of 0.59. Another accurate estimation of r was found using the geometric average of the r.m.s. of the velocities at the surface, u surf = √ u v , and the average size of the convection cells L C , giving K L = 1.1 Du surf /L C as a best fit. The average cell size L C was approximated by 2 L x L y , where L x and L y are the integral length scales based on u in the x direction and v in the y direction, respectively.
Above, the transfer velocity was estimated using information that is available at the water surface. The two estimations were found to have a very good time accuracy in and beyond the quasi-steady state part of the free-fall regime. When only subsurface information is available, the results showed that it is still possible to estimate K L with reasonable time accuracy. Such an estimation can be based, for instance, on the turbulent mass flux w C or on u and the integral length scale L ∞ , where all quantities need to be evaluated somewhere inside the convective mixing region. The estimations will improve in time as the convective mixing layer becomes well developed.