Accounting for surface waves improves gas flux estimation at high wind speed in a large lake

The gas transfer velocity (k) is a major source of uncertainty when assessing the magnitude of lake gas exchange with the atmosphere. For the diversity of existing empirical and process-based k models, the transfer velocity increases with the level of turbulence near the air–water interface. However, predictions for k can vary by a factor of 2 among different models. Near-surface turbulence results from the action of wind shear, surface waves, and buoyancy-driven convection. Wind shear has long been identified as a key driver, but recent lake studies have shifted the focus towards the role of convection, particularly in small lakes. In large lakes, wind fetch can, however, be long enough to generate surface waves and contribute to enhance gas transfer, as widely recognised in oceanographic studies. Here, field values for gas transfer velocity were computed in a large hard-water lake, Lake Geneva, from CO2 fluxes measured with an automated (forced diffusion) flux chamber and CO2 partial pressure measured with high-frequency sensors. k estimates were compared to a set of reference limnological and oceanic k models. Our analysis reveals that accounting for surface waves generated during windy events significantly improves the accuracy of k estimates in this large lake. The improved k model is then used to compute k over a 1-year time period. Results show that episodic extreme events with surface waves (6 % occurrence, significant wave height > 0.4 m) can generate more than 20 % of annual cumulative k and more than 25 % of annual net CO2 fluxes in Lake Geneva. We conclude that for lakes whose fetch can exceed 15 km, k models need to integrate the effect of surface waves.


Introduction
Lakes are universally regarded as significant sources of CO 2 to the atmosphere; however, the accurate quantification of the magnitude of such emissions currently remains challenging (Cole et al., 2007;Tranvik et al., 2009;Raymond et al., 2013). CO 2 fluxes can be directly measured with floating chamber or eddy covariance systems (Vachon et al., 2010;Vesala et al., 2006). However, both approaches have their own constraints. The former suffers from limited time and space integration (from minutes to hours and from centimetres to metres respectively; Klaus and Vachon, 2020), whereas the latter remains technically difficult and can be influenced by non-local processes (entrainment from the shore or advection; Vachon et al., 2010;Esters et al., 2021). Thus, long-term direct flux measurements are mostly restricted to small lakes (Huotari et al., 2011), and fluxes remain mostly estimated with models. CO 2 fluxes at the surface of lakes op-P. Perolo et al.: Accounting for surface waves improves gas flux estimation at high wind speed in a large lake erate through a net diffusive transport, obeying the first Fickian law: where F (mol m −2 s −1 but often expressed as µmol cm −2 h −1 ) is the CO 2 gas flux, α is the CO 2 solubility coefficient (µmol cm −3 µatm −1 ), pCO 2 is the gradient of partial pressure of CO 2 (pCO 2 ) between the water and the atmosphere corrected for altitude (µatm), and k is the gas transfer velocity (cm h −1 ). Therefore, lake carbon emissions are primarily driven by the gradient of partial pressure of CO 2 between the surface lake water and the atmosphere, but the gas transfer velocity controls the rate of CO 2 exchange across the lakeatmosphere interface. Assessing the amount of lake CO 2 emissions to the atmosphere has been a major issue, starting with the Cole and Caraco (1998) seminal paper, with debates regarding both the representativeness of the measurements and the optimal conceptual model for air-water gas transfer (e.g. MacIntyre et al., 2001;Borges et al., 2004). As recent developments in sensor technologies allow continuous and accurate measurements of aqueous CO 2 concentrations, the gas transfer velocity currently remains the main source of uncertainties, which hinders attempts to achieve full carbon budgets (Dugan et al., 2016) or to quantify greenhouse gas emissions by lakes, at local, regional, or worldwide scales (Maberly et al., 2013;Raymond et al., 2013;Engel et al., 2018).
k is inherently tied to turbulent mixing within the surface boundary layer, which enhances the diffusive gas exchange by renewing the surface mass content (Zappa et al., 2007). At the lake-atmosphere interface, turbulent mixing is the product of wind shear (k u ), buoyancy flux (k c ), and wind-driven surface waves, whose effect can be split into wave action (k w ) and wave breaking (k b or k B ), with the latter producing air bubble and water spray ( Fig. 1; Lorke, 2003, Soloviev et al., 2007). Regarding the prominent role of wind action on surface turbulence, first quantitative models have empirically scaled k to wind speed (referenced at a 10 m height; U 10 ), as a proxy for the level of wind-driven turbulence ( Fig. 1; Cole and Caraco, 1998;Crusius and Wanninkhof, 2003). The parameterisations of the k-wind relationships vary between authors (e.g. Klaus and Vachon 2020), as a likely consequence of the local characteristics of the lakes used in the calibration datasets (Table 1). However, all studies suggested a polynomial relationship between U 10 and k with an exponent larger than 1. Further development of empirical models integrated the lake surface area as a second parameter in the k-wind relationships, to account for the role of the fetch length for wind action (Vachon and Prairie, 2013). Generally, empirical wind-based models tend to underestimate fluxes, especially at low wind speed (e.g. Schubert et al., 2012;Heiskanen et al., 2014;Mammarella et al., 2015) where turbulent mixing through buoyancy flux is expected to take over wind shear. Moreover, these empirical models require a proper calibration each time they are applied in a new system with different characteristics, i.e. a new set of lakes and/or meteorological conditions (Klaus and Vachon, 2020), thereby limiting their universal applicability.
In parallel to empirical wind-based models, process-based models attempt to link k directly to near-surface turbulence. The surface renewal model (SRM) is one of the first, and still most widely used, theories (Danckwerts 1951;Lamont and Scott, 1970) with k depending on the product of the turbulent kinetic energy dissipation rate (ε) and the kinematic viscosity of water (υ), both to a power of one-quarter as follows: where a 1 is a calibration constant parameter, and Sc is the Schmidt number. Recently, Lorke and Peeters (2006) and Katul and Liu (2017) demonstrated that this relationship, to which different approaches converge, can be seen as a universal scaling. As opposed to the practical empirical models presented above, process-based models have the potential to predict k using the turbulent dissipation rate over a wide range of environmental conditions extending beyond those encountered in the calibration dataset (Zappa et al., 2007). As for lakes, SRM k models have so far considered the friction velocity at the water side (u *, wat ) and the turbulence created by thermal convection using the buoyancy flux at the surface (B 0 ) ( Fig. 1; Eugster et al., 2003;Mac-Intyre et al., 2010;Read et al., 2012;Tedford et al., 2014;Heiskanen et al., 2014). It is noteworthy that the SRM approach leads to k being related to u *, wat (or U 10 ) to the first order (Wanninkhof, 1992;Lorke and Peeters 2006;see Sect. 2), whereas the empirical models described above predict a higher-order polynomial relationship. This inconsistency is tentatively solved in oceanography by adding another source of gas exchange associated with wind-waves' whitecaps. Early gas flux parameterisation already accounted for wind and buoyancy-driven turbulence as well as surface waves ( Fig. 1; Woolf et al., 1997;Soloviev et al., 2007;Fairall et al., 2011). However, the buoyancy-driven contribution can often be neglected in oceanography, and recent efforts have been dedicated to a better parameterisation of the bubble enhancement term ( Fig. 1; Deike and Melville, 2018). In lakes, wind fetch can be long enough to generate surface waves (Wanninkhof, 1992;Frost and Upstill-Goddard, 2002;Borges et al., 2004;Guérin et al., 2007), implying that surface waves could be a significant driver of k and subsequent CO 2 fluxes (Schilder et al., 2013;Vachon and Prairie, 2013). Thus, the role of surface waves has been essentially empirically accounted for in lake k models, through the polynomial scaling to U 10 in wind-based models, and most often neglected in studies using process-based parameterisations (mainly SRM) (e.g. Read et al., 2012). While this approximation may be appropriate for small-shielded lakes, it is likely to be insufficient in larger, long-fetched lakes.

