Global carbon budget of reservoirs is overturned by the quantification of drawdown areas

Reservoir drawdown areas—where sediment is exposed to the atmosphere due to water-level fluctuations—are hotspots for carbon dioxide (CO2) emissions. However, the global extent of drawdown areas is unknown, precluding an accurate assessment of the carbon budget of reservoirs. Here we show, on the basis of satellite observations of 6,794 reservoirs between 1985 and 2015, that 15% of the global reservoir area was dry. Exposure of drawdown areas was most pronounced in reservoirs close to the tropics and shows a complex dependence on climatic (precipitation, temperature) and anthropogenic (water use) drivers. We re-assessed the global carbon emissions from reservoirs by apportioning CO2 and methane emissions to water surfaces and drawdown areas using published areal emission rates. The new estimate assigns 26.2 (15–40) (95% confidence interval) TgCO2-C yr−1 to drawdown areas, and increases current global CO2 emissions from reservoirs by 53% (60.3 (43.2–79.5) TgCO2-C yr−1). Taking into account drawdown areas, the ratio between carbon emissions and carbon burial in sediments is 2.02 (1.04–4.26). This suggests that reservoirs emit more carbon than they bury, challenging the current understanding that reservoirs are net carbon sinks. Thus, consideration of drawdown areas overturns our conception of the role of reservoirs in the carbon cycle. Globally, reservoirs are net emitters of carbon when drawdown areas are taken into account, according to an analysis of satellite observations of reservoir surface area.

