Effects of surface tension reduction on wind-wave growth and air–water scalar transfer

Abstract Effects of surface tension reduction on wind-wave growth are investigated using direct numerical simulations of air–water two-phase turbulent flows. The incompressible Navier–Stokes equations for air and water sides are solved using an arbitrary Lagrangian–Eulerian method with boundary-fitted moving grids. The wave growth of finite-amplitude and non-breaking gravity–capillary waves, whose wavelength is less than 0.07 m, is simulated for two cases of different surface tensions under a low-wind-speed condition of several metres per second. The results show that significant wave height for the smaller surface tension case increases faster than that for the larger surface tension case. Energy fluxes for gravity and capillary wave scales reveal that when the surface tension is reduced, the energy transfer from the significant gravity waves to capillary waves decreases, and the significant waves accumulate more energy supplied by wind. This results in faster wave growth for the smaller surface tension case. The effect on the scalar transfer across the air–water interface is also investigated. The results show that the scalar transfer coefficient on the water side decreases due to the surface tension reduction. The decrease is caused by suppression of turbulence in the water side. In order to support the conjecture, the surface tension effect is compared with laboratory experiments in a small wind-wave tank.


Introduction
Wind waves are generated by winds blowing over the air-water interface. A great number of studies have been published investigating the effects of wind waves on momentum, heat and mass transfers between the ocean and atmosphere, which strongly affect the local and global weather and climate. The mass transfer across the air-water interface is in particular an important issue for climate change studies because it is related to the airsea CO 2 transfer (Wanninkhof et al. 2009). The study on the growth of wind waves also has a long history. The pioneering works on the theoretical aspect of the wind wave growth were done by Phillips (1957) and Miles (1957), and the understandings have been updated by a large number of analytical, experimental, and numerical studies, including recent studies, e.g., Grare et al. (2013); Paquier et al. (2015Paquier et al. ( , 2016; Zavadsky & Shemer (2017); Perrard et al. (2019); Buckley et al. (2020); Wu & Deike (2021). Sullivan & McWilliams (2010) and Wu & Deike (2021) well reviewed and introduced the literature. While most published studies have considered wind waves on uncontaminated sea water or fresh water with constant surface tension, some studies have discussed the effects of surfactants on the wind-wave structure and have shown that wind waves are suppressed by the surfactants in water (Scott 1972;Mitsuyasu & Honda 1986). However, the detailed mechanism of the wind-wave suppression has not been confirmed. The mechanism of wind-wave suppression due to surfactants has been hypothesised based on the Marangoni effect. More specifically, the interface deformation is suppressed by the surface tension gradient on the interface, which is referred to as the Marangoni stress (e.g., Scott 1972). Non-uniform surfactant concentration on the free surface causes non-uniform surface tension distribution because the surface tension decreases as the surfactant concentration increases. However, this hypothesis is not proven since the Marangoni stress cannot be directly measured.
Furthermore, even in the case where surfactants are uniformly distributed on the free surface, the effect of uniform reduction of surface tension on the growth of wind waves remains unclear. Kawai (1979) discussed the cause of increase of the critical wind speed for soap water, considering the effect of uniform surface tension reduction. Kawai analytically concluded that the critical wind speed does not change significantly when the surface tension is reduced to half the value of water. Tsai & Lin (2004) suggested that uniform surface tension reduction increases the wind-wave growth rate based on linear stability analysis. However, this result can be applied only to the initial wave growth at very early stages. Zonta et al. (2015) investigated the growth of gravitycapillary waves driven by the shear at the air-water interface using three-dimensional numerical simulation and showed that the initial wave growth is faster for a larger Weber number (smaller surface tension). However, they did not clarify the pure surface tension effect because they changed both Froude and Weber numbers. Recently, Wu & Deike (2021) numerically investigated the fundamental mechanism of the growth of gravitycapillary waves with a laminar wind profile, changing the Bond and Reynolds numbers. However, they considered wind waves whose amplitude was smaller than the viscous sublayer thickness. Therefore, their two-dimensional simulation results cannot directly be applied to three-dimensional air-water flows involving turbulent eddies.
In addition, the dependence of uniform surface tension reduction on the mass transfer across the air-water interface is also a subject of interest. Previous studies showed that the mass transfer is dominated by turbulence beneath the air-water interface (Komori et al. 1993a(Komori et al. , 2010. Veron & Melville (2001) showed that the effect of scalar transfer is affected by Langmuir circulations, which are streamwise vortices developing below a wind-sheared air-water wavy interface. In contrast, Takagaki et al. (2015) performed direct numerical simulations (DNSs) and concluded that the scalar transfer across a wind-driven airwater interface is mainly controlled by turbulent eddies, and the effect of the Langmuir circulations is relatively small. Recently, Tejada-Martínez et al. (2020) performed DNS and large-eddy simulation and reported that turbulence on the water side is triggered by the Langmuir forcing, whereas the intensity of the turbulence and the scalar transfer induced by the turbulence are driven by the wind shear. If the uniform surface tension changes the wave growth and the turbulence beneath the interface, the mass transfer across the interface could be significantly affected by the surface tension.
This study aims to clarify the effects of uniform surface tension reduction on wind-wave growth and scalar transfer across the interface under a comparatively low wind speed condition of several meters per second, where intensive wave breaking does not occur. DNS of air-water two-phase flow is performed to investigate the effects of the surface tension reduction on the mechanism of the wave growth and air-water scalar transfer. In addition, the surface tension effects obtained by the DNS are discussed in Appendix A by comparing with the wave-height measurements for aqueous solutions with different surface tensions in a small wind-wave tank.

Direct numerical simulation
2.1. Computational method The governing equations for air and water flows are the continuity and Navier-Stokes equations for incompressible flows: where U i is the velocity component in ith direction, p is the pressure, ρ is the density, ν is the kinematic viscosity, and g is the gravitational acceleration. Different values of ρ and ν were used for air and water flows: ρ = ρ a and ν = ν a for the air flow, while ρ = ρ w and ν = ν w for the water flow. The air-water interface was tracked by the height function η(x 1 , x 2 ), which represents the horizontal distribution of the vertical interface position. The transport equation of η is given by 3) The dynamical balances at the interface in the normal and tangential directions are given by where the subscripts a and w represent the quantities in the air and water, respectively. τ n and τ t are the normal and tangential components of the viscous stress tensor on the interface, respectively. p s is the pressure gap due to the surface tension σ and is given by p s = σκ, where κ is the interface curvature. The definition of the curvature κ is summarised in Appendix B. The transport equation of the non-dimensional passive scalar C is given by where D is the diffusion coefficient of the scalar. The governing equations were solved using an arbitrary Lagrangian-Eulerian (ALE) method with boundary-fitted moving grids (Komori et al. 1993b(Komori et al. , 2010. The grid points were allowed to move in the vertical direction following interface deformation, while they were fixed in horizontal direction. The grid points on the interface moved along with the vertical interface motion, and the vertical positions of the other grid points on the air and water sides were updated at each time step adapting to the interface height change. The governing equations described on the Cartesian coordinate were transformed to the equations on the generalized coordinate and solved in separated domains for air and water flows; i.e., the density and the viscosity for each flow were constant at each grid point. The detailed formulation of the ALE method used here is described in Komori et al. (1993b). The fifth-order upstream and fourth-order central finite difference schemes were adopted for the calculation of the advection and viscous terms, respectively. The velocity and pressure fields were coupled using the marker and cell (MAC) method. The time integration of each equation was calculated by the Euler implicit method. The reliability of the DNS code was confirmed with a comparison with laboratory experiment data (Komori et al. 1993b(Komori et al. , 2010Kurose et al. 2016).
As the interface shape is represented by the height function η(x 1 , x 2 ), the present DNS cannot consider the wave breaking with air entrainment or discontinuity of the surface slope. Occurrence of wave breaking can cause a numerical instability for (2.3). Deike and his colleagues (Deike et al. 2015(Deike et al. , 2016Mostert & Deike 2020;Wu & Deike 2021) performed DNS of two-phase flows with water waves using the volume of fluid (VOF) method with adaptive mesh refinement. Wave breaking can be captured by the VOF method explicitly. Here, we use the ALE method with boundary-fitted moving grids because breaking waves are not targeted. It should be also noted that the DNS of boundary-fitted approach in previous studies (Yang & Shen 2011;Zonta et al. 2015) used spectral methods for the spatial discretization in the horizontal directions. The differentiations by finite difference schemes are less accurate than those by spectral methods for large wavenumbers, and the difference between the spectral and the finite difference methods can be evaluated by the effective wavenumber k eff (e.g., Ferziger et al. 2020). For the case of the finite difference schemes used in the present DNS, we confirmed that the difference between the effective wavenumber k eff and the exact wavenumber k is less that 5% for k 0.47k max , where k max is the maximum wavenumber and k max ≈ 6.28 × 10 3 m −1 for the present DNS. Therefore, the error in the present finite difference schemes is negligibly small for k 3.0 × 10 3 m −1 . Figure 1 shows the schematic diagram of the computational domain. The domain size was set to L 1 = 16δ, L 2 = 3.84δ and L 3 = 3δ in the streamwise (x 1 ), spanwise (x 2 ) and vertical (x 3 ) directions, respectively, where δ was set to 1.25 × 10 −2 m. The initial interface was located 2δ above the bottom, and this height was set to x 3 = 0. The number of grid points was set to 400 in the streamwise direction, 96 in the spanwise direction, and 180 in the vertical direction (60 points on the air side and 120 points on the water side). The uniform grid spacing was used in streamwise and spanwise directions; i,e, the streamwise position of m 1 th grid point is x g (m 1 ) = m 1 ∆x 1 , and the spanwise position of m 2 th grid point is y g (m 2 ) = m 2 ∆x 2 . The horizontal grid spacing is ∆x 1 = ∆x 2 = 5.0 × 10 −4 m. In the vertical direction, the non-uniform grid spacing was adopted to set fine grids near the interface. As explained in the previous subsection, the vertical grid positions change depending on the interface height η; i.e., the vertical position of m 3 th grid point is given by z g (m 1 , m 2 , m 3 ) = −2δ + {η(x g (m 1 ), y g (m 2 )) + 2δ}ζ(M w − m 3 , M w ) for m 3 = 0, · · · , M w (the water side and the interface), and z g (m 1 , m 2 , m 3 ) = δ − {δ − η(x g (m 1 ), y g (m 2 ))}ζ(m 3 − M w , M a ) for m 3 = M w + 1, · · · , M w + M a (the air side), where ζ(m, M ) = tanh{α(1 − m/M )}/ tanh α, and α is the scaling constant (α = 2.8). M w and M a are the number of vertical grid points in the water and air sides, respectively. When the interface is flat, i.e., η(x 1 , x 2 ) = 0, the vertical position of m 3 th grid point is given by

Computational condition
The minimum vertical grid spacing was 9.0 × 10 −6 and 8.8 × 10 −6 m in the air and water sides, respectively. These were sufficiently small to resolve the viscous sublayer thickness δ ′ ≈ 5ν/u * , which was 3.2 × 10 −4 and 6.0 × 10 −4 m for the air and water sides respectively. The periodic boundary condition was applied to the velocity and pressure in the streamwise and spanwise directions, while the Neumann condition was applied in the vertical direction. A wall-bounded turbulent air flow had been developed for a flat and rigid surface until the velocity profile becomes a statistically steady state. The developed flow field was then used for the initial condition in the air side. The friction velocity on the air side was u * a ≈ 0.24 m/s. The corresponding free-stream wind speed was U ∞ ≈ 5.2 m/s based on the empirical relationship of Iwano et al. (2013). The water side was initially in static condition, and the initial displacement of the interface was set to zero. A constant pressure gradient was applied to the air side in the streamwise direction for the driving force. The time integration was calculated for 7 seconds. The air and water densities were set to ρ a = 1.2 and ρ w = 1.0 × 10 3 kg/m 3 , respectively and the kinematic viscosities in the air and water sides were set to ν a = 1.5×10 −5 and ν w = 1.0×10 −6 m 2 /s, respectively. The gravitational acceleration was set to g = 9.8 m/s 2 . This study considers three cases of surface tension σ: (1) the case where σ is equal to that of water (σ = 1.0σ w ), (2) the case where σ is half to that of water (σ = 0.5σ w ), and (3) the case where σ is initially σ = 1.0σ w and suddenly reduced to σ = 0.5σ w at t = 4.0 s. The surface tension of water was set to σ w = 7.2 × 10 −2 N/m. We have confirmed the grid resolution dependence of the DNS results (Appendix C).
The friction Reynolds number is defined for each side of air and water using the friction velocity and the boundary-layer thickness (δ a and δ w for the air and water sides, respectively). The friction Reynolds number on the air side, Re τ,a = u * a δ a /ν a , is approximately 200, and it is sufficiently large to observe the logarithmic region. The friction Reynolds number on the water side, Re τ,w = u * w δ w /ν w , is dependent on the time as the boundary layer develops along with the time. The DNS results show that, for t = 4.0-7.0 s, Re τ,w ranges from 165 to 193.
The scalar C represents the non-dimensional concentration of the imaginary absorbed gas. The scalar transport only in the water side was calculated because the concentration is uniform (C = 1) in the air side. The initial scalar concentration in the water side was set to 0, and the scalar concentration on the interface was fixed to C = 1 uniformly. The scalar transport calculation was started at t = 4.0 s, when wind waves were well developed. The Schmidt number, defined as Sc ≡ ν w /D, was set to Sc = 1 in the water side. Note that real gases usually have much larger Sc; e.g., Sc ≈ 600 for CO 2 . The scalar transfer coefficient k L (which is defined in section 3.2 and also referred to as the transfer velocity) for such large Sc can be estimated from that for Sc = 1 using the relationship k L ∝ Sc −1/2 (Jähne et al. 1979).

Results and discussion
3.1. Effect on wave growth Figure 2 shows the shape of the air-water interface at t = 6.0 s for the water case (i.e., σ = 1.0σ w ) and the reduced surface tension case (i.e., σ = 0.5σ w ). The wind and the propagating waves move in the x 1 direction. For the case of σ = 1.0σ w , the wind waves form ripple waves on the downwind side of the wave crest. For the case of σ = 0.5σ w , the ripple waves are less significant than those for σ = 1.0σ w , and the wave crests are slightly sharper. The above observation can be also confirmed by the streamwise distribution of instantaneous η at the middle position of the spanwise width (x 2 = 2.4 × 10 −2 m) in figure 2(c). For σ = 0.5σ w , the distribution of η in 0.1 x 1 0.2 m does not show clear ripple waves, while that for σ = 1.0σ w shows ripple-like interface fluctuations in the downwind side of the wave crest.
To quantify the effect of surface tension on the wave height, we have calculated the significant wave height, H s , and the significant wavelength, L s (e.g., Toba 1972;Holthuijsen 2007). The significant waves are defined as the waves whose height is the largest one-thirds among all individual waves, and the individual waves are identified using the zero-up crossing method (Pierson Jr. 1954; Takagaki   are the average wave height and wavelength of the significant waves, respectively. Figure 3 shows the temporal variations of H s and L s . Note that, as shown in figure 3(b), the significant waves do not capture wind waves at the initial period (t 1.0 s) because the interface initially oscillates due to the initial impact of imposing the wall-bounded turbulent air flow driven by the pressure gradient. However, in the period of t = 1.0-2.0 s, the wavelength is close to the minimum for both cases, and the wavelength for σ = 0.5σ w is slightly shorter than that for σ = 1.0σ w . During and after this period, the significant wind waves are well captured, and their wavelengths for both σ cases remains less than or comparable to 0.07 m until t = 7.0 s. This means that the wind waves simulated in this study are classified as gravity-capillary waves (0.004 L s 0.07 m) according to the classification in Lin et al. (2008). Figure 3(a) shows that, after the significant waves capture wind waves, waves for σ = 0.5σ w start growing earlier than σ = 1.0σ w , and the wave height for σ = 0.5σ w is remarkably higher than that for σ = 1.0σ w after t = 3.0 s. This result means that the waves are not suppressed by the uniform surface tension reduction in contrast to the effect of surfactants. We have also confirmed the effect of uniform surface tension reduction on the wave height by conducting laboratory experiments using a small wind-wave tank (see Appendix A).
One could argue that the difference in wave height between the cases of σ = 1.0σ w and σ = 0.5σ w in figure 3(a) is caused by hysteresis due to the difference in the initial wave development when the wind shear is loaded on the flat free surface. In order to clarify the dominance of hysteresis, we have performed the additional simulation with the surface tension suddenly changed to σ = 0.5σ w from the wind waves at t = 4.0 s developed using σ = 1.0σ w . The results are shown by blue lines in figure 3. When the surface tension was changed to σ = 0.5σ w , the significant wave height H s becomes higher than the case of σ = 1.0σ w and becomes close to the case with the surface tension σ = 0.5σ w from the initial time (t = 0). This indicates that the hysteresis is not critical, whereas the wave growth speed is increased due to the surface tension reduction.
The surface tension dependence of the significant wavelength in the period of t = 1.0-2.0 s is similar to the analytical result of the Kelvin-Helmholtz (K-H) instability, based on which the initial wavelength is given by L m ≡ 2π γ/g (where γ = σ/ρ w ), which is the wavelength for minimum phase velocity. However, based on the theoretical analyses (Jeffreys 1925;Miles 1957Miles , 1959, the classical K-H instability is not considered as the major mechanism that causes the initial wind-wave formation on an air-water interface. Jeffreys (1925) and Miles (1957Miles ( , 1993 proposed theoretical models for the formation of initial waves under turbulent winds. Their models estimate that the wave energy grows exponentially, and the initial wavelength formed at the minimum wind speed is about 0.06-0.08 m, where the effect of the surface tension is considered to be small. Miles' model is well accepted because the estimate agrees with observation data. However, this model cannot be applied to the present DNS results for the period of t = 1.0-2.0 s. Miles (1957Miles ( , 1993 considered the critical layer, which has the same mean wind speed as the wave speed, in the logarithmic layer of a wall-bounded turbulent flow, whereas, in the the present DNS results, the height for the wind speed equivalent to the wave speed is in the viscous sublayer. Miles (1959) generalized the K-H model to take account the logarithmic layer of a wall-bounded turbulent flow, and concluded that the K-H instability is unlikely for air-water interfaces at commonly observed wind speeds. Apart from the initial wave formation, Bontozoglou & Hanratty (1990) examined the weakly nonlinear K-H instability (Miles 1986) to characterize the instability of finite amplitude waves close to resonance, and showed that the K-H instability becomes strongly subcritical for the waves with wavelength shorter than the resonant wavelength given by L KH ≡ 2π 2γ/g. Their postulate is that, for air-water flows, gravity-capillary waves can be generated by winds and grow by the subcritical K-H instability so that those become sufficiently high to trigger a finite amplitude bifurcation of doubling the wavelength (Chen & Saffman 1979, 1980. The resonant wavelength L KH for σ = 1.0σ w and σ = 0.5σ w are 0.024 and 0.017 m, respectively. In the present DNS results for both surface tension cases, the significant wavelength L s in the period of t = 1.0-2.0 s almost corresponds to the wavelength L KH . Therefore, we can conjecture that the waves in this period are amplified due to the subcritical K-H instability for finite amplitude waves close to resonance. The one-dimensional wave height spectrum S ηη (k 1 ) in the streamwise direction has been calculated to understand the temporal change of H s and L s . The one-dimensional wave height spectrum is defined as is obtained by the Fourier transform in the streamwise direction, and k 1 is the streamwise wavenumber. The asterisk denotes the complex conjugate. Figure 4(a) shows the temporal change of the wave height spectrum for σ = 1.0σ w . Note that the spectrum is temporally averaged for the past 0.5 s with respect to the indicated time instant. In the early growth period (up to t = 2.5 s), as predicted by Bontozoglou & Hanratty (1990), waves with wavenumbers close to and larger than k KH ≡ g/2γ = 2π/L KH grow predominantly. This is because the spectrum for wavenumbers around k KH and 2k KH are amplified due to the subcritical K-H instability for the waves close to the second harmonic resonance at an early stage of wave growth (Bontozoglou & Hanratty 1990). After t = 2.5 s, the peak moves to a lower wavenumber as the wind waves grow, and the spectrum tail broadens on the high wavenumber side. The change of the peak location to lower wavenumbers corresponds to the increase in the significant wavelength in figure 3(b). The decrease of the peak wavenumber is commonly observed for wind waves (e.g., Imasato 1976), and these would be due to the nonlinear wave-wave interactions. The broadening of the spectrum tail on the high wavenumber side is due to a resonant nonlinear wavewave interaction between the fundamental gravity wave and higher harmonics of order N , which forms ripple-like capillary waves (Chen & Saffman 1979;Fedorov et al. 1998;Caulliez 2013). The resonant condition can be approximately satisfied when the phase velocity of the capillary waves matches that of the fundamental wave. The wavenumber k c of the capillary waves that resonate with the significant waves with wavenumber k s is given by k c = k 2 m /k s , where k m = g/γ = 2π/L m , when the non-linearity is negligible. The arrows in figure 4(a) indicate the wavenumber k s = 2π/L s for the time-averaged  significant wavelength L s for 4.0 s < t 7.0 s and the wavenumber k c of the resonant capillary waves. The spectrum at t = 5.0 and 7.0 s shows a gentle slope for k m < k < k c and decreases steeply for k > k c . This confirms that the harmonic resonance causes the broadened spectrum. Figure 4(b) shows similar temporal change of the spectrum for σ = 0.5σ w : Even though k KH for σ = 0.5σ w is larger than that for σ = 1.0σ w by a factor of √ 2, waves with wavenumbers around k KH grow predominantly in the early growth period (up to t = 2.0 s), and then the spectrum broadens on low and high wavenumber sides. Figure 4(c) compares the wave height spectra for the two cases of surface tension. The spectra are temporally averaged over the period of 4.0 s < t 7.0 s. On the low wavenumber side, the peak for σ = 0.5σ w is larger than that for σ = 1.0σ w , corresponding to the larger significant wave height. On the high wavenumber side, the wavenumber k c of the resonant capillary waves increases due to the surface tension reduction because k m is inversely proportional to the square root of the surface tension. The bump on the spectrum also moves to the higher wavenumber following the increase of k c . This indicats t>V@ E g E s >10 −2 -P 2 @ E g E s σ = 1.0σσ σ = 0.5σσ Figure 5. Temporal variations of potential energies due to gravity and surface tension, Eg and Es, for the cases of σ = 1.0σw and σ = 0.5σw.
that the wavelength of the resonant capillary waves becomes shorter. It is also observed that the amplitude at k c for σ = 0.5σ w is about two order of magnitude smaller than that for σ = 1.0σ w . This means that ripple-like capillary waves are suppressed due to the surface tension reduction.
To understand the effect of surface tension reduction on the wind-wave growth mechanism, the potential energies due to gravity E g and surface tension E s were calculated using the following equations: (3.2) Figure 5 shows the temporal variations of E g and E s . The potential energies E g and E s are comparable in the early growth period of 1.5 s < t < 2.5 s for σ = 1.0σ w and 1.0 s < t < 2.0 s for σ = 0.5σ w . According to the small amplitude linear wave theory, E g = E s implies the wavelength is equal to L m . Figures 4(a) and (b) also show that waves with wavenumbers around k KH ∼ O(k m ) grow predominantly in the early growth period. This is also consistent with the fact that the significant wavelength in the early growth period is close to L KH = √ 2L m as shown in figure 3(b). In the early growth period, the potential energy increases exponentially, and the energy growth rate for σ = 0.5σ w is larger than that for σ = 1.0σ w . These results are qualitatively consistent with the results from the linear stability analysis of Tsai & Lin (2004). Lin et al. (2008) reported that, based on their DNS for air-water two-phase flows, initial linear wave growth is followed by the exponential wave growth. In the present DNS results, only the later exponential growth is captured. After the early growth period, E g increases, whereas E s remains almost constant. The difference in the temporal variation between E g and E s corresponds to the broadening of the wave height spectra on both low and high wavenumber sides as shown in figures 4(a) and 4(b). The surface potential energy E s shows larger values for the larger surface tension as expected from the definition of E s . In contrast, the potential energy E g becomes larger for the smaller surface tension. The increase of E g due to surface tension reduction is more significant than the decrease of E s . This also means that E g increases faster for the smaller surface tension. It should be noted that Zavadsky & Shemer (2017) reported the temporal development of the wave height measured by applying nearly impulsive wind forcing. Their results show similar fast initial wave growth followed by the relatively slow wave growth. One would find that the growth of H s becomes even gentler for t > 5.0 s than t = 4.0-5.0 s, while E g for σ = 0.5σ w increases almost linearly along with the time for t > 4.0 s. Zavadsky & Shemer (2017) reported that the quasisteady equilibrium state appears when the waves are long enough to be purely gravity waves, and this stage is considered as 'the principal stage' of Phillips (1957). Therefore, the gentler growth for σ = 0.5σ w at t > 5.0 s could be relevant to this stage because the waves for σ = 0.5σ w are closer to purely gravity waves than the waves for σ = 1.0σ w . However, clear transition to quasi-steady equilibrium state is not observed in the present DNS results.
To clarify the mechanism of faster increase of the potential energy E g after the early growth period for the case of smaller surface tension in figure 5, we focus on the energy input from air to water and the energy dissipation in the water side. The energy fluxes due to the normal and tangential stress at the interface, Q na and Q ta , respectively, are given by where u n is the normal component of the interface velocity, and u t1 and u t2 are the tangential components of the interface velocity. τ t1 and τ t2 are the tangential components of viscous stress in the directions of u t1 and u t2 . Γ represents the interface area, and dS is the infinitesimal area on Γ . Q na and Q ta represent the energy fluxes per unit horizontal area. Q na and Q ta can be considered as the energy fluxes attributed to the form and friction drags, respectively. It should be noted that, in (3.3), we included the term of u n τ na but its contribution was negligibly small. We have defined the energy dissipation rate in the water side per unit horizontal area, ǫ w , as where the local energy dissipation rate ǫ is given by represents the volume in the water side. Figure 6 shows the temporal variations of the energy fluxes from the air side to the water side, Q na and Q ta , as well as the energy dissipation rate on the water side, ǫ w . The dissipation rate ǫ w is shown as negative. A simple centred moving average for the period of 0.2 s is applied. Q na for σ = 0.5σ w is larger than that for σ = 1.0σ w except the period of 2.4 s t 2.9 s, and Q ta for σ = 0.5σ w becomes smaller than that for σ = 1.0σ w after t ≈ 1.9 s. Both Q na and Q ta supply the kinetic energy to the water side, but Q ta mainly accelerates the mean surface velocity and forms a shear layer under the interface. Q na contributes to the wave growth rather than Q ta because the energy flux to gravity-capillary waves is attributed to the form drag (Melville & Fedorov 2015). Thus, the larger values of Q na for σ = 0.5σ w means more energy input to the waves than the case of σ = 1.0σ w . It should be also noted that the energy dissipation could also t>V@ í í Q na Q ta −ε ε >:P 2 @ Q na Q ta −ε ε σ = 1.0σε σ = 0.5σε Figure 6. Temporal variations of the energy fluxes due to the normal and tangential stresses, Qna and Qta, and the energy dissipation rate ǫw in the water side as negative values. A simple centred moving average for the period of 0.2 s is applied.
have a significant effect on the wave growth because steady gravity-capillary waves can be obtained when the energy flux due to the form drag is balanced by the viscous energy dissipation (Melville & Fedorov 2015). Figure 6 indicates that ǫ w is not negligibly small compared to Q na . Therefore, we have examined the relationship between the energy fluxes and the rapid growth of E g : The time derivative of E g (i.e., dEg dt ) has been compared with the energy input to waves and wave energy dissipation. Note that the energy dissipation rate ǫ w also contains the dissipation due to mean shear, which does not contribute to the wave energy dissipation. Hence, we consider the energy dissipation due to fluctuation of the strain rate, i.e., ǫ ′ w ≡ ǫ w − ǫ w,m , where ǫ w,m is the energy dissipation rate due to the mean shear. We evaluated ǫ w,m on the boundary-fitted coordinate as where s ij zg0 is the horizontally averaged strain rate tensor for the same z g0 , i.e., the average over the surface Γ zg0 , which initially locates at x 3 = z g0 and moves along with the boundary-fitted coordinate. Here, we used the approximation of ∂zg ∂zg0 ≈ 1 to calculate ǫ w,m : We assumed the volume change at each grid is negligibly small. We have confirmed that the error of computing ǫ w due to the approximation is negligibly small (less than 4 %). The time derivative dEg dt has been calculated by fitting linear function to the time series of E g for every 0.5 s. Since we focused on the rapid growth of E g , the change rate dEg dt was calculated for 4.0 t < 7.0 s and 3.0 t < 7.0 s for σ = 1.0σ w and σ = 0.5σ w , respectively. When the potential energy E g increases, the kinetic energy associated with the gravity waves would also increases, and that the balance between energy input and dissipation should contribute to the increase of the total wave energy. However, it is difficult to extract the kinetic energy due to the wave motion from the total kinetic energy including the turbulent kinetic energy. Here, we assume that the kinetic energy of the waves is equivalent to the potential energy as it is the case of monochromatic waves. The change rate of the total energy of gravity waves is then given by 2 dEg dt . In figure 7, the change rate 2 dEg dt for every 0.5 s is compared with the time-averaged Q na  and Q na − ǫ ′ w for every 0.5 s in the corresponding period. The dotted line indicates the 2 dEg dt values equivalent to Q na or Q na − ǫ ′ w . In figure 7(a), the markers are located in lower side of the dotted line for both σ = 1.0σ w and σ = 0.5σ w , meaning that the energy flux Q na is larger than 2 dEg dt . In figure 7(b), the markers are also located in lower side but closer to the dotted line. This means that the correlation of 2 dEg dt with Q na − ǫ ′ w is better than that with Q na . Thus, the rapid growth of E g is attributed to the energy flux due to normal stress Q na minus the energy dissipation fluctuation ǫ ′ w . The wave growth is often quantified by using the wave growth rate (Plant 1982;Donelan et al. 2006;Melville & Fedorov 2015), which is given by the change rate of the wave energy. The comparison of the wave growth rate obtained from the DNS results with the previous measurements is described in Appendix D.
The reason of the larger Q na for σ = 0.5σ w is indeed not surprising because Q na is the energy flux attributed to the form drag, and it is reasonable to assume that the form drag becomes larger along with wave growth. The form drag D p is given by where n 1 is the streamwise components of the normal vector to the interface. Melville & Fedorov (2015) evaluated the energy flux due to the normal stress by Q na ≈ CD p , where C is the wave velocity. We have confirmed that Q na and CD p are well correlated also in our results. Here, the obtained form drag D p is compared with the significant wave height H s and the wave slope H s /L s in figure 8. Similar to figure 7, the time-averaged values for every 0.5 s are plotted for 4.0 t < 7.0 s and 3.0 t < 7.0 s for σ = 1.0σ w and σ = 0.5σ w , respectively. In figure 8(a), the form drag D p for the same wave height H s is larger for the smaller surface tension case. The difference in D p between the cases of σ = 1.0σ w and σ = 0.5σ w is robust because D p calculated for the case of sudden change from σ = 1.0σ w to σ = 0.5σ w at t = 4.0 s also becomes larger than D p for the case where the surface tension remained σ = 1.0σ w . Figure 8(b) shows that the form drag D p is rather correlated with the wave slope H s /L s because D p for the three cases of surface tension well collapse. This means that the significant waves become steeper due to the surface tension reduction and that causes the increases of the form drag D p . Therefore, the increase of Q na due to the surface tension reduction can be explained by the increase of the wave slope H s /L s . To investigate the relationship between the energy fluxes and the surface tension, we have evaluated the energy fluxes for the gravity and capillary waves separately. Here, the energy fluxes due to the normal stress at gravity and capillary wave scales (Q na,g and Q na,c respectively) are defined as where k is the magnitude of the two-dimensional horizontal wavenumber vector k = (k 1 , k 2 ), k max is the maximum wavenumber, and S Qna (k) is the Fourier spectrum of the energy flux due to the normal stress at the interface. S Qna (k) can be given by where A(k 1 , k 2 ) indicates the Fourier coefficient of a given horizontal distribution A(x 1 , x 2 ), and R[·] indicates the real part of the complex number. Q na,g and Q na,c can be interpreted as sharp-cut low-and high-pass filtered values of local energy flux −(p a + ρ a gη) ∂η ∂t for the cut-off wavenumber k cut . We set k cut to k m for σ = 1.0σ w . We have confirmed that the above energy fluxes satisfy Q na = Q na,g + Q na,c . Similarly, the energy dissipation at gravity and capillary wave scales, ǫ w,g and ǫ w,c respectively, are defined as where S ǫw (k, z g0 ) is the Fourier spectrum of the energy dissipation rate as a function of initial vertical position z g0 . S ǫw (k, z g0 ) can be given by (3.14) where s ij (k, z g0 (m 3 )) is obtained by the Fourier transform in the x 1 and x 2 direction for s ij (x g (m 1 ), y g (m 2 ), z g (m 1 , m 2 , m 3 )) for each m 3 . We used the approximation of ∂zg ∂zg0 ≈ 1 again to define ǫ w,g and ǫ w,c as (3.12) and (3.13), respectively. We have confirmed that sum of these energy dissipation rates is approximately equal to ǫ ′ w ; i.e., ǫ w,g + ǫ w,c ≈ ǫ ′ w . Figure 9 shows the energy flux due to normal stress and energy dissipation rate at gravity and capillary wave scales. After t = 2.5 s, Q na,g is larger than Q na,c , indicating that the energy flux due to the normal stress at the gravity wave scales has dominant influence on the total energy flux Q na . This is consistent with our intuition that the contribution of the form drag is dominated by the significant gravity waves. Concerning the dissipation rate, ǫ w,c is far larger than ǫ w,g . Therefore, the energy dissipation at the capillary wave scales is dominant for ǫ ′ w . The energy dissipation for capillary wave scales ǫ w,c for σ = 1.0σ w is significantly larger than that for σ = 0.5σ w . It is also observed in figure 9(b) that ǫ w,c is larger than Q na,c . This means that the energy dissipation in the water side exceeds the energy input from the air side at the capillary wave scales, and there must be the energy transfer from the gravity wave scales to the capillary wave scales. The difference between ǫ w,c and Q na,c is smaller for the smaller surface tension case. This indicates that the energy transfer from the gravity wave scales decreases due to the surface tension reduction. There are several mechanisms that can result in the inter-scale energy transfer; the four-wave resonant interaction, harmonic resonance of nonlinear waves, and turbulent energy transfer. Unfortunately, it is difficult to specify the dominant mechanism because of the difficulty of evaluating the energy transfer. However, comparison of Q ta (in figure 6) and ǫ w,c suggests that change in the resonant condition due to the surface tension reduction results in the decrease of the energy transfer from the gravity wave scales because ǫ w,c (i.e., the energy transfer from the gravity wave scales) is larger than Q ta , which produces shear-induced turbulence in the water side, for σ = 1.0σ w .
To fully understand the wave growth mechanism, quantitative evaluation of wave-wave interaction is crucial. However, it is difficult to evaluate the energy flux between waves directly from our simulation data. The present results of the energy fluxes at gravity and capillary wave scales indirectly prove the energy transfer between these scales and the effect of the surface tension reduction on the energy transfer. These results indicate that the significant gravity waves receive kinetic energy from the air side with contributions of form drag, i.e., Q na , and the energy is partially transferred to capillary waves and dissipated by the viscosity of water. When the surface tension decreases, the energy transfer to the capillary waves decreases, and the significant waves retain more energy, resulting in the faster growth of the wave height. The faster wave growth can further cause the increase in the wave slope and, therefore, the increase of the energy received from the air side with contributions of form drag.

Effect on scalar transfer
The effect of the surface tension reduction on the scalar transfer across the air-water interface has been investigated by solving the scalar transport equation (2.6) for t 4.0 s. The scalar transfer was evaluated by the scalar transfer coefficient in the water side k L , defined as (3.15) Here, C i is the scalar concentration on the interface (C i = 1), C b is the bulk scalar concentration on the water side (C b = 0), and F is the horizontally averaged scalar flux   defined as where x n is the normal distance from the interface. Figure 10 shows the scalar transfer coefficient k L for σ = 1.0σ w and 0.5σ w . Immediately after the start of the scalar transport calculation at t = 4.0 s, k L shows large values because large concentration gradient near the air-water interface is formed in this transient period. However, for t > 4.5 s, the scalar concentration boundary layer is sufficiently developed to achieve the local equilibrium of scalar transfer near the interface, and k L for σ = 0.5σ w is smaller than that for σ = 1.0σ w . That is, the scalar transfer decreases as the surface tension decreases. The difference is about 6.5 % for 5.0 s t 7.0 s. One could postulate that the decrease of k L is due to the difference in the surface area of the air-water interface. Figure 11 shows the temporal variation of the surface area expansion ratio, defined as the ratio of the surface area of the interface to the horizontal area. The change in the expansion ratio is less than or approximately equal to 2 %, and the difference of the expansion ratio is even smaller; i.e., the difference of the expansion ratio is too small to explain the difference of k L . These results indicate that the difference in the surface area of the interface does not explain the 6.5 % decrease in k L due to the surface tension reduction. The transfer of gases such as CO 2 across a wind-driven wavy air-water interface is strongly affected by the turbulence under the interface (Komori et al. 1993a). Figure 12(a, b) shows the local scalar flux F local ≡ D ∂C ∂xn on the interface. The colour contours on the side walls shows the scalar concentration in the water side. For the case of σ = 1.0σ w in figure 12(a), there are streamwise streaky structures of the local scalar flux on the interface, which are formed due to the turbulent flow structure in the water side (Komori et al. 2010;Takagaki et al. 2015;Tejada-Martínez et al. 2020). The similar streaky structures of the local scalar flux are also observed for the case of σ = 0.5σ w in figure 12(b), suggesting that the turbulence significantly affects to the scalar transfer even when the surface tension is reduced.
To evaluate the scalar transfer below the moving interface, we calculate the downward scalar flux on the continuous surface Γ zg0 , introduced in Section 3.1. The surface Γ zg0 can be also considered as the surface for a fixed vertical index m 3 , and the surface Γ zg0 for z g0 = 0 corresponds to the air-water interface Γ . The scalar budget is evaluated for the volume below the surface Γ zg0 . The scalar budget equation can be described for each level of z g0 as where χ(z g0 ) is the total scalar per unit horizontal area integrated for the volume below the surface Γ zg0 , and F diff (z g0 ) and F turb (z g0 ) are the downward molecular diffusion and turbulent scalar fluxes across the surface Γ zg0 , respectively. χ(z g0 ), F diff (z g0 ), and F turb (z g0 ) are defined as (3.20) respectively, where z(x 1 , x 2 , z g0 ) is the vertical position of the surface Γ zg0 , n j is jth component of the unit normal vector to Γ zg0 , dS zg0 is infinitesimal area on Γ zg0 , c is the scalar concentration fluctuation (c ≡ C − C zg0 ), w 3 is the vertical fluid velocity relative to the velocity of the boundary-fitted coordinate ( ∂t is the vertical velocity of the boundary-fitted coordinate), and the brackets · zg0 denote the horizontal average over the surface Γ zg0 , i.e., A zg0 = Figure 13. Vertical profiles of downward molecular diffusion and turbulent scalar fluxes, F diff and F turb , respectively, and the time change of stored scalar χ below zg0. Temporal average was taken for the period of 5.0 s t 7.0 s. Horizontal average was taken on the the boundary-fitted coordinate.
1 L1L2 L1 0 L2 0 A(x 1 , x 2 , z(x 1 , x 2 , z g0 ))dx 1 dx 2 for a given function A. The derivation of the scalar budget equation (3.17) is described in Appendix E. Figure 13 shows the vertical profiles of the downward molecular diffusion and turbulent scalar flux, F diff and F turb , respectively, and the time change of χ. The temporal average is taken for the period of 5.0 s t 7.0 s, while the scalar budget equation (3.17) is described for instantaneous scalar distribution. For both cases of σ = 1.0σ w and σ = 0.5σ w , F diff is dominant for the scalar transfer only near the air-water interface.
The profile of F turb shows a peak at z g0 ≈ −2.2 × 10 −3 m, and F turb is dominant below the z g0 . The peak value of F turb for σ = 0.5σ w is approximately 4.7 % smaller than that for σ = 1.0σ w . The magnitude and the difference of the peak values are consistent with the difference in dχ dt near the surface and the difference in k L in figure 10. These results indicate that the decrease in k L due to surface tension reduction is caused by the decrease in the turbulent scalar transfer.
In figure 6, it is observed that Q ta for σ = 0.5σ w is smaller than that for σ = 1.0σ w . This means that the energy input for the shear-induced turbulence weakens due to the surface tension reduction. To confirm the turbulence structures in the water side, the turbulent vortices have been visualized by using the second invariant of the velocity gradient tensor, Q ≡ 1 2 (ω ij ω ij − s ij s ij ), where ω ij is the vorticity tensor defined as ω ij ≡ 1 2 ∂Uj ∂xi − ∂Ui ∂xj , and s ij is the strain-rate tensor defined above. Figures 12(c,  d ) show isosurfaces of Q = 0.002 s −2 , viewed from the water side (bottom-side up). Strong turbulent vortices are distributed near the interface for both surface tension cases: Spanwise vortices are observed in the vicinity of the interface, and streamwise vortices are observed below the wave troughs. The structure similar to the spanwise vorticity was reported by Hung & Tsai (2009). They computed the evolution of two-dimensional gravity-capillary waves and observed that vortices shed from wave troughs distribute in ω 1,rms ω 2,rms >V −1 @ í í í í í z g0 >10 −2 P@ ω1, rms σ =1.0σw σ =0.5σw ω2, rms σ =1.0σw σ =0.5σw Figure 14. Vertical profiles of rms values for streamwise and spanwise vorticity fluctuations, ω1,rms and ω2,rms, respectively. Temporal average was taken for the period of 5.0 s t 7.0 s. the vicinity of the interface. The streaky structure of the scalar flux on the interface is attributed to trains of the streamwise vortices. It is also observed that the streamwise vortices for σ = 0.5σ w are less than those for σ = 1.0σ w .
To compare the turbulent fluctuation, we have calculated the rms values of the streamwise and spanwise vorticities, because it is difficult to separate the effects of wave motion and turbulent flow on the velocity fluctuation. Figure 14 shows the rms values of streamwise and spanwise vorticity fluctuations, ω 1,rms and ω 2,rms , respectively. Large values of ω 1,rms and ω 2,rms are observed in the vicinity of the interface, and these large values correspond to the spanwise vortices in figures 12(c, d ). At the vertical position of z g0 ≈ −2.2 × 10 −3 m, where the peak of F turb is located, the streamwise vorticity fluctuation ω 1,rms exceeds the spanwise vorticity fluctuation ω 2,rms . Thus, the streamwise vortices play major roles for the turbulent scalar transfer quantified by F turb . We can also find that the streamwise vorticity fluctuation is weaker for the case of σ = 0.5σ w . This also indicates that the decrease in the turbulent scalar transfer is due to the decrease in the turbulent fluctuation.
The shear in the water side is produced by the energy flux Q ta as explained above. Thus, the smaller Q ta can result in the weaker shear-induced turbulence. The decrease of Q ta is attributed to the decrease of the friction drag. Note that the friction and form drags satisfy the relationship of D f + D p = ρ a u 2 * a ≈ const., where D f and D p are the friction and form drags, respectively. D p is given by (3.8), and D f is given by where t 1 is the streamwise components of the tangential vector to the interface. Figure  15 shows the friction and form drags, D f and D p . For t 4.0 s, the form drag D p increases with time, and D p for the case of smaller surface tension is larger than that for the larger surface tension. As discussed in the previous subsection, the increase in D p is attributed to the increase in the wave slope H s /L s associated with the wave growth.
In contrast, the friction drag D f decreases with time, and D f for the smaller surface tension case is smaller than that for the larger surface tension. This suggests that the shear-induced water turbulence, which promotes the scalar transfer across the interface, can be suppressed by the decrease in D f due to the surface tension reduction.
These results indicate that, when the surface tension is reduced, the wind wave becomes higher and the surface area becomes slightly larger. However, the scalar transfer across the air-water interface decreases because the turbulence under the interface is suppressed due to the decrease in the friction drag.

Conclusions
Effects of uniform surface tension reduction on the wind-wave growth and the scalar transfer across the air-water interface have been investigated for finite-amplitude and non-breaking gravity-capillary waves using direct numerical simulation (DNS) of a wind-driven air-water two-phase turbulent flow. The incompressible Navier-Stokes equations for both air and water sides have been solved simultaneously using an arbitrary Lagrangian-Eulerian (ALE) method with boundary-fitted moving grids. Finite difference schemes were adopted for the computation of spatial derivatives. The periodic boundary condition was applied in the lateral directions, and the slip condition is applied to the top and bottom boundaries. Initially, the lower water side was static and the free surface was flat, while a wall-bounded turbulent air flow developed for a flat and rigid surface was imposed in the upper air side. The wind speed was comparatively low (approximately 5 m/s) so that no intensive wave breaking occurs. The DNS was performed for the time integration period of 7 seconds to examine growth of gravity-capillary waves with the wavelength of less than about 0.07 m Growth of wind waves was simulated for two values of the surface tension: the value for water and half of that. The results show that the significant waves for the smaller surface tension grow faster, and their wave height becomes higher than that for the larger surface tension. Although wind waves were not well captured during the first 1 second in the present DNS results owing to the initial impact of imposing the wall-bounded turbulent air flow on the air-water interface, the initial wave development is not critical for faster wave growth for the case of smaller surface tension. This is supported by the result that the wave growth speed is increased when the surface tension is suddenly reduced from the wind waves at t = 4.0 s developed with the large surface tension. The effect of surface tension reduction on the wave height has also been confirmed by the wave-height measurements in a small wind-wave tank (Appendix A).
The growth of wind waves have been also analyzed by the wave height spectra. For both surface tension cases, the spectra grow around the wavenumber for the weakly nonlinear Kelvin-Helmholtz instability for finite amplitude waves in the early growth period. After that, the spectra broaden for both low wavenumbers (gravity waves) and high wavenumbers (capillary waves), having a peak on the low wavenumber side. The broadening of the wave height spectra is consistent with the harmonic resonance, which forms ripple-like capillary waves. The spectrum also shows that, when the surface tension is reduced, a bump at the resonant capillary wave scale moves to higher wavenumber and becomes smaller, indicating that the resonant capillary waves become shorter and weaker.
The numerical convergence of the DNS results have been discussed in Appendix C. At the present grid resolution, the temporal variations of the significant wave height and wavelength and the wave height spectrum are weakly dependent on the resolution. However, the resolution dependence is smaller than the difference due to the surface tension reduction. Therefore, the DNS results under the present computational condition are enough to discuss the effects of surface tension reduction.
The temporal variations of wave energy show that the potential energies due to gravity and surface tension are comparable during the early growth period. After this period, the potential energy due to gravity increases rapidly, whereas that due to surface tension remains almost constant. The surface tension reduction results in the increase in the potential energy due to gravity and the decrease in that due to surface tension. The results also show that the growth of the potential energy due to gravity is faster for the smaller surface tension.
The energy flux from air to water has been also investigated to understand the faster growth. The energy flux due to the normal stress, which contributes to wave growth, becomes larger due to the surface tension reduction. The surface tension reduction also causes the decrease in the energy dissipation in the water side. The relationship between the energy flux and the temporal wave-energy change shows that the rapid increase of the wave energy is attributed to the balance between energy flux from the air flow and energy dissipation in the water. The energy flux and energy dissipation filtered for gravity and capillary wave scales show that the energy flux from the air flow is received mostly at the gravity wave scales, and the energy dissipation in the water occurs mostly at the capillary wave scales. This indicates that there is an energy transfer mechanism from the gravity wave scales to the capillary wave scales. The energy dissipation for the capillary wave scales becomes smaller due to the surface tension reduction. This means the energy transfer from the gravity wave scales to the capillary wave scales decreases due to the surface tension reduction. The analyzed results also show that the energy flux due to the tangential stress, i.e., the energy source of shear-induced turbulence, is smaller than the energy dissipation. This suggests that the wave energy due to the normal stress at the gravity wave scales is partially transferred to smaller scales by resonant wave interactions and dissipated at the capillary wave scales. In conclusion, the significant gravity waves receive kinetic energy from the air side with contributions of form drag, and the energy is partially transferred to the capillary waves and dissipated by the viscosity of water. When the surface tension decreases, the energy transfer to the capillary wave scales decreases, and the significant waves retain more energy, resulting in the faster growth of the wave height.
The scalar transfer across the air-water interface has been further discussed based on the DNS results. The scalar transfer coefficient (scalar transfer velocity) for the gas absorption decreases due to the surface tension reduction. The decrease in scalar transfer coefficient cannot be explained by the difference of the interface surface area expansion because influence of the surface tension reduction on the surface area is negligibly small. Since the turbulent scalar flux and the vorticity fluctuation in the water side becomes smaller due to the surface tension reduction, the decrease in the transfer coefficient is caused by the suppression of turbulence under the interface, particularly turbulence associated with streamwise vortices. The suppression of turbulence can be attributed to the decrease of the friction drag, which drives the shear layer beneath the interface. This work was supported by JSPS KAKENHI Grant Number JP19K21937. SK would like to thank the partial support by JST COI Grant Number JPMJCE1316. The authors would like to thank A. Kimura, K. Takane, and R. Hayashi for their help in conducting the experiments. The numerical simulations presented in this paper were carried out on the Earth Simulator supercomputer system in the Japan Agency for Marine-Earth Science and Technology.
Declaration of interests. The authors report no conflict of interest.

Appendix A. Comparison with wind-wave tank experiments
Laboratory experiments have been conducted to examine the DNS results, in which the wave height becomes higher due to the uniform surface tension reduction. Figure  16 shows the schematic diagram of the experimental apparatus. The experiments were conducted using a small recirculation wind-wave tank (length 1.05 m, width 0.20 m and depth 0.05 m), located in a wind tunnel (Nagata et al. 2007) (length 5.0 m, width 0.30 m and height 0.30 m). A straightened air flow was introduced over the wind-wave tank. The water surface was open to the air flow only at the central test section (length 0.90 m and width 0.09 m) of the tank. The surface water flow induced by the wind shear recirculated through the covered side sections so that a counter flow was avoided in the lower part of the central test section. A wave absorber was fixed at the end of the test section to prevent the reflection of wind waves and counter flow. The water-level fluctuation was measured using a capacitance-type displacement meter (Iwatsu Electric Co. Ltd., ST-3572) capable of non-contact measurement. The probe (Iwatsu Electric Co. Ltd., ST-0717A) was set at a fetch length of x = 0.70 m and at an elevation of 0.0020 m or 0.0030 m above the static water level. The sampling time and frequency were 300 s and 1000 Hz respectively. The measurements were conducted for five cases of the free-stream wind speed U ∞ between 1 and 6 m/s. In order to quantify the wave height, significant waves were identified using the zero-up crossing method.
Liquids used in the experiments were filtered tap water and two aqueous solutions of ethanol and glycerine. Table 1 shows the surface tension σ and the viscosity µ of each solution. The surface tension of the 30 wt% ethanol solution was approximately half of the filtered tap water. However, the viscosity of the aqueous ethanol solution was about twice to the filtered tap water. Thus, the concentration of the aqueous glycerine solution was adjusted to approximately 30 wt% so that the viscosity was close to the aqueous ethanol solution while the surface tension remained close to the filtered tap water. Figure 17 shows the measurements of the significant wave height, H s , with the freestream wind speed, U ∞ , for the filtered tap water and two aqueous solutions of ethanol and glycerine. When the wind speed U ∞ is higher than 5 m/s, the significant wave  Table 1. Measured values of surface tension σ and dynamic viscosity µ. height of the aqueous ethanol solution is larger than that of aqueous glycerine solution.
The difference in the significant wave height is due to the difference of surface tension and not to the difference in viscosity because the significant wave height for the aqueous glycerine solution is almost the same as that of the filtered tap water with the half viscosity at U ∞ > 5 m/s. Paquier et al. (2015Paquier et al. ( , 2016 investigated the effect of viscosity on the wind waves experimentally and reported the presence of "the wrinkle regime" under large viscosity or weak wind conditions and the transitions to the regular wave regime under small viscosity or strong wind speed conditions. For ν ≈ 1.0-2.0 × 10 −6 m 2 /s, the winkle regime should be observed for approximately u * a 0.19-0.21 or U ∞ 3.6-4.0 m/s. Thus, the difference in the effect of surface tension reduction between low and high wind speed could be due to the difference in the wave regimes. Note that, apparently from the visualization in figure 2, the wind waves obtained by the DNS are in the regular wave regime. The experimental condition for the filtered tap water for U ∞ ≈ 5 m/s approximately corresponds to the DNS condition of σ = 1.0σ w , where the viscosity was set to ν = ν w and the free-stream wind speed was U ∞ ≈ 5.2 m/s. Here, we have performed the DNS for additional cases for the aqueous ethanol and glycerine solutions. For the case of the aqueous ethanol solution, the surface tension and the kinematic viscosity were set to σ = 0.5σ w and ν = 2.0ν w , respectively, and for the aqueous glycerine solution, σ = 1.0σ w and ν = 2.0ν w , respectively. Other computational conditions were the same as those in Subsection 2.2. Figure 18 shows the temporal variations of the significant wave height, H s , for the cases of the filtered tap water (σ = 1.0σ w , ν = 1.0ν w ), the aqueous ethanol solution (σ = 0.5σ w , ν = 2.0ν w ), and the aqueous glycerine solution (σ = 1.0σ w , ν = 2.0ν w ). To compare the DNS results with the laboratory measurements, we estimate the wave growth time t in the DNS that corresponds to the fetch of x = 0.7 m in the wind wave tank. According to the fetch law (e.g., Wilson 1965), the fetch is formulated as the function of a wave velocity C s , which is defined as C s ≡ gL s /2π + 2πγ/L s for deep-water gravitycapillary waves with a wavelength of L s . This suggests that the significant wave velocity C s or wavelength L s is a suitable parameter for comparison between the DNS results and laboratory measurements. At the fetch of x = 0.7 m, the measured values of the wave velocity and the wavelength are C s = 0.29 m/s and L s = 0.033 m, respectively. In figure 3(b), the significant wavelength obtained by the DNS is close to the measured value for 3 s t 4 s. We have also confirmed that the significant wave velocity averaged over the period of 3 s t 4 s is C s = 0.279 m/s and close to the laboratory measurement. C s for the half surface tension case can be smaller than the filtered tap water case but its effect on C s is smaller than 6 % for the measured wavelength. In the growth stage of 3 s t 4 s, the DNS in figure 18 well supports the experimental results that the homogeneous surface tension reduction enhances the wave height even when the viscosity is doubled, and the viscosity effect on the wave growth is smaller than the surface tension effect. It is also observed that the wave heights in the DNS results are higher than the measured values. This could be due to the difference in the friction velocity, which was not measured in the experiments. Thus, we do not compare the DNS results with the laboratory measurements quantitatively. It would be also interesting to note that the viscosity effects of surface tension reduction on wind-wave growth. However, we omit the detailed discussion about the sensitivity to the viscosity because it is beyond the scope of this paper. t>V@ H s >10 −2 P@ σ = 1.0σ σ ν = 1.0ν σ σ = 0.5σ σ ν = 2.0ν σ σ = 1.0σ σ ν = 2.0ν σ Figure 18. Temporal variations of significant wave height Hs obtained by DNS for the filtered tap water (σ = 1.0σw, ν = 1.0νw), the aqueous ethanol solution (σ = 0.5σw, ν = 2.0νw), and the aqueous glycerine solution (σ = 1.0σw, ν = 2.0νw).
a constant value ξ Γ , and the curvature κ is given by the derivative of n in ξ 1 and ξ 2 directions: where a i ≡ ∂xj ∂ξi e j , and e j is the basis vector of the Cartesian coordinate. The normal vector is given by n = a1×a2 |a1×a2| . The relationship of n · a 1 = n · a 2 = 0 yields κ = h 11 g 22 + h 22 g 11 − 2h 12 g 12 g 11 g 22 − g 12 g 12 (B 2) where h ij ≡ n · ∂aj ∂ξi and g ij ≡ a i · a j .

Appendix C. Numerical convergence
The numerical convergence of the DNS results have been examined by changing the streamwise resolution. The DNS results obtained for M 1 = 400, where M 1 is the number of grid points, were compared with the additional DNS results for M 1 = 200, 300, and 532. The grid spacing was scaled for the fixed domain size. Figure 19 shows the vertical profiles of the mean streamwise air velocity U 1 and the root-mean-square of the streamwise air velocity fluctuations u ′ of the initial wall-bounded turbulent flows obtained for different resolutions. The vertical axis is normalized by u * a , i.e., U + 1 = U 1 /u * a and u ′+ = u ′ /u * a . The horizontal axis is the vertical position in the viscous unit z + = zu * a /ν a . As shown in figure 19, for both U + 1 and u ′+ , the profiles for different resolutions are well converged. Figure 20 shows the temporal variations of significant wave height and wavelength, H s and L s , respectively, for different resolutions. The growth of H s for M 1 = 200 is slower than the growth for higher resolution particularly for t < 3 s, but the wave growth for the other resolutions are well converged. For M 1 300, the difference in H s between the different resolutions increases gradually along with time. The small difference between the different resolution cases for M 1 300 is smaller than the difference due to the surface tension reduction. The convergence is similarly confirmed for temporal variation of L s .  Figure 21 shows the time-averaged wave height spectra for the cases of σ = 1.0σ w and σ = 0.5σ w obtained for different resolutions. As expected from the temporal variations of the significant wave height and wavelength, the wave spectra for low wavenumbers (k 1 < k m ) are well converged. For high wavenumbers (k 1 > k m ), the wave spectra are almost converged. We can still see weak dependence on the resolution but the resolution dependence is sufficiently smaller than the difference due to the surface tension reduction. Thus, the present DNS is sufficiently converged for all over the scales to discuss the surface tension effect. k 1 >P −1 @ í í í í í í í í í S ηη (k 1 )>P 3 @ km =3.7×10 2 km =5.2×10 2 η =1.0ηw η =0.5ηw σ1 =200 σ1 =300 σ1 =400 σ1 =532 Figure 21. Wave height spectra averaged for the period of 4.0 < t 7.0 s, obtained for σ = 1.0σw and σ = 0.5σw with different resolutions. Vertical dotted lines in black and red are the wavenumbers km for σ = 1.0σw and σ = 0.5σw, respectively. d log(Eg +Es) dt . Figure 22 shows β calculated by fitting linear function to log(E g + E s ) for every 0.5 s. The gray symbols are the measurements summarized in figure 2 of Plant (1982). The solid and dashed lines are the relation obtained by Plant (1982), and given by: where u * a is the friction velocity on the air side, ω = 2πf is the radian wave frequency, θ is the angle between wind and wave directions, and C is the phase speed. As shown in figure 22(a), the relationship between β/f and u * a /C agrees with the measurements. Note that ρ a u 2 * a is the sum of form and friction drags, D p and D f , respectively, while u * a in (D 2) is relevant to the contribution of the form drag D p (Melville & Fedorov 2015). Since the energy flux due to the normal stress approximately satisfies Q na ≈ CD p , the wave growth rate β/f should be compared with Q na /ρ a C 3 . In figure 22(b), the wave growth rate obtained from our DNS results is plotted against Q na /ρ a C 3 , and the correlation with Q na /ρ a C 3 is better than that with (u * a /C) 2 .