Main

Reservoirs, like inland aquatic ecosystems in general, are recognized as a globally relevant source of greenhouse gases (GHGs)1. The quantification of this GHG source has been subject to intensive research for at least two decades, and current assessments estimate reservoirs to annually emit 800 Tg carbon dioxide equivalent (CO2e) to the atmosphere1. The exact quantification of this source is also important because GHG emissions from reservoirs may affect the carbon (C) footprint of irrigation2 and hydropower production and thus its perception as C-neutral energy3,4,5. At the same time, reservoirs act as a sediment trap, accumulating organic material from the catchment and in-lake primary production, and they have previously been shown to bury C at higher rates than natural lakes6. Globally, rates of organic carbon (OC) burial in reservoirs have been estimated to exceed rates of C emissions from reservoirs1,6.

Water-level changes are a typical feature of reservoirs and a major discriminator to natural lakes7. They are driven both by natural hydrological dynamics and water management. Water-level fluctuations have several consequences for reservoir limnology and biogeochemical cycling8,9. The most obvious consequence of water-level changes, however, is a changing water surface area (Fig. 1a). Varying portions of reservoirs are temporarily falling dry, forming the so-called drawdown area (Fig. 1b). This has two consequences for upscaling C fluxes from reservoirs: (1) the variable surface area has to be considered in upscaling fluxes based on water surface, and (2) the GHG emissions from the drawdown area have to be included in reservoir C budgets. Although the relevance of these two points has already been acknowledged1, they have not been considered in reservoir C budgets10 and global upscalings1,6 so far.

Fig. 1: The world’s reservoir drawdown areas.
figure 1

a, Theewaterskloof Dam, South Africa, during the water crisis of 2018. The black/white line shows the reservoir margin at maximum capacity31. Drawdown area = 73% of maximum reservoir area. b, Dry part of Hassel pre-dam, Germany, 2014. c, Global distribution of 6,749 reservoir drawdown areas. Colours indicate average drawdown area during the period 1985–2015 as percentage of maximum reservoir area. Landsat-8 image in a courtesy of USGS; photograph in b courtesy of Kurt Friese.

Source data

Full size image

It has been shown that dry aquatic sediments in general, and drawdown areas of reservoirs in particular, emit large quantities of GHGs11,12. GHG emissions from drawdown areas are typically dominated by CO2 while CH4 is of minor importance13. Initial studies reported that the drawdown area has the potential to dominate annual CO2 emissions from reservoirs12,14,15,16 due to areal CO2 emissions, which are significantly higher than those from the water surface11. CO2 emissions from dry inland waters were shown to be controlled mainly by temperature, moisture and organic-matter content of the desiccated sediments11. Extreme water-level drawdown can furthermore create hot moments of GHG emissions and may even offset C burial in the sediments16,17.

To include emissions from drawdown areas into C budgets of reservoirs and, ultimately, into global C inventories, information about the spatial extent and distribution of drawdown areas is crucial. Recently, a new algorithm was developed providing 30-year-long time series of surface-area variations in thousands of reservoirs worldwide18. In this study, we use that dataset to estimate the extent of reservoir drawdown areas globally and to identify drivers influencing the extent of drawdown areas as well as their spatial and temporal patterns. Subsequently, we re-assess global estimates of reservoirs’ C budget (defined as emission-to-burial ratios) by partitioning total emissions into fluxes from the water surface and from the drawdown area, and propagating the uncertainty from the area, and emission and burial rates estimates. The final aim of our study is to reduce the upscaling uncertainty of C fluxes from reservoirs, such as GHG emissions to the atmosphere and OC burial in sediments, relying on accurate estimates of inundated and drawdown areas in reservoirs.

Characterization of global reservoir drawdown

To examine global patterns of reservoir drawdown, we quantified the spatial extent of the monthly drawdown area of each reservoir during the period 1985–2015 (total reservoir-month n = 2,429,640). We subtracted the inundated surface area from the area covered by water at maximum filling level. Reservoirs were completely filled during only 22% of the analysed period (543,353 reservoir-months had drawdown area smaller than 5% of individual maximum reservoir area), implying that most of the time a sizeable portion of the area of reservoirs was dry. On a monthly basis, the drawdown areas covered between 9% (January 2011) and 15% (May 2003) of the total reservoir area (n = 6,749 reservoirs). On average, considering the complete period 1985–2015, 12% (only Global Reservoir Surface Area Dataset (GRSAD)19, 15% when applying Pareto modelling for reservoirs <10 km²) of the global reservoir area has been dry. The drawdown area depended significantly on reservoir size, with smaller reservoirs having larger (relative to reservoir size) drawdown areas (Extended Data Fig. 1). Reservoirs between 100 km² and 1,000 km² contributed most to the total drawdown area (12,349 km², 28% of total drawdown area, n = 349; Extended Data Fig. 2).

