Spatial variability of global lake evaporation regulated by vertical vapor pressure difference

Evaporation (E) from about 300 million lakes worldwide without plant physiological constraints directly reflects hydrological response to atmospheric forcings. However, it remains inadequately understood about what regulate spatial variability of global lake E across seasons. Here we show that vertical vapor pressure difference ( eD ) accounts for 66% of the spatial variability of annual E, followed by wind speed (16%). The eD is also the predominant factor modulating diurnal variability in E and causing greater E at night than during the daytime. As a consequence, spatial variability in nighttime E strongly regulates that in global E across seasons. Therefore, the observed widespread, heterogeneous changes in lake surface temperature that imply spatial variability in eD may have contributed to changes in global E variability.


Introduction
There are more than 300 million lakes with surface areas of greater than 0.001 km 2 in the world which cover about 3% of the terrestrial land (Downing et al 2006). In the climate system, lake evaporation (E) is considered energy limited (Brutsaert and Parlange 1998) and its direct response to climate perturbations independent of varying terrestrial plant ecosystem functionalities and regulations make lakes ideal systems to attribute changes in evaporative processes to atmospheric forcings. However, spatiotemporal variations in global lake E in response to environmental controls remain less studied.
Previous field studies in some individual lakes have indicated the varying effects of atmospheric variables and lake-surface processes on lake E at different time scales (Woolway et al 2020). At diurnal scale, vapor pressure difference (e D ) between lake surface and overlying air, wind-induced turbulence, and atmospheric stability conditions that are largely associated with temperature difference (TD) between the lake surface water temperature (LSWT) and the overlying air primarily affect E variability (Blanken et al 2000, Liu et al 2012. Lake E at sub-seasonal to seasonal time scales is largely regulated by invasions of synoptic weather events and energy inputs (Hostetler andBartlein 1990, Blanken et al 2000). Invasions of synoptic weather events disturb local or mesoscale atmospheric forcings that alter LSWT, leading to changes in water surface energy balance . Lake water bodies store large amount of heat during lake heating periods (e.g. daytime with solar heating at daily time scales and spring and summer at seasonal time scales) and release stored heat to the overlying air during lake cooling periods (e.g. nighttime and autumn and winter) due to the large heat capacity of water , Lenters et al 2005. The heat storage of lake water bodies lessens temporal changes in LSWT that are off phase with temporal changes in the over-lake air temperature (T a ) (Henderson-Sellers 1986, Bonan 1995. Since the saturation vapor pressure (e s ) across the water-air interface is a function of LSWT, such variations in LSWT strongly affect e D variability, leading to changes in lake E across different time scales (Hostetler andBartlein 1990, Woolway et al 2020). Climate warming has caused widespread, but geographically heterogeneous increase in LSWT across the globe (O'Reilly et al 2015). Lakes tend to channel more energy towards latent heat flux (LE) than sensible heat flux (H) in disproportionate rates under warmer climate conditions (i.e. decrease in Bowen ratio) (Wang et al 2018). Thus, the heterogeneous lake warming implies that spatial variability in global lake E might have been changed through various feedback mechanisms.
Lake E and its environmental controls are also indirectly regulated by lake attributes such as geographic locations and depths (Sene et al 1991, Delclaux et al 2007, Rosenberry et al 2007. Lake geographic locations dictate seasonal variations of available energy inputs (i.e. solar radiation). The lake-depth dependent water heat storage controls temporal phase difference between energy inputs and E (Granger and Hedstrom 2011) and therefore plays a critical role in modulating E patterns. In tropical lakes, both shallow and deep, the mean monthly E is generally in-phase with net radiation (R n ), while diurnal variations in E are sensitive to wind speed (Sene et al 1991, Delclaux et al 2007. For mid-and high-latitude shallow lakes, there is a strong correlation between mean monthly E and R n . In mid-and high-latitude deep lakes, however, a phase difference of one to three months (i.e. time delays) affected by lake depths can occur between E and R n , Oswald and Rouse 2004, Li et al 2016, and diurnal variations in E are best explained as a function of wind speeds and e D , Shao et al 2015. While these previous field studies demonstrate varying effects of specific factors among environmental controls and lake attributes on lake E for specific lakes, spatial variability in global lake E as a result of the collective effects from these factors remains less studied.
It is challenging to measure lake E (Lenters et al 2005). Several indirect methods have been used for deriving lake E based on readily measured variables, including pan coefficient (Kohler et al 1955), water budget (Sturrock 1978), Bowen ratio energy balance (Ficke 1972, Sturrock 1978, Penman model (Penman 1948), Priestly and Taylor model (Priestley and Taylor 1972), and Thornthwaite model (Thornthwaite 1948) among many others. More recently, the in-situ lake E measurements using the eddy covariance technique have become increasingly available (Rouse et al 2005, Nordbo et al 2011. Nevertheless, these eddy covariance measurements of lake E remain sparse. To advance understanding of spatiotemporal variations in E and the associated environmental controls, numerical lake models have emerged as a useful approach to simulate global lake E (Wang et al 2018).
This study reports spatial variability of global lake E across seasons and quantifies the impacts of atmospheric forcings on these variations. Using Lake, Ice, Snow, and Sedimentation Simulator (LISSS) to simulate global lake E, we aim to answer the following questions: (a) what is the spatial distribution of global lake E and how available energy is partitioned among surface energy fluxes across seasons; (b) what are the key environmental parameters regulating lake E variability; and (c) how the spatial variability in day and night E regulates global lake E variability and what are their driving parameters.

Model simulation and validation
We use the LISSS in the Community Land Model (CLM v. 4.5) to simulate LSWT and lake surface energy fluxes including LE or E, H, and heat flux into or from lake water (G). For ice/snow covered lakes during cold moths, the ice-surface temperature is used as LSWT. The LISSS is a sub-module within the CLM that uses a nested hierarchy consisting of five land units or tiles (glacier, urban, agricultural, vegetation, and lake) to represent surface heterogeneity at the subgrid level for each CLM grid. In this study, the CLM simulations are conducted at 0.94 • × 1.25 • horizontal resolutions and 30 min temporal resolution over 1930-2016 driven by CRUNECP atmospheric forcing dataset with the first 50 years as the spin-up period. The CLM unifies all lakes within a grid cell by using the mean depth and runs LISSS at subgrid level as one tile to calculate surface energy fluxes. The brief descriptions of LISSS are given in text S1 (available online at stacks.iop.org/ERL/17/054006/mmedia). The outputs from 4591 lake tiles are extracted and then analyzed in this study. The good performance of the LISSS in reproducing fluxes has been demonstrated by several studies (Subin et al 2012a, Hu et al 2017, Wang et al 2018. Our own validation also indicates its good performance (figure S1). The simulated lake E and LSWT are validated against observations (figure S1). The statistical analysis indicates that simulated surface fluxes are significant with r 2 = 0.89, RMSE = 12.02 W m −2 , p < 0.01 for E and r 2 = 0.96, RMSE = 1.8 • C, p < 0.01 for lake surface temperature. Despite of the model good performance in simulating surface fluxes, the LISSS tends to have a positive bias of about 0.69 • C in the simulated LSWT. The implications of this bias for lake E are discussed in section 4.2.

Methods for data analysis
One-year fluxes at the 30 min interval by averaging the fluxes over 2000-2016 are used to examine their spatial variations across seasons. The simulated fluxes at the 30 min resolution were then aggregated to daily and seasonal means. We divided the global lakes in three climate regions and calculated arithmetic means of all values of individual lakes within a climate region (table 1). The time periods with incoming shortwave radiation greater or less than 25 W m −2 are defined as daytime or nighttime, respectively. The unit of mm d −1 is used for both the daytime and nighttime E despite of the different time durations for daytime and nighttime so that their evaporation rates can be directly compared. Fluxes are averaged over single lakes and then over the latitude before the mean zonal ratios of fluxes are calculated. It is important to note that the average fluxes over the latitude are not weighted by lake area and therefore are independent of lake size. The seasons are synchronized for both the hemispheres. For example, spring is defined as a combination of March-May in the Northern Hemisphere and September-November in the Southern Hemisphere. Here and throughout, H and LE from lake surface to the atmosphere are defined as positive, and the downward radiation flux components and G from the water surface to the layers below are defined as positive.
Multiple linear regression (MLR) is performed to analyze the influence of e D , wind speed, T a , specific humidity (q a ), and R n on E. All variables are standardized before analysis is made by removing mean and scaling to unit variance so that the results can be interpreted at comparable scales. The relative importance (i.e. variance explained, β) of each environmental variable as E driver is quantified by a ratio between the regression coefficient of individual variable and the sum of regression coefficient of all variables from the MLR. We use the arithmetic mean of β over latitudes within each climate region to explain the spatial variability of driving variables on lake E. The response of lake E to change in e D and air vapor pressure deficit (VPD) is also compared using E sensitivities to e D and VPD. Such sensitives are calculated using Theil-Sen's linear regression with a 95% confidence interval (Sen 1968). The Theil-Sens is a robust nonparametric approach which calculates all possible slopes and returns the median slope as an estimate of the linear trend. The one-year flux data at the 30 min temporal resolution are used in both the regression analysis and the sensitivity analysis.

Spatial variations of global lake E across seasons
Lake E is intrinsically influenced by the water surface energy balance and partitioning. Globally, the annual mean lake E is 2.3 mm d −1 , and the corresponding annual mean LE (64.8 W m −2 ) accounts for about 82.5% of the annual mean global R n (figure 1(a) and table 1), suggesting that majority of the available energy is partitioned into latent heat fluxes. Annually, the global lakes receive 165.5 W m −2 of incoming shortwave radiation and 300.6 W m −2 of incoming downward longwave radiation (figure S2). The lakes reflect 21.9 W m −2 of the incoming shortwave radiation and absorb remaining 444.3 W m −2 of radiative flux (R abs ). In response, the lakes adjust their surface temperature and radiate longwave radiation (365.8 W m −2 , 82.3%) and fuel LE (64.8 W m −2 , 14.6%) and H fluxes (6.9 W m −2 , 1.5%) ( figure  S2). Here below the surface energy balance equation (R n = LE + H + G) is primarily used to analyze how Each row, from left to right, shows the spatial distribution of E, zonal mean E, and the ratio of zonal mean latent heat flux (LE) and net radiation (Rn), respectively. The shaded areas in zonal mean plots denote 10th and 90th percentile of E. Because of near zero/negative Rn in the latitudes of >40 • N and 40 • S, the LE/Rn ratio is very large/small and therefore has been omitted in autumn and winter zonal plots. The positive values represent flux from lake surface to the atmosphere. the surface absorbed R n is partitioned into LE, H, and G.
Global lake E shows distinct spatial variations across seasons ( figure 1 and table 1). The annual zonal mean E peaks at about 14 • N (5.4 mm d −1 ) and decreases towards the polar regions with a leveloff zone (4.1 mm d −1 ) from 30 • S to the equator and the minimum E at 75 • N ( figure 1(a)). Seasonal variations of the zonal mean E in spring and winter resemble the annual ones, but with the leveloff zone substantially reduced in winter (figures 1(b) and (e)). A second peak occurs around 25 • S in summer and autumn (figures 1(c) and (d)). At the annual and seasonal time scales, latitudinal variations in E are primarily explained by latitudinal variations in R n (figure S3; r 2 > 0.75, p < 0.05), suggesting that latitudinal variations in lake E at seasonal to annual time scales are strongly constrained by latitudinal variations in energy input. However, latitudinal variations in wind speed plays minor roles in affecting latitudinal variations in E at such time scales (figure S3; p > 0.05). LE/R n from 30 • S to 30 • N is about 89.3%, and LE/R n for the latitudes of >30 • N is about 69.3% (figure 1(a) and table 1). The lower LE/R n at latitudes of >30 • N is attributed to the fact that more R n is channeled towards snow/ice melting, as compared with low-latitude lakes. The zonal mean H/R n varies slightly across the latitudes from about 30 • S to 30 • N, but it decreases with increasing latitudes at the latitudes of >30 • N (figure 2). In spring, the negative H and R n at the northern high latitudes of >30 • N (figure 2(b)) are primarily utilized for melting snow/ice, leading to low E ( figure 1(b)). In summer, R n absorbed by high-latitude lakes is partitioned into positive H and LE and heat stored in lakes (figures 1(c) and 2(c)). The negative R n at latitudes of >30 • N in autumn is attributed to the hysteresis of lake temperatures due to heat storage. Such heat storage drives LE and causes positive H in this season (figures 1(d) and 2(d)) as dry cold air masses passes over warmer lakes , Lenters et al 2005. In winter, lakes at the latitudes of >30 • N become frozen, and the negative R n is almost balanced by the negative H, whereas LE approaches zero (figures 1(e) and 2(e)). These results indicate that LE has the large spatial variability; however, it remains unclear as to which atmospheric variables predominantly regulate such variability.

Vapor pressure difference (e D ) regulating lake E spatial variability across seasons
Our statistical analysis indicates that spatial variability in the annual mean E is mainly explained by that in e D (66%) followed by wind speed (16%) ( figure  S4(a)). The remaining 18% of the E variability is explained by a combination of other drivers including T a (figure S5; β = 0.08), q a (figure S5; β = 0.07), and R n (figure S5; β = 0.03). In the latitudes of 30 • S to 30 • N, the regression coefficients between E and e D are negatively correlated with regression coefficients between E and wind speed (figure S6). This negative correlation indicates that e D and wind speed jointly regulate the E variability in this region ( figure S4(a)). However, the relatively high but varying e D in these latitudes acts as the primary driving variable of E under persistently adequate turbulent mixing in the over-lake atmosphere with strong wind (figure 3(a)). In the latitudes of >30 • N and 30 • S, high wind speed provides strong vertical turbulent mixing, but e D is low and keeps decreasing with increasing latitudes ( figure 3(a)). As a result, the low e D acts as a strong constraint on evaporation and regulates E variability, which also explains the high regression coefficients between E and e D . At the seasonal scales, e D still outpaces wind speed in regulating spatial variability in E (figures S4(b)-(e)).

Nighttime E contributes more to spatial variability in global lake E than daytime E
Our analysis shows that global mean nighttime E is about 14% greater than the global mean daytime E (figure 4(a) and table S3). The annually averaged differences in E (△E) between nighttime (E n ) and daytime (E d ) demonstrate substantial spatial variability The left two panels show the spatial and zonal mean distribution of eD, respectively. The positive eD represents higher vapor pressure at air-water interface than air vapor pressure. The right three panels show the spatial distribution of U, the daytime and nighttime zonal mean of U, and the zonal mean of the U difference between daytime and nighttime (△U). The positive △U represents higher wind speeds at nights than during the day. from approximately −0.9-4.2 mm d -1 . Note that the unit of mm d -1 is still used for E n (E d ) which, when dividing by 24 h, results in the nighttime (daytime) evaporation rate in mm h -1 independent of the length of nighttime (daytime). The large positive annual △E primarily occurs across the latitudes from about 30 • S to 30 • N where the zonal mean E n can be 38% higher than the daytime counterpart. The zonally averaged annual △E has two peaks, one (1.4 mm) around 25 • S and the other (1.7 mm) around 14 • N, and becomes negative (i.e. higher E d ) in the latitudes of >30 • S and 30 • N ( figure 4(a)). Seasonal variations in the zonally averaged △E in spring, autumn, and winter resemble the annual △E while approaching zero in the latitudes of >30 • N in winter. In summer, however, △E peaks (2.1 mm) in about 25 • S and remains positive across latitudes ( figure 4(c)).
The spatial variability in △E is predominantly regulated by that of the difference in e D and TD between nighttime and daytime (figures S7-S9). Greater △E between nighttime and daytime corresponds to larger difference in e D and TD. The higher TD at night reflects more unstably stratified atmospheric conditions. More unstable atmospheric conditions promote evaporation thereby contributing to greater evaporation at night than during the daytime. However, the difference in wind speed between nighttime and daytime (△U) plays a minor role in regulating the spatial variability in △E (figure 3). In the latitudes of >30 • S and 30 • N, the negative △E is simply attributed to negative △e D (i.e. the lower e D at night than during the day). The positive △E is large in the latitudes of <30 • S and 30 • N across seasons, being greatest in the winter with the highest positive e D (△E/E = 22%; △e D /e D = 11%), followed by spring (△E/E = 18%; △e D /e D = 9%), autumn (△E/E = 14%; △e D /e D = 8%), and summer (△E/E = 12%; △e D /e D = 8%). Note that the global mean accumulative nighttime evaporation (AE n ) remains higher than accumulative daytime evaporation (AE d ) when day and night length is considered (figure S10). The annual, spring, autumn, and winter AE n is about 9.4%, 12.5%, 28.0%, and 43.6% higher than AE d , respectively. In summer, however, the longer days in the latitudes of >30 • N cause AE n to be 14.6% less than its daytime counterpart (figure S10(c)).

Role of VPD in modulating lake E
Several studies have reported that VPD is a key environmental variable in regulating spatiotemporal variations in terrestrial evapotranspiration (Wang et al 2019, Kimm et al 2020, while its role in the global lake E remains unclear. Although the diurnal variations in lake E are less sensitive to diurnal variations in VPD compared to e D (figure 5), it plays important role in modulating lake E via associated atmospheric dryness. The thermal lag, by which lakes retain their daytime temperature values for longer time due to higher thermal inertia of water (Rouse et al 2005, Liu et al 2009, causes e s to be less sensitive to changes in the nighttime T a (figure 5). Whereas the changes in air saturation vapor pressure are highly (e) Winter. The first and second columns, from left, show spatial distributions of sensitivities of lake evaporation E to vapor pressure deficit (VPD) and vapor pressure difference (eD), respectively. The third and fourth columns show spatial distribution of sensitivities of saturated air vapor pressure and air-water interface vapor pressure (es) to air temperature (Ta), respectively. Each data point represents the sensitivity obtained by linear regression for each lake. sensitive to changes in T a (figure 5), thereby make diurnal variations in lake E less sensitive to VPD. However, the VPD is projected to increase under future warming trends because of faster increase in saturated vapor pressure than actual vapor pressure (Sherwood andFu 2014, Ficklin andNovick 2017). The e s is also increasing due to rapid lakes warming globally (O'Reilly et al 2015, Woolway et al 2020. Higher e s and VPD will lead to elevated e D . Our analysis shows that e D mainly drives lake E across seasons, suggesting that an increase in VPD will have positive feedback on lake E trends under future warming climate.

Uncertainties in simulated lake E
The LISSS simulated lake E is statistically significant while it is acknowledged that uncertainties remain. The LISSS tends to have a mean positive bias of 0.69 • C in simulated LSWT based on our validations ( figure S1(b)). The estimated relation between LSWT and lake E (text S2 and figure S11) implies that the simulate LE could have a 5.8% bias in the annual mean global lake E. However, our test demonstrates that LISSS did not over-predict LE systematically (figure S1(a)). Uncertainties may be associated with (a) aggregating lakes at subgrid level as one tile in each CLM grid using mean depth and calculating light extinction coefficient as a function of depth (Subin et al 2012b); (b) underestimation of wind drag in case of limited fetch for shallow lakes leading to insufficient wind-driven turbulent mixing (Stepanenko et al 2013); (c) inadequate simulations of vertical circulations in deep lakes (Bennington et al 2014). Another uncertainty source is attributed to the quality of the forcing dataset. Ideally, meteorological variables that drive the model should be measured directly over lakes. Given that observation sites over lakes are sparse, all these variables are derived from reanalysis datasets that are primarily based on weather stations over land (Qian et al 2006, Sheffield et al 2006, Saha et al 2010. The use of landbased datasets to force lake models could induce additional uncertainties (Bradley 1968, McJannet et al 2012. The discrepancies among different datasets are another source of uncertainty associated with forcing data (Wang et al 2016). The comparison of terrestrial latent heat flux simulated by CLM v. 4.5 driven by CRUNECP (used in this study) and three other reanalysis products (i.e. MERRA, CFSR, and ERA-Interim) with flux tower measurements demonstrates similar r 2 (0.81-0.83) and RMSE (19.66 W m −2 -21.45 W m −2 ) (Wang et al 2016). The slight difference in RMSE is attributed to the disparities of precipitation, air temperature, and downward solar radiation (Wang et al 2016). However, the independence of lake E on precipitation by maintaining constant water mass (Oleson et al 2013) implies that simulated lake E as driven by different forcing datasets should result in small differences. Nevertheless, our LISSS-derived lake E shows variabilities comparable with measurements in individual lakes across latitudes. Thus, the simulated global E variability should capture its major spatiotemporal features.

Conclusions
Using a state-of-science lake model, our study shows that the annual, spring, summer, autumn, and winter mean global lake E is 2.3, 1.9, 3.4, 2.6, and 1.4 mm d −1 , respectively, with large spatial variability across seasons. The annual zonal mean E peaks around 14 • N where both e D and wind speed are high, and it then decreases towards both high latitudes with a leveled-off zone from the equator to 30 • S where wind speed is low. Seasonal variations of the zonal mean E show similar patterns, but with a second peak in summer and autumn in around 25 • S where both e D and wind speed are high. The e D also predominantly modulates seasonal and diurnal variations in global lake E. Further, our study indicates that the spatial variability in the nighttime E contributes more to that in the global lake E across seasons than that in the daytime E, which is largely attributed to the nighttime spatial variability in e D . The rapid increase in LSWT and VPD under warmer climate leads to elevated e D which will have positive feedback on lake E trends.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.