Greenhouse gas emissions from lakes and impoundments: Upscaling in the face of global change

Lakes and impoundments are important sources of greenhouse gases (GHG: i.e., CO2, CH4, N2O), yet global emission estimates are based on regionally biased averages and elementary upscaling. We assembled the largest global dataset to date on emission rates of all three GHGs and found they covary with lake size and trophic state. Fitted models were upscaled to estimate global emission using global lake size inventories and a remotely sensed global lake productivity distribution. Traditional upscaling approaches overestimated CO2 and N2O emission but underestimated CH4 by half. Our upscaled size‐productivity weighted estimates (1.25–2.30 Pg of CO2‐equivalents annually) are nearly 20% of global CO2 fossil fuel emission with ∼ 75% of the climate impact due to CH4. Moderate global increases in eutrophication could translate to 5–40% increases in the GHG effects in the atmosphere, adding the equivalent effect of another 13% of fossil fuel combustion or an effect equal to GHG emissions from current land use change.

Lakes and other surface waters are globally significant emitters of CO 2, CH 4 , and N 2 O to the atmosphere (Cole et al. 2007;Tranvik et al. 2009;Bastviken et al. 2011;Deemer et al. 2016;Soued et al. 2016). Accurate estimates of greenhouse gas (GHG) fluxes to and from the atmosphere are important to understanding the global carbon (C) budget as well as making projections of the trajectory of climate change. Improvements in sensor technology (e.g., Eugster et al. 2011;Gonzalez-valencia et al. 2014;Maeck et al. 2014;Delwiche et al. 2015) and the use of statistically robust survey designs (Beaulieu et al. 2016;Wik et al. 2016) have improved the accuracy of GHG emission rate measurements in aquatic systems over the last decade. These advancements and the increasing number of observations made around the world have allowed limnologists to identify environmental variables that control GHG emission rates across a broad range of spatial scales. Despite improvements in our understanding of GHG emission rates and factors influencing them, this information has not been combined to inform estimates of global GHG emissions from lakes and impoundments (hereafter referred to as "lakes").
Most traditional analyses of the global role of lakes as sources of GHGs to the atmosphere have scaled local or regional measurements to the globe by multiplying an average emission rate by the global area covered by lentic systems (Cole et al. 2007;Tranvik et al. 2009;Bastviken et al. 2011). For example, Tranvik et al. (2009) multiplied an average rate of CO 2 emission from lakes and impoundments, as determined from a literature review, by the global area of those systems (Downing et al. 2006) to calculate a global emission rate of 0.81 Pg C yr 21 . Bastviken et al. (2011) attempted a similar analysis for CH 4 emission by calculating mean emission rates and lentic area within four latitude zones. This global analysis indicated lake and impoundment CH 4 emissions of 69 Tg yr 21 of CH 4 -C or about 0.85 Pg of C as CO 2 emission equivalents (using a 34-times conversion factor for a 100-yr time-scale; Myhre et al. 2013). N 2 O is another, more potent GHG emitted from aquatic systems. Soued et al. (2016) offered a first-order approximation of the global emission of N 2 O by all inland waters using a similar traditional upscaling approach and suggested that lakes and impoundments could be emitting 0.63 Tg of N 2 O-N per year. In CO 2 equivalents, this is approximately 0.08 Pg of C (using a 298 conversion factor for a 100-yr time-scale; Myhre et al. 2013).
Estimates of global GHG emissions based on this "traditional" average emission-extrapolated approach can be strongly biased, however, if the emission rates are not derived from a representative sample of lakes. For example, estimates of CH 4 flux from northern latitude lakes (i.e., Wik et al. 2016) were substantially reduced when more recent data from large glacial lakes in Alaska were included (Stackpoole et al. 2017). Moreover, studies are increasingly finding that aquatic GHG fluxes are related to certain local drivers, which could be used to improve upscaled emission estimates. For example, Lapierre et al. (2017) found that 31% of summer CO 2 concentration variability in lakes of the United States is explained by chlorophyll a (Chl a), a proxy for system productivity, as well as color and alkalinity. Therefore, traditional regional estimates of CO 2 emissions from these systems would be biased if the average value was not derived from lakes with a distribution of those variables similar to that of all lakes in the U.S. Similarly, Rasilo et al. (2015) report that diffusive CH 4 emission rates are inversely related to lake size in boreal Canada. Estimates of diffusive CH 4 emissions for lakes across this biome would therefore be biased unless the published data were derived from a collection of lakes with a size distribution that mirrors the true size distribution of lakes in the boreal landscape. Ultimately, upscaling with models that explicitly incorporate such environmental drivers (i.e., size and productivity) is likely to produce more accurate regional or global emission estimates (see Downing 2009).
A longstanding challenge to upscaling based on environmental covariates is the limited availability of spatially explicit data sets on lake characteristics at the global scale. Recent advances in satellite imagery are beginning to address this issue, however. For example, Sayers et al. (2015) report the frequency distribution of Chl a in 80,000 lakes around the world, covering the northern and southern hemispheres. Since some metrics of productivity, such as Chl a and nutrients, are potentially good predictors of CO 2 , CH 4 , and N 2 O in lakes and impoundments (Beaulieu et al. 2015;Deemer et al. 2016;DelSontro et al. 2016;Lapierre et al. 2017), the Sayers et al. (2015) data set provide a powerful tool for upscaling GHG emission rates. Similarly, improved satellite imagery and analytical approaches have allowed for several estimates of the distribution of lake and impoundment sizes across the globe (Downing et al. 2006;Verpoorter et al. 2014;Messager et al. 2016).
In this research, we address several gaps in our current knowledge of GHG emissions from lakes and impoundments. First, we estimate GHG emissions from these systems by combining data on the global distribution of lake size and productivity with new empirical models linking GHG emission rates with these variables. Second, we provide global estimates of all three major GHGs (i.e., CO 2 , CH 4 , and N 2 O). This is important because, although they may be reported less frequently than CO 2 fluxes, N 2 O, and CH 4 are more potent GHGs. On a 100-yr time horizon, N 2 O and CH 4 have approximately 265-and 28-times the effect of CO 2 on atmospheric warming, respectively, and these values increase to 298 and 34 when indirect effects of these GHGs on the climate are considered (Myhre et al. 2013). Finally, we use this new analytical approach to assess the sensitivity of lentic GHG emission to expected future changes in lake and impoundment productivity.

