Natural Lakes Are a Minor Global Source of N2O to the Atmosphere

Natural lakes and reservoirs are important yet not well‐constrained sources of greenhouse gasses to the atmosphere. In particular for N2O emissions, a huge variability is observed in the few, observation‐driven flux estimates that have been published so far. Recently, a process‐based, spatially explicit model has been used to estimate global N2O emissions from more than 6,000 reservoirs based on nitrogen (N) and phosphorous inflows and water residence time. Here we extend the model to a data set of 1.4 million standing water bodies comprising natural lakes and reservoirs. For validation, we normalized the simulated N2O emissions by the surface area of each water body and compared them against regional averages of N2O emission rates taken from the literature or estimated based on observed N2O concentrations. We estimate that natural lakes and reservoirs together emit 4.5 ± 2.9 Gmol N2O‐N year−1 globally. Our global‐scale estimate falls in the far lower end of existing, observation‐driven estimates. Natural lakes contribute only about half of this flux, although they contribute 91% of the total surface area of standing water bodies. Hence, the mean N2O emission rates per surface area are substantially lower for natural lakes than for reservoirs with 0.8 ± 0.5 versus 9.6 ± 6.0 mmol N·m−2·year−1, respectively. This finding can be explained by on average lower external N inputs to natural lakes. We conclude that upscaling‐based estimates, which do not distinguish natural lakes from reservoirs, are prone to important biases.


Introduction
Inland waters are important sources of greenhouse gasses to the atmosphere. However, global-scale estimates are afflicted by large uncertainties, in particular for nitrous oxide (N 2 O), which is the third most important long-lived greenhouse gas. In its fifth assessment, the Intergovernmental Panel on Climate Change (IPCC) reports existing estimates of total N 2 O emissions from continental waters with a range of 7.1 to 207.1 Gmol N year −1 , representing a warming potential of 0.1-2.7 Pg CO 2 -eq. year −1 on a 100-year time horizon (Ciais et al., 2013). This estimate for continental waters comprises only rivers, estuaries, and the coastal zone and excludes standing water bodies (SWBs) defined here as the sum of natural lakes and reservoirs (artificial lakes created by damming).
While global estimates of N 2 O emissions from running waters have been presented in a number of studies (Beaulieu et al., 2011;Hu et al., 2016;Kroeze et al., 2005Kroeze et al., , 2010Seitzinger et al., 2000), the first global estimate for SWBs was published in 2015 by Soued et al. They estimated a global N 2 O emission of 45.0 ± 21.1 Gmol N year −1 . Three years later, DelSontro et al. (2018) reassessed these emissions using simple upscaling and different estimates of global SWB surface area, as well as an empirical model relating N 2 O emission rates to lake productivity and lake size, obtaining values ranging from 11.3 to 27.1 Gmol N year −1 . While the studies by Soued et al. and DelSontro et al. did not distinguish natural lakes from reservoirs, Deemer et al. (2016) estimated a global emission of 2.4 Gmol N 2 O-N year −1 from reservoirs only. The large uncertainties in these estimates, however, do not allow to assess with confidence the relative contributions of reservoirs and natural lakes to the overall N 2 O emissions from SWB. For this, a consistent methodological framework that distinguishes both types of water bodies is needed.
The studies by Soued et al. (2015), Deemer et al. (2016), andDelSontro et al. (2018) used an upscaling technique that consists in multiplying an averaged, observation-based N 2 O flux rate by an estimate of the total water surface area. The low number and uneven distribution of SWB N 2 O data explain the large uncertainty in flux estimates and may also cause a bias in the upscaling. Soued et al. (2015) based their estimate on their own observations in Quebec and 157 literature values, of which 137 (88%) were from the temperate zone. Deemer et al. (2016) upscaled their global value from observed fluxes of only 58 reservoirs. DelSontro et al. (2018) have used the largest amount of literature data so far, comprising in total 309 observed N 2 O emission rates. In addition, uncertainties still exist with regard to the total surface area, number, and size distribution of SWBs at global scale (Fluet-Chouinard et al., 2016), and the use of different estimates and data sets leads to very different upscaled results as demonstrated by DelSontro et al. (2018). Moreover, the three global studies so far give lumped estimates of SWB N 2 O emissions without exploring the global-scale spatial patterns. Due to the uneven global coverage of observed N 2 O emission rates, it is not possible to explore the global-scale spatial pattern of N 2 O emissions from SWBs based on empirical upscaling. Global estimates of riverine N 2 O emission (Beaulieu et al., 2011;Hu et al., 2016;Kroeze et al., 2005), on the contrary, rely on empirically estimated emission factors linearly relating N 2 O emission rates to riverine N loads, which were then applied to global-scale assessments of riverine N loads. Existing, spatially explicit representations of riverine N loads (e.g., Mayorga et al., 2010) permitted for similarly spatially explicit estimates of riverine N 2 O emission estimates (e.g., Hu et al., 2016).
Recently, Maavara et al. (2019) developed a process-oriented model to consistently estimate N 2 O emissions from a broad range of inland waters, including reservoirs. Based on a spatially explicit representation of the river-reservoir network and point and nonpoint sources of nitrogen (N) and phosphorous (P), this model made it possible to explore and quantify the global-scale spatial patterns in riverine and reservoir N 2 O emissions in a consistent manner. In this study, we apply the approach of Maavara et al. (2019) to a global data set of 1.4 million SWBs, which are classified into natural lakes and reservoirs . Model results for N 2 O emissions from SWBs are evaluated against regional, observation-driven estimates.
The main objectives of our study are to predict and analyze the global-scale spatial patterns in N 2 O emissions from all standing waters and to quantify the contribution of natural lakes versus reservoirs to the global N 2 O emission flux. We hypothesize that these contributions largely scale with their contribution to the overall water surface area. We thus expect global N 2 O emissions from SWBs to be dominated by the contributions from natural lakes, as they contribute 91% of the global SWB surface area . Further, we hypothesize that global-scale spatial patterns of these emissions follow largely the distribution of SWB surface areas and are therefore different from the spatial patterns of riverine N 2 O emissions. This would be consistent with what was found for riverine versus SWB emissions of CO 2 (Raymond et al., 2013).

