Past anthropogenic activities offset dissolved inorganic phosphorus retention in the Mississippi River basin

The rapid acceleration of anthropogenic phosphorus (P) loadings to watersheds has fuelled massive freshwater and coastal eutrophication and completely changed the global P cycle. Within watersheds, emitted P is transported downstream towards estuaries. Reservoirs can retain a significant proportion of this P. In the long term, this accumulated P can however be re-mobilized, a process lacking in current global P budgets. Here, we include P cycling in a coupled integrated assessment-hydrology-biogeochemistry framework with 0.5 by 0.5-degree spatial resolution and an annual time resolution, and apply it to the Mississippi River basin (MRB). We show that, while reservoirs have aided in the net retention of P, they serve as dissolved inorganic P (DIP) sources due to the transformation of legacy P in sediments. The increasing DIP sourcing in the MRB has been offsetting P retention in streams, especially towards the end of the twentieth century. Due to its bioavailability, DIP is the most likely form to trigger eutrophication. Although P inputs into the MRB have decreased since the 1970s, legacy effects are delaying positive outcomes of remediation measures.


Introduction
Phosphorus (P) is a critical nutrient for all living species, and its application in croplands has become essential to sustain global food production (Smil 2000). From an agricultural perspective, P is a valuable resource that is becoming increasingly scarce, with major deposits located in a handful of countries (Elser and Bennett 2011). Since the beginning of the industrial era, the intensification of agriculture and extensive use of fertilizer, together with fast urbanization and increase in wastewater effluent, have drastically amplified and accelerated Abstract The rapid acceleration of anthropogenic phosphorus (P) loadings to watersheds has fuelled massive freshwater and coastal eutrophication and completely changed the global P cycle. Within watersheds, emitted P is transported downstream towards estuaries. Reservoirs can retain a significant proportion of this P. In the long term, this accumulated P can however be re-mobilized, a process lacking in current global P budgets. Here, we include P cycling in a coupled integrated assessment-hydrology-biogeochemistry framework with 0.5 by 0.5-degree spatial P transport to global surface freshwaters. From an environmental perspective, these additional P inputs to watersheds have subsequently increased both the P export to estuaries, and its retention in river networks globally (Beusen et al. 2016). Excess P in the hydrosphere can lead to widespread eutrophication and increase the frequency of hypoxia events (Rabalais et al. 2014;Jenny et al. 2016), constituting threats to human and ecosystem health.
Nevertheless, not all P forms have the same effects on aquatic environments (Mainstone and Parr 2002;Baker et al. 2014). Dissolved inorganic phosphorus (DIP, e.g., orthophosphate) is the most readily available P form for biota, but it only constitutes a small fraction of the overall P inputs into the hydrosphere. A major scientific challenge for sustainable P management thus resides in quantifying the transformations that various P forms undergo as they move toward the mouth of the river. This includes pinpointing major accumulation points and the stock change of each P species in the face of this net anthropogenic acceleration of the global P cycle. This knowledge can eventually aid in the development of strategies to minimize the presence of P in the freshwater environment, and negative effects on its ecosystems. P point sources can generally lead to high risk of eutrophication, due to bioavailability (mainly DIP) and potential dominance at low flow, in periods of high biological activity (Jarvie et al. 2006). Interactions of DIP with river sediments in streams can mask the effect of these sources (Jarvie et al. 2006(Jarvie et al. , 2012. This retained DIP may however be re-mobilized, for instance due to reduction of inputs, or constitute legacy pools in downstream low-flowing systems, i.e., lakes and reservoirs. Growing evidence suggests that P accumulated in the landscape [legacy P (Haygarth et al. 2014)] can be re-mobilized on the long term, thereby potentially masking the results of reduction measures (Meals et al. 2010;Sharpley et al. 2013). Reservoirs, accounting for an estimated 22% of P retention in river networks worldwide (Beusen et al. 2016), cannot indefinitely accumulate nutrients within bottom sediments (Powers et al. 2015). Several studies have shown that reservoirs can constitute net P sources, due to legacy effects (e.g. Teodoru and Wehrli 2005;Powers et al. 2015). P accumulated in bottom sediments can for instance be desorbed from bed inorganic particles or released via the mineralization of accumulated organic matter (Søndergaard et al. 2003).
Spatially explicit, dynamic models are needed to better understand all the processes involved in the development of P legacy in a river basin, as both long-term (climate change, land use changes, longterm P budgets in agriculture, accumulation in lakes and reservoirs) and shorter-term changes (new dam construction, changes in agricultural P runoff and erosion) need to be accounted for. In order to further understand long-term human-induced effects on the in-stream P cycle, we apply a new process-based coupled hydrological-biogeochemical model [IMAGE-DGNM (Vilmin et al. 2020)] to a major watershed, the Mississippi River basin (MRB). For this study, a new representation of P cycling was implemented into IMAGE-DGNM's biogeochemical module (Dynamic Instream Chemistry, a.k.a. DISC module) to quantify the P cycle within the MRB from 1970 to 2000, based on biogeochemical transformations throughout the watershed.
The MRB drains ca 40% of the continental USA into the Gulf of Mexico and has been subjected to a number of human pressures, including reservoir building, increased agricultural land use, and increased urban population. Currently, P inputs mainly originate from agricultural lands (USEPA 2007;Jacobson et al. 2011). Nutrient export from the MRB enhances primary production, causing oxygen depletion in the lower water column of the northern Gulf of Mexico, and leading to the largest "dead zone" of the entire western Atlantic Ocean (Rabalais et al. 2002;Diaz and Rosenberg 2008). P is often the limiting nutrient for primary production in freshwaters, and was shown to have an important control over phytoplankton development during spring and summer in the lower salinity, near-shore regions of the Gulf of Mexico, thereby affecting the size of the hypoxic zone (USEPA 2007). Water, sediment, and nutrient transfers to the Gulf of Mexico have also been deeply modified by reservoir trapping, with nowadays thousands of dams dispersed throughout the basin (Lehner et al. 2011), serving as temporary storage for large amounts of P (Beusen et al. 2016).
However, DIP retention in these low-flowing systems is most likely lower than that of particulate forms, which could modify downstream P speciation. Our spatially-explicit simulation of P cycling in the MRB, including re-mobilization processes, allows for the investigation of the controlling processes of P transfers and speciation throughout the basin, notably legacy effects. We focus on the impact of increasing human infrastructure on DIP transfers throughout the river network, and the export of this most bioavailable P form to the estuary.

Methodology
P cycling within the IMAGE-DGNM framework IMAGE-DGNM allows for the dynamic simulation of biogeochemical cycling in river networks at a half-degree spatial resolution (Vilmin et al. 2020 (Stehfest et al. 2014)], simulating interactions between the Human and Earth systems, and a modular, process-based representation of sediment dynamics and biogeochemical cycling [Dynamic In-Stream Chemistry module -DISC (Vilmin et al. 2020)].
PCR-GLOBWB simulates water fluxes and volumes in streams of Strahler order ≥ 6 (herein referred to as "large" or "high order" streams) and connected lakes and reservoirs. Smaller streams are characterized using the formulations from Wollheim et al. (2008), who defined the characteristics of the subgrid river network for the non-lake portion of each cell using accepted geomorphic principles. The number, drainage area, and stream length are calculated for each stream order 2-6 on the basis of specified characteristics of first order streams, whereby geomorphic parameters are assumed to be globally uniform (for details see Supplementary Information S.1.1, Table  S.1). In the present study, yearly derived water fluxes and volumes are used.
The IMAGE-DGNM model describes long-term (twentieth century) changes in land use, agricultural production for different crops and grasslands, P budgets and P runoff from agricultural land (Bouwman et al. 2017) and land under natural vegetation, P released from weathering, aquaculture and in-stream P inputs from vegetation in flooded areas at a yearly time step (Beusen et al. 2015). The delivery of DIP, PIP and OPDET to the river network is estimated with the approach presented by Vilmin et al. (2018), who defined speciation (based on literature) for each simulated nutrient delivery source. Furthermore, all runoff water reaching the river network is assumed to contain 1 µgChla L −1 , converted to OPPP using a C:Chla ratio of 0.035 mgC µgChla −1 and the Redfield C:P molar ratio (Garnier et al. 1998). Assessment of TSS inputs, simulation of sediment dynamics, and advective transfers are described in detail by Vilmin et al. (2020).
The new P-cycling implementation of DISC simulates four P forms: DIP, particulate inorganic phosphorus (PIP), and two forms of organic phosphorus, i.e. detrital OP (OP DET ), and OP constituting primary producers (OP PP ) (Fig. 1a). We keep track of PIP and OP DET accumulated in the benthic layer (denoted PIP b and OP DETb , respectively) for the calculation of DIP sorption to bed sediments, or release by benthic biogeochemical processes. Total suspended sediments (TSS) and total bottom sediments (TBS) are also computed, as they can have a significant impact on P processing (i.e., influence on sorption equilibrium and light penetration). The model output time step is 1 year, but processes are solved implicitly using an adaptive time step (≤ 0.5 year). DISC represents processes linked to primary production in the water column, as well as organic matter degradation (i.e., OP mineralization) and sorption interactions both in the water column and in the benthic layer (Fig. 1a). Particulate species can settle in the benthic layer and be re-suspended by the flow. Over a certain threshold of accumulated bed sediments, compaction occurs, and part of benthic particulate species become unavailable for erosion or biogeochemical transformations. The fraction of DIP transferred through the sediment-water interface due to benthic dynamics (f DIFF [−]) is described with an empirical formalism defined by Billen et al. (2015). The change rates of the simulated P species due to biogeochemical transformations and sediment-water exchanges are expressed as: All notations are provided in Table 1. Details on the equations and parameters used to calculate the different process rates can be found in the Supplementary Information (S.1.2, Tables S.2 and S.3). No model calibration was performed. To represent biogeochemical processes as realistically as possible, we used parameters issued from the literature that are relevant to the MRB conditions (Froelich 1988;Billen et al. 1994;Garnier et al. 2000;Chao et al. 2010).

Application to the MRB
The model is applied to the MRB, 4th largest river basin worldwide (in terms of surface area), over the  Net flux of particulate species i from the benthic layer due to sedimentation and erosion processes twentieth century, assuming equilibrium in 1900. Prior to 1900, fertilizer use in the MRB was negligible (Bouwman et al. 2013) and the number of dams and their volume was small and only located in upstream areas (Lehner et al. 2011). Results are analyzed for the period 1970-2000, over which extensive validation data is available (USGS 2013). According to a review from Jeppesen et al. (2005), internal P loading from internal sources can delay eutrophication mitigation strategies by 10-15 years. The effect of any legacy P from before 1900 should therefore not impact results for the study period. The model's accuracy for the simulation of yearly sediment and P dynamics is assessed by comparing its results to measured concentrations at 11 United States Geological Survey (USGS) monitoring stations for the period 1970-2000 (Fig. 1b, see Table S.4 for description of the stations). These stations are representative 9 subbasins of the MRB, which is documented in detail in USGS (2007). Observed annual concentrations of TSS, TP and DIP (orthophosphate measurements), weighted with the reported discharge, are calculated and compared to model results: i) visually ( Fig. 2 and Figs. S.1-S.4), and ii) statistically (Table S.5). For reasons of representativity, results are compared only for years with more than 6 measurement dates.
The model is then used to investigate changes in P export to the estuary over these three decades, as well as to identify the controlling biogeochemical processes of P retention in the MRB and the role of different waterbody types, notably reservoirs. Therefore, the export and retention are aggregated over the whole basin, or for these different systems only. Retention of each P form is calculated for every waterbody type as the difference between total input of the concerned form and its total export from the system (expressed in GgP, or in % of the total input). This retention is the net result of water residence times and in-stream biogeochemical reactions. Since the model is run from 1900, pre-1970 accumulated P loads may be contributing to the P release at the end of the twentieth century. Note that we focus here only on in-stream retention mechanisms, nutrient retention in other compartments of the river basin (e.g., soils) is accounted for in the model inputs that result from previous studies.
To identify potential sources of uncertainty in the model, we analyzed the sensitivity of major model outputs (TP export to the mouth, TP retention in the MRB and DIP retention in lakes and reservoirs) to the variation of 44 model inputs and parameters using the Latin Hypercube Sampling method (Saltelli et al. 2000). Therefore, 750 sensitivity runs were carried out for the period 1995-2000, using as initial state the 1995 conditions simulated in the default run. Each of these runs have the same spatial scale and set-up as the default model run, but the values of the 44 investigated input variables and parameters are randomly picked within a ± 5% range of their reference values. Based on the750 simulations, the sensitivity is expressed using standardized regression coefficients (SRCs), which represent the relative change of one output due to the relative change of one input (see S.1.3 for mathematical expressions).

Model validation and sensitivity analysis
The model properly reproduces the yearly TP and DIP concentrations close to the outlet of the basin, with root mean square errors (RMSE) of 17 and 29%, and average differences between observations and simulation of − 1 and 13%, respectively ( Fig. 2 and Table S.5). On average, DIP concentrations are overestimated at MRB USGS stations, but there is underestimation at 4 out of 10 stations over the 1970-2000 period (Table S.5). DIP and TP underestimations do not always coincide, which could mean that the discrepancy comes from uncertainties in point sources principally affecting DIP, or in in-stream process rates affecting the distribution between the different P forms in the water column.
RMSE and bias results for TP and DIP for all upstream stations in the main branch of the Mississippi are good (stations 3, 4, 5 and 8 -see Table S.4 for stations' descriptions). Upstream stations in the Ohio, Missouri, and Red River show poorer model performance for TP and DIP, with RMSE values ranging from 26 to 78%, but with small biases, except for station 1 for DIP (Ohio, Cannelton dam), station 6 for DIP (Missouri in Omaha) and station 9 for TP and DIP (Arkansas near Little Rock).
The sensitivity analysis of model outputs to 44 input variables and parameters shows that TP export at the river mouth is strongly influenced by the diffuse delivery of PIP, OP and DIP to the river network (SRCs of 0.66, 0.46 and 0.21, respectively), and by TSS loads (inputs from soil loss, SRC = − 0.38) ( Table 2). All other inputs and model parameters have an influence smaller than 4% on TP export and retention. DIP retention in lakes and reservoirs is significantly negatively influenced by diffuse inputs of PIP (SRC = − 0.34) and OP (SRC = − 0.26), since these forms tend to be desorbed or mineralized in these water bodies, and positively influenced by TSS loads (SRC = 0.55), additional soil loss providing more adsorption sites for DIP to be retained. The Langmuir isotherm maximum adsorption (SRC = 0.48) and DIP half-saturation concentration (SRC = − 0.40) are the only process parameters in the model that influence DIP retention by more than 4% (Table 2), while the results are little sensitive to other biogeochemical transformations.

P dynamics in MRB
Between 1970 and 2000, DIP, OP, and PIP inputs into the MRB have decreased by 10, 5 and 6 GgP (23, 8 and 6%), respectively (Fig. 3). The decrease in DIP inputs is quite homogeneous over most of the MRB (Fig. 4). DIP delivery to the river network has however increased in the central region of the watershed (eastern part of the Missouri basin), where agriculture has been identified as the dominant P source by Beusen et al. (2016), and in localized, densely populated areas (in the north-western part of the MRB and close to the river outlet), where sewage constitutes the major input. Likewise, outputs have decreased by 5, 4 and 6 Gg (15, 28 and 10%), respectively.
Due to deposition (ca 60% contribution to OP retention related to in-stream biogeochemical processes) and re-mineralization (ca 40% contribution),  Fig. S.6) and remains the dominant P form along the transfer from land to estuary in the MRB. PIP is produced in the water column via sorption of DIP onto suspended sediments ( Fig. S.5). However, large amounts settle in the benthic layer, leading to a net retention of 43% within the whole river network in 2000 (Fig. 3c). DIP represents 18-22% of the TP inputs to the basin's river network during 1970-2000 (Fig. 3a) and is 30-43% of the TP export (Fig. 3b). Retention of DIP is 11% over the whole MRB in 2000 (Fig. 3c). The budget of the MRB indicates that retained DIP has decreased by 5 Gg (from 20% of the inputs in 1970 to 11% in 2000). This decrease is mainly linked to expanding areas of "negative retention", i.e., in-stream DIP production, in the southern part of the basin (Fig. 4). In contrast, OP, and PIP retention have increased, from 79 to 83% and from 41 to 43% of the inputs, respectively. Overall, we estimate that the MRB accumulated a total of 3024 Gg of TP in its water/sediment network during the three decades period 1970-2000.
Low order streams contribute little to the overall P retention in the MRB (Fig. 5). 17% of all the OP entering the MRB is retained within high order streams between 1970 and 2000 (Fig. 5). Even though large amounts of PIP are retained within the whole MRB river network over the 1970-2000 period (1428 GgP), sorption of DIP onto particles in high order streams exceeds PIP settling, leading to a net production (i.e. negative retention) of PIP in these systems (358 GgP). According to our results, 473 Gg of DIP was retained in high order streams over these three decades, almost fully due to watercolumn sorption processes.
Results show that, within the river network, lakes and reservoirs play the largest role in the retention of different P forms. Particle P (i.e., PIP and OP) retention in the MRB is indeed mainly driven by processes within these systems (Fig. 5). On the contrary, our estimates show that lakes and reservoirs from the MRB constitute a net source of DIP over the 1970-2000 period. According to the model, TP retention in reservoirs was 75 Gg in 2000, with 35 Gg of OP retention through mineralization in the water column and benthic accumulation. In the same year, 11 Gg of DIP was released, mainly via benthic processes (i.e. re-mobilization of legacy P), so that net OP + DIP retention is 24 Gg.

Discussion
Model performance IMAGE-DGNM captures the measured total suspended solids (TSS), TP, and DIP trends at various MRB USGS stations (Section "Model validation and sensitivity analysis", Fig. 2 Moriasi et al. (2015). This process-based biogeochemical framework improves the estimates of Fig. 3 a Total delivery of P forms to the MRB river network, b export at the river mouth, and c retention due to in-stream biogeochemical reactions TP concentrations previously estimated with the IMAGE-GNM model, which simulates in-stream retention using the nutrient spiraling approach (Beusen et al. 2015). IMAGE-GNM results had a RMSE of 51% for the outlet station (against 29% in the present study), and consistently underestimated TP at the USGS monitoring stations.
However, simulations for TP concentrations still diverge from the measured data by a mean difference of only 15%. This discrepancy is likely due to an overestimation of both discharge [as suggested by Beusen et al. (2015)], and the soil loss estimates from IMAGE (soil loss being the dominant particulate P source). Overestimation of TSS concentrations confirm that these estimates are most likely slightly on the higher side in the MRB.
Looking at the various monitoring stations in the MRB, we see that P concentrations at the upstream stations are not consistently underestimated or overestimated (Figs. S3 and S4). RMSE values and biases for stations in the Mississippi main branch and most tributary stations indicate a good model performance. Disagreements in some of the tributaries may have several reasons. First, large fluctuations in the measured discharge and TSS/TP concentrations at some stations, with peaks during periods of high rainfall, may have a large influence on the derived yearly aggregated value, especially years with little sampling. For example, a peak value that is assumed to represent 2 months for a year with only 6 measurements also yields a peak in the aggregated yearly value. It is not known if this peak is in reality representative of conditions during 1 day (which would Retention of the different P forms in the water column and in the benthic layer of small streams (Strahler orders < 6), high orders (≥ 6), and for lakes and reservoirs in the MRB over the 1970-2000 period. The black bars indicate the net retention, resulting from both benthic and pelagic biogeochemical transformations yield a much lower aggregated annual value) or 2 months. Another possible reason for discrepancy, at this level of comparison, is that the spatial data for land use and wastewater discharge locations in urban areas may not be fully realistic. IMAGE land use distributions become less reliable in subbasins, where agriculture area is a small share of the total area and wastewater discharge occurs in all grid cells with urban population, while, in reality, wastewater emissions occur at discrete locations (wastewater treatment plants).
The results of the thorough sensitivity analysis show that the model reliability largely depends on diffuse inputs, in particular soil loss, and on the representation of adsorption. The negative SRC values for TP export by the MRB relative to TSS inputs (Table 2) indicate the enhanced retention due to P absorption to sediment particles. In the present work, we used sorption parameters from the literature that were determined within the MRB (Chao et al. 2010). These could however be spatially-variable due to, e.g., changes in sediment characteristics along the aquatic continuum. These results indicate that future research aiming at validating/refining soil loss estimates and the representation of P adsorption dynamics will help further strengthen assessments of the dynamics of bioavailable P in large, heavily dammed rivers, where nutrient inputs are dominantly diffuse, such as the MRB.

P dynamics in the MRB
Many studies quantifying nutrient dynamics at the basin scale focus on TP. From 1970 to 2000, we estimate that the MRB accumulated a total of 3024 Gg of TP in its water/sediment network. During this same timeframe, Beusen et al. (2016) calculated a higher TP retention in the MRB of 4479 Gg. They estimated that 9225, 4203 and 4545 GgP accumulated in other major basins such as the Amazon, the Nile, and the Yangtze basins, the latter facing a fourfold increase in P retention rate between 1970 and 2000 due to escalating anthropogenic inputs. Unlike these other basins, inputs into the freshwater network of the MRB have been reduced since the 1970s (Fig. 3 and Fig. S.6), which is mimicked by a decrease in the P export by the river basin to the estuary. While DIP retention has decreased, OP and PIP retention has risen, presumably as a response to dam building in the basin (increase in the total reservoir surface area by 14% between 1970 and2000, see Fig. 1b). P retention processes (including both physical retention and biogeochemical transformations) therefore not only reduce the amount of TP reaching the estuary, but also change its speciation at the mouth of the river (Fig. 3). This retained P remains mainly as sedimentary deposits.
Environmental impacts of P loading to watersheds depend on specific forms of P in the river network (i.e., its bioavailability). A striking example is the one of Lake Erie, where, despite a decrease in TP delivery during the past decades, the DIP load has concomitantly risen, fueling recent harmful algal blooms (Jarvie et al. 2017). Our detailed biogeochemical framework, by explicitly modelling DIP, allows for a more thorough understanding of transfers of this P form in relation to the other forms and hydrological characteristics (i.e., lake and reservoir versus small or large streams). Even though DIP represents only 18-22% of the TP inputs to the basin's river network for the study timeframe, it constitutes 30-43% of the TP reaching the estuary. Retention of DIP (11% over the whole MRB in 2000) is much lower than that of OP and PIP (Fig. 3c). DIP sorption in the water column (forming PIP) is counterbalanced by the release of DIP through OP mineralization in the water column and via benthic processes (mineralization of deposited OP or desorption from benthic particles) (Fig.  S.5).
The significant sediment loads in the MRB, high with respect to other temperate basins (Ludwig and Probst 1998), originate mainly from the muddy parts of the Missouri basin (Meade and Moody 2010), and exert a sizable control on DIP dynamics. For one, they provide ample sites for DIP to sorb, whilst also reducing the light penetration and limiting primary production, a sink for DIP. These effects have been directly observed even near the mouth of the MRB (Chao et al. 2010). Headwater streams, accounting for most of the length of river networks, typically transport small amounts of sediments, due to low flows, high roughness, etc., thereby potentially acting as temporary sediment reservoirs (Benda et al. 2005). At the basin scale, they however have little effect on the retention of TP or of individual P forms (Fig. 5). In larger streams, sediments are kept in suspension by the flow, leading to yearly concentrations as high as 2000 mg L −1 (see measurements at USGS stations in Fig. S.2). Sediment dynamics can hence significantly influence P retention within these systems.
Our results confirm the strong effect of sediment dynamics and P transfers in high order streams, where P fluxes and speciation are strongly affected by sorption processes. Large DIP amounts are delivered to high order streams, with high dilution potential, through point source effluents (see Fig. S.6). This DIP is rapidly sorbed onto to the suspended sediments, mainly originating from soil loss from adjacent agricultural lands, which is unsaturated in P. DIP concentrations in the mainstream are more than one order of magnitude lower than the maximal sorption capacity defined by Chao et al. (2010) (for sediments of a Mississippi Delta lake), and used in the model. Therefore, PIP concentrations tend to be linearly related to DIP concentrations. In this way, high order streams constituted a net sink for DIP during the 1970-2000 period and a net source for PIP.

Effects of lakes and reservoirs on DIP dynamics
Particle P (i.e., PIP and OP) retention in the MRB is mainly driven by processes in lakes and reservoirs (Fig. 5). Accumulation in bottom sediments is high due to low flow velocities in these systems. According to Maavara et al. (2015), who simulated the long-term retention of different forms of P in global reservoirs, the MRB remained the top P-retaining catchment in the world for the period 1970-2000. They estimated, for the year 2000, a retention of 67 Gg of TP (53% of the total load) and 30 Gg of reactive P (i.e., total dissolved P, exchangeable P, and particulate OP, 52% of the total load). Our results for reservoirs are very similar: we estimate that 75 Gg of TP is retained in reservoirs in 2000, with an OP + DIP retention of 24 Gg. Yet, their model does not account for spatial variability of P delivery and processes throughout the basin, and does not represent the potential re-mobilization of accumulated P. According to our results, not all "reactive P" is retained in reservoirs: while OP is retained through mineralization and benthic accumulation, DIP is overall produced in reservoirs. We find that, while DIP retention is highest in the north-western part of the MRB (Missouri basin), receiving large sediment loads, negative DIP retention (i.e., release) occurs mainly in the most southern part (Fig. 4), where the density of dams is high (Fig. 1b).
Even if lakes and reservoirs tend to trap large amounts of the P load, sediment P can then constitute a long-term source for the aquatic biota (Carignan and Kalff 1980). James and Barko (2004), for example, showed that Lake Pepin, in the Upper Mississippi River, exports 30% more soluble reactive P than it receives in the summer. This re-mobilization process, also referred to as "internal loading" source, has been identified as significant in different lakes and reservoirs over the globe. It can potentially mask the effects of external loading reduction measures (Sharpley et al. 2013). The recovery period following P loading reduction depends on the loading history and the accumulation of P in the sediment (Sondergaard et al. 2003). In some cases, a negative P retention can occur for decades. In the Iron Gate I Reservoir, nutrients deposited between 1970 and 1990, when loads to the Danube were 30-40% higher, seemed to be re-mobilized decades later due to early diagenesis (Teodoru and Wehrli 2005). Sharpley et al. (2013) also provide the example of Loch Leven, where in the years following efforts to restore lake water quality, summer P concentrations increased, with a dominance of internal to external loading. This internal loading source can depend on many environmental factors; it can notably be enhanced by anoxia of bottom water (Nürnberg 2009). This source is not commonly included in large scale P budgets. However, it typically releases DIP, and thus, when investigating effects of P transfers on ecosystem functioning, our results show that it should be taken into account.
Besides potentially releasing DIP to the water column, trapping of sediments by dams can also alter the downstream propensity of sorbing DIP inflows. This altering of the ability of suspended sediments to retain river bioavailable P has been highlighted downstream the Three Gorges Dam (Cao et al. 2011) and in the Yellow River (Pan et al. 2013), thereby increasing eutrophication risk. River ecosystems downstream from impoundments might therefore be more sensitive to point sources, due to lower buffering capacity. The risk of re-suspension of accumulated particle P might also be greater, if the sediment yield is below carrying capacity of the flow at the outlet of the reservoir [increased erosion by sediment-starved flows (Kondolf 1997)].

Conclusions
Damming can play a pivotal role in the shift in P speciation during its transfer in whole river networks. In the case of the MRB, we show that DIP retention in the basin is significantly reduced due to net release in lakes and reservoirs, leading to an increase in the proportion of this mainly bioavailable form downstream.
The sensitivity analysis points out that future research aiming at validating/refining soil loss estimates and the representation of P absorption will help to further strengthen assessments of the fate of bioavailable P in large, heavily dammed rivers, where P emissions predominantly originate from diffuse sources, such as the MRB.
The analysis with IMAGE-DGNM results demonstrates that distinguishing between different P forms, particularly DIP, in large-scale transfer estimates is crucial to better assess potential ecological risks in aquatic environments from anthropogenic activities. The relative amounts of these forms depend not only on inputs, but also on in-stream transformations. Quantifying them at a fine spatio-temporal resolution requires a process-based representation of the physical and biogeochemical reactions involved. This includes P internal loading processes that have been identified in earlier lake studies, but typically lack in large-scale watershed P budgets.
The implementation of P dynamics in IMAGE-DGNM presented here, highlights the importance of including sediment dynamics and benthic transformations in large-scale biogeochemical models. In the MRB basin, release of previously accumulated P in lakes and reservoirs in the highly bioavailable DIP form threatens to exceed DIP uptake along the aquatic continuum, increasing its export to downstream ecosystems. Accounting for such legacy effects is crucial to better evaluate the effect of future pollution mitigation strategies over time. This has significant implications for the eutrophication management in aquatic ecosystems in the context of booming dam construction worldwide (see e.g. Zarfl et al. 2015).