Drawdown-area extent was not evenly distributed across the globe (Fig. 1), and the mean annual drawdown area of reservoirs was related to latitude (Fig. 2a). Small relative drawdown areas were typical in the boreal zone whereas drawdown-area extent was largest close to the tropics, where the climatic regime is characterized by a pronounced seasonality (Fig. 2b). Smaller drawdown areas were also typical towards the Equator, where diurnal variations in climate patterns generally exceed seasonal variations. The drawdown area depended also on reservoir type (Extended Data Fig. 3). Large drawdown areas were typical for reservoirs used for irrigation (median = 23%) and flood control (median = 15%). Small drawdown areas were particularly typical for hydropower reservoirs (median = 9%). The latitudinal pattern of reservoir drawdown inversely related to the latitudinal share of the fraction of reservoirs devoted to hydropower (Fig. 2c). Because hydropower reservoirs tended to have smaller drawdown areas, the larger share of hydropower reservoirs at higher latitudes also contributed to the smaller drawdown areas at high latitudes. We used stepwise multiple linear regression (MLR) to identify the influence of environmental and anthropogenic factors on drawdown-area extent (Methods). Models containing climatic variables (precipitation seasonality, mean annual temperature, water stress) and variables representing anthropogenic factors (maximum area, main use) explained 27% and 25% of the variability, respectively. Thus, the global distribution of reservoir drawdown areas was affected by both climatic and anthropogenic drivers (global MLR model, r² = 0.41, P < 0.001, n = 6,749; Supplementary Table 1 and Extended Data Fig. 4), confirming recent findings based on a one-year dataset of global reservoir levels20.

Fig. 2: Latitudinal distribution of drawdown area and influential parameters.
figure 2

a, Average drawdown area of reservoirs. b, Precipitation seasonality. c, Share of hydropower reservoirs to total number of reservoirs. Data were binned into 10° latitude intervals. Data are means ± s.d. (a,b). Dashed lines indicate the northern tropic, Equator and southern tropic. Numbers in the box on the right indicate the total number of reservoirs per latitudinal bin.

Source data

Full size image

The global monthly average drawdown area did not show a significant temporal trend during the analysed period (Extended Data Fig. 5). However, increases of drawdown area were significant in 21% of all reservoirs (P < 0.05 for 1,404 of 6,749 reservoirs). A similar proportion of reservoirs had significantly decreasing drawdown areas (26%, P < 0.05 for 1,760 of 6,749 reservoirs). However, these findings should be handled with care because the launch of Landsat 7 in 1999 affected the temporal homogeneity of the dataset (more-frequent observations); thus, trends from individual reservoirs can be altered. Depending on latitude, total drawdown area showed a pronounced seasonality (Extended Data Fig. 5), with annual maxima occurring in fall for reservoirs located between 40° N and 90° N and either summer or spring for reservoirs located between 10° N and 40° N and between 10° S and 60° S.

Extreme drawdown (drawdown areas exceeding a given threshold in respect to maximum reservoir extent) was observed in 2,666 ± 791 reservoirs built before 1985 (40% of total reservoirs). Younger reservoirs (791 reservoirs, 12% of the dataset) were excluded here to prevent effects of reservoir filling. This estimate is based on averaging over a range of thresholds (40% to 70% of maximum reservoir area; Methods and Extended Data Fig. 6) because there is no universal definition of ‘extreme drawdown’ as what is ‘extreme’ depends on local conditions. Cape Town, South Africa, for example, prepared for ‘Day Zero’ in 2018, when the drawdown area of Theewaterskloof Dam was close to 70% of its maximum area (Fig. 1a)21.

Implications for C cycling in reservoirs