Materials and Methods
In this study, we make use of a methodological framework recently established to estimate N 2 O emissions from inland waters (Maavara et al., 2019). The model structure was kept simple to be applicable to a broad range of inland water types, and indeed, the global-scale results from Maavara et al. (2019) using this model are comparable to the results of recent, measurement-based upscaling methods used for reservoirs (Deemer et al., 2016) and rivers (Hu et al., 2016). To date, this methodology has been applied to rivers, reservoirs, and estuaries (Maavara et al., 2019). Here we apply it to a data set of 1.4 million SWBs, which include both natural lakes and reservoirs . In the following, we provide an overview about the model construction (section 2.1), before we describe the preparation of input data sets and the global-scale application of the model (section 2.2) and the evaluation of model outputs against observation-driven estimates based on literature data (section 2.3).

N Model
For their global estimate, Maavara et al. (2019) used a set of simple equations that predict the N budget and N 2 O emission of each inland water body based on inputs of total N (TN in ) and total P (TP in ) from the watershed and the water residence time τ r . It was, however, impossible to train these equations directly on observed data due to the limited amount of empirical work available to constrain these specific relations. To overcome that limitation, Maavara et al. (2019) developed a generalized, mechanistic mass balance model, which represents all important biogeochemical processes involved in aquatic N cycling and which is valid for a broad range of inland water types including rivers, reservoirs, and estuaries. The major biogeochemical processes represented by the model include primary production, nitrification, denitrification, mineralization, solubilization, N fixation, and burial of N in sediments as well as N 2 O production and N 2 O emission (see Figure S1 in the supporting information). That model itself could not directly be applied for a spatially explicit estimation of N 2 O emissions from inland waters at the global scale due to the lack of geodata representing all necessary drivers. Instead, the model was fed with hypothetical but realistic combinations of parameters defining inland water properties (water body volume, discharge, surface area, temperature, and elevation), parameters defining quantity and quality of nutrient inputs, and parameters defining the kinetics of all biogeochemical processes involved. For all these parameters, Maavara et al. (2019) defined probability density functions (PDFs) approximating the global statistical distribution of those parameters as derived from literature values and databases. Then, the model was run in a Monte Carlo simulation with 8,000 iterations randomly choosing parameter values from these PDFs and generating a database of annual N cycling, N 2 O production, and N 2 O emission rates in/from 8,000 hypothetical water bodies spanning the entire parameter space of physical characteristics. This database was then used to fit the simpler equations relating inland water N and N 2 O budgets to TN in , TP in , and τ r , which were then used for the global upscaling. A similar strategy was previously applied to quantify global-scale silica, phosphorous, and organic carbon cycling in reservoirs (Maavara et al., 2014(Maavara et al., , 2015(Maavara et al., , 2017. In the following, we briefly summarize the fitted equations used for upscaling the nitrification, denitrification, N fixation, and burial fluxes defining the N and P budget of each water body. For more details on the mechanistic model used to produce that database, we would like to refer to the original publication by Maavara et al. (2019). Next, we give more detailed explanations on the assumptions and equations used to simulate N 2 O production and emission, which were slightly adapted here and which are of central interest for our study.
The nitrification (Nitrif) and denitrification (Denit) rates as well as the burial flux of N (TN burial ) were fitted to equations using as predictors the total input of N, that is, the sum of TN in and the nitrogen fixation flux (Fix), and the water residence time τ r , assuming that the dependence follows an error function (erf, equations (1)-(3)), that is, a sigmoid-shaped, monotonically increasing function that asymptotically approaches a maximum value for high τ r . The mean and median τ r of the 1.4 million SWBs in the HydroLAKES data set used here are 4.9 and 1.2 years, respectively, with 1% of the water bodies having a τ r of 50 years and longer. If the molar ratio TN in :TP in is lower than 30, Fix is calculated according to equations (4) and (5), while for higher TN in :TP in ratios, Fix is assumed to be zero. First, the fraction of Fix on the total N input to the water body (N fix ) is estimated from TN in :TP in and τ r (equation (4)). Then, Fix is calculated from TN in and N fix (equation (5)). The outflow of N is then calculated by subtracting the losses due to Denit and TN burial from the total N inputs (TN in + Fix). Phosphorus is lost from the water column only by burial in sediments (TP burial ), which is as well estimated from τ r (Taylor Maavara et al., 2015;equation (6)).
In the above equations, τ r is given in years, N fix is expressed in %, and all other variables are in moles of N or P per year. Numerical values in equations (1)-(3) and (6) were obtained from fitting the mechanistic model outputs.
As reviewed in Maavara et al. (2019), N 2 O production and/or emissions have in earlier studies been estimated from nitrification/denitrification rates or from N loads based on empirically derived emission factors (EFs). The wide range of published EFs is one of the main reasons for the wide range of existing 10.1029/2019GB006261

Global Biogeochemical Cycles
N 2 O emission estimates. When scaled to process rates, the EFs correspond to the production of a relatively small fraction of N in the form of N 2 O during incomplete nitrification and denitrification. When these processes are complete, ammonium is entirely oxidized to nitrate (nitrification), and nitrate is entirely reduced to dinitrogen gas (denitrification). Following Maavara et al. (2019), the EFs used in our mechanistic model allow to estimate the production of N 2 O (N 2 O prod ) from simulated Denit and Nitrif rates (equation (7)). Further in agreement with that study, we adopted an EF of 0.9% (for both denitrification and nitrification) as best estimate, which was derived as arithmetic mean of 40 literature values reviewed by Beaulieu et al. (2011).
In order to estimate the N 2 O emission flux (N 2 O em ), we applied the two "default simulation" setups DS1 and DS2 used in Maavara et al. (2019). DS1 represents a simplified approach that assumes N 2 O em to equal N 2 O prod , thus using equation (7) for the global upscaling of N 2 O em . For DS2, in contrast, it was taken into account that in the process of denitrification, the initially produced N 2 O can further be reduced to N 2 if it is not evading sufficiently rapidly from the water column. Inclusion of this loss term required the explicit representation of the N 2 O pool in the mechanistic model ( Figure S2). To account for the progressive consumption of N 2 O produced via denitrification in water bodies with long τ r , the EF (0.9%, 0.3%, or 1.5%; see above) was modulated by the inverse function of the term in square brackets in equation (2) in such a way that N 2 O production now reads Note that for nitrification, we are considering that N 2 O is not consumed, and the EF in DS2 is the same as in DS1. N 2 O emission in DS2 was then calculated assuming that only the fraction of N 2 O that is supersaturated with respect to the equilibrium atmospheric N 2 O concentration is emitted. Since N 2 O production and emission differ in DS2, a new Monte Carlo simulation was performed for this scenario, and from the output, a single equation for the total (nitrification + denitrification) N 2 O emissions was fitted (equation (9)). Note that this empirical equation relates the N 2 O emissions directly to TN in and τ r and thus represents denitrification and nitrification as N 2 O sources in an implicit manner.
with a = 2.28e-3, 0.79e-3, or 3.79e-3 and b = 1.63, 1.96, or 1.62 for an EF = 0.9%, 0.3%, or 1.5%, respectively. Figure 1 illustrates the difference between DS1 and DS2 in simulated N 2 O em relative to N inputs as function of τ r and the EF chosen. As upper-bound (EF = 1.5%) and lower-bound (EF = 0.3%) estimates deviate from the best estimate (EF = 0.9%) in a close to symmetric way, we report our bounded estimates as best estimate ± average difference from the bounds. Note that up to a τ r of 7 months, both scenarios give relatively similar emission rates. For higher τ r , DS2 gives consistently lower N 2 O em estimates, and the maximum ratio N 2 O em to N inputs is reached already after about 1 year, while under DS1 that ratio keeps increasing up to a τ r of about 5 years.

Preparation of Standing Water Body Data Set
The global HydroLAKES database  contains information on 1.4 million SWBs with a minimum surface area of 0.1 km 2 , including their surface polygons and associated area, their upstream watershed area, and estimates of water volume and water residence time τ r . HydroLAKES distinguishes three different types of SWBs: (1) natural lakes without dam constructions, (2) reservoirs, and (3) natural lakes hydrologically regulated by dam constructions. The latter two types of water bodies cover only 6,796 cases and were taken from the Global Reservoir and Dam database (Lehner et al., 2011), which was used in the study of Maavara et al. (2019). For our study, we use the classification of SWBs as defined in HydroLAKES but, as long as not indicated otherwise, combine classes (1) and (3) as "natural lakes" (i.e., natural origin). We chose the HydroLAKES database over remote sensing derived maps in raster format (e.g., Verpoorter et al., 2014) or statistical approaches (e.g., Downing et al., 2012), as it uniquely provides individual SWB surface polygons that can be topologically connected within the river network based on their outlet coordinates, which are associated with the drainage network of the global HydroSHEDS database at 15 arc second (~500 m) spatial resolution (Lehner et al., 2008).
In order to estimate the TN in and TP in loads to each SWB, we combined the information of HydroLAKES with a spatially explicit representation of TN and TP mobilization from the watershed into the river network (after Bouwman et al., 2009;Van Drecht et al., 2009;Mayorga et al., 2010; see Maavara et al., 2019, for details). In short, we redistributed the yields of dissolved inorganic and organic, as well as particulate, N and P from point and nonpoint terrestrial sources, as they are represented in the GlobalNEWS v2 model, over a 0.5°raster. Then, we overlaid that raster with the HydroLAKES database. To find the N and P terrestrial contributions to each natural lake or reservoir, we multiplied the TN and TP yields by the direct upstream watershed area (net watershed), that is, excluding those areas that are intercepted by upstream lying SWBs. Furthermore, we added the TN and TP outflows of upstream SWBs as TN and TP inputs to downstream lying SWBs, following the topologic information provided by the HydroSHEDS drainage network. Following the methodology described in Maavara et al. (2019), we used the watershed area of and the distance between SWBs to estimate the average riverine travel time of N and P from the terrestrial sources to the SWB and from the outlet of one SWB to the next downstream lying SWB, respectively. This travel time is used as τ r to estimate gains of N through fixation (equations (4) and (5)) and losses of N through denitrification (equation (2)) during riverine transport. Losses of N and P due to burial are assumed to be zero during riverine transport. In order to apply this model to the global data set, we calculated the N and P budgets of SWBs in sequence from upstream to downstream so that outflows from SWBs can easily be cascaded to the next downstream lying SWB. Further, we estimated the TN concentration in each SWB by dividing the simulated TN outflow by the discharge reported in HydroLAKES, assuming concentrations in the outflow to equal average concentration in the water column.

Validation Against Observation-Based Flux Estimates
We compared our simulation results of N 2 O emissions against regional estimates based on literature values of N 2 O emission rates or N 2 O concentrations. In case only N 2 O concentrations were reported, we calculated emission rates based on an a gas exchange velocity k 600 estimated from water body size (Read et al., 2012), water temperature and an assumed atmospheric partial pressure of N 2 O (p N2Oatm ) of 0.32 μatm for studies after 2000, and 0.31 μatm for studies before 2000, respectively (see below; equations (10) and (11) (7), panel a) and DS2 (equation (9), panel b). For each scenario, we plotted the function using the EFs of 0.9%, 0.3%, and 1.5% (see text). Simulated N 2 O em scale linearly with N input from upstream (TN in ; DS2) or the sum of TN in and N fixation (Fix; DS1). In this graph, we normalized N 2 O em accordingly.
area, which is reflecting the stronger wind shear effect over a larger fetch (Read et al., 2012). Water temperature was used to calculate the solubility constant K N2O after Weiss and Price (1980;equation (11)) and the Schmidt number (Sc) for N 2 O after Wanninkhof (1992), assuming salinity to be zero (equation (12)). In case water temperature was not reported, we used the empirical equation from Lauerwald et al. (2015;equation 13) to estimate the mean annual water temperature from mean annual air temperature at the location of the SWB, which was in turn derived from the WorldClim2 database at 1-km spatial resolution (Fick & Hijmans, 2017). We chose to compare simulated N 2 O emissions versus observation-based estimates per regions, that is, for a group of SWBs, rather than performing an evaluation of simulation results for individual SWBs for three main reasons: 1. The HydroLAKES database used for our model gives names only for the largest SWBs, and it is thus difficult to identify individual medium sized or small SWBs described in the literature. 2. Our global-scale modeling approach is too coarse to give reliable estimates for individual SWBs, in particular as TN and TP mobilization from the watershed to individual SWBs are derived from global scale, empirical models at 0.5°resolution (see section 2.2). 3. The number of observations per SWB is in most cases rather low, and the uncertainties on whether these observations are representative of average annual emission fluxes are high and unconstrained considering the high temporal and small-scale spatial variations within each water body.
In total, we identified seven regions from the literature review for which we could compare our simulation results to observation-driven estimates. Table 1 lists for each of these regions the sources and methodology behind the observation-based estimate and the selection of simulation results for comparison. In addition, we compared simulation results against direct observations from three large SWBs: Lake Kivu in East Africa (Roland et al., 2017), as well as Lake Poyang (Liu et al., 2013) and Lake Taihu (Wang et al., 2009), which are the two largest freshwater lakes in China. In order to compare our simulation results to observation-driven estimates of N 2 O emission rates per area of water surface (in mmol N·m −2 ·year −1 ), we divided the simulated N 2 O emission flux by the water surface area reported in HydroLAKES for each SWB and for both simulations DS1 and DS2.  Table 1. Model results correspond to the simulated means using the best estimate for the EF (0.9%) and to the range obtained with the lower-bound (0.3%) and upper-bound (1.5%) EF values. Table 2 provides the numerical values for the simulated and observed mean and standard deviations of N 2 O em /A water rates for each region or large SWB. Overall, we find that our model is able to reproduce the large-scale spatial trends in average  (Huttunen et al., 2002(Huttunen et al., , 2004 Results from all natural lakes in Finland (N = 7,891) Switzerland Calculated from reported N 2 O concentrations water surface area from 14 SWBs (Mengis et al., 1997) and average T water estimated from T air Results from the same SWBs identified in HydroLAKES by name and location (N = 14)

Model Validation
Tianjin region (China) Observed N 2 O emission rates from 10 natural lakes (Liu et al., 2015) All All SWBs within administrative boundaries of Quebec (N = 195,000) Great Lakes Reported N 2 O emission rates from Lake Huron, Lake Erie, and Lake Ontario (Lemon & Lemon, 1981) Results from the same large SWBs identified in HydroLAKES by name and location (N = 3) Figure 2. Regional averages of observed versus simulated N 2 O emission rates for the two simulation runs DS1 and DS2. The round points correspond to the best model estimates based on an EF = 0.9%; the whiskers represent the range defined by alternative lower and upper bound estimates based on EFs of 0.3% and 1.5%, respectively. Blue: regions with groups of SWBs; light blue: alternative results for two regions that were constrained by a maximum TN concentration of 26 μM (#5′) and 280 μM (#3′) as reported in the reference papers (see Table 1); red: single large SWBs. The dashed line represents the 1:1 line. For exact averages and standard deviations in flux rates per region, see Table 2.

10.1029/2019GB006261
Global Biogeochemical Cycles N 2 O em /A water rates across middle to high latitudes. Moreover, simulated average N 2 O emission rates per water surface area for DS2 are relatively close to the observation-driven estimates, while DS1 appears to overestimate average flux rates. Exceptions are the regional estimates for Tianjin region in China and the Republic of Ireland, which are substantially overestimated with both DS1 and DS2. The Tianjin region is highly urbanized, and thus, the map of estimated N inputs to the inland water network, which we used to calculate TN inflows to lakes, gives particularly high yields. As Liu et al. (2015) report TN concentrations for their 10 observed lakes in this region, we can further reduce the set of SWBs used for the model-data comparison to those SWBs with a simulated TN concentration that does not exceed their maximum reported TN concentration of 290 μM. Applying this maximum TN concentration, we retain only 141 of the 365 SWBs in the Tianjin region. The simulated N 2 O em /A water from these 141 SWBs are in fact quite close to those reported by Liu et al. (Figure 2 and Table 2). Similarly, the Irish lakes studied by Whitfield et al. (2011) are described as oligotrophic, which implies a maximum TN concentration of 26 μM (Carlson, 1977;Kratzer & Brezonik, 1981). Applying this value as upper threshold, we only retain nine of the 350 previously selected SWBs, and the simulated average N 2 O em /A water rate is again quite close to the observations by Whitfield et al. (Figure 2 and Table 2).
As stated in the methods section, a SWB-by-SWB comparison of simulated versus observation-driven N 2 O em /A water rates is a questionable approach. However, for two of the three large SWBs for which we performed this comparison, we find a good agreement between simulated and observation-driven estimates, while N 2 O emission from Lake Kivu appear to be overestimated by our model.

Patterns of Simulated N 2 O Emissions Across Space and Per Type of Standing Water Body
In what follows, we restrict our analysis to results for the more elaborated scenario DS2, which is identified in section 3.1 as the better performing model setup. At the global scale, we estimate a N 2 O emission from SWBs of 4.5 ± 2.9 Gmol N year −1 . The simulated spatial distribution of SWB N 2 O emission per continental surface area (N 2 O em ; Figure 3a) is mainly following the geographic distribution of these water bodies ( Figure S2a). The simulated spatial patterns of N 2 O em /A water rates (Figure 3b), on the contrary, mainly follow those of TN inputs to the inland water network ( Figure S2b). According to the HydroLAKES database used here, about 44% of the area of SWBs is found at latitudes higher than 50°(here and in the following, latitudinal zonation includes both hemispheres). However, this latitudinal zone contributes only to about 20% of the global SWB N 2 O emissions due to low average N 2 O emission rates per water surface area resulting from  Figure 3). Accordingly, the average N 2 O em /A water rates on both continents are rather low compared to those estimated for Asia, South America and Africa (Table 3). Classifying SWBs according to the definitions used in the HydroLAKES database, we estimate global N 2 O emissions from natural lakes without dams, natural lakes that are regulated by dams, and reservoirs of 1.8 ± 1.2, 0.3 ± 0.2, and 2.4 ± 1.5 Gmol N year −1 , respectively (Table 3). Surprisingly, reservoirs contribute 53% of the global lake N 2 O emissions although they represent only about 9% of the global SWB area. On the other hand, natural lakes without dam constructions contribute only 40% to global emissions, while they represent 84% of the global SWB surface area. Accordingly, reservoirs show N 2 O em /A water rates that are 1 order of magnitude larger than those for natural lakes with and without dam constructions. As N 2 O em / A water from both types of natural lakes are rather similar and as natural lakes with dam constructions play globally a minor role, we aggregate both types of water bodies in the following simply as natural lakes. The relative contribution of reservoirs to N 2 O emissions is highest in the low latitudes (<25°) with 67% and lowest in the high latitudes (>50°) with only 30%. Table 4 lists the correlations between simulated N 2 O em /A water rates, TN concentrations, and the environmental drivers used in our modeling approach, which largely explain these spatial patterns. Simulated N 2 O em /A water rates show substantial positive correlations to TN inputs (TN in ), which is one of the main drivers of our model. This strong correlation explains the high N 2 O em /A water rates in the densely populated and/or agriculturally used areas of the temperate zone ( Figure 3b) where inland water TN loads are increased by fertilizer application and sewage water inputs . We also find a high correlation of simulated N 2 O em /A water rates to TP inputs and the watershed area per SWB, which are due to the high correlations between these parameters and TN in . In contrast, we find a significantly negative correlation between N 2 O em /A water rates and the water residence time τ r , although τ r has by definition a positive effect on the ratio N 2 O em /TN in (equation (9); Figure 1b). This can be explained by the strong negative correlation of τ r to TN in as the main driver of simulated N 2 O emissions (Table 4). Simulated TN concentrations show a very high correlation to simulated N 2 O em /A water rates. This relationship also becomes apparent when  Note. All parameters have been log-transformed. All reported correlations are statistically significant at a p = 0.05 level. A wshd = watershed area, A water = surface area of the standing water body, τ r = water residence time, TP in = inflow of total phosphorous from upstream, TN in = inflow of total nitrogen from upstream, N 2 O em /A water = N 2 O emission rate per water surface area.

10.1029/2019GB006261
Global Biogeochemical Cycles comparing the regional differences in both parameters, or when comparing global estimates for natural lakes versus reservoirs (Table 3). Maavara et al. (2019) gave estimates of N 2 O emissions from all water bodies in the Global Reservoir and Dam database comprising reservoirs as well as natural lakes regulated by dam construction, which represent the lake type classes 2 and 3 in the HydroLAKES database used here. For the ensemble of these two classes, our global best estimate of 2.7 Gmol N year −1 is slightly lower than the 3.0 Gmol N year −1 given by Maavara et al. (2019) though using a similar model framework. This can be explained by the reduced simulated N inputs to these water bodies resulting from the additional N losses through denitrification and burial in the numerous natural lakes without dam constructions, which are accounted for in the present study. In contrast to Maavara et al. (2019), we group natural lakes with dam construction and natural lakes together and analyze their behavior in opposition to reservoirs as man-made water bodies.

Contribution of Natural Lakes to Inland Water N 2 O Emissions
Nevertheless, global N 2 O emissions from reservoirs are slightly higher than those from natural lakes because of higher simulated TN concentrations and N 2 O emission rates per lake surface area (N 2 O em /A water ). These results can be explained by the location of reservoirs, which are often found in more populated or agriculturally impacted areas where the anthropogenic contribution to riverine N loads is substantial (Deemer et al., 2016), while a large proportion of natural lakes is found in the more pristine boreal zone (Hastie et al., 2018;Messager et al., 2016;Verpoorter et al., 2014). Also, most natural lakes are regionally concentrated and rather disconnected from the main river network or found in mountainous headwater regions, while many reservoirs are built lower in the river network to serve their particular purposes of storing water from larger rivers for irrigation or hydropower demands. Accordingly, reservoirs have on average larger watershed areas than natural lakes, as reflected by the global averages of 12,853 versus 617 km 2 , respectively, and thus receive on average higher loads of TN.
When additionally taking into account the riverine N 2 O emissions of 3.3 Gmol N year −1 (Maavara et al., 2019), we estimate that natural lakes contribute only about ¼ to the global N 2 O emission flux from inland waters. Using the recent estimate of global river surface area of 773 ± 79 10 3 km 2 (Allen & Pavelsky, 2018), the estimated river N 2 O emission flux would translate into an area-weighted average rate of 4.3 mmol N m −2 year −1 , which is again substantially higher than that for natural lakes (0.8 ± 0.5 mmol N·m −2 ·year −1 ; see Table 3) but lower than the rate of 9.6 ± 6.0 mmol N·m −2 ·year −1 estimated here for reservoirs. This again suggests that reservoirs are predominantly found in areas where riverine N loads are anthropogenically increased.
The global-scale spatial patterns in N 2 O emissions from SWBs found here (Figure 3) are quite different from those of rivers simulated by Maavara et al. (2019). The spatial patterns depend mainly on the areal distribution of the respective inland water surface but also on the TN loads ( Figure S2). Latitudes >50°contribute about 20% to global SWB N 2 O emissions and 44% to their surface area, while the contribution of these regions to global river N 2 O emission and river surface area are only about 10% and 27%, respectively (Lauerwald et al., 2015;Maavara et al., 2019). Thus, in comparison to the rest of the globe, the contribution of high latitudes to inland water N 2 O emission is disproportionally low compared to their contribution to inland water surface area, a result that can be explained by generally low inland water TN loads ( Figure S2). On the contrary, low latitudes <25°account for about 70% of global river N 2 O emissions (Maavara et al., 2019) and 50% of river surface area (Lauerwald et al., 2015) but only about 30% and 15% of SWB N 2 O emissions and area, respectively (Table 3). Here the contribution to global inland water N 2 O emissions are disproportionally high compared to their share in water surface area, which is related to higher TN loads ( Figure S2). Moreover, the N 2 O emission from SWBs in low latitudes is dominated by reservoirs, which contribute about two thirds of the total flux and one fifth of the SWB surface area in that latitudinal band.

Comparison With Previous Studies
To our knowledge, our study is the first to give a global, spatially explicit estimate of N 2 O emissions from SWBs, and further the first study to give separate estimates for natural lakes and reservoirs in a consistent 10.1029/2019GB006261 Global Biogeochemical Cycles manner. However, a few global estimates have been published in recent years, either on N 2 O emission from standing waters without distinguishing natural lakes from reservoirs (DelSontro et al., 2018;Soued et al., 2015) or from reservoirs only (Deemer et al., 2016). In Table 5, we compare those estimates to our own simulation results. The results vary due to different methods of estimation but also due to different estimates or data sets of standing water surface area (A water ) used. The simplest method to estimate global fluxes is the upscaling based on averaged, observation-based N 2 O em /A water rates, which are then multiplied by estimates of total A water (Deemer et al., 2016;DelSontro et al., 2018;Soued et al., 2015). Using this method, DelSontro et al. (2018) estimated a N 2 O emission flux from SWBs of 11.4 to 19.3 Gmol N year −1 based on three different estimate of water surface area. These estimates are a factor 2 to 4 higher than our estimate of 4.5 (1.7-7.4) Gmol N year −1 . Soued et al. (2015) performed an upscaling for three latitudinal zones and then summed up the N 2 O emission flux estimates to a global value of 45.0 ± 21.1 Gmol N year −1 , which is even about 10 times higher than our estimate. Note however that Soued et al. (2015) attribute more than 80% of this global estimate to the latitudinal band <25° (Table 5). For this latitudinal band, they derived an average N 2 O em /A water rate of 47.2 mmol N·m −2 ·year −1 , which is more than 20 times higher than the average flux rates found for the other two latitudinal bands (Table 5) and which is based on observations from only six reservoirs in Latin America (4xBrazil, 1xPanama, and 1xFrench Guiana) taken from the study of Guérin et al. (2008). Natural tropical lakes have thus not been taken into account in the Soued et al. study. Note further that for three of these reservoirs a N 2 O em /A water rate of 5.1 mmol N·m −2 ·year −1 or lower is reported, while for the other three the reported N 2 O em /A water rate is 67.9 mmol N·m −2 ·year −1 or higher. This reveals that not all reservoirs in that area show similarly elevated N 2 O em /A water rates and highlights the huge uncertainty related to an empirical approach based on very few local measurements. High N 2 O emission from tropical reservoirs can be due to the high organic matter content in the submerged soils (Guérin et al., 2008), and decreasing greenhouse gas emission with increasing age  Downing et al. (2006) to estimate area of smaller water bodies not represented in HydroLAKES.

10.1029/2019GB006261
Global Biogeochemical Cycles of reservoirs has been reported (Barros et al., 2011) as concentration and reactivity of the submerged organic matter decrease with time (Maavara et al., 2017). In a recent study on small reservoirs in tropical Africa (Okuku et al., 2019), much lower N 2 O em /A water rates of 0.7-1.2 mmol N·m −2 ·year −1 have been reported, referring to a low organic matter content of submerged soils as explanation. Submerged soils as sources of N 2 O do not play a role for natural lakes. Moreover, the low N 2 O em /A water rate for Lake Kivu (0.3 mmol N·m −2 ·year −1 ; Table 2), which is based on an intensive sampling campaign (Roland et al., 2017), suggests that tropical lakes are generally not strong sources of N 2 O.
Interestingly, our global simulation results for N 2 O emission from reservoirs only of 2.4 ± 1.5 Gmol N year −1 are actually very close to the upscaling-based estimate by Deemer et al. (2016) of 2.4 Gmol N year −1 ( Table 5).
That finding, together with the regionalized validation of our simulation results against observation-based flux estimates for natural lakes (section 3.1), is supportive of our model-based approach. It thus represents an important methodological alternative to empirically based upscaling studies. DelSontro et al. (2018) explored an intermediate approach by applying a regression equation to predict N 2 O emission rates from lake/reservoir size and chlorophyll-a concentrations. The predicted global emissions are roughly 40-60% higher than those derived from their simple upscaling approach (Table 5). However, given the low predictive power of that equation (R 2 = 0.09), it is difficult to conclude which of the two methodologies is the most reliable for upscaling. Moreover, even if chlorophyll-a concentrations may be related to organic matter availability and N availability via N fixation by cyanobacteria while lake size may increase the organic matter degradation time, these two parameters are not directly linked with nitrification and denitrification rates.
On the contrary, some findings suggest that nitrification and denitrification are actually negatively affected by the presence of algae (Enrich-Prast et al., 2016;Risgaard-Petersen, 2003). This lack of mechanistic connection could explain the low predictive power of the empirical equation by DelSontro et al. (2018).

Model Limitations and Challenges for Future Developments
In our study, we estimated N 2 O emissions from global SWBs based on simulated TN and TP flows through each SWB and the residence time τ r within each SWB, using the model framework developed by Maavara et al. (2019). This model framework allows to estimate annual nitrification and denitrification fluxes and uses EFs to estimate N 2 O production from these fluxes, a methodology so far mainly applied to rivers. To estimate N 2 O emissions, we followed two scenarios. The first rather simplistic scenario (DS1) assumes N 2 O emission to equal N 2 O production. The other scenario (DS2) accounts for a reduced net production of N 2 O in the process of denitrification due to the reduction of N 2 O to N 2 , which becomes significant at water residence times longer than 7 months. For the 37% of the 1.4 million SWBs, which have τ r of 7 months or less, the global emissions estimated by model configurations DS1 and DS2 are nearly identical at 2.0 Gmol N year −1 . Maavara et al. (2019) also found that DS1 and DS2 gave very similar results for rivers, where τ r is generally low and gas exchange is quick For the remaining SWBs with a higher τ r , the scenario DS1 would give a total N 2 O emission flux, which is more than twice than that simulated under DS2. In our model validation against observation-based flux estimates, we found that simulation of N 2 O emissions under DS2 performs better and concluded that the reduced net production of N 2 O from denitrification at higher residence times is a nonnegligible process in SWBs, which precludes the use of constant EFs as done for rivers (e.g., Beaulieu et al., 2011;Hu et al., 2016;Seitzinger et al., 2000).
A considerable source of uncertainty in our global model results is related to the three drivers used: TN and TP inputs to the lakes and reservoirs and τ r of the lakes and reservoirs. We thus performed additional simulations quantifying the sensitivity of our global SWB N 2 O emission estimate to those three drivers by systematically varying each driver individually. When increasing or decreasing TN inputs from the watershed by 50%, the simulated global N 2 O emission is increased or decreased by 44%, respectively. Increasing or decreasing TP inputs from the watershed by 50% increases or decreases simulated emissions by 6%, respectively, which is due to the effect of TP/TN ratios in the water column on the simulated fixation rates of atmospheric N (equation (4)) and is thus a secondary effect on the total inputs of TN to SWBs. Finally, increasing τ r by 50% leads to an increase in simulated N 2 O emissions of 11%, while a decrease of τ r by 50% leads to a decrease in simulated N 2 O emission by 16%. The relatively low sensitivity of global N 2 O emission estimates to τ r can be explained by the fact that the ratio between N 2 O emission to TN inputs reaches its maximum already at a τ r of about 1 year (Figure 1b) As drivers of our model, TN inputs from upstream and τ r are not easy to determine and thus usually not reported in empirical studies. On the contrary, the concentrations of TN or important N species like nitrate are often reported and are important variables as Deemer et al. (2016), for instance, report a strong correlation between observed N 2 O flux rates and nitrate concentrations (R 2 = 0.49). Similar correlations have been reported at regional scale for natural lakes (McCrackin & Elser, 2011). Consistent with these observations, our simulation results also reproduce a high correlation between TN concentrations and N 2 O emission rates (R 2 = 0.85; Table 4). For the entirety of SWBs in our study, we calculate an mean (±standard deviation) TN concentration of 42 ± 75 μM (based on values between 1th and 99th percentile to exclude the effects of outliers; Table 3), which is actually very close to the observation-based estimate of global natural lake and reservoir TN concentrations of 38 ± 186 μM by Chen et al. (2015). That means that on global average, our simulated TN concentrations are reasonable and not the source for a general underestimation or overestimation of N 2 O emissions.
Our use of TN as substrates for N 2 O production, as opposed to using dissolved inorganic nitrogen only (Kroeze et al., 2005;Seitzinger et al., 2000), is justified by the fact that both organic and inorganic N play an important role in aquatic N cycling. Inorganic N (e.g., NH 4 + , NO 2 − , and NO 3 − ) concentrations represent only the direct availability of substrates for nitrification and denitrification. The aerobic or anaerobic degradation of organic nitrogen will result in the in situ release of NH 4 + that can be oxidized to NO 2 − and NO 3 − via nitrification and will potentially release N 2 O, depending on oxygen concentrations. Both NO 2 − and NO 3 − produced aerobically by nitrification can subsequently be reduced to N 2 O via denitrification, at low oxic conditions (Fenchel et al., 1998).
At regional scale, however, large uncertainties may result from the global-scale modeling approach. The examples of the Irish headwater lakes and lakes in the Tianjin region (China) discussed in section 3.1 illustrate that while the model is able to reproduce observed N 2 O emission rates given correct TN loads, the coarse representation of TN inputs from point and nonpoint sources (0.5°≈ 50 km) is a considerable source of uncertainty in our simulations. In particular for small lakes with small watershed areas, this may lead to a wrong attribution of TN loads. While strong point sources of TN associated with urban areas are important for the TN loads of larger rivers, they are in reality likely not linked to small lakes. Similarly, small oligotrophic head water lakes in an otherwise agriculturally impacted region will not be represented at this spatial resolution. Based on these considerations, our modeling approach may overestimate TN loads and thus N 2 O emissions from small lakes in agricultural and urban areas.
Another source of uncertainty is the representation of natural lakes and reservoirs at global scale. As seen in Table 5, different global estimates of SWB N 2 O emissions depend on different estimates of the number and total area of SWBs. The HydroLAKES database used here considers only natural lakes and reservoirs with a surface area ≥0.1 km 2 and thus gives lower estimates of SWB surface area than the estimates by Downing et al. (2006) and Verpoorter et al. (2014). Downing et al. (2006) used statistics to estimate the number of small lakes and reservoirs not represented in global databases. Verpoorter et al. (2014) used high-resolution satellite data to map areas of standing water bodies as small as 0.002 km 2 . However, the Verpoorter et al. study relied on an automated classification algorithm and ground truthing was only performed for Sweden, thus questioning its reliability for global-scale applications. Nevertheless, we may underestimate the global N 2 O emissions from SWBs, as we do not account for water bodies smaller than 0.1 km 2 . The underestimation is likely highest for high latitudes where small glacial lakes are numerous. This is supported by the total SWB surface area for latitudes >54°reported by Soued et al. (2015) on the basis of the work by Downing et al. (2006), which is nearly twice our estimate (Table 5).
Finally, an additional source of uncertainty is the necessary simplicity of our model, which simulates SWB N cycling and N 2 O emission rates at the annual time scale and which describes each SWB only through the variable τ r for the spatially explicit upscaling step. The water column and aquatic sediments are not explicitly represented, but simulated N transformation rates, including the production of N 2 O, integrate processes in both media. SWBs, however, are heterogeneous ecosystems with an intrinsic spatial and temporal variability in the different processes of N cycling. Mengis et al. (1997) and Beaulieu et al. (2015) have described sources and sinks of N 2 O for stratified SWBs varying in morphology and trophic status. Generally, in the well-mixed and oxic upper layer, the epilimnion but also in the oxic parts of the lower layer, the hypolimnion, nitrification is the dominant source of N 2 O. In anoxic parts of the hypolimnion, denitrification is the dominant process, which can be source as well as sink for N 2 O. Mengis et al. (1997) found for deep Swiss lakes that this anoxic part of the water column is often undersaturated with respect to atmospheric N 2 O, which clearly indicates a net consumption of N 2 O through denitrification. Moreover, Beaulieu et al. (2015) found for a group of reservoirs in Ohio that denitrification in the anoxic hypolimnion is a sink for N 2 O at nitrate concentrations below 3.6 μM, while at higher nitrate concentrations denitrification is a source of N 2 O. Beaulieu et al. explain that observation by the low energy yield for the reduction of N 2 O to N 2 compared to the reduction of nitrate, thus favoring full reduction only when the availability of nitrate is limited. High peaks of N 2 O in the water column have also been reported at the oxic-anoxic interface, mainly resulting from high nitrification rates sustained by underlying ammonia-rich anoxic waters (Beaulieu et al., 2015;Mengis et al., 1997;Roland et al., 2017). Conversely, nitrification of ammonia in the oxic top layers can be the main source of nitrate for denitrification in the lower layer, as observed in aquatic sediments Zhu et al., 2015). Aquatic sediments as sources for lake N 2 O emission have mainly been reported for eutrophic and shallow SWBs, and in particular for shallower parts of the SWB Zhu et al., 2015). In terms of temporal dynamics, only the epilimnion exchanges gas with the atmosphere during periods of thermal stratification. Therefore, periods of lake turnover have been reported to represent hot moments of nitrification, N 2 O production and emission as anoxic water rich in ammonium mix with oxic waters (Beaulieu et al., 2015;Roland et al., 2017). Temporal (and spatial) variations in N 2 O emissions have also been reported to be related to algae blooms and the associate increase in nitrification rates (Liu et al., 2018). However, while the global coverage with direct observations of SWB N 2 O emission is already poor, studies systematically investigating their seasonal pattern, that is, during all phases of lake mixing and aquatic growing season, are even scarcer.
More systematic fieldwork is thus needed to better understand the annual and seasonal N 2 O budgets of natural lakes and reservoirs and to better constrain their contribution to atmospheric N 2 O budgets at global scale. This scientific progress should be aided by the development and application of more sophisticated models, in particular for the global-scale upscaling of N 2 O emission estimates. Future models should represent these N processes and the individual N-species (N 2 O, N 2 , ammonium, nitrite, nitrate, dissolved organic N, and particulate organic N) in a spatially resolved framework, distinguishing processes in benthic sediments and in the water column, and further processes in the littoral and pelagic zones of SWBs, while taking into account stratification and the geometry of the SWB bed. Such a model should further represent interactions of N-processes with the organic C cycling and lake pH (Bajwa et al., 2006;Kowalchuk & Stephen, 2001). In particular, the dynamics of N 2 O sources and sinks have to be represented more explicitly to replace the use of empirical EFs. Moreover, the models will have to resolve the seasonal cycle, representing the growing season of aquatic vegetation and distinguishing periods of SWB stratification and turnover. To realize that, the representation of biogeochemical processes should be coupled to models of SWB physics driven by temporally resolved climate forcing data, like for instance the FLake model (Thiery et al., 2014). Validation of these models will require empirical studies describing the physical and biogeochemical mechanisms of N and N 2 O cycling for a broad range of SWBs of different morphology, trophic status, and climate-induced mixing regimes.

Conclusions
Using a two-step modeling approach based on a mechanistic N model and predictive equations for spatially explicit upscaling, we reestimated global N 2 O emissions from standing water bodies (natural lakes and reservoirs) at 4.5 ± 2.9 Gmol N year −1 , at the far lower end of existing estimates. Moreover, our results show that reservoirs as man-made lakes contribute more than half of that flux, although they represent less than 9% of the total SWB area. We conclude that natural lakes are a relatively small contributor to global inland water N 2 O emissions, even though their global surface area is~10 times larger than that of reservoirs and about four times larger than that of streams and rivers. While our approach is admittedly coarse, it still helps to improve our understanding of global inland water N 2 O emissions beyond what can be achieved by direct upscaling of averaged observational data or by applying empirically derived relationships such as the one relating N 2 O emissions to chlorophyll-a concentrations and SWB size. Our results, despite bearing significant uncertainties resulting from limited available observations and fragmentary knowledge about N cycling in SWBs from a global perspective, are important step forward, as it provides the first spatially explicit estimate of N 2 O emissions for 1.4 million SWBs. Further progress in this research field would benefit greatly from improved inland water databases containing smaller water bodies, from a more detailed data sets of N sources to the inland water network at higher spatial resolution and from more fieldwork on SWB N-and N 2 O cycling to improve our mechanistic understanding of these processes on regional to global scales.