Methods
The general approach was to perform a secondary analysis (sensu Kiecolt and Nathan 1985) of GHG fluxes and the factors influencing them (i.e., lake size and trophic status) and to upscale these relationships to global estimates following Downing (2009), while using a global, joint lake size and productivity distribution. Although several global size distributions have been published (Downing et al. 2006;Verpoorter et al. 2014;Messager et al. 2016), global productivity distributions are less well constrained so we derived a distribution using remotely sensed lake Chl a data from the best available source (Sayers et al. 2015).

Data collection
We collected previously published GHG flux and ancillary data from 223 different studies for a total of 8233 system measurements from at least 54 countries (Supporting Information Table S1). About 7637 systems were natural lakes and 596 were considered impoundments of some type. More data were available for CO 2 than for CH 4 and N 2 O (7824 unique systems for CO 2 vs. 561, 144, and 166 with diffusive, ebullitive, and total CH 4 estimates, respectively, and 309 for N 2 O). Because we aimed to upscale based on the premise that productivity and lake-size impact emissions, we also collected nutrient (total phosphorus [TP] and total nitrogen [TN]), Chl a, and surface area (SA) data when available. The vast majority of collected water quality and gas flux measurements were made during the warm, open water season and fluxes were in daily areal units of mg m 22 d 21 . Flux data, sources, and ancillary variables are available in the Supporting Information or at figshare.com/s/c6a4133f3595b67a9816. Our compilation dataset should not be considered exhaustive as we chose to use data from studies with large and/or easily accessible datasets as well as containing the ancillary data necessary for our approach (Supporting Information  Table S1).

