Bulk Transfer Coefficients Estimated From Eddy-Covariance Measurements Over Lakes and Reservoirs

The drag coefficient, Stanton number and Dalton number are of particular importance for estimating the surface turbulent fluxes of momentum, heat and water vapor using bulk parameterization. Although these bulk transfer coefficients have been extensively studied over the past several decades in marine and large-lake environments, there are no studies analyzing their variability for smaller lakes. Here, we evaluated these coefficients through directly measured surface fluxes using the eddy-covariance technique over more than 30 lakes and reservoirs of different sizes and depths. Our analysis showed that the transfer coefficients (adjusted to neutral atmospheric stability) were generally within the range reported in previous studies for large lakes and oceans. All transfer coefficients exhibit a substantial increase at low wind speeds (<3 m s −1 ), which was found to be associated with the presence of gusts and capillary waves (except Dalton number). Stanton number was found to be on average a factor of 1.3 higher than Dalton number, likely affecting the Bowen ratio method. At high wind speeds, the transfer coefficients remained relatively constant at values of 1.6·10 −3 , 1.4·10 −3 , 1.0·10 −3 , respectively. We found that the variability of the transfer coefficients among the lakes could be associated with lake surface area. In flux parameterizations at lake surfaces, it is recommended to consider variations in the drag coefficient and Stanton number due to wind gustiness and capillary wave roughness while Dalton number could be considered as constant at all wind speeds.


Introduction
The major process that governs the interaction between the atmosphere and surface waters is the turbulent exchange of momentum, heat and gases at the air-water interface.Although lakes and reservoirs occupy only about 3% of the land surface area (Downing et al., 2006), they are known to have an impact on local weather and climate.For example, lakes affect the stability of the atmosphere above (Sun et al., 1997), leading to the formation of clouds and precipitation on the shores (Changnon & Jones, 1972;Eerola et al., 2014;Kato & Takahashi, 1981;Thiery et al., 2016).Furthermore, lakes and reservoirs are recognized as significant contributors to the global carbon cycle by emitting significant amounts of carbon dioxide and methane (DelSontro et al., 2018;Rosentreter et al., 2021).
The past three decades have seen a rapid development of lake models (Stepanenko et al., 2014) and their incorporation into numerical weather and climate prediction models (Ljungemyr et al., 1996;Mironov et al., 2010;Salgado & Le Mogne, 2010).Experiments on the coupling of lakes and the atmospheric model revealed their beneficial impact on the weather prediction quality (Balsamo et al., 2012).A number of case studies have demonstrated the importance of lakes for extreme local weather phenomena, such as lake-effect snow over Great American lakes (Fujisaki-Manome et al., 2020), deep hazardous convection over Great African lakes (Thiery et al., 2016), wind speeds over Lake Superior (Desai et al., 2009), or stratiform cloudiness in winter over Lake Ladoga (Eerola et al., 2014).Thus, an accurate representation of the exchange of momentum, heat and water vapor at the air-water interface in water bodies is essential to understand those processes.
In state-of-the-art, momentum, sensible and latent heat fluxes are usually determined based on gradient approaches utilizing transfer coefficients (also known as bulk transfer coefficients) and easy to measure meteorological and limnological variables, that is, wind speed, air temperature, air humidity and water surface temperature (Stull, 1988).The exchange at the air-water interface and therewith the bulk coefficients are controlled by boundary-layer turbulence.The bulk exchange coefficient of momentum, known as the drag coefficient (  , DN ) (Garratt, 1977), is of particular importance for all air-water fluxes.The coefficients of heat (   , HN ) and water vapor exchange (  , EN ) are also known as Stanton and Dalton numbers, respectively.Here, "N" stands for "neutral" transfer coefficients, corresponding to the neutral thermal stability of the atmosphere.The transfer coefficients depend on the measurement height of the mean wind speed, air temperature and humidity, respectively, and for this reason, they are usually reported for the reference meteorological height of 10 m.
A considerable amount of studies has been published on the momentum flux and the drag coefficient starting from the early 1950s when the fundamental work, presenting the theory later on named as Monin-Obukhov similarity theory, was published (Monin & Obukhov, 1954;Obukhov, 1971).The theory aims at describing the structure of turbulence in the atmospheric surface layer about several tens of meters thick with the assumption of the fluxes being constant and independent of height.Similarity laws introduce functional relations to derive the universal shapes for the vertical profiles of different quantities for atmospheric thermal stability other than neutral.During the past decades, considerable effort has been devoted to define the exact form of these similarity functions (Businger et al., 1971;Högström, 1988;Paulson, 1970;Zilitinkevich & Calanca, 2010).
As the drag coefficient is one of the key parameters in atmospheric and lake models, the errors in its parameterization lead to errors in the bulk flux estimates.Therefore, numerous early studies focused on exploring different parameterizations of the drag coefficient over the land and oceans in terms of wind speed, atmospheric stability, and surface roughness, which are a function of the surface wavefield (for oceans) (Garratt, 1977;Kantha & Clayson, 2000).Most of the extensive field measurement campaigns over the oceans have been conducted during the last 30 years of the twentieth century (Fairall et al., 1996;Godfrey & Beljaars, 1991;Large & Pond, 1981;Smith et al., 1996).Several of these studies agreed that the drag coefficient linearly increases with increasing wind speed ignoring the state of the wavefield.More recent parameterizations of the drag coefficient (e.g., the COARE algorithm by Fairall et al. (2003); Edson et al. (2013)), however, include a wave dependence.There is still an ongoing scientific discussion concerning the importance of waves and how their impact could be included in the models (Wu et al., 2019).
Along with the studies in the marine environment, the research started to focus on the drag coefficient estimated from measurements over large and medium-sized lakes (e.g., Donelan, 1982;Graf et al., 1984;Hicks, 1972;Simon, 1997).To date, in total, about two dozen studies focusing on lakes have been published since the beginning of the 1970s.In reviewing these studies, we separated them by the wind speed regime they were interested in.