Figure 1.
Conceptual scheme of the four main processes driving gas transfer velocity (k) in a large lake induced by wind and cooling events. These four processes are split into two types of k: k-bubble for the bubble formation (k B = k b ) and k-no bubble for the convective mixing, wind shear, and wave action term which are added (k NB = k c + k u + k w ). Below this scheme, a non-exhaustive review of the conceptual approaches of k models used in first Fickian law is given. From left to right, the increase in the complexity level of k models and their study site (limnological to oceanic case) is visible. All of these variables are described in Sect. 2.4 and Table 1.
Herein, we aim to identify the most adequate k model for Lake Geneva, a large, clear, hard-water lake in the Swiss Alps, to assess k values over a full annual cycle. We compare the performances of different models of gas transfer velocity, in their original or slightly modified published formulations from the limnological and oceanic literature. This set of models includes different levels of complexity, ranging from empirical models integrating wind speed and lake size to process-based models including wind shear, convection, and surface waves. Continuous pCO 2 measurements by in situ automated sensors and CO 2 fluxes, obtained from a new generation of automated (forced diffusion) flux chamber, were collected during specific periods of intensive field survey covering a wide range of natural conditions. Empirical k values computed from chamber data are then compared to outputs from the different k models. Owing to the size of Lake Geneva, we anticipate that models accounting, implicitly or explicitly, for the four key exchange drivers (i.e. wind shear, convective mixing, wave action, and bubble formation) will show the highest accuracy and precision in their estimation of k and that a precise integration of surface wave effects in such a large system should enhance model predictions. Thereafter, the relative distribution of these components is computed over a full year and analysed in the scope of the temporal variability of the gas transfer velocity. Finally, we expect that extreme wind and associated wave events should contribute disproportionately to accumulated k values over the year. In such a case, episodic weather events could generate large CO 2 fluxes over very short timescales that should be accounted for when computing annual CO 2 emission budgets.