Data treatment
All CH 4 and N 2 O flux data we compiled from the literature were used in our analysis. The compiled CO 2 data set contained an order of magnitude more observations than the other GHG data sets, however, and was pruned prior to analysis to reduce the effect of outliers. There are two main concerns that must be addressed when analyzing large data sets: (1) that "large patterns" (sensu Ramaswamy et al. 2000) be identified by analyzing the bulk of the dataset, and (2) outliers that represent aberrant observations (Yu et al. 2002; e.g., detect credit card fraud or irregularities in gene expression) be identified. Because we were interested in determining characteristic patterns exhibited in a significant portion of the CO 2 emission rates, we pruned this large dataset using the approach of Ramaswamy et al. (2000). This approach defines outliers not by the distance from the mean but based on the distance defined by the k th nearest neighbor to the mean, where k is defined as a fraction of the dataset, above and below the mean observation. We thus limited the CO 2 data to those observations within the 5 th and 95 th percentiles of the distribution. The efficacy of the Ramaswamy et al. (2000) pruning approach is illustrated by the fact that fitting to the CO 2 outliers (i.e., using unpruned data) resulted in a biased regression analysis that predicted the three largest lake size bins to be net CO 2 sinks, which is clearly unrealistic (cf., Alin and Johnson 2007). To reduce heteroskedasticity and normalize model residuals, all data were log 10 -transformed prior to analysis and a constant was added to data sets containing flux values less than or equal to zero prior to log transformation.
We analyzed the relationship between emission rates and predictor variables using simple and multiple linear regression. Candidate predictor variables were lake size, Chl a, TP, and TN. For each gas and emission mechanism, multivariate models were constructed for all combinations of size and one other predictor variable, including their interaction. Variable selection was conducted using the partial t-statistic at a significance level of 0.05. The model with the greatest explanatory power (i.e., r 2 ) and smallest model error (i.e., mean absolute error, root mean squared error) for each gas and emission mechanism was selected as the "best" model and used for all subsequent analysis and upscaling.

Global lake size and productivity distribution
In order to upscale gas fluxes to global rates, it is necessary to know the worldwide joint distribution of lake size and productivity. Three size analyses are available that quantify the global size distribution and abundance of lakes and impoundments (Downing et al. 2006;Verpoorter et al. 2014;Messager et al. 2016). These analyses were done with different methods and yield slightly different global SA and size distributions. Downing et al. (2006) estimated the global size and abundance distribution of lakes using data on the largest systems and a Pareto function to model the remaining systems. This analysis covered nine logarithmic size categories from 0.001 km 2 to > 100,000 km 2 (Supporting Information  Table S2). Verpoorter et al. (2014) estimate was created by applying an automated algorithm to high resolution satellite imagery of lakes and impoundments, covered eight logarithmic size categories, did not include the Caspian Sea, and could not resolve the smallest (< 0.002 km 2 ) waterbodies. We determined the total lake SA within each of the eight size categories by digitizing Fig. 2b of Verpoorter et al. (2014). We then added a new category for the Caspian Sea using data from Downing et al. (2006) (Supporting Information  Table S3). Messager et al. (2016) estimated the number and area of lakes and impoundments in six logarithmic sizecategories (0.1 km 2 to > 10,000 km 2 ) from their "HydroLAKES" database created by unifying several GIS databases and calculating areas of smaller lakes from Pareto distributions. We calculated the area of lakes < 0.1 km 2 to match the previous two distributions by subtracting the area listed for lakes from 0.1 km 2 to > 10,000 km 2 in Table 1 of Messager et al. (2016) (2,672,900 km 2 ) from their estimate of the total area of lakes > 0.01 km 2 (3,232,200 km 2 ). We further modified the Messager et al. (2016) distribution by substituting the data in the largest size category (> 10,000 km 2 ) with the canonical area of the largest systems reported in Table 2 of Downing et al. (2006) (Supporting Information Table S4).
We estimated the global distribution of lake productivity using satellite-based remote sensing measurements of Chl a on 80,000 lakes around the world, covering the northern and southern hemispheres (Sayers et al. 2015). The data divide the world's lakes into productivity bins of 5 lg L 21 width (e.g., 0-5, 5-10, etc. up to 1001). We used these data to calculate the proportion of global lake and impoundment SA within each Chl a bin. We then propagated this Chl a distribution across the Downing et al. (2006), Verpoorter et al. (2014), and Messager et al. (2016) global size distributions mentioned above, except for the two largest size bins. We therefore assumed that each size bin would have the same Chl a distribution as derived for the world's lakes, which is supported by the lack of correlation globally between lake size and both Chl a and TP in our collected dataset (Supporting Information Fig. S1). For the two largest size bins, however, we used literature data to determine the size and Chl a of the world's 17 largest systems (Supporting Information Table S1). The joint size-productivity distribution is insensitive to error in individual lake Chl a estimates because the bins are of 5 lg L 21 breadth. We defined the joint distribution of TP and system size by converting Chl a bins to TP bins using a regression relationship between TP and Chl a calculated using the global data of McCauley et al. Table 2. Global carbon emissions (Pg C-CO 2 eq yr 21 for each GHG) from lakes and impoundments and individual contributions to total radiative forcing from aquatic waterbodies.  (1989). The resulting joint size-productivity distributions for the Downing et al. (2006), Verpoorter et al. (2014), and Messager et al. (2016) lake and impoundment size distributions are presented in Supporting Information Tables S2-S4.