of 20
It is usually assumed that surface wave development starts when the wind speed exceeds 3-4 m s −1 (Ataktürk & Katsaros, 1999;Kantha & Clayson, 2000).This is also supported by wave measurements in several lakes (Guseva et al., 2021;Simon, 1997).Therefore, we intend to separate the two wind speed regimes using this threshold.
At the "high" wind speed regime (wind speed exceeds 3 m s −1 ), in the most simplified way, the surface waves are assumed to be fully developed, and the surface roughness length is described as a function of wind stress, which is commonly known as Charnock relationship (Charnock, 1955).However, this assumption might not hold for lakes with limited wind fetch (Donelan, 1990;Geernaert, 1990).Thus, some research has been made to study the drag coefficient as a function of the surface wave state, for example, taking into account wave characteristics such as the wave age (Ataktürk & Katsaros, 1999;Donelan, 1982).Vickers and Mahrt (1997) reported that for a given wind speed the drag coefficient tends to be larger for younger steeper waves representative of short wind fetches than for longer fetches.Ataktürk and Katsaros (1999) could significantly reduce the scatter in the estimated drag coefficients by considering waves in the parameterization of the surface roughness length.However, these studies mainly examined large lakes and only a few were performed in lakes with short fetch and young wave states (Babanin & Makin, 2008;Lükő et al., 2020).Given the fact that the surface wave measurements in lakes are often not available, their effect still could be investigated via analyzing the relationship between the drag coefficient and fetch length.
At the "low" wind speed regime, several studies found that the neutral drag coefficient in lakes and oceans tended to increase by an approximate factor of two up to ten compared to the value of 1.3•10 −3 (corresponding to a typical value of open water surface roughness, Foken, 2008) (Woolway et al., 2017;Wüest & Lorke, 2003).Although the wind speed dependence is obvious, many numerical and empirical studies employ a constant value for the drag coefficient, which is often considered as a model tuning parameter (Stepanenko et al., 2014).Despite the fact that there have been many attempts to address the reasons of such increase, there is still no consensus in the scientific community.The low wind speed regime was first described as the aerodynamically smooth flow, when the surface waves are buried within the viscous sublayer and the surface roughness is described as a function of the thickness of this layer (Schlichting, 1968).On the contrary, Wu (1988) proposed that the flow is aerodynamically rough and that capillary gravity waves play a key role at low wind speeds.Surface roughness length was described as a function of the water surface tension.As an additional reason for the increase of the drag coefficient at low wind speed, Godfrey and Beljaars (1991) and Grachev et al. (1998) considered the concept of gustiness, which assumes that at "zero" wind speeds there are dry random convective motions-gusts-in the convective boundary layer (CBL).Thus, the "traditional" formulation of the drag coefficient has been modified using the scalar-averaged wind speed (not the vector-averaged wind speed) to account for gusts.All the possible mechanisms mentioned above were addressed in the recent work by Wei et al. (2016).They concluded that none of them explained the increase of the drag coefficient at low wind speeds.However, they found it can be explained by the increase in the turbulent kinetic energy due to buoyant energy production.Similar to Grachev et al. (1998), Sahlée et al. (2014) and Liu et al. (2020) related the increase of the drag coefficient with nonlocal effects, such as the penetration of large convective eddies into the surface layer from the atmosphere above.Liu et al. (2020) introduced the factor describing this effect and estimated it from two-level measurements of wind speed (however, over the land surface and only for neutral conditions).Another formulation of the drag coefficient at low wind speeds was done by Zhu and Furst (2013) relating the drag coefficient to the turbulent kinetic energy budget.However, their fitting coefficients for the drag coefficient formula were found to be site-specific (Liu et al., 2020).
Other studies on the bulk transfer coefficients in lakes branched off from the main direction-potential physical mechanisms-with a focus on the possible correlation between the bulk transfer coefficients and some lake characteristics.Among them are lake depth at the measurement location (Panin et al., 2006), lake surface area (Read et al., 2012;Woolway et al., 2017), wind fetch at the measurement location (Lükő et al., 2020) and lake biota, for example, submerged macrophytes (Xiao et al., 2013).All studies showed a strong dependence of the transfer coefficients on these lake characteristics.The drag coefficient tends to decrease with increasing water depth, lake area, fetch and in the presence of water plants at the water surface.It is important to note that although Panin et al. (2006) and Woolway et al. (2017) revealed the correlation between the transfer coefficients and the lake parameters, the estimation of the transfer coefficients was based either on bulk parameterization (Woolway et al., 2017), or was compared to other studies where there were no direct flux measurements (Panin et al., 2006).
Fewer studies have been published on the Stanton and Dalton numbers.Although the measurements in the oceans showed their obvious increase at low wind speeds, both transfer coefficients were considered as fairly constant GUSEVA ET AL.

10.1029/2022JD037219
4 of 20 with a value of 1.1•10 −3 (review of these measurements in Kantha and Clayson (2000)).First measurements conducted in lakes revealed this value being higher and equal to ∼1.5•10 −3 (Harbeck, 1962;Hicks, 1972) or 1.9•10 −3 (Strub & Powell, 1987).Harbeck (1962) and Brutsaert and Yeh (1970) reported a dependence of the Dalton number on the lake surface area.Heikinheimo et al. (1999) summarized that the Dalton number is generally known to be less dependent on the wind speed.From the most recent studies (Dias & Vissotto, 2017;Li et al., 2016;Wei et al., 2016;Xiao et al., 2013), there is evidence that both coefficients depend on the wind speed and that the Stanton number is higher than the Dalton number by approximately a factor of 1.3.This indicates that the earlier assumption of the equality of both coefficients may not be valid for lakes.
The eddy-covariance (EC) technique is a micrometeorological method to directly measure momentum, heat, water vapor and greenhouse gas fluxes (Foken, 2008).It is based on the correlation between turbulent fluctuations of vertical wind speed and scalar air properties.Using this technique, one can obtain the spatial and temporal average of turbulent fluxes originating from an area called footprint and a period of meteorological stationarity (Burba & Anderson, 2010;Foken et al., 2012;Lenschow et al., 1994;Sun et al., 2006).Nowadays, the EC technique is commonly used over lakes (Blanken et al., 2000;Golub et al., 2021;Lee et al., 2014;Mammarella et al., 2015;Nordbo et al., 2011;Spank et al., 2020;Vesala et al., 2006).However, several studies reported difficulties in measuring the wind stress at weak winds, which resulted in large uncertainties (Kantha & Clayson, 2000).Low wind speed conditions are more relevant for lakes and specifically small lakes that are the most abundant inland water bodies (Downing et al., 2006).
In this study, we evaluate the first multiple water body estimates of bulk transfer coefficients and their dependencies on wind speed and water body characteristics using EC data measured above lakes.The analysis aimed at answering the following research questions: (a) What are the typical values for the bulk transfer coefficients and their variability among lakes and reservoirs?(b) How do the values compare with the reported transfer coefficients for oceans and other lakes?(c) Can the mechanistic approaches mentioned above describe the transfer coefficients at low wind speed regime?(d) Is there a consistent dependence of the transfer coefficients on lake characteristics, such as water depth, lake area and wind fetch?In the sections below, we examine possible answers.

