Modelling approach

In-reservoir OC transformations and burial rates are simulated using a simple, biogeochemical mass balance (or box) model (Supplementary Fig. 1). The parameters of the kinetic expressions representing the in-reservoir processes are assigned probability density functions (PDF) that take into account each parameter’s inter- and intra-reservoir variability at the annual timescale. The PDFs are derived from available data on photosynthetic C fixation, OC mineralization and OC burial in lentic systems. A similar approach was previously applied to predict global-scale phosphorus and nutrient silicon retention in dam reservoirs27,28.

A virtual database of model reservoirs is generated by performing 6,000 MC simulations with rate parameters selected randomly from their corresponding PDFs. From this virtual database, global relationships are derived that predict burial, mineralization and photosynthesis rates from the reservoir’s hydraulic residence time. The global relationships are applied to a real-world database that includes existing reservoirs (GRanD)9 in 1970 and 2000, and reservoirs under construction or planned to be completed by 2030 (ref. 10). Together with the DOC and POC loads to reservoirs obtained from the Global-NEWS model, the burial, mineralization and photosynthesis fluxes in the reservoirs of all major watersheds worldwide are then simulated for four selected time points (1970, 2000, 2030 and 2050). The MC analysis also yields estimates of the uncertainties on the OC fluxes associated with parameter variability.

The Global-NEWS model is well suited for the proposed modelling approach: it differentiates between POC and DOC, it implicitly accounts for in-stream OC losses, and it has been used to hindcast nutrient loads to watersheds in the years 1970 and 2000, and to forecast the 2030 and 2050 loads according to the MA scenarios22,35. Global-NEWS predicts sediment, OC and nutrient yields based on land use and land cover (for example, wetlands, cropland and grasslands), climate variables (including temperature and precipitation), geomorphological parameters (including slope and lithology) and anthropogenic alterations (including consumptive water usage). The MA scenarios are storylines for a future world that will become either more globalized (Global Orchestration (GO) and TechnoGarden, TG) or regionalized (Adapting Mosaic, AM, and Order from Strength, OS), and take either a proactive (TG and AM) or reactive approach to environmental management (GO and order from strength, OS)35. Each of the MA scenarios assumes changing land uses, climate regimes and anthropogenic perturbations, which in turn modify the fluxes of sediments, OC and nutrients delivered to river systems.

A caveat of Global-NEWS is the lack of representation of extreme hydrological events with a low recurring probability but potentially large contributions to long-term riverine fluxes. These events are not captured by the Global-NEWS model itself or by the observed data it was trained on. These events can involve geomorphologic processes such as landslides, which mobilize the carbon stored in the entire soil horizon plus the standing biomass. Such events, which can be related to extreme weather phenomena (for example, tropical cyclones) or earthquakes, have been described for steep catchments in the sub-tropics (for example, Taiwan43) and temperate regions (for example, New Zealand44). Their contributions to OC mobilization at the global scale remain to be quantified.

In our modelling approach, reservoir infilling due to sedimentation is not taken into account. Vörösmarty et al.45 estimate that sedimentation has reduced the volume of reservoirs in the United States by up to 2 km3, that is, only ∼0.2% of the total US reservoir volume. While reservoir infilling may vary significantly from one reservoir to another, the effect of sediment accumulation on water residence time likely represents a relatively minor source of uncertainty on the impact of dams on the global riverine OC flux.

Mass balance model

The mass balance model for in-reservoir allochthonous and autochthonous OC cycling is shown in Supplementary Fig. 1. Riverine inflow fluxes of allochthonous POC and DOC are calculated by multiplying inflow concentrations (mol km−3) by the river discharge (Q, km3 per year) through the dam. Assuming inflow of water to a reservoir equals the outflow through the dam, values for Q are retrieved from the GRanD database9, augmented with the database of hydroelectric dams of Zarfl et al.10 for the 2030 and 2050 projections. Outflow fluxes of OC through the dam are obtained by multiplying the in-reservoir POC and DOC masses (mol) by the flushing rate, ρ (per year), which is equal to Q/V where V is the reservoir volume in km3. The key role of ρ (or its inverse, the water residence time) in material mass balances of reservoirs is well-established6,24,27,46.