Upscaling procedures
Size-productivity weighting approach We used the size and/or productivity model that best explained the emission rates of the individual GHGs in our datasets to predict an areal emission rate and 95% confidence interval (CI) for each unique combination of lake size and productivity represented in the three joint size-productivity distributions (i.e., Supporting Information Tables S2-S4). Next, we multiplied the predicted areal emission rate for each sizeproductivity bin by the total SA of the bin, then summed across productivity bins to find the total amount of gas emitted daily from each size class. Finally, summing daily emissions from all size classes and multiplying by 365 provided an annual global estimate of GHG emissions from lakes and impoundments according to the three different size distributions. These calculations were repeated using the lower and upper 95% CI of the predicted emission rates, resulting in a 95% CI for the emission estimate for each gas and lake size distribution. This 95% CI reflects uncertainty in the statistical models used to predict emission rates, but does not account for uncertainty in the flux estimates reported in the literature. Spatiotemporal variability (Tranvik et al. 2009;Bastviken et al. 2011) and the use of different measurement approaches (Deemer et al. 2016) can add multiple levels of uncertainty to flux estimates. This uncertainty is rarely quantified, however, and is not reflected in the 5% and 95% CI of the statistical models used in this secondary analysis.

Traditional approach
The traditional approach to upscaling entails multiplying an average areal emission rate by a global lentic SA estimate (e.g., Cole et al. 2007;Tranvik et al. 2009;Bastviken et al. 2011). We applied this approach to each GHG using the mean emission rate calculated from the collection of measurements used to derive the models described above and scaled that rate to the globe using each of the three global lake and impoundment size distributions described above (Downing et al. 2006;Verpoorter et al. 2014;Messager et al. 2016).

Results and discussion
Nutrients and lake size drive aquatic CO 2 emissions Out of the diverse candidate models (Supporting Information Table S5), the model that best explained CO 2 flux included lake size and TP as well as their interaction (Supporting Information Table S6). CO 2 emission rates decreased with size across all productivity levels (Fig. 1A). CO 2 flux was positively correlated with TP in the smallest systems and negatively correlated with TP in medium to large-sized systems, but Chl a was not included in the best CO 2 model (Supporting Information Table S5). TP is likely serving as a proxy for watershed inputs of terrestrial materials. High TP concentrations in smaller lakes typically reflect large watershed inputs of water and solutes (Arbuckle and Downing 2001) including terrestrially sourced CO 2 and/or dissolved organic carbon (DOC) (Palviainen et al. 2016). Indeed, others have found that TP tends to covary with DOC in lakes (Lapierre and del Giorgio 2012) and that direct inputs of CO 2 (Maberly et al. 2013;Weyhenmeyer et al. 2015) from the surrounding watershed, as well as inputs of biologically and photochemically degradable DOC (Lapierre et al. 2013), drives CO 2 concentration in lakes. Thus, the positive correlation between CO 2 flux and TP that we observed in the smallest lakes and impoundments could merely represent the impact of nutrient and C runoff from the watershed on C dynamics in the smaller systems, and the negative effect of TP in larger systems reflects enhanced primary production in pelagic waters that tends to be more important in those systems. Our results support previous findings that CO 2 dynamics in the world's lentic waters are directly impacted by activities in their watersheds (Tranvik et al. 2009;Regnier et al. 2013) and that future changes in watershed properties, such as climate-induced changes to runoff (Milly et al. 2005), can alter the C balance of the world's lakes and impoundments.

Productivity-driven aquatic CH 4 emissions
The only factor retained in the best models for ebullitive and total (i.e., diffusive 1 ebullitive) CH 4 emission rates was a positive effect of Chl a ( Fig. 1B; Supporting Information Table S6). Chl a was also positively related to the diffusive CH 4 emission rate, though the strength of the effect was stronger for large rather than small systems. While previous studies have reported a positive relationship between Chl a and CH 4 emissions for lakes (Bastviken et al. 2004) and impoundments (Deemer et al. 2016), this is the first study to generalize the finding to lentic systems across the globe. This relationship could arise through several mechanisms, including the correlation between Chl a and the oxygenpoor and C-rich environments that favor CH 4 production and often develop in eutrophic systems. There is also evidence that the labile C derived from autochthonous production is more readily converted to CH 4 than more recalcitrant allochthonous C sources (West et al. 2012). Our findings suggest that the eutrophication of the world's lentic ecosystems could increase global CH 4 emissions.