Eddy-Covariance Data Set
For this analysis, most of the existing EC data measured by various researchers over lakes and reservoirs were extracted from open access databases and repositories of published papers.The fluxes that are reported in the datasets were calculated using different software (e.g., EddyPro, LI-COR, Inc, 2021; TK3, Mauder & Foken, 2015;EddyUH, Mammarella et al., 2016).In total, we obtained data for 23 lakes and 8 reservoirs located in the arctic, subarctic, temperate and subtropical zones (Figure 1, Table S1 in Supporting Information S1).The water bodies are located in different landscapes, including mountains (e.g., Lake Lunz, Austria or Lake Klöntal, Switzerland), forests (e.g., Lake Vanajavesi, Finland), and arctic landscapes.The EC mast at each lake or reservoir was installed either on a floating or bottom-fixed platform, on shore, or on small islands.The measurement height ranged between 1.3 and 16.1 m with 2 m being the most frequent height among all data sets.Elongated shapes of the lakes or shore/island locations were the subject of wind direction filtering to ensure that the measured surface fluxes were originating from water.Approximately half of the water bodies in this study had a surface area (A s [km 2 ]) smaller than 10 km 2 with an average wind fetch (F ave [m]) ranging from 160 to 1,550 m.The fetch grid was estimated from the map as the distance from the measurement location to the shore with the corresponding wind direction.Then, the time series of the fetch was interpolated from this grid using the measured wind directions.The average fetch was calculated as the mean distance for the filtered wind directions.The rest of the lakes and reservoirs were larger: the maximum surface area of 2.6•10 4 km 2 and the maximum mean fetch of 2.6•10 4 m refer to one of the North-American Great Lakes-Lake Erie in the USA.The maximum depth (D max [m]) varied between 1.3 m (Lake Villasjön, Sweden) and 89 m (Rappbode Reservoir, Germany).Each EC data set contained the estimated variables averaged over 30 min intervals.).The skin temperature was observed with an infrared thermometer or calculated from outgoing longwave radiation, both corrected for the reflectance of incoming longwave radiation.Some lakes or reservoirs had only momentum flux data, resulting in fewer estimates of heat and water vapor transfer coefficients.Precipitation (parameter which was considered as a factor for filtering the data) was not available for all data sets.The duration of the EC measurements ranged from 11 days (Lake Wohlen, Switzerland) to 2,243 days (or ∼6.1 years, Lake Dagow, Germany) with a median duration of 155 days.

Data Filtering and Averaging
The individual data sets used in the analysis were subject to filtering with the following different criteria.
1. filtering based on stationarity and integral turbulence test quality flags; 2. restriction of the wind directions to ensure >90% of footprint was originated from water; 3. removing periods with ice cover; 4. removing periods with precipitation (if data on precipitation was available); 5. removing periods when the difference between water surface temperature (surface water specific humidity) and air temperature (air specific humidity) is less than 0.2°C (1.5•10 −3 kg kg −1 ); 6. removing periods with floating vegetation on the water surface (only for Lake Suwa, Japan).
Quality screening of EC data is known to be site-and instrument-specific (Burba & Anderson, 2010).The data were either available in filtered form, or they contained the quality flags provided by the software.Non-filtered data sets included quality flags for each flux value (momentum, sensible and latent heat fluxes) to ensure the stationarity of the time series (homogeneity of the flow) and developed turbulent conditions (Foken et al., 2004;Foken & Wichura, 1996).
Removing wind directions was site-specific and we carefully studied each individual site.We only accepted the data from periods when wind was blowing from the lake with sufficient fetch.We specified the accepted wind directions for each site in Table S1 in Supporting Information S1.We focused on open-water conditions and discarded ice-covered periods either using the water temperature time series or interval camera data.For Lake Suwa we removed the approximate periods when floating vegetation appeared on the water surface using interval camera data, however, for other sites this kind of data was not available.For some sites, all erroneous data due to rain interference and site maintenance were filtered by data providers, or we removed periods with precipitation (if data were available).Based on the standard accuracy of temperature and humidity sensors, we  2020).(b) The maximum depth versus surface area for each lake or reservoir.The circles in blue and red color show lakes with average wind fetch less than or more than 1,000 m, respectively, with the lake name next to it.The diameter of the circle represents the relative average fetch: the larger diameter, the larger the average fetch is.

10.1029/2022JD037219
6 of 20 also applied a minimum threshold for temperature and specific humidity difference between water surface and air.We describe the effect of these filters on the data in Text S1 in Supporting Information S1.Note, that in SI and in the sections below all transfer coefficients are shown as median values following suggestions from former studies (DeCosmo et al., 1996;Fairall et al., 2003).

Transfer Coefficients
Turbulent fluxes of momentum (τ), sensible heat (H) and latent heat (L v E) at the water surface are expressed as: where   ′ ,  ′ ,  ′ are the horizontal and vertical wind velocity fluctuations, respectively; are temperature and specific humidity scales, respectively.The standard sign convention is that the momentum flux is defined as positive downward, while sensible and latent heat fluxes as positive upward (Kaimal & Finnigan, 1994).Equations 1a-1c are not the only way to describe the transfer coefficients.The widely used COARE 3.0 bulk algorithm (Fairall et al., 2003) includes scalar-averaged wind speed instead of  10 (vector-averaged wind speed) in Equations 1a-1c, which includes the gustiness (will be discussed below in Section 2.4.4).
Using measured flux data from the obtained EC datasets, the transfer coefficients can be derived from Equations 1a-1c as follows: Wind speed, air temperature (   ) and specific humidity (   ) measured at a certain height   were converted to a standard height of 10 m considering stability of the atmosphere following the equations: where ) -the same for sensible and latent heat (Businger et al., 1971).In the literature,  ∕ is usually denoted as the non-dimensional stability parameter   .To remove the effect of atmospheric stability on the magnitude of the transfer coefficients,  ,  ,  are GUSEVA ET AL.

𝐴𝐴
DN, HN, EN were estimated for 31, 25, 23 water bodies under study, respectively, depending on the flux data availability (see details about each lake or reservoir in Table S1 in Supporting Information S1 and in data repository).
In the scientific community, there has been an ongoing discussion on the form of the transfer coefficients to be presented.For example, some studies focused only on neutral values of the drag coefficient (Li et al., 2016) or some considered the drag coefficient non-adjusted to their neutral counterpart (   ).Other studies addressed the so-called "effective" drag coefficient, which was derived as the slope coefficient for the linear relationship between   2 * and   2 10 (Xiao et al., 2013).We examine the difference between  ,  ,  and  DN, HN, EN in Section 3.1.

Smooth Flow
Previous studies focused on the parameterizations of surface roughness length  0 (see Equation 4a) to assess wind speed dependence of the drag coefficient (e.g., Ataktürk & Katsaros, 1999).In our study, we compared  DN estimated from measured momentum fluxes with the existing approaches.One of the approaches is based on the smooth flow regime at low wind speed (<3 m s −1 ), where the thickness of the viscous sublayer (   ) determines the aerodynamic roughness of the interface (Schlichting, 1968), and not the physical roughness of the water surface: where  = 0.11 [-] and   = 1.6 × 10 −5 [m 2 s −1 ] is kinematic viscosity of air. 0 can be derived from Equation 4a as:

Capillary Waves
As an alternative method to parameterize  DN at low wind speeds, we considered the approach proposed by Wu (1994).He suggested that the wind shear stress in the absence of large gravity waves is related to the ripples (capillary waves).For the capillary waves, the roughness length is related to surface tension (σ) as: where α Wu = 0.18 [-] is an empirical constant and   [kg m −3 ] is water density.Surface tension at a temperature of 20°C is σ = 7.28  ⋅ 10 −2 N m −1 .