Except for primary production, the fluxes (F, mol per year) describing the in-reservoir OC processes assume first-order kinetics with respect to the mass (M, mol) of the OC pool from which the flux originates:

where k is the first-order rate constant (per year). Primary production of autochthonous POC (P, mol per year) is assumed to be phosphorus limited31,32,47, according to

where P max is the maximum value of P under nutrient saturated conditions, [TDP] the concentration of total dissolved phosphorus in the reservoir and K s the half-saturation TDP concentration. For each reservoir and time point, [TDP] is extracted from the previously published global reservoir phosphorus model28. The mass balance equations for the various OC pools are solved using the Runge–Kutta 4 integration scheme with 0.01-year time steps, and run for the length of time since dam closure (that is, if the dam is 10 years old, the model is run for 1,000 time steps).

Model parameterization

For each parameter value a PDF is derived from data reported in the literature (Supplementary Table 1). For the river inflow DOC concentration we impose the gamma distribution proposed by Sobek et al.48 based on data for 7,500 lakes. This is justified as the mean and range of DOC concentrations in rivers are comparable to those observed in lakes49,50. For the inflow POC concentrations, a generalized Pareto distribution yields the highest log-likelihood when fitted to the POC concentrations derived by applying the Global-NEWS loads and discharge values to the GRanD reservoirs9,22. Similarly, the PDFs for reservoir volume, discharge, latitude, altitude and age (that is, years since dam closure) are the functions that were found to best fit the observed distributions in the GRanD database (Supplementary Table 1). Reservoirs in GRanD further yield the following relationship between surface area (SA, km2) and volume:

In the MC analysis, SA values from equation (3) are multiplied by a random number between 0.1 and 10, thus allowing deviations in SA by±one order of magnitude from equation (3). We compute SA from V, rather than vice versa, because V is the primary size variable used to calculate the water residence times in the up-scaling procedure.

The temperature dependence of the allochthonous POC and DOC mineralization (biological respiration and photochemical degradation) rate constants, k min,allo , is assumed to follow the same expression as that derived by Hanson et al.29,30 from a compilation of ecosystem-scale mineralization rates of allochthonous DOC in lakes:

where T is the average water temperature (°C), k 20 is the rate constant at 20 °C, and θ is a dimensionless coefficient equal to 1.07 (ref. 29). The PDF of k 20 in Supplementary Table 1 is directly derived from the data of Hanson et al.29,30. While the same k 20 PDF is applied to both allochthonous DOC and POC29, the k 20 values for POC and DOC are allowed to vary independently in the MC analysis. This approach allows for reactivities of allochthonous POC that in some cases are higher, in others lower, than those of allochthonous DOC, as has been reported for lentic systems51,52.

In-reservoir produced autochthonous OC tends to be more labile than allochthonous OC53,54. The k 20 value for autochthonous OC is therefore multiplied by a variable scaling factor. An average three-times higher autochthonous k 20 yields the best fit of our model-predicted total carbon mineralization rate constants to the rate constants compiled by Catalan et al.17 across a wide range of aquatic ecosystems and climate zones. To capture the range of the Catalan et al.17 rate constant data (with the exception of four outliers, three of which are from the same small streams), the autochthonous scaling factor is assumed to follow a normal distribution between 1 and 6 (Supplementary Table 1). With this range of the scaling factor, the model-generated total OC mineralization rate constants reproduce the range of observed rate constants of Catalan et al.17 (Supplementary Fig. 7). Because of its greater lability, autochthonous OC exported through a dam is assumed to be mineralized before reaching a downstream reservoir.

The reservoir water temperature T is computed as a function of latitude and altitude following Hart and Rayner55. For 2030 and 2050, the average air temperature increases predicted by Fekete et al.33 for each individual MA scenario are converted to reservoir water temperature increases using the relationship of Lauerwald et al.20, and added to the temperature predicted with the Hart and Rayner equation (Supplementary Table 6).