Productivity and lake size influence N 2 O emissions
In aquatic ecosystems, N 2 O is produced primarily by the oxidation and reduction of nitrogen (N) via microbial nitrification and denitrification, respectively. Nitrous oxide is one of several N-species that can be produced by these processes and rates of N 2 O production are determined by the overall rates of de/nitrification and the fraction of de/nitrified N converted to N 2 O (i.e., the N 2 O yield; Davidson et al. 2000). In this study, we found that N 2 O emission rates exhibited a positive relationship with lake size and Chl a ( Fig. 1C; Supporting Information Table S6). Chl a is a proxy for algal biomass which can serve as a labile carbon source for denitrification and consequent N 2 O production (Sirivedhin and Gray 2006;McMillan et al. 2010;Chen et al. 2012). High algal biomass can also result in hypoxia in the hypolimnia of stratified lakes, further stimulating denitrification, an anaerobic process. Finally, Chl a may serve as an indicator of N loading, which can stimulate both algal growth (McCauley et al. 1989) and denitrification (Beaulieu et al. 2011).
The positive relationship between N 2 O emission rate and lake size has not been previously reported and may be related to differences in the relative abundance of littoral and pelagic habitats along the lake size continuum. In general, well oxygenated pelagic habitats will favor nitrification, the aerobic oxidation of ammonium, and will be relatively more abundant in large lakes, whereas carbon rich littoral zones will be relatively more abundant in small lakes. Systematic differences in the rates of de/nitrification and/or the N 2 O yield among these habitats could account for the effect of lake size on N 2 O emissions. For example, the denitrification N 2 O yield decreases as the relative availability of C to N increases (Miller et al. 2008); therefore, littoral habitats, which represent a larger area in smaller lakes, may yield less N 2 O than pelagic habitats, which dominate larger lakes. Regardless of the mechanisms driving the relationship with lake size and Chl a, our analysis illustrates that these covariates can be used to improve global estimates of N 2 O emissions from lakes and impoundments.