Emerging drawdown areas are inevitably coupled to shrinking water surfaces of reservoirs. Thus, all global estimates of C cycling in reservoirs that are based on their global surface area are directly affected by reservoir drawdown. Currently, global GHG emissions from reservoirs are estimated to be about 800 (500–1,200) TgCO2e yr–1 (ref. 1) by multiplying global average GHG emission rates times the global reservoir area. However, this number is based on aquatic GHG emissions only, assuming that all reservoirs are completely filled1,10, whereas dry aquatic sediments can emit large quantities of CO211,12. We re-assessed the most recent estimate of C emissions from reservoirs1 by multiplying reported C emission rates from wet and dry reservoir areas times inundated or dry areas calculated here. We also estimated the uncertainty of our calculations by propagating the uncertainty from areas and emission rates estimates using Taylor series expansion and Monte Carlo methods (Methods and Supplementary Discussion). Multiplication of the mean CO2 flux from dry sediments (mean = 611 gC m−2 yr−1, range = 9–2,301 gC m−2 yr−1, n = 47; Extended Data Fig. 7 and Supplementary Table 2) times the total extent of drawdown areas (43,539 (26,939–60,130)) km2; Methods and Supplementary Discussion) results in a total CO2 carbon flux of 26.2 (15–40) TgCO2-C yr–1. Adding CO2 emissions from the drawdown area to the emissions from inundated areas as quantified in this study implies a global emission of 60.3 (43.2–79.5) TgCO2-C yr–1. This increases the re-scaled current estimate of global CO2 emissions from reservoirs (39.5 (27.1–53.8) TgCO2-C yr–1 (ref. 1)) by 53% (Cohen’s D effect size = 2.6; Table 1 and Fig. 3).

Table 1 Summary of global estimates of CO2 and CH4 emissions and OC burial in reservoirs as derived from recent emissions and burial rates reported in the literature (‘Without drawdown areas’) and updated values after the inclusion of drawdown areas (‘With drawdown areas’)
Full size table
Fig. 3: Gaseous C emissions and OC burial in global reservoirs.
figure 3

a, Emission and burial expressed as annual teragrams of C. b, Emission and burial expressed as CO2e fluxes. Green ribbon indicates the potential range of buried CO2e assuming conversion of OC between 100% CO2 (lower bound) and 50% both CO2 and CH4 (upper bound). a,b, CO2 and CH4 emission: emissions from reservoirs comprising both drawdown area and water surface. Total emission: sum of CO2 and CH4 emissions. Vertical dot–dashed line indicates current global estimate of reservoir drawdown (15%). Shaded areas represent 95% CIs. c–e, Probability density functions as derived by Taylor series expansion and Monte Carlo uncertainty propagation. c, CO2-C emission from global reservoirs (Cohen’s D effect size = 2.6). d, GHG emission expressed as CO2e fluxes (Cohen’s D effect size = 0.1). e, Ratio between C emission and OC burial in reservoirs (Cohen’s D effect size = 1.2). C emission is the sum of CO2-C and CH4-C emission. The dashed line indicates ratio = 1. c–e, The estimate without drawdown area (red curve) is assuming the total reservoir area to be inundated while the blue curve is considering the desiccation of the drawdown area. Densities are derived by 100,000 Monte Carlo simulations (Extended Data Fig. 10).

Source data

Full size image

Total reservoir C emissions are dominated by CH4 in terms of climate forcing1. CH4 emissions from drawdown areas, however, have been shown to be low compared with CO2 emissions, but data on CH4 emissions from drawdown areas are scarce12. Assuming zero CH4 emissions from drawdown areas would reduce current global CH4 emissions from reservoirs (12.7 (6.8–23.2) TgCH4-C yr–1 (ref. 1)) by 1.8 TgCH4-C yr−1 (Table 1 and Fig. 3). However, when expressed in terms of CO2e of the three main GHGs emitted from reservoirs (that is, considering CO2, CH4 and N2O, assuming a global warming potential for CH4 of 34 (ref. 22) and for N2O of 298 (ref. 22)), no substantial change of total CO2e emissions from reservoirs can be detected given the uncertainty of the calculation (758.4 (484.3–1,235.3) to 745.4 (491.9–1,174.3) TgCO2-C yr−1, Cohen’s D effect size = 0.9; Table 1 and Fig. 3). CO2e emissions are dominated by CH4 for drawdown areas smaller than ~30% of total reservoir surface. Because this finding is based on assuming zero CH4 emissions from drawdown areas, the dominant role of CH4 will be even increased if future studies reveal substantial CH4 emissions linked to drawdown areas, potentially leading to a total increase in CO2e emissions. In fact, exposure of hypolimnetic sediments at extreme drawdown may be a substantial methane source at least in the short term16,17. Furthermore, substantial decrease of water levels in reservoirs is expected to increase CH4 emissions via ebullition due to lowering of the hydrostatic pressure23,24. Our analysis shows that 40% of reservoirs built before 1985 experienced extreme drawdown at least once. Indeed, water withdrawals from reservoirs globally have more than doubled since the 1960s due to increased demand25. Our findings show that drawdown areas have to be considered for the global C budget of reservoirs but that we need more data for refining global upscalings of GHG emissions from reservoirs.

