Increased extreme rains intensify erosional nitrogen and phosphorus fluxes to the northern Gulf of Mexico in recent decades

Soil erosion delivers enormous amounts of macro-nutrients including nitrogen (N) and phosphorus (P) from land to rivers, potentially sustaining water column bioavailable nutrient levels for decades. In this study, we represent erosional N and P fluxes in the Energy Exascale Earth System Model (E3SM) and apply the model to the continental United States. We estimate that during 1991–2019 soil erosion delivers 775 Gg yr−1 (1 Gg = 109 g) of particulate N (PN) and 328 Gg yr−1 of particulate P (PP) on average to the drainage basins of the northern Gulf of Mexico, including the Mississippi/Atchafalaya River and other rivers draining to the Texas Gulf and the Eastern Gulf. Our model simulation shows that in these rivers PP is the dominant P constituent and over 55% of P exported by erosion comes from soil P pools that could become bioavailable within decades. More importantly, we find that during 1991–2019 erosional N and P fluxes increase at rates of about 15 Gg N yr−1 and 6 Gg P yr−1, respectively, due to increased extreme rains in the Mississippi/Atchafalaya river basin, and this intensification of erosional N and P fluxes drive the significant increase of riverine PN and PP yields to the northern Gulf of Mexico. With extreme rains projected to increase with warming, erosional nutrient fluxes in the region would likely continue to rise in the future, thus complicating the effort of reducing eutrophication in the inland and coastal waters.