Upscaling global GHG emissions with size-productivity weighting
The global lake and impoundment size distribution chosen for upscaling impacted the global emission estimates for all three GHGs. With the traditional upscaling approach, the Verpoorter et al. (2014) distribution produced global emission rates 22% higher for each gas than when using the Downing et al. (2006) distribution, which in turn produced rates 35% higher than when using the Messager et al. (2016) distribution (Table 1). This pattern is consistent with the total SAs predicted by each distribution. We found a similar pattern when comparing the size-productivity weighting (SPW) approach across the three size distributions, except for CO 2 emissions. Despite the fact that Verpoorter et al. (2014) reports 22% more total lake and impoundment area than Downing et al. (2006), there is only a 4% difference in the global CO 2 emissions estimated by SPW from these two distributions. This is because the Verpoorter et al. (2014) distribution estimates less SA in the smallest class (< 0.01 km 2 ) than the Downing et al. (2006) distribution, and these small systems support the highest CO 2 emission rates (Fig. 1A). The SA discrepancy in the small size bins between Verpoorter et al. (2014) (Verpoorter et al. 2014;Feng et al. 2016;Messager et al. 2016); therefore, as the Downing et al. (2006) size distribution provides estimates falling between the two other distributions, we used it in the following sections for comparison between our upscaling approaches and with literature values that also used the Downing et al. (2006) distribution.
We found differing trends for each GHG when comparing the traditional upscaling approach to the SPW approach ( Table 1). The traditional approach, when applied to our data, produced annual emissions 50% and 80% higher than the SPW approach for CO 2 and CH 4 , respectively. Even the upper 95% CI of the SPW approach was lower than the corresponding traditional approach estimate for these two gases, except for diffusive emission of CH 4 (Table 1). For N 2 O, the SPW approach predicted 35% higher annual emissions than the traditional approach. Differences between the SPW and traditional approaches can be attributed to the effect of lake size and productivity on GHG emission rates. These relationships are explicitly incorporated into the SPW approach, whereas they are ignored by traditional published calculations. Consequently, we would expect that a simple multiplication of an average emission rate by global lake and impoundment area (i.e., the traditional approach) would yield an inaccurate GHG emission estimate unless the systems sampled to calculate the average emission rate exhibited a size and productivity distribution identical to that of the true distribution that exists worldwide.
Lack of a representative mean emission rate is likely what caused previous global estimates to diverge from that found using the SPW approach. Tranvik et al.'s (2009) estimate of 810 Tg yr 21 CO 2 is 62% higher than our mean CO 2 estimate, suggesting that lakes in their dataset were smaller and more productive than those in our global dataset. Soued et al.'s (2016) N 2 O emission estimate of 0.63 Tg yr 21 is nearly twice the mean rate we calculated, suggesting that the lakes in their dataset must have been larger and more productive than those in our global dataset. Contrary to the other gases, our SPW approach resulted in a mean total CH 4 emission rate double that of the literature-derived CH 4 emission rate of 69 Tg yr 21 from Bastviken et al. (2011), although our lower 95% CI estimate is equivalent to their estimate. The difference in estimates thus suggests that systems included in the Bastviken et al. (2011) emission rate were much less productive than those in our global dataset. Our new estimate of CH 4 emissions from lakes exacerbates the discrepancy between top-down estimates of natural global CH 4 emissions, derived from atmospheric sampling, and bottomup estimates that are derived from in situ flux measurements, with bottom-up exceeding top-down estimates by 60% (Saunois et al. 2016). Saunois et al. (2016) suggests that much of the discrepancy could be attributed to uncertainty in the areal extent of inland waters, as demonstrated by the difference in SPW emission estimates from the three lakesize distribution models (Table 2). Furthermore, difficulty in differentiating wetlands from small/shallow lakes may lead to those systems being double counted in bottom-up estimates. Finally, the accuracy of reported CH 4 emission rates for lakes is sensitive to how ebullition, the dominant CH 4 emission pathway in many systems, is measured and integrated into a system-scale estimate. System-scale ebullition rates can be greatly overestimated if investigators focus on emission hot spots without accounting for lake areas with low ebullition rates. This highlights a major limitation of using the traditional approach for global estimates in which point measurements are merely averaged and scaled upward. In general, CO 2 and CH 4 fluxes can vary by a few orders of magnitude with size or productivity (Fig. 1) so that a few divergent fluxes could easily skew a predicted global rate using the traditional approach (Table 1). Upscaling emissions based on predictors of GHG emission rates, as we did with the SPW approach, should result in a more accurate estimate of global fluxes and one that can be upscaled to diverse regions with divergent lake size and productivity distributions. Our analysis implies an evolving view of the "active pipe" model (sensu Cole et al. 2007), which describes how lakes and impoundments store terrestrially carbon in sediments and emit CO 2 to the atmosphere (Fig. 3). As traditional approaches used improved databases for the calculation of CO 2 emissions, the global role of lakes was amplified substantially (i.e., Cole et al 2007vs. Tranvik et al. 2009). Application of the SPW approach used here modulates the estimated role in emissions from 0.81 Pg yr 21 (Tranvik et al. 2009) to 0.5 Pg yr 21 when weighting global estimates by the size-and productivity-distributions of the world's lentic systems. Although the difference between the SPW estimate and the estimate of Tranvik et al (2009) seems small, it is 1.5-times the amount of C buried annually by the world's oceans and 1 =3 of the C delivered annually to the world's oceans from the continents. Accounting for important drivers of GHG fluxes can make a substantial difference in the evaluation of the role of lakes and impoundments in global GHG flux and atmospheric effects.
The relative contribution of different lake size classes to GHG emissions differed among the three gases (Fig. 2, Supporting Information Figs. S2, S3). The smaller size classes had the greatest global CO 2 emissions ( Fig. 2A), which can be attributed to the large cumulative SA of small systems and the negative relationship between CO 2 emission rates and lake size. Areal CH 4 emissions were predicted by Chl a, which has a similar distribution across most lake sizes (Supporting Information Fig. S1), therefore global CH 4 emissions scale with the total SA in each size class (Fig. 2B). Areal N 2 O emission rates, on the other hand, increased with both productivity and size, resulting in a lack of correlation between global N 2 O emission and size class (Fig. 2C). Finally, in all GHG emission distributions, the size bin in which the US Great Lakes fall is easily visible by the spike in emissions in the second to largest size class (Fig. 2). Previous studies have reported that GHG emissions are greatest for large systems (Alin and Johnson 2007) or small ones (Holgerson and Raymond 2016), yet the nature of these relationships had not been heretofore thoroughly resolved. Our analysis shows that, in terms of GHG emissions, the relative importance of large and small aquatic systems depends upon which gas is being considered.