Study site
Lake Geneva is a peri-Alpine lake defining part of the Swiss-French border, at 372 m a.s.l. (metres above sea level) (46 • 26 N, 6 • 33 E). Its surface area (582 km 2 ) and its maximum depth (309 m) make it the largest freshwater body in western Europe, with a volume of 89 km 3 (Fig. 2). Lake Geneva is monomictic. The two prevailing winds are almost diametrically opposed and come from the southwest and northeast respectively (Fig. 2). The lake water has been surveyed monthly or fortnightly since the late 1950s (OLA-IS, AnaEE-France, INRAE of Thonon-les-Bains, and CIPEL; Table 1. Summary of the characteristics of k Sc models for predicting the air-water gas transfer velocity based on wind speed (CC98 and CW03) and lake size (VP13), surface renewal model (T14, R12, and S07), COARSE approach (DM18), and both adapted models, namely SD21 and SD21-fit, from a combination of S07 and DM18. Mass balance by gas tracer Lake Area (0.128 km 2 ) U 10 (< 6 m s −1 ) VP13 k 600 = 2.51 + (1.48 · U 10 ) + (0.39 · U 10 · log 10(Lake size)

Field data at LéXPLORE
All field data were collected from the LéXPLORE platform, a 10 m by 10 m pontoon equipped with high-tech instrumentation and installed on Lake Geneva in 2019 (Wüest et al., 2021). LéXPLORE is moored at a 110 m depth, 570 m off the northern lake shore (Fig. 2).
On LéXPLORE, local weather conditions (air temperature, wind speed and direction, relative humidity, shortwave radiation, and atmospheric pressure) were continuously recorded (at 10 min intervals) by a Campbell Scientific auto-matic weather station. Lake surface temperature was measured every minute at 50 cm depth using a Minilog II-T (VEMCO, resolution 0.01 • C). The partial pressure of water surface CO 2 (pCO 2 ) was also measured at 50 cm depth during specific surveys (see Sect. 2.3) using a mini CO 2 sensor (Pro-Oceanus Systems Inc.) with an accuracy of ± 30 ppm. Values of pCO 2 in parts per million (ppm) were converted into micro-atmospheres (µatm) following the basic equation correcting for altitude (Russell and Denn, 1972). Therefore, we assume that the concentration and the temperature are homogeneous over the first 50 cm.
Fetch distance (m) from LéXPLORE to the lake shores considering wind direction was computed using data from the Federal Office of Topography online portal (Swisstopo geoportal: https://www.geo.admin.ch/, last access: 14 October 2020). The position of LéXPLORE is particularly relevant for this study as the fetch ranges from ∼ 0.5 to ∼ 30 km for the two prevailing winds. Significant wave height, H s (in m), was computed after Hasselmann et al. (1973) according to the following equation: where g is the gravitational constant. This variable H s is defined as the average height of the highest one-third of the waves (crest to trough) corresponding to the thickness over which the wind can push laterally (Wüest and Lorke, 2003). This equation is equivalent to the formulation by Carter (1982) that is more widely used in the oceanic literature. Simon (1997) tested the model for significant wave heights in Lake Neuchâtel (a lake close to Lake Geneva) with a fetch distance of 9 km. These results showed that the significant wave height in this lake was consistent with this oceanic formulation. However, Simon (1997) highlighted that the Joint North Sea Wave Project (JONSWAP) wave breaking parameterisation did not hold for winds greater than 5 m s −1 , producing faster wave breaking and with a higher probability in the case of surface waves that were not fully developed. Such lake waves are characterised by steeper slopes that favour their wave breaking and wave action (Wüest and Lorke, 2003). The net CO 2 flux at the lake-atmosphere interface, F , was directly measured with an automated (forced diffusion) floating CO 2 flux chamber (eosFD, Eosense: environmental gas monitoring; Fig. A1; Risk et al., 2011), which was originally developed for soil flux studies. The flux chamber had a detection limit close to 0.05 µmol cm −2 h −1 and measured F every 15 min in summer and 30 min in winter for battery-saving purposes. The standard floating chambers require quiet surface conditions (e.g. Cole et al., 2010;Vachon et al., 2010;Bastviken et al., 2015), thereby limiting studies to low to moderate wind speed conditions. One typical problem with floating chambers arises from the possible atmospheric leakage under rough surface conditions (Fig. A2a). To work around this problem, Vachon et al. (2010) recommend the creation of 10 cm long-edges entering the water (Fig. A2b), and this design also reduces artificial turbulence generated by the chamber's walls at surface. A second typical issue with this method is potential flux enhancement by artificial (chamber-generated) turbulence. This was also studied by Vachon et al. (2010), who demonstrated that the overestimations due to this effect can be as high as 1000 % at low wind but less than 50 % when the wind speed exceeds 4 m s −1 in large lakes. At even higher wind speed, this overestimation should decrease further because the surface water turbulence becomes much greater than that produced by the floating chamber. Thus, our flux chamber was specifically conceived to increase stability under calm and windy conditions and to limit artificial turbulence, but we do not exclude a bias under low and moderate wind conditions (Figs. A1, A2c).
Regarding the operation of the eosFD, it has two independent cavities: one for the chamber and one for the atmosphere (Fig. A1). These are connected to the same CO 2 sensor by a pump which sends either the chamber gas or the air gas to the sensor at regular intervals (about 20 s) and then completely flushes the chamber cavity according to the programmed measurement time step (15 or 30 min). Therefore, the advantage of this new instrument is to have a constant monitor-ing of the chamber's variation but also of the atmosphere. In addition, the use of the same CO 2 sensor for the two measurements limits the need for intercalibration between CO 2 sensors. We tested the performance of the floating chamber by comparing the standard deviation of the CO 2 concentrations of the atmosphere and in the chamber estimated from two separated cavities ( Fig. A1; Risk et al., 2011). We did not observe any difference in the standard deviation between high and low wind conditions (Fig. A3), suggesting that the measured fluxes remained reliable at high wind speed without leakage of the chamber.
We assessed the performances of our flux chamber during five specific periods over the annual cycle: 13-14 June 2019, 27-28 August 2019, 1-5 October 2019, 18-20 December 2019, and 20-26 February 2020. To select the most robust dataset for comparison with k estimates derived from models, we discarded flux data that were below the detection limit as well as CO 2 gradients that ranged within the uncertainty of sensors (i.e. ± 20 ppm for air and ± 30 ppm for water, resulting in a ± 50 ppm uncertainty) (Fig. A2). Accordingly, we were able to retain the most robust data points during the following deployment periods: 18-20 December 2019 and 20-26 February 2020. Finally, all these field data were standardised at a 1 h time step.

Computed k values from field data
k values (cm h −1 ) from field observations (k obs ) were computed from the gas transfer velocity equation: where F is the measured CO 2 flux (µmol cm −2 h −1 ); α is the gas solubility coefficient (µmol cm −3 µatm −1 ), which depends on the measured water temperature (Wanninkhof, 1992); and pCO 2 is the differential of pCO 2 measured at 0.5 m below the surface (pCO 2-water ) and pCO 2 at saturation (pCO 2-sat ; in ppm) measured from the flux chamber corrected by altitude (µatm). It is noteworthy that the chemical enhancement factor (Wanninkhof and Knox, 1996) was not considered in this equation, as the fluxes retained corresponded to conditions of moderate pH (i.e. < 8). k obs was then standardised in k 600 using the dimensionless Schmidt number (Sc) of CO 2 : k 600−obs = k obs · (600/Sc) −1/2 (600 for freshwater standardised at 20 • C).

Models for air-water gas transfer velocity
After years of debate, a consensus has begun to emerge on the relationship linking k, the intensity of turbulence, and Sc (Eq. 2), even when starting from different physical assumptions (see Katul et al., 2018). In this study, we selected six parameterisations widely used in limnology and oceanography, combining specific calibration characteristics (Table 1). We first show that they can all be expressed following Eq. (2) for wind shear and convection, despite their different formulations. We then develop the effects of surface waves from oceanic models and adapt the wave action for a large lake. The final lake model integrating the wave effect is ultimately calibrated using our field data (Table 1).

Wind shear stress
We start with the case where near-surface dynamics are driven by a weak to moderate wind, in the absence of heat exchange. In this case, the contribution of surface waves can be neglected and the wind stress (τ 0 = ρ air · C 10 · U 2 10 , where ρ air is air density and C 10 is the drag coefficient at 10 m) is equal to the tangential shear stress (τ t = ρ wat · u * ,wat 2 , where ρ wat is water density). The relationship between ε and the sheared velocity on the water side, u * ,wat , is then derived from a lawof-the-wall scaling for the velocity profile: ε = u 3 * ,wat /κz (0), where κ is the von Kármán constant (= 0.41), and z(0) is the thickness of the diffusive boundary layer. This relationship leads to The challenge is then to define z(0). Tedford et al. (2014) followed an ad hoc observational approach and chose z(0) = 0.15 m as the shallower depth where ε was measured. In contrast, theoretical studies have linked z(0) to the thickness of the diffusive or viscous sublayer (∼ 0.1-1 cm). In line with theory, we scale this layer as z(0) = cν/u * ,wat (Wüest and Lorke, 2003;Lorke and Peeters, 2006) with c as a constant value. Taking c = 114 (Soloviev et al., 2007), the thickness of this layer typically ranges from 0.04 to 0.14 m under a wind regime of 10 to 1 m s −1 . Therefore, we modify Eq. (5) to compute the interfacial (no bubble, NB) exchange coefficient: or k NB = a 1 (ρ air /ρ wat )C 10 U 10 (1/κc) 1/4 Sc −1/2 .
These equations show that the SRM formulation (Table 1, Fig. B1; Soloviev et al., 2007;Read et al., 2012) is analogous to the Coupled Ocean-Atmosphere Response Experiment Gas transfer algorithm (COAREG) flux (Fairall et al., 2011) and to the formulation used in Deike and Melville (2018) with the sheared velocity on the atmosphere side, u * ,atm (Table 1: DM18). Indeed, when equating the expression by Deike and Melville (2018), With Eq. (6a), we find that the coefficient a 1 = 0.29 in Soloviev et al. (2007) and Read et al. (2012) is essentially equivalent to the coefficient A NB = a 1 (κc) −1/4 (1/600) 1/2 (ρ wat /ρ air ) −1/2 ≈ 1.5 × 10 −4 in Deike and Melville (2018) (Fig. B1), which, in turn, was found to be equal to the coefficient of A = 1.5 in Fairall et al. (2011). These results agree with Lorke and Peeters (2006), who derived a unified relation for interfacial fluxes (air-water and water-sediment) through a linear relationship of u *, wat with k, especially at the bottom interface where shear is the only relevant process. Furthermore, Eq. (6) has a similar (i.e. quasi-linear) k-wind relationship to the data-driven parameterisation from VP13 but cannot explain the higher-order polynomial relationship reported in CW03 and CC98.

Convection
A second source of dissipation at the surface is the convection (ε c ) resulting from surface cooling. In the SRM formulation, only the negative buoyancy flux is considered when this term directly enters into the turbulent kinetic equation as a production term. The combination of wind shear and free convection near a boundary is described by the Monin-Obukhov similarity theory (MOST) with a general form derived from a turbulent kinetic energy balance (Lombardo and Gregg, 1989;Tedford et al. 2014): where L MO is the Monin-Obukhov length scale defined as L MO = u 3 *, wat /κB 0 , including ε(z) = c u ·ε u +c c ·ε c in Eq.
(2). The latter expression can be rearranged as follows: where a 1 ranges in the literature from 0.2 to 1.2 (Soloviev et al., 2007;MacIntyre et al., 2010;Tedford et al., 2014;Heiskanen et al. 2014;Winslow et al., 2016), c u ranges from 0.84 to 1 , and c c ranges from 0.37 to 2.5 (Wyngaard and Coté, 1971;Tedford et al., 2014). Hereafter, we use the following set of values: c u = 1 and c c = 1. Fairall et al. (2011) used an essentially equivalent approach but formulated in terms of a Richardson number to describe the partitioning between dissipation from convection and wind shear, expressing the wind shear in terms of the air-side friction velocity: R f = B 0 υ/u 4 *, atm , which can be integrated into Eq. (9) as where . The details of this demonstration can be seen in Soloviev and Schlüssel (1994).

Wave action
The effect of surface waves is commonly implemented in oceanography but barely considered in limnology. All process-based models rely on the same parameterisation of energy dissipation by wind shear and convection. However, they differ in how they parameterise energy dissipation by wave action and wave breaking.
The contribution of the wave action (ε w ) is, accordingly, added as a third source of turbulence (Fig. 1). In the presence of surface waves, the balance between τ t and τ 0 no longer holds. Therefore, Soloviev and Schlüssel (1994) added a corrective factor, ϕ, using the Keulegan number (Ke = u 3 *, wat /(gν)), in order to decrease the component τ t = τ 0 · ϕ (where ϕ = 1/(1 + Ke/Ke c )), with the critical Keulegan number (Ke c ) defined in Soloviev and Lukas (2006). As a result, the equation for shear-driven dissipation ε u (z) is as follows: Following this step, the turbulent kinetic energy dissipation rate from wave action (ε w ) is added and defined with the Keulegan number by Soloviev et al. (2007) as follows: where C T = (z 0 /H s ). z 0 is the surface roughness scale from the water side, and the C T value is set as a constant at 0.6 (more details are given in Soloviev et al., 2007). This definition does not hold for closed basins because, in the case of incompletely developed waves, the dissipation of energy from wind shear transmitted to the waves is not fully redistributed in the water body (Simon, 1997). Hence, for the application in Lake Geneva, we followed Terray et al. (1996), who defined a varying C T : where C p is the peak speed of the wave spectrum defined in Deike and Melville (2018) according to Toba (1972Toba ( , 1978. This leads to C T 1. This allows one to increase the effect of ε w (inversely proportional to C T in Eq. 12) on k. Here, we used this formulation to adapt the S07 ocean model for a large lake (closed basin). Henceforth we refer to the adapted uncalibrated and calibrated models as SD21 and SD21-fit respectively (Table 1). Finally, these three terms of ε (ε u , ε c , and ε w ) can be added before computing the SRM (Eq. 2) for determining k-no bubble (k NB ).

Bubble enhancement
Additional deviations from the linear relationship to U 10 are explained by the gas transfer resulting from bubbles and sprays during wave breaking. This mechanism is accounted for by adding a k-bubble (k B ) term to the previously mentioned k NB . Soloviev et al. (2007) used the empirical kbubble parameterisation from Woolf et al. (1997): where W is the fractional whitecap coverage only expressed as a function of wind (3.84 × 10 −6 · U 3.41 10 ), and O s is the Ostwald gas solubility. This formulation does not take wave height into account (Fig. B1). Nevertheless, a recent study (Deike and Melville, 2018) performed a new numerical process-based parameterisation for gas transfer velocity from bubble enhancement considering H s through the following equation: where A B is an empirical factor with dimension (= 10 −5 m −2 s 2 ), and O s is defined by the ideal gas constant (R), the surface water temperature (T 0 ), and the CO 2 solubility coefficient in freshwater (α) (Reichl and Deike, 2020). The gas transfer velocity is the expressed as a sum of the no bubble k NB and bubble k B components (Table 1: S07, DM18, SD21, and SD21-fit) following Keeling (1993) and Woolf et al. (1997Woolf et al. ( , 2005. Our adapted model for lake includes a refined parameterisation of the wave action term ε w from S07 along with the bubble term from DM18. For these reasons, the model will be called SD21 in the rest of the paper. In addition, for the appellation SD21-fit, the a 1 parameter of Eq. (2) and the A B parameter of Eq. (15) were fitted to the k 600 observations (a 1 = 0.33 and A B = 3 × 10 −5 m −2 s 2 ). With this review of existing and adapting parameterisations, we show that (i) there is a discrepancy between an SRM-based model with shear stress as the only energy source and empirical parameterisations with a polynomial (order > 1) wind-based relationship. Such a discrepancy is tentatively resolved by adding the effect of convection and surface waves. (ii) We further highlight that most fitting parameters from the different SRM-based models are in good agreement. (iii) We finally recall that it is possible to provide a unifying parameterisation of k with an SRM model including wind shear, wind-induced waves, and convection with only a few input parameters such as U 10 , B 0 , and Fetch.

Observed and predicted k
After quality check, our dataset contains 94 discrete CO 2 flux observations. We first assess the representativeness of our sampling by comparing the survey-specific and annual distributions of the three main inputs for k models: U 10 (all models), B 0 during convective periods (T14, S07, SD21, and SD21-fit), and H s (S07, DM18, SD21, and SD21-fit) (Fig. 3; the temporal evolution of these three terms is shown in Fig. C1). From 13 June 2019 to 12 June 2020, the average wind speed over Lake Geneva is 2.9 m s −1 with a mode at 2.5 m s −1 ; very low wind speeds (< 1 m s −1 ) are encountered 12 % of the year, whereas high (> 5 m s −1 ) to very high (> 10 m s −1 ) wind events represent 15 % and 2 % of the year respectively. The sampling surveys covered the full annual range of U 10 . Average and modal values of B 0 over the year are close to 0.25 × 10 −7 m 2 s −3 . However, the sampling covered only the lowest 50 % of the annual distribution and undersampled conditions of potentially strong convection. Considering that the dissipation by buoyancy flux, as parameterised in the process-based models, is already well known in the literature and that it is not the central point of our study, we posit that the undersampling of B 0 is not expected to significantly affect our analysis. The predicted modal H s value is 0.15 m over the year. Events of high H s (> 0.4 m) represent 6 % of the year, with a maximum H s of 1.1 m. As for U 10 , the surveys covered the full range of annual H s .
k 600 values based on observations are shown in Fig. 4a with their error bars corresponding to the uncertainties of the pCO 2 in air and in water (± 50 ppm). We notice that all of the measurements with a wave height > 0.4 m were observed for wind speeds > 5 m s −1 , and the corresponding k 600 values are located above the linear function (i.e. from a linear regression against wind shear velocity; Fig. B1) scaling k 600 to u * (i.e. first-order relationship). We then compare the k 600 observed during the specific surveys to the values computed with all k 600 models throughout the annual cycle, in relation to U 10 (Fig. 4b, c). Table 2 provides the rootmean-square errors (RMSEs) for all model estimates compared to k obs during the flux surveys, (i) for the full dataset (All Wind), and split (ii) for low wind (< 5 m s −1 , LW) and (iii) strong wind conditions (≥ 5 m s −1 , SW). The three empirical wind-based models only depend on wind (Fig. 4b).
Both CC98 and CW03 were originally calibrated for small lakes, using a mass balance calibration method (Table 1). However, they lead to divergent gas transfer velocities, particularly above 5 m s −1 , illustrated by a RMSE for SW as high as 22.8 cm h −1 for CC98, whereas CW03 performs better (RMSE SW = 12.8 cm h −1 ). Furthermore, both models underestimate k 600 at low wind (Fig. 4a), with a higher deviation for CW03 (Table 2). The k values predicted by VP13 are closer to those of the process-based models that explicitly  (c) k 600 process-based models (T14, S07, DM18, SD21, and SD21-fit) computed with annual data. Observed k 600 derived from CO 2 flux chamber measurements is shown using the "+" symbol.
integrate wave actions (S07 and DM18) (Fig. 4b, c), demonstrating that lake size integration in the empirical model captures at least part of the wave action on k. Performances of VP13 at strong winds (RMSE SW = 12.7 cm h −1 ) were better than those of the ocean-derived models integrating surface waves (RMSE SW = 13-15.9 cm h −1 ). However, VP13 shows a positive offset during calm periods, along with the highest RMSE of the set of models at low wind speed (Fig. 4b, Table 2).
The process-based models (Fig. 4c) provide different k 600 values for a given wind speed, owing to the integration of additional environmental components (i.e. the varying drag coefficient, the convective mixing in R12 and T14, and the effect of waves in S07, DM18, SD21, and SD21fit). All process-based models are similar at low winds, as they share a common physical basis for parameterisation of wind shear and convection. Therefore, they lead to similar RMSEs (2.9-3.5 cm h −1 ) under such conditions, where surface waves are negligible (Table 2). Divergences occur at higher wind speeds. T14, initially developed for small lakes with limited wind exposure, performed the worst (RMSE = 19.8 cm h −1 ). This increased k underestimation at high winds can be attributed to (i) dissipation by wave action and bubble formation not being considered (in R12 and T14) and (ii) to the use of a constant z(0) = 0.15 m in the T14 model (Eq. 5). This approximation of the diffusive layer is consistent with low wind speed but is almost 1 order of magnitude too large under strong wind speed. Other processbased models, designed for greater wind range (> 10 m s −1 ), integrate surface waves and, as a result, lead to better estimates than R12 and T14 (RMSE = 10.4-15.9 cm h −1 ). However, the ocean wave model of DM18 shows lower performances at strong winds than CW03 and VP13 (Fig. 4b, c). Finally, the specific fit parameterisation of the SD21-fit model improves the performance at high wind speeds by ∼ 30 % (RMSE = 10.5 cm h −1 ), outperforming all the other methods.  (d) SD21 is similar to S21 but the k-bubble term of DM18 is added; (e, f) SD21-fit is similar to SD21 but with a 1 and A B fitted to observed k.

Surface wave integration
Herein, we scrutinise how those varying parameterisations ultimately alter the shape of relationship between k 600 and U 10 . In R12 (Fig. 5a), wind is only included through wind shear, resulting in a linear relationship between k 600 and U 10 , as already anticipated. Adding the wave action (no bubble) through the S07 parameterisation (Fig. 5b) does not lead to any significant departure from the minimal R12 model. Adapting wave action by decreasing C T (for the no bubble term) leads to a departure from the wind shear linear relationship for H s > 0.4 m (Fig. 5c).
Adding the k-bubble term related to wave breaking of DM18 further increases this deviation from the linear k 600 -U 10 relationship but also scatters k 600 estimates for a given U 10 (Fig. 5d). Finally, the fitting with observationally based k 600 improves the estimation for strong wind (Fig. 5e, Ta-ble 2). Given the range of wind fetch from ∼ 0.5 to ∼ 30 km, the contribution of waves varies for a given wind speed depending on the fetch, as evidenced by the scattering of the parameterised k 600 for a given U 10 (Fig. 5f). A significant modification of k 600 by wave action and wave breaking occurs for a fetch length > 15 km and U 10 > 5 m s −1 (Fig. 5f), in the case of Lake Geneva, generating wave of H s > 0.4 m (Fig. 5e).
Compared with k 600 estimated by R12 (Fig. 5a), the SD21 and SD21-fit models provide k estimates that are 20 %-50 % higher for U 10 = 10 m s −1 respectively and 40 %-70 % higher for U 10 = 15 m s −1 respectively. Therefore, adapting the surface waves, through the change in the wave action for incompletely developed waves and the fitting to observed data encountered in local lake conditions, leads to better performances of the SD21 models. SD21-fit reached the lowest RMSE at all wind speeds and was thereafter used as a reference for the modelling of the annual gas transfer velocities.

Annual cumulative gas transfer velocity and the effect of extreme conditions
We show above that models accounting for surface waves better represent the non-linear increase in k 600 at high winds. Because high-wind events remain rare, we test whether a better representation of k 600 during rare, high-wind events affects the local estimates of k 600 over a full year. To this end, cumulative sums of hourly k 600 were computed over a full annual cycle (13 June 2019-12 June 2020) for all k models (Fig. 6). The annual dynamics, such as the annually averaged k 600 , were compared using SD21-fit as a new reference model. Cumulated k 600 computed for SD21-fit shows some episodic steep increases between December and March, due to the winterly prevalence of high-wind events (the winter average wind speed, 3.25 m s −1 , was greater than the summer mean, 2.55 m s −1 , by 25 %) and greater significant wave height (the winter average wave height, 0.15 m, was greater than the summer value, 0.10 m, by 50 %) (Fig. C1). The average hourly k 600 by the SD21-fit model is 7.3 ± 7.4 cm h −1 (mean ± SE, Fig. 6a). Periods of high winds, although accounting for only 15 % of data points, contribute 44 % of annually cumulative k 600 in the SD21 and SD21-fit models (Fig. 6b), whereas periods of high waves (H s ≥ 0.4 m), accounting for only 6 % of data points, contribute to more than 20 % of annually cumulative k 600 . The wind-based models are those for which cumulative k 600 diverges the most from the SD21-fit reference model, with the lowest annual averaged k 600 for CC98 and CW03 (3.9 ± 2.7 and 4.8 ± 9.3 respectively) and the highest for VP13 (9.9 ± 6.1). These divergences arise from the low performances of these models at low wind regimes (Figs. 4a, 6b; Table 2), which represent 85 % of annual data points. All of the other process-based models have relative dynamics of cumulative k 600 similar to that of the SD21-fit model and end up with annually averaged k 600 that are 15 % lower than for the SD21-fit. The representation of k 600 at low wind speeds is similar for all processbased models, and the divergence arises from the representation of the rarer high-wind-speed episodes, which contribute to 43 %-46 % of annual cumulative k 600 (Fig. 6b).

Discussion
The history of k models, simulating the gas transfer velocity for surface waters, dates back from the early 1990s. k models have been developed and tested in small lakes sheltered from winds (e.g. Crusius and Wanninkhof, 2003;Tedford et al., 2014), large lakes under low to moderate wind speed (Vachon and Prairie, 2013), and oceans (e.g. Soloviev et al., 2007;Fairall et al., 2011;Esters et al., 2017;Deike and Melville, 2018). While the effects of surface waves on k can be neglected in small lakes, we question whether this assumption holds for large lakes such as Lake Geneva, in which surface waves are frequently observed (Figs. 2, C1). We evaluated the performance of different experimental-based and process-based models to estimate k 600 in the large Lake Geneva. We show that integrating the effect of wave formation at high wind speeds and long fetch better represents the sharp increase in the k 600 values during such episodic windy events.

Choice of k models
Wind-based models have long been known to misestimate k 600 at low wind speeds (Eugster et al., 2003;MacIntyre et al., 2010;Erkkilä et al., 2018). Consistently, wind-based models showed the lowest performances for Lake Geneva, especially at low wind speeds (CW03 and VP13), which resulted in large discrepancies in annually averaged and cumulative k 600 over the full year. They are, however, easy to 1180 P. Perolo et al.: Accounting for surface waves improves gas flux estimation at high wind speed in a large lake compute, require few inputs (only U 10 ), and remain by far the most used to estimate lakes' CO 2 emissions worldwide (e.g. Raymond et al. 2013). Another possibility is to broadly adopt process-based models. The presented process-based models require input data that are currently more easily accessible and are routinely acquired at high frequency in many lakes: wind speed, heat flux, and wind fetch (i.e. distance from the shore). The development of R packages, such as "LakeMetabolizer" , in which the calculations of process-based models are implemented also alleviates their computational difficulty. Both increased data availability and computational tools should foster the use of process-based k models, which hold great potential to obtain more accurate global k 600 estimates.
The analysis of the models adapted from the existing literature to account for the effect of lake surface waves, SD21 and SD21-fit (Fig. 5d, e, f), showed that the wave contribution to k becomes significant for H s > 0.4 m, corresponding, for Lake Geneva, to winds blowing at 5 m s −1 from the southwest where the fetch length is maximal (> 15 km) with respect to the measurement site. A significant contribution to the gas transfer velocity by surface waves is expected in lakes where H s > 0.4 m is not infrequent. Wave heights beyond this threshold value of H s are frequently encountered in lakes that are larger than or a similar size to Lake Geneva (6 % of annual time at the LéXPLORE platform). In the Great Lakes of North America, Hubertz et al. (1991) showed that the mean wave heights of all of these lakes were > 0.4 m in summer and close to 1 m in winter with a maximum of up to 5 m. H s > 0.4 m can also form over elongated lakes of smaller size, such as smaller Swiss lakes (e.g. Lake Neuchâtel and Lake Bienne) (Amini et al., 2017). As SD21-fit is a process-based model integrating the four main processes in a mathematically coherent way, we would expect that it can be applied to lakes experiencing H s > 0.4 m and can improve the accuracy of k estimates. Because waves can physically damage inshore and offshore infrastructure, many large lakes benefit from wave forecasts. H s data from those forecasting systems (e.g. National Data Buoy Centre -NOAA and Wave Atlas from SwissLakes.net; Amini et al., 2017) could allow one to test whether the SD21-fit models can be applied to those lakes and whether k NB and k B through a 1 and A B need to be recalibrated or fitted to the local context if flux measurement data are available, as for this study. Energy dissipation during high-wave events increases the gas transfer velocity well beyond the linear relationship derived for wind shear alone. Therefore, we expect that computed gas fluxes at the air-water interface should be significantly improved by the integration of surface waves into the k models.
4.2 Implication of four components on the annual k estimation and the annual CO 2 fluxes

Seasonal and hourly distribution of k CO 2
Converting k 600 to k CO 2 using the Schmidt number (Wanninkhof, 1992) highlights the importance of water temperature in gas exchange dynamics. Indeed, the seasonal distribution of the cumulative k 600 is ∼ 20 % and ∼ 30 % for the warm (spring and summer) and cold (fall and winter) seasons respectively. Once the temperature effect is accounted for, this distribution increases to 26.1 % for summer and decreases to 24.9 % for winter but remains unchanged for spring and fall. While R12 only uses wind shear and convective terms, the selected process-based model (SD21-fit) allows a decomposition into the four main drivers of the gas transfer velocity, thereby paving the way to a better understanding of the implication of these processes throughout an annual cycle. Wind shear remains the dominant component of the gas exchange velocity over the different seasons (Fig. 7a). The annual contribution of surface waves (wave action and bubble formation) is limited to 9 %-10 % of the cumulative k in fall and winter. The contribution of the buoyancy flux at the surface to k is even smaller for both models (R12 and SD21-fit) at this seasonal scale. However, both the buoyancy flux and the surface waves can significantly increase k during episodic events, during which they can contribute disproportionately to k at hourly (up to 80 % for convection) and daily (up to 25 % for surface waves) timescales (Fig. 7b). Several studies have emphasised the disproportionate contribution of episodic mixing events to the annual flux, bringing CO 2 back to lake surfaces such as after ice break in dimictic lakes (Karlsson et al., 2013;Finlay et al., 2019) or during fall mixing in a eutrophic deep lake (Reed et al., 2018). Processbased k models integrating both the buoyancy flux and the wind-induced waves offer the opportunity to mechanistically investigate how much these episodic events contribute to annual emissions through short-term modifications of the gas exchange velocity.

Consequences of the choice of k model on the seasonal to annual CO 2 flux estimation
We produced coarse estimates of monthly CO 2 fluxes with the objective of scaling the effects of wave integration at seasonal and annual scales. Monthly fluxes were computed based on k estimates at LéXPLORE from the different models at an hourly time step as well as the monthly average of water temperature and recorded pCO 2 at the lake surface (OLA-IS, AnaEE-France, INRAE of Thonon-les-Bains, CIPEL; Rimet et al., 2020;Perga et al., 2016) as well as a constant pCO 2 in the atmosphere (400 µatm). For months for which surface pH > 8.4, k values were computed with and without considering the chemical enhancement (CE; Wanninkhof and Knox, 1996) (Table 3). The dependency of the Figure 7. (a) Distribution of k CO 2 generated by two main processes (k u and k c ) in R12 and four main processes (k u , k c , k w , and k b ) in SD21-fit for each season: spring (April-May-June), summer (July-August-September), fall (October-November-December), and winter (January-February-March). The height of the bar represents the cumulative k CO 2 by season for both models (R12 and SD21-fit). (b) Distribution of four k generated by wind shear, convection, wave action, and bubble enhancement (k u , k c , k w , and k b respectively) along the annual cycle.
For the SD21-fit model, Table 3. Seasonal to annual CO 2 flux estimation (mmol C m −2 d −1 ) from k models (with, denoted using "-CE" (italic font), and without (roman font) chemical enhancement), the monthly CO 2 average (µatm), and the models' deviation from SD21-fit. CE was only considered for seasons when pH > 8.  (Wanninkhof and Knox, 1996) might generate a further uncertainty in estimated CO 2 fluxes (e.g. with a greater k value being related to a lower chemical enhancement factor). As predicted by Fick's law, the highest outgassing fluxes occur in fall and winter, when water mixing brings CO 2 up to the lake surface, whereas low gas uptake fluxes occur in spring and summer, when primary production depletes surface CO 2 below saturation. However, annual estimates of net CO 2 outgassing vary from 14.7 to 37.1 mmol C m −2 d −1 (Table 3) depending on the k model used for computation. Consistently, differences between model estimates are relatively low in summer, as both the pCO 2 gradient (100-200 µatm) and wave occurrence are limited. Adding CE during the spring and summer months causes an increase in influx of the order of 5 %-128 % depending on the k model used. For example, using CC98 in summer leads to a 93 % increase in the estimated influx (compared with fluxes without CE), whereas the increase is only 39 % for SD21-fit. This chemical enhancement factor deserves more attention in future studies, especially for high pH lakes and summer seasons. However, the variability introduced by the CE at the annual scale remains low (∼10 %) compared with that introduced by the choice of the k model. Estimated fluxes are also strongly dependent on the chosen k model in winter when both pCO 2 (475 µatm) and surface wave occurrence are higher (Fig.7, Table 3). Therefore, while high-wave events represent only 6 % of the total surface wave occurrence (H s > 0.4 m), an incomplete consideration and description of their contribution may lead to an annual flux underestimation of about 20 %-25 %.
The weak contribution of convection is at odds with observations in small lakes, although not unexpected, as large lakes are exposed to stronger winds, such that wind-sheardriven ε u often outpaces convectively driven ε c (Read et al., 2012). However, the limited impact of the buoyancy flux on k does not rule out its contribution to CO 2 exchange. Indeed, convective mixing plays a central role in the deepening of the mixed layer, allowing the export of the CO 2 stored in the hypolimnion towards the surface during the cold period and, thus, controlling the pCO 2 gradient (Zimmermann et al., 2019) and the observed wintertime outgassing. Altogether, both surface oversaturated CO 2 concentrations (as a result of convective mixing) and wind-induced waves are more relevant in fall and wintertime for the monomictic Lake Geneva, leading to most of the annual outgassing during this season (Table 3). As for many monomictic lakes, these seasons drive most of the annual CO 2 budget of Lake Geneva (Perga et al., 2016), whereas they usually correspond to those where direct measurements are the scarcest. An improved quantification of k values through SRM models including wind-induced waves should contribute to refining the overall estimation of large lakes' contribution to regional CO 2 emissions. Wind fields from COSMO-1 for two episodes of northeast and southwest wind directions; (c, d) wave fields on Lake Geneva considering the two prevailing winds (northeast and southwest); (e, f) gas transfer velocity from SD21-fit; (g) boxplots of the spatial variability, at the lake scale, of k values computed from CC98, R12, and SD21-fit under both meteorological conditions. Diamonds represents the spatial mean, and the cross (+) represents the k value computed from the averaged fetch distance (NE: 9.5 km; SW: 9.3 km).

Wind and wave field on Lake Geneva and their
impact on the spatial integration of k 600 SD21 and SD21-fit were built on the basis of a single measurement point on the lake, just as for most of the existing k models. Therefore, the question of the extrapolation of the model to the whole lake remains essential. Herein, we showcase two snapshot situations of high wind (i) to illustrate how process-based models could enable spatially resolved estimates of k values and (ii) to exemplify how much k can vary as a result of the spatially variable wave and wind fields during a single episode. Two events of high and similar wind speed (11 m s −1 ) but different directions (NE on 30 March 2020 at 08:00 UTC and SW on 10 February 2020 at 06:00 UTC; Fig. 8a, b) were extracted from the 0.01 • hourly resolved numerical weather model of the Swiss Federal Office of Meteorology and Climatology (COSMO-1, MeteoSwiss). The two fetch distances from both prevailing wind directions (NE, SW) were then measured at each pixel of the grid (n = 583), from which the wave height field was mapped (Fig. 8c, d) considering Eq. (3) and the two wind grids. The maps were qualitatively consistent with previous studies on wind waves for Lake Geneva using the spectral wave model (SWAN) for wave height (Amini et al., 2016). Spatially resolved k values were computed from the fetch and wind grids, using the wind-based model CC98, the processbased model without lake wave implementation R12, and the SD21-fit model containing the lake wave parameterisation.
Taking H s > 0.4 m as a threshold for the significant effect of waves on k, the two prevailing winds show opposite wave responses, with long fetch and higher wave heights affecting either the northern or the southern shores for the SW and NE winds respectively. Under both conditions, more than 60 % of the full lake surface area experiences H s values > 0.4 m, leading to k values as high as 68 cm h −1 as computed by SD21-fit (Fig. 8g). The eastern part of lake experiences the lowest wind speeds and wave heights in both situations, as a consequence of the orographic effect of the Alps surrounding the Grand Lake.
Both the range and the mean of estimated k values increase with increasing model complexity. Accounting only for wind speed, through the wind-based model CC98, leads to a spatially integrated k value of 15.7 cm h −1 (range of 1-22 cm h −1 ). k values computed from R12, accounting for both the wind shear and buoyancy flux, are on average 55 % greater (24.3 cm h −1 ) with a moderate effect on the range of spatial variability. Finally, including surface waves results in a spatially averaged k value that is more than double the k value computed from wind speed only (35.5 cm h −1 ), with a variability that is almost 3 times greater (range of 0-68). It is noteworthy that the spatial average of the k values computed by the wind grids (Fig. 8g, diamonds) is equivalent to the average of k values computed from the spatial mean fetch (9.5 km for NE and 9.3 km for SW; Fig. 8g, crosses) under these specific weather conditions. Thus, the application of an average fetch would be relevant to estimate a spatially averaged k value.
For all models, the integration of spatially resolved wind fields may improve the accuracy of k at the lake scale, but accounting for wind only would underestimate both the average gas piston velocity and its spatial variability. Therefore, a better understanding of wave behaviour in large lakes, using different approaches such as field and laboratory measurements, new physical models, and technical development, would improve the accuracy of gas exchange estimates at the lake-water interface, at both temporal and spatial scales. Further, estimates of lake-scale CO 2 fluxes would also require one to account for the spatial variability of CO 2 . The question of the spatial variability of the CO 2 is still open and remains difficult to analyse at high frequency in large lakes.

Conclusions
Investigations of the four main processes generating the gas transfer velocity in the large Lake Geneva demonstrated the importance of considering surface waves during episodic windy events responsible for more than 44 % of annual cumulative k 600 . The in-depth study of the behaviour of the process-based models has enabled us to underscore their consistent predictions under low and strong wind conditions, especially considering the new combination and adaptation model, SD21-fit. This last model significantly improves the estimation of the CO 2 flux when these three thresholds appear in the field, U 10 > 5 m s −1 , Fetch > 15 km, and H s > 0.4 m, making it applicable to a wide range of lake sizes. Furthermore, SD21-fit is assembled on solid theoretical bases coming from limnological and oceanic literature and allows one to analyse the distribution of these four main terms (k u , k c , k w , and k b ) across a variety of timescales depending on the kind of study.
To conclude, the development of high-resolution gridded atmospheric models such as COSMO-1 is an asset for future estimates of the gas transfer velocity and the spatial heterogeneity of lake biogeochemical cycles. Moreover, this study sheds light on the complexity of large lakes located at the interface between small, sheltered lakes and the open oceans, experiencing a combination of processes relevant for both small and large systems. The possibility of using processbased models in a fairly simple way with few inputs to improve the precision of the gas transfer velocity and, therefore, the gas flux should be supported in future research. In addition, this approach is very promising with respect to uncovering long-term trends of CO 2 emissions from lakes as well as for finer estimation of fluxes during more intense episodic events.

1184
P. Perolo et al.: Accounting for surface waves improves gas flux estimation at high wind speed in a large lake Appendix A Figure A1. Schematics of eosFD operation (https://eosense.com, last access: 18 August 2020), its mini-platform construction, and its positioning for measurements in the field (Lake Geneva at LéXPLORE platform). The raft design also complies with recommendations to minimise artificial turbulence induced by the chamber's walls, with 10 cm long-edges entering the water (Vachon et al., 2010).   Appendix B Figure B1. Panel (a) shows the comparison of Soloviev et al. (2007) and Deike and Melville (2018) for the first-order function of friction velocity at the water side (u * ,wat ) (blue points) and at the atmosphere side (u * ,atm ) (green points) with their linear regression (black line), the linear function of Vachon and Prairie (2013) for a lake size of 582 km 2 (yellow points), and the linear regression from u * ,wat or u * ,atm . Panel (b) is a visualisation of S07 with empirical parameterisation of the bubble term (Woolf, 1997) regardless of wave height as a function of wind speed at 10 m. Panel (c) is a visualisation of DM18 as a function of wind speed, only using the effect of the bubble term from 10 m s −1 .
Financial support. This research has been supported by the Schweizerischer Nationalfonds zur Förderung der Wissenschaftlichen Forschung (grant no. SNF 200021_175530; CARBOGEN). Review statement. This paper was edited by Tom Shatwell and reviewed by two anonymous referees.