Charnock Relationship
With increasing wind speed, the thickness of the viscous sublayer becomes smaller, and the aerodynamic roughness of the water surface (  0 ) becomes minimal, before surface gravity waves evolve.At wind speeds exceeding GUSEVA ET AL.

10.1029/2022JD037219
8 of 20 3 m s −1 , waves protrude from the viscous sublayer and surface roughness length increases with increasing wind speed, indicating the transition from a smooth to a rough flow regime.Charnock (1955) proposed the following equation for surface roughness length over fully developed surface waves, which account for typical oceanic conditions: where   ranges from 0.011 to 0.0185 [-] (Garratt, 1994),   [m s −2 ] is the gravitational acceleration.

The Concept of Gustiness
Under strong convective conditions, the wind stress at the water surface is governed by random convective motions-gusts-in the convective boundary layer (CBL), whereas the mean wind speed vector can even become zero (Godfrey & Beljaars, 1991).These large convective eddies embrace the entire CBL and affect the turbulence regime in the atmospheric surface layer.Grachev et al. (1998) formulated an approach to estimate the drag coefficient using this concept.In their study, the gustiness could explain the apparent increase of the drag coefficient estimated using the traditional equation (Equations 2a, 4a) at low wind speeds.The estimated drag coefficient accounting for gusts was a factor of 1.5-6 smaller at wind speeds below 2 m s −1 in comparison to the drag coefficient calculated from Equations 2a, 4a.The gustiness concept is widely accepted and used in the COARE algorithm to estimate air-sea fluxes (Fairall et al., 2003).
The gustiness factor,   corresponds to the ratio of the scalar-averaged (  Ũ10 , for the definition see (Grachev et al., 1998) to vector-averaged wind speed: Taking into account that ) Since the scalar-averaged wind speed is not a standard output parameter in the EC-software, the estimation of   from velocity standard deviations allowed calculation of the transfer coefficients for more lake data sets.Alternatively,   can be parameterized in terms of the convective velocity scale   * (Fairall et al., 2003;Grachev et al., 1998).We denote the gustiness factors derived from measured wind speed and from CBL scaling as  wind and  conv , respectively: where   = 1.2 [-] is an empirical constant (Beljaars, 1995) and   * is expressed as: where,   [K] is the virtual temperature,  CBL [m] is the CBL height, defined as the height of the lowest inversion.Previous studies used the fixed height of the CBL equal to 600 m (Fairall et al., 2003) or 1,000 m (Beljaars, 1995).The neutral gustiness drag coefficient  DNG and its gustiness counterpart  CDG is: and other gustiness transfer coefficients,  HNG and  ENG , can be derived as: 10.1029/2022JD037219 9 of 20 where ) 2 and    * is the scalar-averaged friction velocity.The majority of the available EC datasets from lakes did not contain the scalar-averaged friction velocity.Akylas et al. ( 2003) investigated the combinations with different averaging procedures and suggested that vector-averaged friction velocity   * is more appropriate to use with scalar-averaged wind speed for all wind speed classes.Thus, we used the vector-averaged friction velocity and therefore In a first step, we applied the gustiness approach for the drag coefficient with the cases corresponding to unstable atmospheric condition (    0 ).We estimated two types of the gustiness drag coefficient, using as  wind and conv and we considered the former one as a reference.These calculations were possible only for a subset of the data sets (11 lakes) for which both scalar-averaged wind speed (or the standard deviations of the wind speed components) and virtual sensible heat flux (   ′  ′  ) were available.Then we fitted the coefficient   in Equation 11 in order to get best agreement between the parametrized  DNG and the referenced one.As the second step, we applied the parametrization of gustiness with the fitted   (Equation 11) to larger subset of the data (26 lakes) for which the virtual sensible heat flux was available.Thus, we estimated all gustiness transfer coefficients (Equations 13a-13c).
We compared the gustiness transfer coefficients with the standard formulations (Equations 4a-4c) and with the results from the COARE 3.0 bulk algorithm (Fairall et al., 2003), which includes the parameterization of gustiness (Equations 11 and 12).In addition, the COARE algorithm considers the effect of changing water surface roughness by using parameterizations of  0 that combines the smooth flow approach (Equation 5) and Charnock's relationship (Equation 8).