Introduction
Soil erosion is one of the most important pathways for the lateral fluxes of nitrogen (N) and phosphorus (P) from land to rivers (Pimentel and Kounang 1998, Beusen et al 2005, Liu et al 2010, Berhe et al 2018. Unlike soil leaching that dominates the lateral fluxes of dissolved N and P, erosional N and P fluxes are mostly in particulate forms and can contribute as much as 34% and 78% of riverine N and P yields to oceans, respectively, across the globe (Seitzinger et al 2010). With carbon (C) and nutrient cycling closely coupled in terrestrial and aquatic ecosystems (Peñuelas et al 2013), this massive transfer of N and P has immense influences on the global and regional ecology and biogeochemistry (Ward et al 2017). On land, the loss of N and P would degrade soil quality and limit atmospheric CO 2 fixation by vegetations (particularly crops) (Berhe et al 2018). In surface waterbodies, erosional particulate N and P (hereafter referred to as PN and PP, respectively) gained from land can be remobilized in the water column through chemical transformations (e.g. desorption of inorganic N and P from reduced iron, manganese and other redox-sensitive compounds, and mineralization of organic N and P) (Sutula et al 2004, Ward et al 2017. These mobilized nutrients will then be easily accessible to aquatic organisms. Further, data showed that PP is frequently more abundant than dissolved forms of P in inland waters (Meybeck 1982, Froelich 1988, Goolsby et al 1999, Seitzinger et al 2010.
Recent surveys showed that excessive N and P loadings mainly from anthropogenic sources are responsible for high levels of N and P eutrophication in over a quarter of streams and one-fifth of lakes in the United States (US)  and for development of the world's second largest coastal hypoxic zone during summertime in the northern Gulf of Mexico (Rabalais et al 1999, Sylvan et al 2006, Dale et al 2007. Although the US government has invested billions of dollars in managing N and P inputs from anthropogenic sources (Jeppesen et al 2005, Mississippi River/Gulf of Mexico Watershed Nutrient Task Force 2015, the results are not always as good as expected (Finlay et al 2013, Van Meter et al 2018. These unsatisfactory outcomes could be partly attributed to the nutrient legacy effect, i.e. N and P that entered watersheds a long time ago continue to sustain soluble reactive N and P levels in inland and coastal waters through groundwater discharge, inorganic matter desorption and organic matter decomposition (Froelich 1988, Finlay et al 2013, Van Meter et al 2018, Kreiling et al 2019. As such, realistic estimates of erosional N and P fluxes are essential to understand and manage the N and P enrichment in the northern Gulf of Mexico and over the continental US. Furthermore, erosional N and P fluxes (i.e. soil erosion driven PN and PP loads to inland waters) are tightly connected to both water and biogeochemical cycles. It is thus important to represent such processes in Earth system models to understand their evolution in the context of climate change (Gonzalez-Hidalgo et al 2010, Tan et al 2017 and their interactions with other environmental changes (Scavia et al 2002, Heisler et al 2008. For this purpose, we develop a new process-based model to represent erosional N and P fluxes within the framework of the E3SM (Golaz et al 2019, Burrows et al 2020. Upon successful validation of the model in the continental US, we use it for two scientific questions: (a) how may observed climate variability in recent decades impact erosional N and P fluxes in the continental US? (b) How may changes in erosional N and P fluxes influence eutrophication, defined as over-enrichment in bioavailable nutrients, in inland and coastal waters? Particularly, we focus on the Mississippi/Atchafalaya river basin (MARB) that has a strong influence on hypoxia and other environmental issues in the northern Gulf of Mexico.

Model description
We develop a new erosional N and P model building on an existing soil erosion module within the E3SM land model (ELM) (hereafter referred to as ELM-Erosion) (Tan et al 2020). ELM-Erosion is essentially an improved version of the Morgan-Morgan-Finney soil erosion model (Morgan 2001) and runs with a daily time step (Tan et al 2018). It simulates soil erosion and erosion-induced lateral fluxes of sediment and particulate organic carbon (POC) based on ELM simulated hydrological, ecological and biogeochemical conditions (Tan et al 2020). We validated ELM-Erosion at small catchment scales (Tan et al 2018), and further demonstrated its ability to capture the spatial variability of soil erosion and sediment and POC yield in large river basins of the continental US, including the MARB and its sub-basins (Tan et al 2020). The effectiveness of ELM in simulating the spatial and temporal variability of land hydrology and soil C, N and P pools have been well documented in previous studies (Ricciuto et al 2018, Golaz et al 2019, Zhu et al 2019, Burrows et al 2020. In this work, we extend ELM-Erosion to include erosional N and P fluxes defined as the product of erosional POC flux and soil particle C/N and C/P mass ratios in surface soils. ELM has two representations of soil biogeochemistry based on the Equilibrium Chemistry Approximation (ECA) and Convergent Trophic Cascade (CTC) routines (Burrows et al 2020). Soil C, N and P cycle among different forms (i.e. organic and inorganic forms), with inputs from litterfall, fertilization and atmospheric deposition, and outputs through decomposition, plant uptake, nitrification/denitrification, erosion and runoff/leaching (Yang et al 2014, Zhu et al 2016. Specifically, the particulate forms of soil C, N and P in ELM are represented by three soil organic matter pools with different turnover times in the ECA routine versus four soil organic matter pools with different turnover times in the CTC routine (Yang et al 2014, Zhu et al 2016. In addition, both routines represent the particulate forms of soil P with four mineral pools (inorganic labile P, parent material P, secondary mineral P and occluded P) (Yang et al 2014, Zhu et al 2016. Correspondingly, the ELM-Erosion model assumes that erosional PN originates dominantly from soil organic N pools, whilst erosional PP originates from both soil organic P pools and mineral P pools, which is also consistent with Berhe et al (2018).
In this study, instead of relying on a single soil biogeochemistry routine, we use an ensemble approach to estimate erosional PN and PP fluxes. Using an ensemble approach is important because (a) soil particle C/N and C/P ratios simulated by soil biogeochemistry routines usually have large uncertainties (supplementary figure S1 (available online at stacks.iop.org/ERL/16/054080/mmedia)), (b) soil biogeochemistry routines usually do not well resolve the vertical variation of C:N:P stoichiometry in the soil column (Rostad et al 1997, Tipping et al 2016, and (c) the C/N and C/P ratios in eroded soils are usually much lower than those in surface soils due to the preferential erosion of N-and P-rich fine soil particles (Meybeck 1982, Lindström et al 2010. In the ensemble approach, the C/N and C/P ratios in eroded soils are estimated using three methods: the CTC routine, the ECA routine and the empirical method of Beusen et al (2005). The empirical equations of the C/N and C/P ratios in Beusen et al (2005) were compiled from a large observational dataset of riverine N and P: where PN c and POC c are the contents of PN and POC in percentage in eroded soils, respectively, and PP m and POC m are the mass fluxes (kg yr −1 ) of eroded PP and POC. POC c is calculated as the ratio of soil organic C density to soil bulk density. Noticeably, because the exponent of POC m is close to 1, the C/P ratio calculated from equation (2) is mostly at a constant value of 22 for rivers worldwide. This C/P ratio of 22:1 is consistent with the reported value of Meybeck (1982). The change of soil N and P pools after erosion are modeled in ELM-Erosion following the method of Naipal et al (2018). In this method, the loss of C, N and P by erosion from the surface soil layer causes the vertical displacements of C, N and P from the deeper soil layers. At the bottom of the soil column, where no soil layer exists underneath, we assume that only parent material P can be moved upward from the bottom.
Because E3SM does not yet represent river biogeochemistry, we cannot directly estimate how much PN and PP are lost in river systems during transport through either deposition or degradation. Instead, we estimate riverine PN and PP yields from river basins based on the nutrient spiraling theory of inland waters which assumes that the nutrient retention in rivers is an exponential function of the travel time of nutrient molecules in the river network (Runkel 2007): where PN i and PP i are erosional PN and PP fluxes from a grid cell in a river basin (kg yr −1 ), respectively, PN o and PP o are erosional PN and PP yields from the river basin (kg yr −1 ), respectively, τ i is the water residence time for a grid cell in the river basin (yr), λ PN and λ PP are PN and PP retention rates in the channels of the river basin (yr −1 ), respectively. In E3SM, the river water resident time is estimated by the Model for Scale Adaptive River Transport (MOSART) river model which simulates runoff routing using the kinematic wave method, considering channel geometry and roughness . MOSART has been coupled with a water management model that represents reservoir operation and irrigation and their impacts on streamflow (Voisin et al 2013). The basin-scale PN and PP retention rates λ PN and λ PP are calibrated by minimizing the difference between estimated and observed riverine PN and PP yields, respectively.

Simulation protocol
Following Tan et al (2020), numerical experiments are conducted over the North American Land Data Assimilation System (NLDAS) 1/8-th degree resolution latitude/longitude grid for the period of 1991-2019 driven by the NLDAS-2 atmospheric forcing data (Mitchell et al 2004). The ELM land biogeochemistry is first spun up for 200 years with accelerated decomposition spin-up and then another 600 years with regular spin-up, using repeated atmospheric forcing from an E3SM preindustrial simulation (Golaz et al 2019). More details about the ELM spin-up process are described in Burrows et al (2020). After the total ecosystem C has reached an equilibrium, the model is then run from 1850 to 1990 using atmospheric forcing from an E3SM historical simulation to generate initial conditions for the land states. Cropland erosion management (e.g. no tillage) is a critical factor governing soil erosion over agricultural land (Montgomery 2007, Meade andMoody 2010). In theory, because the application of conservation practices affects the canopy cover and ground cover, these conservation practices should be directly parameterized in ELM-Erosion. However, as ELM-Erosion and Revised Universal Soil Loss Equation (RUSLE) have very different model structures and run on different spatial resolutions, it is difficult for ELM-Erosion to effectively parameterize conservation practices. For example, the support practice factor and prior-land-use factor used by RUSLE need high-resolution data of conservation measures and soil properties to calculate, which is difficult for the coarse resolution ELM. Further, conservation practices usually affect soil erosion in a combined manner, which means that simply implementing one or two RUSLE subfactors in ELM may not improve the model performance substantially. Thus, ELM-Erosion adopts a data-driven approach to mimic its effect. First, we extracted the statelevel soil erosion intensity in cropland from the USDA National Resources Inventory (NRI) benchmark estimates (USDA 2018). The NRI estimates were based on the RUSLE approach which includes three dynamic factors (Renard et al 1997): the rainfall erosivity factor R, the cover-management factor C, and the support practice factor P. During our simulation period of 1991-2019, the effect of the R factor on the erosion dynamics is minor because there is no discernible trend in mean annual rainfall (see our discussions below) and this factor is less sensitive to extreme rainfall when calculated based on coarse resolution precipitation data (Tan et al 2018). As such, the temporal variability of the estimated soil erosion intensity should mainly reflect the temporal variability of cropland management (the C and P factors). In the next step, we interpolated the 5 year interval NRI estimates linearly in time to construct annual erosion intensity time series at each grid cell during 1991-2019, which were then normalized by the erosion intensity of the grid cell in 1991 (supplementary figure S2). We then apply the normalized soil erosion intensity as an adjustment multiplier to the equations of rainfall-and runoff-driven erosion (Tan et al 2020 and supplementary text S1) to account for the change in crop management effect relative to 1991.
We conducted three ensembles of simulations in this study, with each ensemble corresponding to one scenario. The first scenario S0 uses the adjustment factor for cropland erosion management described above. The second scenario S1 does not use the adjustment factor. The third scenario S2 is similar to S0 except that the rainfall forcing is manipulated. Each ensemble consists of three simulations with the C/N and C/P ratios in eroded soils estimated using the CTC routine, the ECA routine and the empirical method. Between S0 and S1, another difference is that S0 includes the change of land use based on the Land Use Harmonized version 2 transient dataset (Hurtt et al 2011, supplementary figure S3) but S1 does not. Because our study focuses on the erosional loss of PN and PP in soils while N and P sourced from fertilizer and manure applications are mostly lost through runoff and leaching (Berhe et al 2018), we do not expect fertilizer and manure applications to affect our simulations significantly. However, as nutrient limitation is important for plant growth, we apply N fertilization in the S0 and S1 simulations to better quantify the impact of plant canopy cover and ground cover on erosion. In ELM, N fertilization is applied to the crop land units based on prescribed crop-specific fertilization rates (Oleson et al 2013). P fertilization is not considered in our study because it has not been implemented in ELM (Oleson et al 2013), and it is perceived to mainly impact plant growth in tropical regions (Zhu et al 2019). Model validation is based on S0 which is a more realistic configuration accounting for land use change and temporal variations in crop management effect, but analysis of how climate impacts erosional nutrient fluxes is based on the climate-only scenario S1.
The rainfall manipulation scenario S2 was conducted to isolate the effect of extreme rains on erosional PN and PP fluxes. The major difference between S2 and S0 is that the interannual variations of non-extreme rainfall in the NLDAS-2 atmospheric forcing data are removed for S2. Specifically, for each year of 1991-2019, we first divide all the days in each grid cell into two groups: the days with extreme rains and the days without extreme rains. Days with extreme rains are determined by the 95th percentile of all daily rains of the year. We then repeat the rainfall data of 1991 for 28 years to construct the rainfall data of 1992-2019. Finally, we replace the extreme rainfall in the constructed rainfall data for each year of 1992-2019 with the extreme rainfall from the real rainfall data of 1992-2019.

Model validation and statistical analysis
To validate the simulated C/N and C/P ratios of detached sediment that enters rivers, we collected published C/N and C/P ratios of suspended sediment for selected large US rivers (figure 1(a)) from a literature review (supplementary table S1). Noticeably, due to biogeochemical processes, the C/N and C/P ratios of suspended sediment can vary during transport, which might introduce uncertainties to the validation. To validate the dynamics of estimated riverine PP flux, we used the observed riverine PP (total P minus orthophosphate) yield data for the MARB and its sub-basins from USGS (Lee and Reutter 2019). Our validation focuses on the temporal variability of riverine PP flux but not PN flux because we do not have observed PN flux for such a comparison. Noticeably, using a similar method (total N minus dissolved N) to estimate riverine PN flux from the USGS total and dissolved N data is not appropriate. It is because terrestrial PN only contributes to a small fraction of total riverine N (Meybeck 1982, Goolsby et al 2000), and dissolved N in rivers is strongly controlled by N fixation, nitrification and denitrification. Model performance is assessed using three metrics: relative error, the coefficient of determination (R 2 ) and Nash-Sutcliffe coefficient.
Here for each grid cell, we define extreme rains as daily rainfall that is above the 95th percentile of all daily rains during 1991-2019. Correspondingly, the surface runoff and erosional N and P fluxes during the days with extreme rains are marked as extreme rainfall driven. For both hydrological and nutrient fluxes, their trends of annual total and extremes are calculated using linear regressions and verified by the Mann-Kendall test with the significant value of 0.05.

Model validation
Model validation shows that ELM-Erosion can reproduce the C/N ratios of suspended sediment in the large rivers of the continental US (figure 1). Because ELM-Erosion already captures the spatial variability of POC fluxes in the continental US (Tan et al 2020, supplementary figure S4), the realistic simulation of C/N ratios means that the model can also reproduce the spatial variability of erosional N flux. In fact, in of simulated eroded organic matter C/N ratios (diamonds) with previously estimated C/N ratios of suspended organic matter (circles). The MARB subbasins are outlined with grey boundaries. The mean and standard deviation of previously estimated C/N ratios can be found in supplementary table S1. The mean and standard deviation of simulated C/N ratios are from the ensemble of three S0 runs. most cases (even including the Arkansas, Sacramento and Rio Grande rivers), the estimated C/N ratios from equation (1) have relative errors less than 10%. Noticeably, the rivers with large positive model biases of C/N ratios are all in the regions of warm climate, suggesting that the positive biases could partly be caused by in-river biogeochemical processes, such as primary production, which are not represented in E3SM. As illustrated, the simulated C/P ratios of suspended sediment by equation (2) are very close to 22 across the continental US, which is consistent with observations documented by Meybeck (1982). Encouragingly, the interannual variability of estimated riverine PP yield is consistent with observations in the large rivers of the MARB (figure 2). In some cases, our method fails to capture the extreme high or low PP yield in the observations, such as the observed extremes of PP yields in the Missouri and Tennessee rivers (figure 2), which could be caused by the oversimplification of river PP retention in our method. For example, river PP deposition may deviate from the normal condition during extremely wet and dry years. Both model estimates and observations show that the riverine PP yield in the MARB increases during the water year of 1992-2018: the observed trend is 1.3 Gg P yr −1 (1 Gg= 1 × 10 9 g) (Mann-Kendall test, p < 0.05) and the estimated trend is 1.7 Gg P yr −1 (Mann-Kendall test, p < 0.05). Here a water year is defined as the 12 month period starting from October 1 prior to a given year. Due to in-river deposition and biogeochemistry, only a small fraction of erosional PP flux eventually reaches river outlets (supplementary table S2): on average 38 ± 5% in the MARB.
Results from simulations S0 show that soil erosion can drive enormous amounts of PN and PP transferring from land to inland waters in the continental US during 1991-2019. On average, soil erosion delivered about 1047 ± 269 Gg N of PN and 445 ± 215 Gg P of PP, respectively, to rivers each year (figure 3). In the continental US, the MARB contributed the largest erosional N and P fluxes: about 718 ± 102 Gg N yr −1 and 302 ± 105 Gg P yr −1 , respectively. Within the MARB, most of the erosional N and P fluxes originate from the Missouri, Ohio and Upper Mississippi river basins (Goolsby et al 1999, 2000, Goolsby and Battaglin 2001 because of their vast cropland areas (Tan et al 2020, figure 3). Noticeably, as more than 60% of eroded soils are redeposited on land and never entered the river networks in the region (Tan et al 2020), the amounts of soil N and P disturbed by erosion are much larger than the amounts of erosional PN and PP loads to rivers ( figure 3 and supplementary figure S5). Thus, the impact of soil erosion on soil N and P cycling is likely much larger than what we presented here.
The modeled erosional N and P fluxes cannot be directly compared with riverine PN and PP yields, because of retention of PN and PP in river networks by reservoir trapping and chemical transformations (Maavara et al 2020). However, the modeled erosional N and P fluxes are expected to be within the bounds of the measured riverine post-dam PN and PP yields (the lower bound) and pre-dam PN and PP yields (the upper bound). Using gauge data, Goolsby and Battaglin (2001) estimated the mean annual postdam PN yield from the MARB during 1980-1996 to be about 210 Gg N yr −1 . Meanwhile, the pre-dam PN yield from the MARB is about 600 Gg N yr −1 , which is based on the assumption of 0.15% N content in suspended sediment (Goolsby et al 2000) and 4 × 10 5 Gg of suspended sediment yield in the pre-dam condition (Milliman and Syvitski 1992). Our modeled erosional N flux of 718 ± 102 Gg yr −1 in the MARB is at the higher end of the above range. It is also larger than the modeled mean annual PN flux from the MARB (320 ± 160 Gg N yr −1 ) estimated by Beusen et al (2005). A possible reason for our higher estimate is that soil C and N contents simulated by the ELM soil biogeochemistry are higher than the prescribed values of Goolsby et al (2000) and Beusen et al (2005). Another possible reason is that as shown in our study, the erosional N flux has continued increasing since 2000s. Our estimate of PN flux from the Saint Lawrence river basin (31 ± 7 Gg N yr −1 ) is of the same order of magnitude as the river measurements (22-42 Gg N yr −1 ) during 1981-1985(Pocklington and Tan 1987. For the MARB, our modeled erosional P flux is within the range of the observed post-dam riverine PP yield (on average 95 Gg P yr −1 during [1980][1981][1982][1983][1984][1985][1986][1987][1988][1989][1990][1991][1992][1993][1994][1995][1996] (Goolsby et al 1999) and the estimated pre-dam riverine PP yield (about 415 Gg P yr −1 ) (Meybeck 1982, Milliman andSyvitski 1992).

Impacts of climate and land management changes
The impact of cropland erosion management on erosional N and P fluxes during 1991-2019 is moderate. Compared with S0, erosional N and P fluxes in the continental US in S1 (ignoring cropland erosion management and land use change) are only about 76 Gg N and 32 Gg P higher, respectively, on average (table 1). Similarly, compared with S0, erosional N and P fluxes from the MARB in S1 are only about 70 Gg N and 29 Gg P higher, respectively, on average (table 1). The small difference in erosional N and P fluxes between simulations S0 and S1 can be explained by (a) the gradual maturity of cropland erosion management methods since 1990s that slowed down erosion reduction (supplementary figure S2) and (b) the small fraction of eroded soils that can enter the river networks due to limited sediment transport capability by surface runoff (Tan et al 2020). However, without including cropland erosion management and land use change, erosional N and P fluxes would have increased more quickly in the continental US (particularly in the MARB) during 1991-2019. For example, the increasing trend of erosional P flux in the MARB in S1 is 8.2 Gg P yr −1 (Mann-Kendall test, p < 0.05), which is 2.1 Gg P yr −1 (34%) higher than that in S0 (table 1). Based on the ensemble of three S0 runs, we estimate that soil erosion generates 775 Gg yr −1 of PN flux and 328 Gg yr −1 of PP flux on average to the drainage basins of the northern Gulf of Mexico during 1991-2019 (table 1). Since the MARB is the largest source of PN and PP for the northern Gulf of Mexico (on average 93% and 92%, respectively), increases in erosional N and P fluxes from the MARB drove significant positive trends (Mann-Kendall test, p < 0.05) of erosional N and P fluxes from land to the drainage basins of the northern Gulf of Mexico (figure 4): about 15 Gg N yr −1 and 6 Gg P yr −1 , respectively during 1991-2019 (table 1). Far behind the MARB, the Texas-Gulf region (supplementary figure S4) is ranked second in contributing erosional N and P fluxes to the northern Gulf of Mexico: on average 3% of N and 3% of P, respectively. As indicated by the riverine N and P data in the MARB, the fraction of PN in the total riverine N is only 30%, but the fraction of PP in the total riverine P is as large as 70% (Goolsby et al 2000, Goolsby and Battaglin 2001, Lee and Reutter 2019. The relative contributions of dissolved and PN and PP to the total riverine N and P are probably determined by the difference in land N and P export pathways (i.e. leaching, runoff and soil erosion) rather than in-stream processes (supplementary figures S7-S8). As such, it is reasonable to conclude that the increase of erosional nutrient flux should have a much more significant impact on P enrichment than N enrichment in the northern Gulf of Mexico. Moreover, erosional PP in the region consists of a large fraction of organic and inorganic P that are bioavailable to aquatic organisms (supplementary figure S8). Over 55% of eroded P consists of soil organic P and mineral P that can be degraded or desorbed within decades (active and slow soil organic matter, labile PP and secondary mineral PP), hence contributing to the nutrient legacy effect. Because P is an important limiting nutrient for algal growth during spring and early summer in the northern Gulf of Mexico (Sylvan et al 2006, Dale et al 2007, the increase of erosional P flux may enhance P enrichment and exacerbate other environmental issues in the coastal waters (Cai et al 2011, Laurent andFennel 2017).
Without considering cropland erosion management and land use change, erosional N and P fluxes to the drainage basins of the northern Gulf of Mexico could have increased at faster rates (Mann-Kendall test, p < 0.05) of 20.2 Gg N yr −1 and 8.4 Gg P yr −1 , respectively during 1991-2019 (table 1). Soil erosion in ELM-Erosion is jointly controlled by the change of climate and vegetation cover (supplementary text S1). As the change of vegetation cover only plays a minor role in soil erosion evolution during the period (Tan et al 2020), the positive trends of erosional N and P fluxes are thus mostly caused by the change of rainfall and runoff. Furthermore, as shown in figure 5(a), the increase of rainfall and runoff mainly occurs in the Ohio, Tennessee and Upper Mississippi river basins. Noticeably, erosional N and P fluxes share similar spatio-temporal variabilities because they are both dominated by the variability of sediment flux instead of that of soil PN and PP pools (Berhe et al 2018). We differentiate the simulated erosional P flux into two groups: those occurring on days with extreme rain events, and those occurring on days without extreme rain events. We find that as much as 80% of the positive trend of erosional P flux in the three river basins is dominated by the increase of erosional P flux in the first group (figures 5(b)-(d)). The weaker P flux trend in the second group is due to the negligible increase of mean annual rainfall and runoff (figures 5(e) and (f)). Further analysis shows that the total extreme rainfall volume per year actually only increases mildly during 1991-2019 (figure 5(g)). However, due to the nonlinear relationship between rain intensity and As documented in our observation-based analysis, soil erosion is highly responsive to rain and runoff extremes (Tan et al 2017). Hence, the increase of extreme rainfall and runoff in the three river basins drives significant increases of erosional N and P fluxes to the northern Gulf of Mexico. With the same river retention ratios as simulations S0 but allowing only extreme rains to vary temporally while interannual variability of non-extreme rains is removed, simulations S2 shows a similar increase of riverine PP yield in the MARB in recent decades (supplementary figure S9): 2.3 Gg P yr −1 (Mann-Kendall test, p < 0.05). Hence the rainfall manipulation test S2 provides further support of the role of extreme rains in the increased erosional N and P fluxes in recent decades. Furthermore, we compared the histograms of simulated erosional P fluxes at the two time periods of 1991-1995 and 2015-2019 to highlight the trend (figure 4). Consistently, this analysis also suggests that large erosional P production occurs more frequently during 2015-2019 than during 1991-1995 in the three river basins (figure 6). Particularly in the Ohio river basin, daily erosional P flux larger than 9 mg P m −2 d −1 was absent during 1991-1995 but starts emerging during 2015-2019.

Limitations and future development
Previous studies showed that sediment yields in the MARB and many of its sub-basins have been decreasing in recent decades due to reduction of soil erosion and increase of sediment trapping in reservoirs (Meade and Moody 2010, Mize et al 2018, Murphy 2020. There are two possible explanations for the contradicting trends between the observed sediment and PP yields in the MARB. First, despite the reduction of soil erosion, the increase of overland sediment transport (supplementary figure S10) may actually have increased sediment loads to rivers . Simulated (a) decadal trend of erosional P flux, (b) decadal trend of erosional P flux associated with 95th percentile daily rain extremes, (c) the percentage of decadal erosional P flux trend associated with rain extremes, (d) the percentage of decadal erosional P flux trend associated with non-extreme rains, (e) decadal trend of annual rainfall, (f) decadal trend of annual surface runoff, (g) decadal trend of annual rainfall associated with 95th percentile daily rain extremes, and (h) decadal trend of annual surface runoff associated with 95th percentile daily rain extremes during the period of 1991-2019. In (a)-(h), colors are only drawn on grid cells with significant trends (Mann-Kendall test, p < 0.05). In (c) and (d), negative trend is not shown. (a)-(d) are based on the ensemble of S1 runs.
in recent decades. However, because sediment yield in the MARB was strongly regulated by the river sediment budget which has been decreasing in recent decades (Meade and Moody 2010), the increase in sediment loads was not potent enough to increase the sediment yield. Second, despite the reduction of total sediment loads, loads of fine sediments that contain more PP may actually have increased. Notably, the Figure 6. Comparison of histograms of daily erosional P flux in the first five simulation years (1991-1995, red) with those in the last five simulation years (2015-2019, blue) in the large river basins of the northern Gulf of Mexico based on the ensemble of S1 runs.
observations also confirm the increase of sediment yield in the Ohio river basin during 1991-2012 (Murphy 2020) where our estimated PP yield trend is the strongest (figure 4). But we acknowledge that the possibility of the second explanation cannot be fully excluded by our results because our model does not represent management actions (such as the use of riparian buffer strips) to attenuate sediment loads to rivers.
Our findings that erosional N and P fluxes are enhanced significantly in the MARB by the increase of extreme rainfall events have two important implications. First, because extreme rains are projected to further increase in general under climate warming and the increase in extreme rains will be larger than the increase of mean annual rainfall (Allen andIngram 2002, Du et al 2019), both the total amount and mean concentration of riverine PN and PP may further increase. As such, the impact of soil erosion on N and P enrichment in the northern Gulf of Mexico is expected to become larger in the future. Further, because soil erosion is the dominant pathway for riverine P loading, the increase of erosional flux will influence P enrichment more substantially. Second, the mechanism described above could also affect other hazardous materials (e.g. heavy metals) that can be detached by erosion and transported to the MARB and the northern Gulf of Mexico (Presley et al 1980). Therefore, we suspect that fluxes of these hazardous materials could also increase in the future.
We note that a few limitations in the model may introduce large uncertainties in our simulations. These include the lack of representation of river nutrient retention, cropland erosion management, riparian management, and the interactive effect of climate and land use. For example, river retention of PN and PP is essential to determine their fluxes to coastal waters but related to complex physical and biogeochemical processes that cannot be fully represented by a single basin-scale retention rate. Practices of cropland and riparian managements for soil erosion and sediment loads can alter the response of erosional PN and PP fluxes to climate. Thus, the prediction of PN and PP loads to rivers must include the spatiotemporal evolution of these practices. Also, as climate and land use changes are interdependent, it is important to study their interactive effect on erosional N and P fluxes. Our future work will prioritize the representation of PN and PP transformations in the biogeochemistry component of the E3SM river model . These processes have already been included in a few river biogeochemical models and are shown to influence fluxes of bioavailable nutrients and phytoplankton growth in the water column significantly (Minaudo et al 2018, Vilmin et al 2020. We also intend to include a realistic representation of cropland erosion management in ELM-Erosion, such as the management of plant residues (Smets et al 2008), to improve understanding of landriver-ocean interactions in regions outside of the continental US.

Conclusion
This study presents a new erosional N and P model building on an existing continuous soil erosion module within ELM. The model is validated for the large rivers in the continental US. Using the model, we estimate that during 1991-2019 soil erosion yields 775 Gg yr −1 of PN and 328 Gg yr −1 of PP on average to the MARB and other rivers draining to the northern Gulf of Mexico. In particular, the MARB sub-basins such as the Ohio, Missouri and Upper Mississippi river basins that have intense agricultural activities contribute the most of these lateral PN and PP fluxes. We find that due to intensified extreme rains in the MARB erosional N and P fluxes in the drainage basins of the northern Gulf of Mexico increase at significant rates of 15 Gg N yr −1 and 6 Gg P yr −1 during 1991-2019. Based on the nutrient spiraling model with the simulated river water residence time from MOSART, we indicate that the intensification of erosional N and P fluxes has driven the significant increase of riverine PN and PP yields to the northern Gulf of Mexico in recent decades. Because extreme rains are projected to rise with warming and the majority of erosional P comes from bioavailable soil P pools, erosional nutrient fluxes in the region would likely continue to rise in the future and contribute more to riverine/estuarine nutrient enrichment. Therefore, it is important to enhance soil erosion control measures (e.g. reduced tillage, cover management and riparian buffer strips) to reduce the transport of P and N to the northern Gulf of Mexico.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).