Global C emissions from lakes and impoundments and their atmospheric potential
We combined the traditional and SPW emission results of the three GHGs, as well as those from the literature, into total CO 2 -equivalent (CO 2 eq) emissions (using 100-yr timeframes for CH 4 and N 2 O conversions, including indirect effects; Myhre et al. 2013) in order to evaluate the role that lakes and impoundments play in the global C cycle and their potential impact on the atmosphere and climate. While the traditional approach suggested a doubling of the literature value for aquatic CO 2 eq emissions, global CO 2 eq emissions based on our SPW approach were remarkably similar to the previously reported estimates ( Table 2). The combined literature total of 1.74 Pg C yr 21 falls within the 95% CI range of the SPW estimate for each lake size distribution, but is closest to (just 10% lower than) the Downing et al. (2006) estimate of 1.92 Pg C yr 21 (Table 2). Moreover, the 95% CI for all three lake size distributions overlap suggesting (1) that we are close to constraining the true aquatic C emission value, and (2) that using the middle estimate found with the Downing et al. (2006) distribution is justified for comparative purposes.
Interestingly, even though the upscaling estimates using the Verpoorter et al. (2014) and Messager et al. (2016) distributions were somewhat higher and lower than the Downing et al. (2006) estimate, respectively, all three distributions revealed that CH 4 has the greatest atmospheric influence of the GHGs being emitted from lakes and impoundments. On average, CH 4 accounts for 75% of total CO 2 eq emissions, while CO 2 and N 2 O account for 23% and 2%, respectively (Table 2). Despite CO 2 being emitted from lakes and impoundments at rates about four times higher than CH 4 (Table 1), the enhanced climate change potential of CH 4 implies that it is responsible for three times more of the radiative forcing by lentic systems globally than CO 2 . Although N 2 O is about 300-times more potent in the atmosphere than CO 2 (on a 100-yr time scale), aquatic N 2 O emissions from lakes and impoundments appear to have a small effect.
Ultimately, we found a global CO 2 eq emission estimate similar to that reported in the literature, but the traditional approach used to calculate those literature values misattributed the majority of the CO 2 eq emission to CO 2 . By using an approach based on the drivers of GHG emission rates, we found that CH 4 may be the dominant climate changerelated gas being emitted from lentic systems. Similar results were indicated in a recent global synthesis of emissions of all three GHGs from impoundments (Deemer et al. 2016) and together these findings support the conclusion that the productivity of lakes and impoundments plays an important role in the alteration of Earth's climate. However, because productivity enhances GHG emissions but the global Chl a distribution used in our analysis did not exceed 100 lg L 21 (Supporting Information Tables S2-S4), our estimate of the global role of lakes and impoundments is conservative. Furthermore, as a result of ongoing global change, the warming of lakes globally has already been observed (O'Reilly et al. 2015), as has the continued eutrophication of inland waters, despite concerted efforts to thwart it (Smith et al. 2014). If eutrophication leads to shifts in the GHG balance of lentic systems toward the more potent CH 4 , then we will consequently see an increase in the impact that these systems have on the climate.
To estimate the potential impact of continued eutrophication on the GHG emissions of lakes and impoundments, we incremented our joint size-Chl a distributions (Supporting Information Tables S2-S4) by 1 lg L 21 , 5 lg L 21 , and 10 lg L 21 and calculated the CO 2 eq emission of these simulated conditions using the SPW models. This analysis suggests that these moderate levels of enhanced eutrophication could increase the atmospheric effect of GHGs emitted from lakes and impoundments by 5%, 26%, or 42%, respectively (Supporting Information Table S7). This increased emission would be equivalent to around 1 Pg CO 2 eq yr 21 or about 13% of the effect of the current global emission of CO 2 by the combustion of fossil fuels, and about equal to the excess CO 2 emissions to the atmosphere from global land use change (Ciais et al. 2013).