Transfer Coefficients Over Lakes
Bulk transfer coefficients for neutral atmospheric stability  DN ,  HN and  EN (Equations 4a-4c) were estimated using data from 23 lakes and 8 reservoirs (see data availability details in Table in the data repository).The transfer coefficients varied between the water bodies and differed on average by a factor of 2-3 for wind speeds exceeding 3 m s −1 .However, we identified three water bodies for which the estimated drag coefficients (  DN ) were exceptionally large at all wind speeds (up to a factor of five, Lake Quinghai, China, Nam Theun 2 Reservoir, Laos), or exceptionally low (factor of four, Bol'shoi Vilyui Lake, Russia), when compared to other water bodies with similar surface area.These three water bodies contributed largely to the variability among systems (Figure S5a in Supporting Information S1 shows the estimates for individual water bodies).We did not find possible sources of errors and considered these data as outliers.In the overall estimates and in the range of variability shown in Figure 2a, we included the complete data set.
All transfer coefficients showed a similar wind speed dependence (Figures 2a-2c).At high wind speeds (>3 m s −1 ),  DN ,  HN, EN had relatively constant values of 1.8•10 −3 , 1.4•10 −3 , 10•10 −3 , respectively.All transfer coefficients increased towards the lowest wind speeds.The strongest increase was found for  DN , which was one order of magnitude higher (1.1•10 −2 ) at the lowest wind speed (0.5 m s −1 , the first bin) compared to values at higher wind speeds.A similar, but less pronounced increase was observed for  HN and  EN : their values at the lowest wind speed were 3.1•10 −3 and 2.2•10 −3 , respectively.

Unstable atmospheric conditions (
0 ) prevailed over all water bodies, particular during the evening and at night time, when >80% of all data were obtained under unstable conditions (Figure 2d).Stable atmospheric conditions occurred most frequent during the day (12-19 hr).In addition, we estimated the percentage of time when the wind speed was less than 3 m s −1 (Figure 2e).Low wind speed conditions prevailed slightly during the evening and at night, when the atmosphere was mostly unstable.This means that the significant increase of the transfer coefficients at low wind speeds frequently coincides with unstable atmospheric conditions, when the water is still warm and the atmosphere starts cooling at the end of the day.
To analyze the effect of atmospheric stability on the transfer coefficients, we compared the transfer coefficients

10.1029/2022JD037219
10 of 20 it is evident that at low wind speeds (0-2 m s −1 ) these transfer coefficients under in-situ conditions were systematically higher (up to a factor of 2-3) than their neutral counterparts  DN ,  HN and  EN .
Estimation of   and   (Equations 2b, 2c and Equations 4b, 4c) involves water surface temperature, for which the skin temperature is the most appropriate measure.However, these measurements were not available for some sites.Instead, we used water temperature measured at some depth (often varying between 0 and 0.5 m between datasets).We compared three types of calculations of   using two subsets which use: (a) only skin temperature (b) only water temperature and (c) the combined lake data set which includes either skin or water temperature (Figure S7 in Supporting Information S1).
HN estimated with water temperature tends to be slightly lower than the estimates using skin temperature (the mean percentage difference is approximately 15%).As a result, we presented   and   (Figures 2b and 2c) calculated using all available data, independent of how water surface temperature measured.When both the skin and water temperatures were available for one site, the skin temperature was used.

Parametrizations of the Drag Coefficient
As described in Section 2.4.4,we estimated drag coefficients accounting for gustiness (  DNG ) using gustiness factors derived from measured scalar-averaged wind speed (  wind ), and from the parametrization of convective velocities (  conv ).To test the applicability of the parametrization, we compared  DNG estimated using both Colored lines show the results from previous studies: LK, brown line-Lake Kasumigaura, Japan (eddy covariance, Wei et al., 2016); LN, red line-Lake Neuchâtel, Switzerland (dissipation method, Simon, 1997); LG, pink line-nearshore site at Lake Geneva, Switzerland (wind profile method, Graf et al., 1984); RG, light orange color-Reservoir Gorkiy, Russia (wind profile method, Kuznetsova et al., 2016) approaches during unstable atmospheric conditions for a subset containing 11 lake data sets with all required data (Figure S8 in Supporting Information S1).The estimated  DNG based on  conv was slightly higher than the one based on  wind .We obtained a coefficient β = 1.4 (compared to a value of 1.2 suggested by Beljaars (1995)) for the best agreement between both drag coefficients (minimum of the mean absolute value of their difference).Note, that for this subset of lake data set the gustiness drag coefficient was on a factor of two lower than  DN (Figure S8 in Supporting Information S1).
Using the fitted parameterization, we calculated  DNG ,  HNG and  ENG based on  conv for a larger subset of the data-26 and 23 lakes (Figure 3). DNG and  ENG decreased approximately by a factor of two at wind speeds less than 1.5 m s −1 in comparison to the standard formulation  DN and  EN . ENG was almost constant with a value of 1.0•10 −3 throughout all ranges of wind speed.
decreases only (∼10%) with  HN at wind speeds.For wind speeds 3 m s −1 , all gustiness transfer coefficients almost with their standard  DNG was slightly smaller: ∼1.6•10 −3 .examined the possible mechanisms (Equations 5 and 7, Section 2.4) that could explain the increase of the drag coefficient at low wind speed and the Charnock relationship (Equation 8), which wind speed dependence at high wind speed.We applied the COARE 3.0 algorithm to the same subset of the data (26 and 23 lakes) and compared the results with our estimates of the gustiness transfer coefficients.It was evident that the smooth flow approach (together with Charnock relationship) could not explain the increase of  DNG at low wind speeds.However, replacing this model with the approach which considers capillary waves (Equation 7with  Wu = 0.8) led to reasonable agreement between the bulk parameterization and  DNG estimated from measured fluxes.In contrast,  ENG calculated using this approach overestimated values of the EC-derived  ENG .There, smooth flow was a more appropriate parameterization.
When analyzing the standard formulation of the drag coefficient, we found that the function describing the wind speed dependence of  DN proposed by Liu et al. (2020) based on EC measurements over terrestrial surfaces ) , can also be applied to observations over all lakes (Figure 3a).In a similar way, we applied this empirically derived function to   and   estimates (Figure S7).The fitted coefficients for our data are provided in Table 1.proposed by Liu et al. (2020) with fitted coefficients (Table 1).The green and blue lines show transfer coefficients as they are used in the COARE 3.0 algorithm for bulk parameterizations.The green line shows the original COARE parameterization, which combines the effects of thesmooth flow approach (Equation 5) and Charnock relationship (Equation 8).The blue line shows a modified parameterization, which includescapillary wave roughness (Equation 7with  Wu = 0.8) and Charnock relationship.
Note.The coefficients were obtained from least-square fits of the Function to the bin-averaged data shown in Figures 3a-3c.

10.1029/2022JD037219
12 of 20 The gustiness approach, together with increasing surface roughness due to capillary waves allowed to explain the increase of the transfer coefficients at low wind speeds, suggesting that these formulations of the transfer coefficients can provide most accurate flux parameterizations at lake surfaces.
The mean ratio of  HNG to  ENG is 1.3 and has its maximum value of 2 low wind speeds (Figure 4).wind speeds from 2 to 9 m s −1 the ratio at relatively constant of 1.2.

of the Bulk Coefficients on the Lake Characteristics
We examined the dependencies of the bulk transfer coefficients accounting for gustiness on lake characteristics, including the maximum and average water depth, water depth at the measurement site, maximum and average wind fetch, and water surface area.As the transfer coefficients at high wind speeds were relatively constant, we first analyzed effects of lake characteristics on the median values of the transfer coefficients for wind speeds exceeding 3 m s −1 estimated for each individual water body.
We found that  DNG and  HNG decreased significantly (Pearson correlation coefficient r = −0.5, p-value <0.05) with increasing lake surface area (Figures 5a and 5d).These relationships could be expressed as power law dependencies (  =   exp( ln 10) , where   and   are the slope and intercept of the linear regression  log 10  =  log 10  +  ) with exponent of −0.06.Most variability in  DNG was found to be explained by the lake surface area (for log-transformed data the coefficient of determination was R 2 = 0.3).
The correlation between  DNG ,  HNG and  ENG and mean or maximum fetch was low (r ∼ −0.2, −0.3, Figures 5b, 5e and 5h, Figures S9a, S9d, S9g in Supporting Information S1).A principal component analysis revealed that lake surface area has a largest predictive power (Figure S10 in Supporting Information S1).We did not find a significant correlation (r ∼ −0.3, p-value >0.05) between  ENG and surface area, however, a similar trend as for  DNG and  HNG could be observed.
Using the principal component analysis, we identified that there was no significant correlation of the transfer coefficients at high wind speeds with maximum, average or local water depth (Figures 5c,5f and 5i,Figures S9 and S10 in Supporting Information S1).We used the exponential dependence from Panin et al. (2006) to compare with our results.However, we did not have sufficient sites with larger depth to confirm the depth dependence reported in their study.
At low wind speeds (<3 m s −1 ), the transfer coefficients were strongly wind speed dependent (Figures 2a-2c) and their relationships with lake characteristics are examined separately for each different wind speed interval.Here we found that  DNG significantly increased with increasing water surface area for wind speeds between 0.5 m s −1 and 2 m s −1 .At higher wind speeds these correlations become negative, as in the analysis for wind speed >3 m s −1 presented above.As an example, we show the transfer coefficients for a wind speed of 1 m s −1 in Figure S11 in Supporting Information S1.At the same time, we found significant correlation of  DNG with measurement height at low wind speeds, which was not present at high wind speeds (Figure S11d in Supporting Information S1).No correlation with measurement height was found for  HNG and  ENG .
Additionally, we looked for a possible relationship between the averaged wind speed (estimated over entire time series for each individual water body) and surface area.We found a significant correlation between them in a double-logarithmic domain (r = 0.5, p-value <0.05, Figure S12 in Supporting Information S1), resulting in increasing mean wind speed with increasing lake size following a power-law dependence  10 =  0.05  exp(0.5 ln 10) .