Future research also needs to constrain remaining uncertainties. Although a recent study reported no systematic difference of CO2 emissions from dry inland waters among climate zones, the variability of CO2 fluxes across locations has been shown to be substantial11. Thus, data with better spatial coverage are required to allow for a more comprehensive geostatistical upscaling considering regional differences in both CO2 fluxes and drawdown-area extent. Analogously, more knowledge of the temporal variability of C emissions from drawdown areas is required to constrain the uncertainty induced by seasonal and long-term dynamics of C fluxes26. It was, for example, shown that CO2 fluxes from dry sediments increase with temperature11. Hence, emissions could be enhanced in regions where seasonal variability of the water level exposes an above-average extent of drawdown areas in seasons with conditions boosting CO2 emissions. Furthermore, drawdown areas with pronounced seasonal cycles of desiccation and reflooding could have a disproportionate effect on C emissions due to the seasonal resupply of fresh sediments containing more labile C. In addition, it is well known that the trophic state of reservoirs affects their water-bound GHG emissions27. One may argue that the trophic state also has an impact on emissions from the drawdown area, for example, due to sedimentation of labile OC from autotrophic production. While we combined global estimates of drawdown-area extent, C emissions and OC burial for a global reassessment of these estimates, refining our results to regional or local scales would require comprehensive studies measuring all these processes in single reservoirs. Unfortunately, the amount of data currently available precludes such an analysis. Finally, desiccation of sediments allows the emergence of terrestrial vegetation12, a typical feature of many reservoirs. The fate of the C temporarily stored in this vegetation upon reflooding (burial versus mineralization to CO2 and/or CH4) may affect the C budget of reservoirs, but again, current data available are too scarce to meaningfully incorporate this process in our calculations.

Carbon gases emitted from reservoirs to the atmosphere are recycled in the biosphere on contemporary time scales, while C stored in sediments enters the long-term geological cycle. Reservoirs have previously been shown to bury C at higher rates than do natural lakes6. Total OC burial in reservoirs is estimated to be 41.7 (21.6–73.0) TgC yr−1 by multiplying area-specific burial rates6 times global reservoir area (Methods and Supplementary Discussion). Thus, globally, reservoirs are assumed to emit 1.26 (0.66–2.58) times the C they bury, an estimate that would not allow stating emissions are significantly larger than OC burial. However, considering the effect of drawdown areas on both emissions (including dry areas) and burial (subtracting dry areas from calculations) substantially increases the global emission-to-burial ratio to 2.02 (1.04–4.26) (Table 1 and Fig. 3). This is substantially different from the previous value (Cohen’s D effect size = 1.2) and implies an overturn of the C budget of reservoirs tilted now towards C emissions. It may be argued that it is unclear how the global OC burial may vary with changing reservoir surface area due to its complex system-specific dependencies (for example, sediment focusing, dependence on allochtonous production of organic matter, and particle transport), However, even assuming a non-reduced OC burial in the C budget, consideration of drawdown areas on C emissions increases the emission-to-burial ratio to 1.72 (0.91–3.49) (not shown in Table 1). Hence, including drawdown areas in the global C budget of reservoirs is expected to shift the C mass balance towards C emission with emissions substantially exceeding OC burial, especially if global drawdown areas increase in the future (Fig. 3). In this respect, reservoirs become more similar to natural lakes, where OC burial represents ~20% of C emissions6.

Implications for further environmental services

