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

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:

$$\begin{equation}{\text{P}}{{\text{N}}_{\text{c}}} = 0.116\,{\text{PO}}{{\text{C}}_{\text{c}}} - 0.019,\end{equation} \tag{ 1 }$$

$$\begin{equation}{\text{P}}{{\text{P}}_{\text{m}}} = \frac{{{\text{POC}}_{\text{m}}^{1.002}}}{{22.1536}},\end{equation} \tag{ 2 }$$

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⁻¹) 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):

$$\begin{equation}{\text{P}}{{\text{N}}_o} = \sum {{\text{P}}{{\text{N}}_i}{{\text{e}}^{ - {\lambda _{{\text{PN}}}}{\tau _i}}},} \end{equation} \tag{ 3 }$$

$$\begin{equation}{\text{P}}{{\text{P}}_{\text{o}}} = \sum {{\text{P}}{{\text{P}}_i}{{\text{e}}^{ - {\lambda _{{\text{PP}}}}{\tau _i}}},} \end{equation} \tag{ 4 }$$

where PN_i and PP_i are erosional PN and PP fluxes from a grid cell in a river basin (kg yr⁻¹), respectively, PN_o and PP_o are erosional PN and PP yields from the river basin (kg yr⁻¹), 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⁻¹), 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 (Li et al 2013). 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.

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 and Moody 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 state-level 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.

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²) and Nash–Sutcliffe coefficient.

Zoom In Zoom Out Reset image size

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 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 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 Gg= 1 × 10⁹ g) (Mann–Kendall test, p < 0.05) and the estimated trend is 1.7 Gg P yr⁻¹ (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.

Zoom In Zoom Out Reset image size

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⁻¹ and 302 ± 105 Gg P yr⁻¹, 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.

Zoom In Zoom Out Reset image size

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 post-dam PN yield from the MARB during 1980–1996 to be about 210 Gg N yr⁻¹. Meanwhile, the pre-dam PN yield from the MARB is about 600 Gg N yr⁻¹, which is based on the assumption of 0.15% N content in suspended sediment (Goolsby et al 2000) and 4 × 10⁵ Gg of suspended sediment yield in the pre-dam condition (Milliman and Syvitski 1992). Our modeled erosional N flux of 718 ± 102 Gg yr⁻¹ 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⁻¹) 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⁻¹) is of the same order of magnitude as the river measurements (22–42 Gg N yr⁻¹) 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⁻¹ during 1980–1996) (Goolsby et al 1999) and the estimated pre-dam riverine PP yield (about 415 Gg P yr⁻¹) (Meybeck 1982, Milliman and Syvitski 1992).

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⁻¹ (Mann–Kendall test, p < 0.05), which is 2.1 Gg P yr⁻¹ (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⁻¹ of PN flux and 328 Gg yr⁻¹ 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⁻¹ and 6 Gg P yr⁻¹, 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 and Fennel 2017).

Zoom In Zoom Out Reset image size

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⁻¹ and 8.4 Gg P yr⁻¹, 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 runoff generation (Gonzalez-Hidalgo et al 2010), the mild increase of extreme rains causes a much larger increase of surface runoff extremes (figure 5(h)).

Zoom In Zoom Out Reset image size

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⁻¹ (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⁻² d⁻¹ was absent during 1991–1995 but starts emerging during 2015–2019.

Zoom In Zoom Out Reset image size