Bulk Transfer Coefficients Estimated for Lakes and Reservoirs
We examined the bulk transfer coefficients describing the transport of momentum, heat and water vapor at the water surface estimated based on EC data collected at 23 lakes and 8 reservoirs of different size, depth, and location.At first, we used their standard formulations commonly used in literature (Equations 4a-4c).All GUSEVA ET AL.

10.1029/2022JD037219
13 of 20 transfer coefficients tended to increase toward low wind speeds and remained relatively constant at wind speeds exceeding 3 m s −1 .This increase reported in previous studies for lakes (see, e.g., Wei et al., 2016;Xiao et al., 2013) and for the land surface (Grachev et al., 2011;Liu et al., 2020) and has been extensively investigated, yet remained unexplained up to now.Authors of Grachev et al. (2011) referred to lakes as an "intermediate" GUSEVA ET AL. 10.1029/2022JD037219 14 of 20 case between over-sea and over-land locations in terms of turbulent transfer.The lower bound for  DN ,  HN, EN among the water bodies at high wind speeds were within the range reported by previous studies, including for large lakes (>200 km 2 (Kuznetsova et al., 2016;Wei et al., 2016), classical open ocean measurements (Fairall et al., 2003;Large & Pond, 1981) and coastal sites under fetch-limited conditions (Lin et al., 2002).Indeed, we also considered large lakes (Figure 1b) that were expected to have the smallest drag coefficient as they had the largest fetch (e.g., Lake Erie, Lake Taihu, Lake Balaton).The mean  DN for winds exceeding 3 m s −1 was equal to 1.8•10 −3 and this value corresponded to an upper bound for the water surface roughness (0.001 m) reported by Foken (2008), but was a factor of two higher than the values reported for oceans and large lakes or reservoirs (Fairall et al., 2003;Large & Pond, 1981).

Values of
DN varied considerably depending on the type of measurements used for its estimation.For example, in Simon (1997),  DN was calculated from the dissipation rate of turbulent kinetic energy measured at the water side of the air-water interface in the relatively large Lake Neuchâtel (218 km 2 , Switzerland).
DN was significantly lower than our estimates (factor of ten) and the estimates from lakes or marine measurements (factor of five).However, these estimates also confirmed the increase of  DN at low wind speeds. DN at high wind speeds calculated from the wind profile method at the nearshore site in Lake Geneva (Graf et al., 1984) was in close agreement with our estimates.The strong increase of  DN ,  HN, EN at low wind speeds was similar to the one observed for a large lake using surface fluxes measured by EC method (Wei et al., 2016), but it was not supported by measurements in the marine environment.

Bulk Transfer Coefficients at Low Wind Speed
Low wind speeds are typical conditions for lakes (Woolway et al., 2018), especially for smaller ones (Figure S10 in Supporting Information S1), which are most abundant by number (Downing et al., 2006).The most pronounced increase in bulk transfer coefficients at low wind speed was observed for  DN , which was up to one order of magnitude higher at low wind speeds compared to its value at high wind speeds.We found less pronounced increases of  EN and  EN , yet their values at low wind speed can be larger compared to the constant values at high wind speed by up to factor of three and two, respectively.Periods with low wind speeds mostly corresponded to periods with unstable atmospheric conditions or enhanced convective transport, which is the most prevailing condition at all studied lakes during the ice-free period.It confirms former findings of Read et al. (2012) and Woolway et al. (2017).
While Wei et al. (2016) suggested that the contribution of gusts (different formulation of the  DN ) was not significant in their measurements obtained over a large lake, we found that the strong increase can partially be attributed to missing consideration of gustiness in the definition of  DN .Correction involving the gustiness factor (  conv ) could reduce the values of  DN and  EN by up to a factor of two at wind speeds around 0.5 m s −1 .Thus, the gustiness approach is recommended to be applied in bulk parametrizations of both fluxes in lakes.
It is important to note, that the CBL scaling requires knowledge about the height of the unstable boundary layer, for which we used a fixed value of 600 m, as in other studies, including the COARE algorithm for bulk parametrization of air-sea fluxes (Fairall et al., 2003).Over lakes, convective motions can be associated with lake-land breeze circulation (Crosman & Horel, 2012), which additionally depends on lake size and on the atmospheric conditions in the surrounding landscape.Nevertheless, we found that the gustiness parameterization based on the convective velocity scale (  conv ) agreed reasonably well with observed gustiness in wind speed fluctuations (  wind ).Best agreement was obtained by using a slightly higher value of the empirical constant   (Equation 11), for which we estimated a value of 1.4 (compared   = 1.2 used in the COARE algorithm.The close agreement between parameterized and measured guestiness suggests that the parameterization can improve bulk parametrizations also under more complex CBL dynamics in lakes in comparison to the open ocean. Comparison of our measurements with the bulk parametrizations in the commonly used COARE 3.0 algorithm, which includes the gustiness correction (Fairall et al., 2003), showed a remaining underestimation of the parameterized  DNG at low wind speed.Obviously, the remaining increase in the gustiness-corrected drag coefficient at low wind could not be described by the smooth flow parametrization (Equation 5), which is applied in the COARE 3.0.Instead, we applied the capillary wave approach for parametrizing the surface roughness at low wind speed (Equation 7) and found reasonable agreement between observed and parameterized gustiness drag coefficients.For  ENG the smooth flow model performed better, as it became a nearly constant value after applying the gustiness approach.The effect of gustiness was not significant for  HNG and similar to  DNG the capillary wave parametrization performed better than smooth flow approach.
The wind speed dependence of the bulk transfer coefficients in their standard formulation (especially at low winds), could be well described by an empirical function that was originally proposed for the land surface (Liu et al., 2020).This suggests that this function can be used to describe the transfer coefficients without consideration of gustiness.
Unexpectedly, we found that the gustiness drag coefficient significantly increase with increasing lake surface areaat wind speeds less than 2 m s −1 .This result is counterintuitive, because at low winds we did not expect a dependence on lake surface area or fetch, as it should only be important for the development of surface waves, which appears only at wind speeds exceeding 3 m s −1 (Guseva et al., 2021;Simon, 1997).At the same time, we found a significant positive correlation between  DNG and the measurement height, which was also unexpected.This finding could be the result of measurement limitations and require a separate detailed investigation.

Bulk Transfer Coefficients at High Wind Speed
The estimated  DN and  DNG agreed closely at high wind speeds, thus, can be considered identically.In comparison to the results from the COARE 3.0 algorithm,  DNG at wind speeds larger than 3 m s −1 was higher than predicted by Charnock relationship (Equation 8).This result was expected as Charnock relationship is based on the assumption that the water surface roughness is controlled by fully developed surface gravity waves.When comparing with the results from the modified COARE 3.0 algorithm (capillary wave instead of smooth flow approach), we found that the consideration of Charnock relationship had a weak influence only, so that it could probably be omitted without large changes in bulk parameterizations.
We found a significant correlation between  DNG ,  HNG and lake surface area.For large lakes, the transfer coefficients at high wind speed tended to be lower.At these higher wind speeds, the surface gravity waves could potentially reach the fully developed in large water For  ENG , the correlation with lake surface area was not significant, but a similar trend could be observed.The values of  ENG our were mately a factor of two lower than previous studies (Brutsaert Yeh, 1970; Harbeck, Our could not confirm a bilinear decrease  DNG with increasing lake size with a weaker dependence for large lakes, as suggested by Woolway et al. (2017) (Figure 5a).This can potentially be attributed to the parameterization of the transfer coefficients used in Woolway et al. (2017) in lack of flux measurements.In contrast to the results reported in (Panin et al., 2006), we did not find evidence for the existence of an influence of water depth on the bulk transfer coefficients.

Study Limitations
The estimated bulk transfer coefficients show large scatter, even after filtering the data.The scatter is particularly high at light winds, that is, in the first three to four wind speed intervals (0.5-2 m s −1 ).It could be associated with limitations of the EC measurements, namely, the validity of the underlying assumptions, including the homogeneity and stationarity of the flow, as well as by increasing random errors.These effects may not have been removed completely by data filtering.Moreover, our estimates of  DNG differ from former results obtained using different types of measurements, such as water-side energy dissipation rates (Figures 3a, e g., Simon, 1997).Thus, the combination of water-and air side measurements could be beneficial for further investigation of the bulk transfer coefficients.
The parameterization of gustiness used in the COARE 3.0 algorithm assumes the CBL height as a constant value of 600 m, however, this model has been used mostly for open ocean conditions, not for lakes.Considering the fact, that the thermal internal boundary layer can develop above lakes due to temperature difference between land and water (Glazunov & Stepanenko, 2015), the exact layer structure of the atmosphere (corresponding to the cases of the high drag coefficients-low wind speed and convection) should be studied.
In our analysis we considered 30-min intervals for averaging, however, the results will likely depend on it.For example, selection of longer time interval leads to larger gustiness velocity.But since the 30-min time interval is the standard averaging interval for the EC data processing, we did not have an opportunity to test the effect of averaging time scale.Hwang (2004) suggested that the standard height of 10 m at which the transfer coefficients are reported is inappropriate for analyzing  DN and its dependence on surface roughness under wave conditions.They argue that the only relevant parameter that could serve as a reference height is the wavelength that describes the decay rate of the waves with the distance from the water surface.The adjustment of the transfer coefficients to 10 m height may not be very relevant for lakes and reservoirs and the flux measurements at two different should be considered in future These measurements would additionally provide confirmation for the existence of a constant flux layer, which is another important assumption underlying EC measurements.

Broader Implications
Bulk transfer coefficients are usually applied in numerical models for the atmospheric boundary layer, as well as in hydrodynamic models of lakes and reservoirs.Currently, the global modeling studies focusing on the lake mixing and phytoplankton blooms for climate change predictions use constant values for  DN (Grant et al., 2021;Jöhnk et al., 2008;Read et al., 2014;Woolway & Merchant, 2019), or consider  DN as a tuning parameter of the models (Stepanenko et al., 2014).Inadequate values of  DN result in biased estimates of the current velocities in lake models (Chen et al., 2020).The increase of the transfer coefficients at low wind speeds observed in our analysis can therefore lead to significant errors, as these conditions are the most prevailing conditions for lakes.We found that instead of the standard formulation of the transfer coefficients, their gustiness counterparts  DNG ,  HNG , and  ENG can improve flux parameterizations in lake models.In the absence of data for the gustiness approach, the empirical parameterizations of the wind-speed dependence of the bulk transfer coefficients provided in Table 1 can potentially be applied in modeling lake-atmosphere interactions.The observed dependence on the lake surface area is more complicated to implement, as we observed contrary dependencies at low and high wind speeds.
While the Stanton and Dalton numbers are commonly assumed to be equal, we found that  HNG was on average by a factor of 1.3 higher than  ENG (averaged over all wind speeds and all water bodies under study).The finding of  HNG being higher than  HNG confirmed the results reported by, for example, (Dias & Vissotto, 2017;Wei et al., 2016).The mean value of  HNG for high wind speeds (1.0•10 −3 ) was found to be the same as in (Kantha & Clayson, 2000), but  HNG was larger (1.4•10 −3 ) as in (Harbeck, 1962;Hicks, 1972).The fact that  HN > EN may have significant implications, because it results in biased estimates of lake evaporation based on the energy-budget (Bowen ratio method).In particular, it would result in larger (smaller) sensible heat (latent heat) fluxes than those predicted under that assumption.Both the physical mechanisms underlying their difference and the extent of the differences in the predicted sensible and latent heat fluxes require further investigation.
In state-of-the-art weather and earth system models, lakes are included as separate tiles in the model cells, where the surface fluxes over the tiles are computed via Monin-Obukhov similarity scaling.The models provide constant meteorological variables for each grid cell, which is a so-called blending height concept (von Salzen et al., 1996).To use the bulk transfer coefficients derived in this study to compute fluxes, specific values of wind, temperature and humidity over lakes should be used, which can be obtained in generalization of the tile approach, involving the parameterization of internal boundary layers over contrasting surfaces (Arola, 1999;de Vrese et al., 2016;Molod et al., 2003).MacKay (2019) presents a specific example of such an approach developed for lakes and wind speed only.

Conclusions
We were the first to analyze the bulk transfer coefficients of momentum, sensible and latent heat from directly measured surface fluxes above various lakes and reservoirs.We observed a pronounced increase of the transfer coefficients at low wind speeds (<3 m s −1 ) and relatively constant values at high wind speeds (>3 m s −1 ).The strong increase in the transfer coefficients at low wind speed was found to be associated with the existence of gusts.It is recommended to use the gustiness approach in calculation of the transfer coefficients.The reduced, yet still increasing values of the gustiness drag coefficient and Stanton number at low wind speeds can be explained by capillary wave roughness.The gustiness Dalton number, however, remains constant at all wind speeds and is better described by the smooth flow roughness.The Stanton number was systematically higher than the Dalton number by a factor of 1.3, which has implications for the Bowen ratio method.At high wind speed, the drag coefficient and the Stanton number decreased with increasing surface area of the water body and with The variables included wind speed (U z [m s −1 ]), wind direction (WD [°]), friction velocity (   * [m s −1 ]) as a quantity characterizing the momentum flux (   =  2 * [kg m −1 s −2 ],   -air density [kg m −3 ]), air temperature (T a [°C]), turbulent fluxes of sensible heat (H [W m −2 ]), and latent heat (L v E [W m −2 ]), the latter referred to in this paper as GUSEVA ET AL.
flux as well.Water temperature was provided either as skin surface temperature (T s [°C]) or bulk water temperature, measured at 0-0.5 m water depth (T w[°C]

Figure 1 .
Figure 1.(a) Geographical distribution of the eddy-covariance measurements over the lakes and reservoirs used in this study (red circles).Map was created using the software by Pawlowicz (2020).(b) The maximum depth versus surface area for each lake or reservoir.The circles in blue and red color show lakes with average wind fetch less than or more than 1,000 m, respectively, with the lake name next to it.The diameter of the circle represents the relative average fetch: the larger diameter, the larger the average fetch is.

Figure 2 .
Figure 2. Neutral (a) drag coefficient (  DN ), (b) Stanton number (heat transfer coefficient,  HN ), (c) Dalton number (water vapor transfer coefficient,  EN ) versus wind speed at 10 m height.The shaded gray area indicates the variability among water bodies by marking the range from the 5th to the 95th percentiles of the median values estimated for each wind speed bin (0.5 m s −1 ).The black line with circles, marked as L&R denoting "lakes and reservoirs," represents the median values for all lakes and reservoirs.The small inset in (a) shows the data beyond the scale.Vertical and horizontal black dashed lines mark a constant wind speed of 3 m s −1 and typical values of  DN ,  HN = EN 1.3•10 −3 , 1.1•10 −3 , respectively.Colored lines show the results from previous studies: LK, brown line-Lake Kasumigaura, Japan (eddy covariance,Wei et al., 2016); LN, red line-Lake Neuchâtel, Switzerland (dissipation method, Simon, 1997); LG, pink line-nearshore site at Lake Geneva, Switzerland (wind profile method,Graf et al., 1984); RG, light orange color-Reservoir Gorkiy, Russia (wind profile method,Kuznetsova et al., 2016); OO, dark green color-open ocean (eddy covariance, Large & Pond, 1981); OO, light green color-open ocean (eddy covariance, Fairall et al., 2003); CO, dark yellow color-coastal ocean at limited fetch conditions (eddy covariance, Lin et al., 2002).(d) Mean diel pattern of the percentage of time periods with unstable atmospheric conditions (stability parameter    0 ).(e) Mean diel pattern of the percentage of time with low wind speed (<3 m s −1 ).The red and black lines in (d) and (e) show the mean value and the shaded area shows ± standard deviation of all data.
Figure 2. Neutral (a) drag coefficient (  DN ), (b) Stanton number (heat transfer coefficient,  HN ), (c) Dalton number (water vapor transfer coefficient,  EN ) versus wind speed at 10 m height.The shaded gray area indicates the variability among water bodies by marking the range from the 5th to the 95th percentiles of the median values estimated for each wind speed bin (0.5 m s −1 ).The black line with circles, marked as L&R denoting "lakes and reservoirs," represents the median values for all lakes and reservoirs.The small inset in (a) shows the data beyond the scale.Vertical and horizontal black dashed lines mark a constant wind speed of 3 m s −1 and typical values of  DN ,  HN = EN 1.3•10 −3 , 1.1•10 −3 , respectively.Colored lines show the results from previous studies: LK, brown line-Lake Kasumigaura, Japan (eddy covariance,Wei et al., 2016); LN, red line-Lake Neuchâtel, Switzerland (dissipation method, Simon, 1997); LG, pink line-nearshore site at Lake Geneva, Switzerland (wind profile method,Graf et al., 1984); RG, light orange color-Reservoir Gorkiy, Russia (wind profile method,Kuznetsova et al., 2016); OO, dark green color-open ocean (eddy covariance, Large & Pond, 1981); OO, light green color-open ocean (eddy covariance, Fairall et al., 2003); CO, dark yellow color-coastal ocean at limited fetch conditions (eddy covariance, Lin et al., 2002).(d) Mean diel pattern of the percentage of time periods with unstable atmospheric conditions (stability parameter    0 ).(e) Mean diel pattern of the percentage of time with low wind speed (<3 m s −1 ).The red and black lines in (d) and (e) show the mean value and the shaded area shows ± standard deviation of all data.

Figure 3 .
Figure 3. Neutral gustiness transfer coefficients (a)  DNG , (b)  HNG and (c)  ENG (calculated using  conv , see Section 2.4.4) versus  10 (red line) estimated for 26, 23 and 23 lakes, respectively.The red shaded area marks the range between the 5th and 95th percentiles.The black line with symbols shows estimates of  DN ,  HN and  EN from Figure 2. The dark yellow line in all plots shows the function  DN = 1 [ 1 + 2 exp(310) ] proposed by Liu et al. (2020) with fitted coefficients (Table1).The green and blue lines show transfer coefficients as they are used in the COARE 3.0 algorithm for bulk parameterizations.The green line shows the original COARE parameterization, which combines the effects of thesmooth flow approach (Equation5) and Charnock relationship (Equation8).The blue line shows a modified parameterization, which includescapillary wave roughness (Equation7with  Wu = 0.8) and Charnock relationship.

Figure 4 .
Figure 4. Ratio of  HNG to  ENG estimated for each individual data set (21 water bodies, shown by gray lines).Black line with circles shows the median values for all data sets.The horizontal dashed black line shows a ratio of 1:1 and the horizontal dashed red line indicates the overall mean value of the ratio (1.3).

Figure 5 .
Figure 5. Neutral transfer coefficients accounting for gustiness (a, b, c)  DNG ; (d, e, f)  HNG ; (g, h, i)  versus surface area of the water body, mean fetch length, and water depth at the measurement site.Panels show the median values of the transfer coefficients for wind speeds exceeding 3 m s −1 for each individual water body.The red lines show the linear regressions of log-transformed transfer coefficients.The corresponding power laws as well as the Pearson correlation coefficient (r) and p-value are provided in the upper left corner of each panel.Significant correlation is marked by bold labels.The single red symbol in (a) marks the data from Nam Theun 2 Reservoir, which was not considered in the linear regression analysis and in the Pearson correlation as an obvious outlier.Blue, green and dark yellow lines show results from previous studies by Woolway et al. (2017), Panin et al. (2006), Harbeck (1962), respectively.
are the covariances of vertical wind velocity and air temperature (   ′ ) and specific humidity (   ′ ) fluctuations. 10 [m s −1 ] is wind speed at 10 m height,   and  10 [K] are the surface water temperature and the air temperature at 10 m height, respectively,   and  10 [kg kg −1 ] are the specific humidity at the air-water interface (estimated from surface temperature) and at 10 m height, respectively.
[J kg −1 K −1 ] is the specific heat of air at constant pressure, and