The half-saturation constant K s in equation (2) is represented by a uniform PDF ranging from 2.0 × 107 to 7.0 × 108 mol km−3 in units of total dissolved phosphorus56,57,58,59,60,61,62,63,64. The maximum primary production P max in equation (2) is estimated with:

where B is the average annual depth-integrated chlorophyll concentration in the reservoir (mg Chl-a km−3), P Chl is the maximum, annually integrated, chlorophyll-specific carbon fixation rate (mol C (mg Chl-a)−1 per year), M is a dimensionless metabolic correction factor for water temperature and D f is the yearly proportion of ice-free days. The value of P Chl is fixed at 0.15 mol C (mg Chl-a)−1 per year, which is equivalent to 2.5 g C (g Chl-a)−1 h−1, assuming an annual average of 12 h daylight per 24 h (refs 65, 66, 67, 68); M is equal to 1 at water temperatures ≥28 °C, and decreases at lower temperatures with a Q 10 of 2 (ref. 69); the duration of ice cover required to estimate D f is calculated from latitude following Williams et al.70; B is updated in each model iteration using the relationship provided by Reynolds71:

where k c is the absorbance of photosynthetically active radiation (PAR) per unit of chlorophyll (fixed at 0.014 m2 (mg Chl-a)−1 globally69), PP/RP is the ratio of maximum gross photosynthesis to algal respiration per unit chlorophyll, fixed at 15 (ref. 72), DL is the hours of daylight, fixed annually at 12 h per day71, I o,max is the maximum site-specific PAR (μmol m−2 s−1), I k is the PAR at the onset of photosaturation, fixed at 120 μmol m−2 s−1)71, Z mix is the reservoir mixing depth (m) and K dw +K dp +K dg is the nonalgal PAR attenuation (m−1). The value of I o,max is calculated for each reservoir based on the annually averaged latitude-specific values provided by Lewis69 assuming linear interpolation between the latitudinal bands provided; Z mix is estimated based on the empirical relationship with lake fetch69, where fetch is assumed equal to the diameter of a circle with the same area as the reservoir; PAR attenuation in pure water, K dw , is fixed at 0.13 m−1, PAR attenuation by inorganic suspended particulate matter (tripton), K dp , is fixed at 0.06 m−1 , and that by dissolved organic matter (gilvin), K dg , is calculated for each reservoir with the relationship proposed by Lewis69:

where [DOC] in p.p.m. is generated using the same PDF as for inflow [DOC]48. A t-test shows that B-values generated in the MC analysis are statistically indistinguishable (P<0.05) from values found in the literature65,69,71.

The burial rate constant, k bur , is an effective parameter describing the long-term retention of POC with the sediments accumulating in the reservoir, that is, the POC that is not remineralized or otherwise remobilized and exported over the reservoir’s lifetime. The value of k bur aggregates all the factors controlling the POC burial efficiency other than the water residence time, including sedimentation rate, reservoir morphology (which controls deposition patterns), oxygen exposure time, temperature and sediment resuspension events53. The upper and lower bounds for the uniform k bur distribution, 1 and 15 per year respectively, are based on published rate constants73 and values back-calculated from global and regional data sets of OC accumulation rates in lakes and reservoirs38,53,74. More studies quantifying burial rate constants in a diversity of reservoir settings will be needed to further delineate the form of the k bur PDF. The rate constant for allochthonous POC solubilization to DOC is kept fixed at 0.1 per year, as the model is insensitive to this parameter (see Model sensitivity and uncertainty section)75.

Global upscaling

Global regression relationships are obtained from the MC analysis by carrying out 6,000 iterations with the 1970 and 2000 parameter values each, plus 6,000 additional iterations for each of the MA scenarios in 2030 and 2050 taking into account the projected changes in air temperature. For all time points and MA scenarios considered, the water residence time, τ r , provides the best predictor for the fate of allochthonous OC entering a reservoir. The in-reservoir burial and mineralization fluxes, F i,bur and F i,min, produced by the MC simulations are fitted to