Further improvement of aquatic C emission upscaling
The accuracy of our predicted GHG emissions from lentic systems is limited in part by the availability of data on the global distribution of key limnological variables. For example, while the Sayers et al. (2015) Chl a distribution represents an important advancement in the mapping of limnological variables at the global scale, the distribution has a maximum Chl a value of 100 lg L 21 , whereas GHG flux rates have been reported from hypereutrophic systems with Chl a concentrations in excess of 700 lg L 21 (Supporting Information Fig. S1). Therefore, we are likely underestimating GHG emissions because the most eutrophic lentic systems are not accounted for. This bias could be significant for some regions, such as the 1.1 million km 2 Southern Plains ecoregion of the United States where 25% of the lakes have Chl a in excess of 100 lg L 21 (USEPA 2017). This effect is also likely to be more important for N 2 O and CH 4 than CO 2 emissions because Chl a was a stronger covariable for those GHGs. In this study, we estimated TP from the Sayers et al. (2015) Chl a distribution, thereby allowing us to upscale CO 2 emissions using the CO 2 TP model, which explained more of the variation in CO 2 emission rates than the CO 2 Chl a model. However, CO 2 emission estimates could be improved by better global mapping of TP in lakes and impoundments. Another source of variability in modeled global emission estimates arises because of the disagreement among published global lake and impoundment distributions. We used three distributions for our analyses, but others can be found in the literature (e.g., Feng et al. 2016). Further, these distributions tally only permanently wetted environments and ignore the likely large influence of seasonally flooded ecosystems. Until a consensus is reached on the global area and size distribution of lakes and impoundments, the global lentic GHG emission estimate will remain sensitive to the chosen SA distribution model.
The accuracy of the emission estimate is also limited by methodological constraints, particularly for CH 4 . Ebullition can be the dominant CH 4 emission pathway, especially in smaller systems, yet our ability to accurately measure this flux is hampered by its spatiotemporal variability and difficulty in assessing ebullition rates at the system scale (e.g., Beaulieu et al. 2016). Much work is still needed to sufficiently address this issue. The geographic distribution and temporal resolution of reported GHG emission rates also limits the representativeness of global annual estimates. A substantial fraction of the global data is from studies in the boreal biome, and more measurements are needed across a range of lake sizes and trophic levels from poorly studied regions of the world. In addition, the global flux models we found may not be the best ones for lakes or impoundments on a regional basis. For example, differing emission drivers have been found for CO 2 emissions from lakes across various regions of Canada (Lapierre and del Giorgio 2012) and the United States (Lapierre et al. 2017). Another important gap is the relative lack of data that represent annual emissions (cf., Jones et al. 2016). Most available data in our analysis represent the open water/summer season and do not include large emissions that can occur during ice-out (Michmerhuizen et al. 1996;Striegl et al. 2001;Karlsson et al. 2013) or turnover (Kankaala et al. 2007;Schubert et al. 2012;Beaulieu et al. 2015). For example, the ventilation of accumulated gases during ice-out constitutes on average 17% and 27% of annual emissions of CO 2 and CH 4 , respectively (Denfeld et al. unpubl.). While not directly considered here, our approach for scaling daily rates to annual rates may somewhat compensate for this. As previous publications have done, we multiplied daily rates by 365 to arrive at annual rates, which implicitly assumes that emissions continue to occur during periods of ice cover. While this assumption is likely inaccurate, ice-out emissions represent GHGs that were produced and stored during periods of ice cover; therefore, our assumption that emission occurs year-round may compensate for the poor coverage of ice-out emission in the literature. Nevertheless, improved temporal resolution of measurement campaigns will reduce uncertainty in upscaled emission estimates. Finally, predictive models may be further improved through the inclusion of other environmental covariables that predict GHG emission rates from lakes and impoundments at the global scale. For example, global maps of climate, land use, and soil variables are available and could be used for upscaling if empirical relationships with GHG emissions could be demonstrated.

Conclusions
Here, we offer the first joint global analyses of all GHGs emitted from both lakes and impoundments. We found that lake size and productivity are important drivers of emission rates and used these drivers in creating new, global estimates of GHG fluxes from lakes and impoundments. This analysis indicates that lentic systems are important sources of GHGs driving climate change, but concludes that CH 4 emissions may be of disproportionate importance due to their link with lake/impoundment productivity and the high potency of CH 4 as a GHG. While we expect that as data availability evolves that our understanding of the role of lentic waters in the global GHG budget will as well, we suggest that linking emission rates with important drivers is a path toward a more accurate understanding of their global GHG role. Our study shows that GHG emissions from lakes and impoundments are equivalent to 20% of global fossil fuel CO 2 emission (9.3 Pg C-CO 2 yr 21 ; Le Qu er e et al. 2016) and that emissions will rise even further with the continued eutrophication of Earth's lentic ecosystems.