R eservoirs, 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 (CO 2 e) to the atmosphere 1 . The exact quantification of this source is also important because GHG emissions from reservoirs may affect the carbon (C) footprint of irrigation 2 and hydropower production and thus its perception as C-neutral energy [3][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 lakes 6 . Globally, rates of organic carbon (OC) burial in reservoirs have been estimated to exceed rates of C emissions from reservoirs 1,6 .
Water-level changes are a typical feature of reservoirs and a major discriminator to natural lakes 7 . They are driven both by natural hydrological dynamics and water management. Water-level fluctuations have several consequences for reservoir limnology and biogeochemical cycling 8,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 acknowledged 1 , they have not been considered in reservoir C budgets 10 and global upscalings 1,6 so far.
It has been shown that dry aquatic sediments in general, and drawdown areas of reservoirs in particular, emit large quantities of GHGs 11,12 . GHG emissions from drawdown areas are typically dominated by CO 2 while CH 4 is of minor importance 13 . Initial studies reported that the drawdown area has the potential to dominate annual CO 2 emissions from reservoirs 12,14-16 due to areal CO 2 emissions, which are significantly higher than those from the water surface 11 . CO 2 emissions from dry inland waters were shown to be controlled mainly by temperature, moisture and organic-matter content of the desiccated sediments 11 . Extreme water-level drawdown can furthermore create hot moments of GHG emissions and may even offset C burial in the sediments 16,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 worldwide 18 . 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 Global carbon budget of reservoirs is overturned by the quantification of drawdown areas Philipp S. Keller 1 ✉ , Rafael Marcé 2,3 , Biel Obrador 4

and Matthias Koschorreck 1
Reservoir drawdown areas-where sediment is exposed to the atmosphere due to water-level fluctuations-are hotspots for carbon dioxide (CO 2 ) emissions. However, the global extent of drawdown areas is unknown, precluding an accurate assessment of the carbon budget of reservoirs. Here we show, on the basis of satellite observations of 6,794 reservoirs between 1985 and 2015, that 15% of the global reservoir area was dry. Exposure of drawdown areas was most pronounced in reservoirs close to the tropics and shows a complex dependence on climatic (precipitation, temperature) and anthropogenic (water use) drivers. We re-assessed the global carbon emissions from reservoirs by apportioning CO 2 and methane emissions to water surfaces and drawdown areas using published areal emission rates. The new estimate assigns 26.2 (15-40) (95% confidence interval) TgCO 2 -C yr −1 to drawdown areas, and increases current global CO 2 emissions from reservoirs by 53% (60.3 (43.2-79.5) TgCO 2 -C yr −1 ). Taking into account drawdown areas, the ratio between carbon emissions and carbon burial in sediments is 2.02 (1.04-4.26). This suggests that reservoirs emit more carbon than they bury, challenging the current understanding that reservoirs are net carbon sinks. Thus, consideration of drawdown areas overturns our conception of the role of reservoirs in the carbon cycle.
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 levels 20 . 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) TgCO 2 e 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 filled 1,10 , whereas dry aquatic sediments can emit large quantities of CO 2 11,12 . We re-assessed the most recent estimate of C emissions from reservoirs 1 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 Table 1 and Fig. 3).
Total reservoir C emissions are dominated by CH 4 in terms of climate forcing 1 . CH 4 emissions from drawdown areas, however, have been shown to be low compared with CO 2 emissions, but data on CH 4 emissions from drawdown areas are scarce 12 . Assuming zero CH 4 emissions from drawdown areas would reduce current global CH 4 emissions from reservoirs (12.7 (6.8-23.2) TgCH 4 -C yr -1 (ref. 1 )) by 1.8 TgCH 4 -C yr −1 (Table 1 and Fig. 3). However, when expressed in terms of CO 2 e of the three main GHGs emitted from reservoirs (that is, considering CO 2 , CH 4 and N 2 O, assuming a global warming potential for CH 4 of 34 (ref. 22 ) and for N 2 O of 298 (ref. 22 )), no substantial change of total CO 2 e emissions from reservoirs can be detected given the uncertainty of the calculation  Cohen's D effect size = 0.9; Table 1 and Fig. 3). CO 2 e emissions are dominated by CH 4 for drawdown areas smaller than ~30% of total reservoir surface. Because this finding is based on assuming zero CH 4 emissions from drawdown areas, the dominant role of CH 4 will be even increased if future studies reveal substantial CH 4 emissions linked to drawdown areas, potentially leading to a total increase in CO 2 e emissions. In fact, exposure of hypolimnetic sediments at extreme drawdown may be a substantial methane source at least in the short term 16,17 . Furthermore, substantial decrease of water levels in reservoirs is expected to increase CH 4 emissions via ebullition due to lowering of the hydrostatic pressure 23,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 demand 25 . 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 CO 2 emissions from dry inland waters among climate zones, the variability of CO 2 fluxes across locations has been shown to be substantial 11 . Thus, data with better spatial coverage are required to allow for a more comprehensive geostatistical upscaling considering regional differences in both CO 2 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 fluxes 26 . It was, for example, shown that CO 2 fluxes from dry sediments increase with temperature 11 . 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 CO 2 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 Green ribbon indicates the potential range of buried CO 2 e assuming conversion of OC between 100% CO 2 (lower bound) and 50% both CO 2 and CH 4 (upper bound). a,b, CO 2 and CH 4 emission: emissions from reservoirs comprising both drawdown area and water surface. Total emission: sum of CO 2 and CH 4 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, CO 2 -C emission from global reservoirs (Cohen's D effect size = 2.6). d, GHG emission expressed as CO 2 e 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 CO 2 -C and CH 4 -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).
known that the trophic state of reservoirs affects their water-bound GHG emissions 27 . 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 vegetation 12 , a typical feature of many reservoirs. The fate of the C temporarily stored in this vegetation upon reflooding (burial versus mineralization to CO 2 and/or CH 4 ) 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 lakes 6 . Total OC burial in reservoirs is estimated to be 41.7 (21.6-73.0) TgC yr −1 by multiplying area-specific burial rates 6 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 emissions 6 .

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 summer 28 . Extreme drawdown events during summer also pose a threat to fish communities and water quality in general 29,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.

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/ s41561-021-00734-z.

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 GRSAD 18 . It comprises all reservoirs from the Global Reservoir and Dam dataset 19 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 dataset 31 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: 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 (Area max ). 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 Area max . Because there is no uniform definition of 'extreme drawdown' , we used the Cape Town water crisis 2018 as a reference 21 . 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 Area max 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 2018 21 (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: For the calculation of total global drawdown area, used for the upscaling of GHG emissions, we combined data for reservoirs larger than 10 km 2 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-2013 32 . 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-2010 33,34 . Data on baseline water stress were extracted from Aqueduct 3.0 25 . 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 model 35,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.4 37 .

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 burial 6 in and GHG emissions 1 from reservoirs to allow for comparison (Extended Data Fig. 10). We fitted empirical distributions to CO 2 emission data from drawdown areas (Supplementary Table 2 and Extended Data Fig. 7) as well as the published OC burial rates 6 and published GHG emission data 1 from water surfaces of reservoirs. For CO 2 emissions from drawdown areas, we used a gamma distribution to account for non-normality of the data (Extended Data Fig. 7). For CO 2 and N 2 O emissions from the water surface, we fitted a skewed normal distribution because of the occurrence of negative values (Extended Data Fig. 7). For CH 4 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 rate 6 (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 dataset 18 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 propagate 38 .

Data availability
The GRSAD dataset containing the surface-area time series of 6,818 global reservoirs is available online at https://doi.org/10.18738/T8/DF80WG. The CHELSA climate dataset, which gives high-resolution climate information for global land areas is available online at https://chelsa-climate.org/. Data on climate zones according to the Köppen-Geiger classification are available at http:// koeppen-geiger.vu-wien.ac.at/. Data on baseline water stress extracted from Aqueduct 3.0 are available at https://www.wri.org/publication/aqueduct-30. Source data are provided with this paper.

Code availability
Scripts used for the analyses described in this study can be obtained from the corresponding author upon request.  Table 1. a-c, histograms of modelled continuous variables. c, Please note that maximum reservoir area is presented on log 10 scale. d, normal Q-Q plot of model's residuals. Fig. 9 | Modelling of relationship between drawdown area and reservoir area. a linear model (LM), b Square root model (Sqrt), c Logarithmic model (Log), d Polynomial model (Poly). a-d Black and grey dots show ratios as observed in the GRSAD dataset. Black dots were used for model fitting, red dots show extrapolated values. Reservoirs larger 10 km² were used to fit models a-c, as coverage of smaller reservoirs is considered to be increasingly incomplete in the underlying data set. In contrast, all available data was used for fitting the polynomial model (d) in order to achieve a hypothetical best fit.