The relevance of drawdown areas on the C footprint of reservoirs suggests that water-level management (for example, avoiding drawdown in the dry season) could be a promising tool to control C emissions from reservoirs and minimize this source of GHGs. However, reservoir drawdown has further effects on the services these systems provide to humans that should be taken into account. For example, exposure of sediments can have additional effects on water quality beyond GHG emissions, such as enhanced phosphorus release from exposed sediments on reflooding, which lead to recommendations for avoiding drawdown events during summer28. Extreme drawdown events during summer also pose a threat to fish communities and water quality in general29,30, and many leisure activities (boating, angling) and aesthetic values from reservoirs would also benefit from a management strategy avoiding drawdown events during summer. However, the provision of drinking water, irrigation or hydropower may not be compatible with keeping high water levels during the dry period. In any case, decisions on water-level management in reservoirs usually require finding trade-offs between multiple uses and environmental considerations; thus, general recommendations should not be given lightly.

Methods

Data for estimating drawdown areas

The calculation of drawdown areas was based on monthly time series of surface-area values for 6,818 reservoirs provided by GRSAD18. It comprises all reservoirs from the Global Reservoir and Dam dataset19 except of 45 reservoirs without reported geometric information. In accordance with ref. 1, we further removed 24 reservoirs classified as natural lakes that have been modified with water regulation structures (this includes lakes Victoria, Baikal and Ontario). The GRSAD dataset comprised entries from March 1984 to October 2015. To have a constant number of data points per year, we restricted our analysis to the period from January 1985 to December 2014.

GRSAD was created by correcting the Global Surface Water dataset31 for images contaminated with clouds, cloud shadows and terrain shadows. With this correction, the number of effective images that can be used in each time series has been increased by 81% on average. Substantial improvements have been achieved for reservoirs located in regions with frequent cloud cover and high-latitude reservoirs in the Northern Hemisphere, where low illumination has previously resulted in missing area values during winter months.

Calculation of drawdown areas

We calculated monthly drawdown areas for all reservoirs contained in GRSAD according to:

$${\rm{DA}}=\left({{\rm{Area}}}_{{\rm{\max }}}-{\rm{Area}}\right)/{{\rm{Area}}}_{{\rm{\max }}}$$
(1)

where DA is the relative extent of the drawdown area for a given reservoir considering the current monthly surface area (Area) and the maximum area recorded during the period 1985–2015 (Areamax). We assumed that the maximum area of each reservoir recorded during the 30-year period is a valid representation of its nominal surface area (the area of the reservoir at maximum filling level).

Complete filling of reservoirs was defined by a drawdown area smaller than 5% of Areamax. Because there is no uniform definition of ‘extreme drawdown’, we used the Cape Town water crisis 2018 as a reference21. The number of reservoirs experiencing extreme drawdown was estimated by averaging the number of reservoirs with drawdown areas exceeding 40%, 50%, 60% or 70% of Areamax at least once. To prevent initial filling of reservoirs being identified as extreme drawdown, 791 reservoirs built during the analysed period (year built ≥ 1985) were excluded from this analysis. The upper bound (70%) corresponds to the drawdown-area extent during the Cape Town water crisis 201821 (Fig. 1a). The lower bound (40%) corresponds to a reservoir capacity (storage water volume) of approximately 35%, as remained available during that water crisis, assuming an idealized, triangular reservoir shape (Extended Data Fig. 8). This was estimated according to:

$$0.36=\frac{{\left(0.6\times \sqrt{2}\right)}^{2}}{2}$$
(2)

For the calculation of total global drawdown area, used for the upscaling of GHG emissions, we combined data for reservoirs larger than 10 km2 with values derived from a Pareto model for smaller reservoirs. First, we estimated total reservoir surface area for nine size classes following a Pareto distribution. Subsequently, we estimated total drawdown area for each size class by multiplying the size-class-specific relative drawdown-area extent by the total reservoir surface area of each size class (Supplementary Table 3). Because the relative drawdown-area extent for reservoirs smaller than 0.001 km² is unknown and furthermore considered as being imprecise for reservoirs smaller than 10 km², we derived estimates for these size classes on the basis of four different statistical models (linear, square root, logarithmic, polynomial; Extended Data Fig. 9). Reservoirs larger than 10 km² were used to fit linear, square root and logarithmic models, whereas all available data were used for fitting a second-degree polynomial model to achieve a best representation of the available data. The four models all have a constant (linear model) or decreasing (square root, logarithmic, polynomial) slope. We have refrained from using models with increasing slopes (for example, exponential) to not overestimate the drawdown extent of small reservoirs and, thus, consider these estimates as conservative.

Data analysis