and

where the subscript i stands for allochthonous POC or DOC (POC only in the case of burial), in stands for inflow, T is average annual reservoir water temperature (°C), α i and β i are first-order rate coefficients describing loss due to burial and mineralization, respectively, and a i and b i are dimensionless coefficients. The best-fit values of α i , β i , a i and b i are given in Supplementary Table 2. Note that equations (8) and (9) are formally similar to the equation derived by Vollenweider46 to describe phosphorus cycling in lakes, though scaled with a i and b i to account for the separation of loss fluxes into mineralization and burial, rather than lumped together as in Vollenweider’s derivation.

The inflows of allochthonous DOC and POC, F i,in , are estimated by spatially overlaying the GRanD reservoirs9 onto the Global-NEWS watersheds. Global-NEWS is used because it generates spatially explicit riverine OC and phosphorus yields for the four MA scenarios. The projected MA global average air temperature increases are 0.91–1.09 °C in 2030 and 1.29–2.11 °C in 2050, relative to 2000. These projections are similar to the temperature increases of the Representative Concentration Pathways scenarios RCP4.0 and RCP6.5 (ref. 76). Reservoir water temperature, T, is calculated for each reservoir in GRanD using the relationship with latitude and altitude derived by Hart and Rayner55. For the 2030 and 2050 scenarios, an increase in water temperature is added to each reservoir based on the average latitude-specific temperature increases predicted by Fekete et al.33.

All reservoirs are spatially routed to account for dams that occur longitudinally in series, and dams on tributaries as illustrated in Supplementary Fig. 2. The total influx (mol per year) of allochthonous POC or DOC to a reservoir k is then given by:

where is the sum of the fluxes of allochthonous POC or DOC leaving all dams immediately upstream of reservoir k on tributaries 1 to n, W k is the total upstream watershed area (km2), is the sum of the upstream watershed areas of dams k-1, and Y k is the POC or DOC yield (mol km−2 per year) in the undammed, , watershed area downstream of dams k-1. In other words, the OC loads leaving all upstream reservoirs are added to the OC load entering the river from its undammed upstream watershed area.

In Global-NEWS, carbon and nutrient loads from landscapes to rivers are empirically adjusted to match the fluvial DOC and POC export fluxes measured at the mouths of large rivers. To reconcile our estimates of OC elimination in reservoirs with the observed OC export fluxes to the coastal ocean, the Global-NEWS allochthonous OC loads to rivers are recalibrated. As initial hypothesis we impose the original Global-NEWS DOC and POC loads to calculate in-reservoir elimination fluxes and export fluxes to the coastal ocean. The loads are then iteratively adjusted until they reproduce the Global-NEWS coastal OC export fluxes for 1970 and 2000. For the 2030 and 2050 scenarios, which are based on projections rather than calibrated to data, we apply the same correction factor to any given watershed as that for 2000. On average, the revised OC loads to watersheds differ only by ±1.3% from the original Global-NEWS estimates.

The combined outputs of the MC analyses of reservoir OC and phosphorus models yield the following relationship between in-reservoir primary production, P, and the inflow flux of total dissolved phosphorus, TDP in (mol per year):

where the uncertainty associated with each parameter corresponds to ±1 s.d. (see Model sensitivity and uncertainty section). Equation (11) is then used in relationships of the form of equations (8) and (9) to describe the fractions of in-reservoir produced autochthonous OC that are buried and mineralized:

and

where the subscript j stands for burial and mineralization of autochthonous OC. The values of the parameters α j , β j , a j and b j are listed in Supplementary Table 2.

The initial amount of degradable soil and biomass OC flooded after closure of a dam, TOC 0 , is estimated from global soil77 and biomass carbon78 maps, scaled to the surface area of the reservoir. The decay of the flooded OC follows:

where t is the number of years since dam closure. We assume the same temperature function (that is, θ(T–20)) for k min,flood as in equation (4) and adjust the pre-function coefficient k 20 so as to reproduce the decay timescale of flooded soil organic matter typically observed after dam closure (10–15 years)16,79. The resulting expression is then:

The Strahler stream orders of the GRanD reservoirs (Supplementary Table 3) are estimated with the empirical scaling law relating stream order to catchment area derived by Lauerwald et al.20. This scaling law is derived from the HydroSHEDS stream network for third order stream and higher80, and extrapolated for the two lowest stream orders that are not represented in HydroSHEDS20. To test the validity of the scaling law, we have compared calculated stream orders with actual stream orders from the HydroSHEDS stream network for all dammed streams and rivers in Europe (n=2,192) and South America (n=1,602), yielding statistically significant R2 values of 0.90 and 0.87, respectively.

Model sensitivity and uncertainty

Sensitivity analyses are performed separately for in-reservoir allochthonous and autochthonous OC cycling. The effects on burial and mineralization fluxes of changing one parameter value at the time are summarized in Supplementary Tables 4 and 5. In most cases, the parameter is varied ±10% from the assigned ‘default’ value in 100 iterations, the results of which are compared to the burial and mineralization fluxes obtained using the default value. Reservoir volume and discharge (and hence water residence time) are varied according to the distributions in Supplementary Table 1. Sensitivity to initial conditions is assessed by comparing the results of model runs with the initial reservoir DOC, POC and TOC masses set equal to either 0 or 1 × 106 mol. Two reservoir ages are tested: 10 and 40 years. Sensitivity to a 1 °C increase in global air temperature is determined, which is equivalent to a water temperature increase of 0.82 °C. Not surprisingly, allochthonous and autochthonous OC burial fluxes are most sensitive to the rate constant of burial k bur while mineralization fluxes are most sensitive to k 20 and latitude, given the dependence of k min on temperature, which is in turn related to latitude, and to a lesser degree, altitude (Supplementary Tables 4 and 5).

A bootstrap analysis is used to estimate the uncertainties associated with primary productivity P calculated with equation (11): sampling with replacement was conducted drawing 5,000 samples from the 6,000 model runs generated in the MC analysis, and fitted to a power law as in equation (11). After 5,000 iterations, the bootstrap analysis yields the s.d. estimates for each parameter in equation (11). When scaled up globally, the uncertainty on P translates into uncertainties in the mineralization and burial fluxes of autochthonous OC of ±15%.

The total uncertainties in the global reservoir OC burial and mineralization are assessed as follows. Burial and mineralization fluxes predicted by the 6,000 MC iterations are binned according to water residence times as shown in Supplementary Fig. 3. Gamma functions are fitted to the distributions of fluxes in each bin. A second MC analysis is then carried out in which, for each GRanD reservoir, the burial and mineralization fluxes are randomly selected from the gamma functions. A total of 20 simulations yields ±8% s.d. on the global OC fluxes for allochthonous DOC mineralization, ±15% for allochthonous POC mineralization and ±12% for allochthonous burial. For autochthonous OC fluxes, this analysis yields ±2% and ±3% s.d. for mineralization and burial fluxes, respectively. Combined with the ±15% error associated with the Global-NEWS OC loads to rivers, we estimate uncertainties on the order of ±23%, ±30% and ±27%, for allochthonous DOC mineralization, and allochthonous POC mineralization and burial, respectively. For the global values of autochthonous OC burial and mineralization, we estimate an uncertainty on the global fluxes of ±32% for mineralization and ±33% for burial. The higher uncertainties for autochthonous OC reflect the uncertainty in estimates of P (see above). Note that the estimated uncertainties in the burial and mineralization fluxes are those associated with the mass balance modelling approach, and do not account for inaccuracies and omissions in the GRanD or future dam databases.

Data availability

Reservoir burial and mineralization fluxes of OC for individual river basins, in 1970, 2000 and 2030 (GO scenario), are given in Supplementary Data 1. For access to computer code, please contact Taylor Maavara.