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

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¹⁸. It comprises all reservoirs from the Global Reservoir and Dam dataset¹⁹ except of 45 reservoirs without reported geometric information. In accordance with ref. ¹, 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³¹ 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²¹. 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²¹ (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² 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³². 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²⁵. 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³⁷.

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⁶ in and GHG emissions¹ from reservoirs to allow for comparison (Extended Data Fig. 10). We fitted empirical distributions to CO₂ emission data from drawdown areas (Supplementary Table 2 and Extended Data Fig. 7) as well as the published OC burial rates⁶ and published GHG emission data¹ from water surfaces of reservoirs. For CO₂ emissions from drawdown areas, we used a gamma distribution to account for non-normality of the data (Extended Data Fig. 7). For CO₂ and N₂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₄ 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⁶ (mean ± s.d. = 144 ± 75.83 gC m⁻² yr⁻¹) 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¹⁸ 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³⁸.

Source: Resources - nature.com