Statistical models to predict drawdown-area extent for each reservoir were developed using stepwise MLR. Climatic data (mean annual temperature, precipitation seasonality) for all reservoir locations were extracted from the Climatologies at High Resolution for the Earth’s Land Surface Areas climate dataset, which gives high-resolution (0.5 arcmin) climate information for global land areas over the period 1979–201332. Climate zones in the Köppen–Geiger system were determined from the high-resolution (5 arcmin) global climate map derived from long-term monthly precipitation and temperature time series representative for the period 1986–201033,34. Data on baseline water stress were extracted from Aqueduct 3.025. Baseline water stress measures the ratio of total water withdrawals to available renewable surface and groundwater supplies and is derived from high-resolution (5 arcmin) hydrological model outputs using the PCR-GLOBWB 2 model35,36.

Dates were categorized into four seasons on the basis of their meteorological definition depending on hemisphere. Therefore, for the Northern Hemisphere, spring begins on 1 March, summer on 1 June, autumn on 1 September and winter on 1 December. For the Southern Hemisphere, spring begins on 1 September, summer on 1 December, autumn on 1 March and winter on 1 June.

For the analyses of reservoir use types, we used the information provided in the column ‘MAIN_USE’ of the Global Reservoir and Dam dataset. Reservoirs where the main use was not specified (n = 1,554) were combined with those having MAIN_USE = ‘Other’ (n = 205).

To identify the magnitude of trends in time series, we used the non-parametric Theil–Sen estimator and the Mann–Kendall test because they do not require prior assumptions of statistical distribution for the data and are resistant to outliers. The Theil–Sen estimator was used to compute the linear rate of change, and the Mann–Kendall test was used to determine the level of significance. We analysed differences between groups using the Kruskal–Wallis test and Dunn’s post hoc test. The threshold to assess statistical significance was 0.05 for all analyses, The statistical analyses were performed using R 3.4.437.

Upscaling of GHG emissions and OC burial

Because the global reservoir area derived in this study differed from the area used in previous studies, we recalculated the published global estimates for both OC burial6 in and GHG emissions1 from reservoirs to allow for comparison (Extended Data Fig. 10). We fitted empirical distributions to CO2 emission data from drawdown areas (Supplementary Table 2 and Extended Data Fig. 7) as well as the published OC burial rates6 and published GHG emission data1 from water surfaces of reservoirs. For CO2 emissions from drawdown areas, we used a gamma distribution to account for non-normality of the data (Extended Data Fig. 7). For CO2 and N2O emissions from the water surface, we fitted a skewed normal distribution because of the occurrence of negative values (Extended Data Fig. 7). For CH4 emissions from the water surface, we fitted a log-normal distribution (Extended Data Fig. 7). Because the global estimate of OC burial was derived using geostatistical modelling, we fitted a gamma distribution to the published moments of OC burial rate6 (mean ± s.d. = 144 ± 75.83 gC m−2 yr−1) where the s.d. is calculated as the s.d. of the four scenarios used in that study. The final global empirical distributions for all fluxes were estimated by multiplying average emission and burial rates derived from resampling the preceding distributions times the total water surface area and drawdown area of reservoirs, resulting also from resampling their distributions after uncertainty propagation (see Treatment of uncertainty).

Treatment of uncertainty

As in all upscaling exercises, the global analysis conducted in this study is subject to substantial uncertainty. In our case, the uncertainty results from both the quantification of water surface and drawdown area of reservoirs and the estimation of global rates for GHG emission and OC burial. To comprehensively take all sources of uncertainty into account, we propagated all uncertainty throughout the whole analysis using a combination of Taylor series expansion and Monte Carlo simulations (Extended Data Fig. 10). In brief, we applied customary equations for uncertainty propagation derived from the Taylor series expansion method when propagating uncertainty of moments (for example, mean) or simple arithmetic calculations (for example, multiplication). For more-complex situations or when non-normality was conspicuous, we used Monte Carlo propagation. To obtain global estimates and standard error of water surface and drawdown area of reservoirs, both the systematic (bias) and random uncertainties of the remote-sensing-derived dataset18 as well as the uncertainty induced by our Pareto modelling for reservoirs <10 km² were propagated through all arithmetic calculations. These estimates and their uncertainty were again propagated to the calculation of global GHG flux and burial rates by combining them with rates and uncertainties derived from empirical distributions as described in the preceding. In a final step, Monte Carlo propagations were used to calculate the emission-to-burial ratios. All instances of Monte Carlo propagation were performed with 100,000 simulations by using the R package propagate38.