Modification of exploration of long‐term nutrient trajectories for nitrogen (ELEMeNT-N) model to quantify legacy nitrogen dynamics in a typical watershed of eastern China

As a typical rainfall-dominated watershed, the Yongan watershed is located in the subtropical region of eastern China (120.23–121.01°E/28.47–29.04°N; elevation ∼15–1000 m; figure 1). Dominant land-use types are forest (∼68%) and agricultural (∼15%) lands. The Yongan River is the third largest river of Zhejiang province draining 2474 km² and flowing into the East China Sea with an average discharge of 72.9 m³ s⁻¹ at the downstream sampling site. Climate is subtropical monsoon with average annual precipitation and temperature of 1400 mm and 17.4 °C; rainfall mainly occurs in May–October. During 1980–2010, total population and poultry production (chicken and duck) increased ∼81% and 4.8-fold, respectively, whereas other domestic livestock production decreased by ∼25% (Chen et al 2015). Over the study period (1980–2010), paddy field (two crops per year) and dryland (tea, vegetables and fruits) area contributed 43%–70%, and 18%–48% of total agricultural area, respectively (Chen et al 2015). Since 2000, the watershed experienced a remarkable reduction (∼49%) in N fertilizer application (Hu et al 2018). Renovations to the old agricultural drainage system (stone/mud construction in 1950s) and construction of new drainage systems (cement channels/pipes) increased the area of drained agricultural land by ∼2-fold since 2000 (4800 ha in 2000 and 10 300 ha in 2010; Chen et al 2014a).

Zoom In Zoom Out Reset image size

ELEMeNT-N utilizes a coupled framework that integrates source-zone dynamics depicting SON accumulation/depletion processes within the source zone, and a travel-time model that describes hydrological processes for groundwater (figure 2). The model conceptualizes the landscape as a bundle of stream-tubes, each with a unique travel time to the watershed outlet, thus establishing a watershed travel-time distribution function f(τ). The travel-time distribution function controls N load trajectories (M_out(t)) at the watershed outlet:

$$\begin{align}{M_{{\text{out}}}}\left( t \right) & = \int\limits_{\text{0}}^{{\text{inf}}} {{J_{{\text{s,wshd}}}}} \left( {t - \tau } \right)f\left( \tau \right){{\text{e}}^{ - \gamma \tau }}{\text{d}}\tau \nonumber\\ & \quad + {\text{DMST}}{{\text{C}}_{{\text{sw}}}}(t)\end{align} \tag{ 1 }$$

Zoom In Zoom Out Reset image size

where J_s,wshd is annual nitrate flux leaching from source zone to groundwater; f(τ) represents the travel-time distribution (f(τ) = 1/μ × e^(−τ/μ); where τ is travel time and μ is mean travel time); γ is first-order rate constant describing N removal via denitrification along hydrologic pathways; and DMSTC_sw(t) is wastewater-N discharge to stream. More details of ELEMeNT-N are described in supporting information (SI, Part A).

SON mineralization rate is mainly mediated by microbial reactions (Gruber and Galloway 2008), which are affected by various factors. Temperature usually plays the most important role in well-drained, humid ecosystems (Francis et al 2007, Brookshire et al 2011). Laboratory/field experiments demonstrate that temperature alters soil microbial metabolic processes and rates, thus significantly altering SON dynamics (Dai et al 2020). Moreover, temperature data are relatively easy to access for most watersheds. Therefore, we employed mean annual temperature (MAT) to describe variations of SON mineralization rate for the active N pool (figure 2) in the form of an exponential function (Rey and Jarvis 2006, Brookshire et al 2011):

$$\begin{align}{k_{\text{a}}} = {k_{\min }}{{\text{e}}^{bT}},\end{align} \tag{ 2 }$$

where k_a is mineralization rate of the active organic N pool, ${k_{\min }}$ is basal N mineralization rate at 0 °C, T is MAT, and b is the sensitivity of k_a to temperature.

Data inputs for ELEMeNT-N include annual riverine N flux, SON content, water discharge, land use and the components contributing to the N surplus. Annual riverine N fluxes were estimated by the Weighted Regressions on Time, Discharge, and Season (WRTDS) model based on the 31 year record of total N concentrations and water discharge at the Yongan watershed outlet as calculated using EGRET package in the R framework (Hirsch et al 2010, Wu et al 2022). The WRTDS model excludes the influence of random temporal variability in streamflow to evaluate trends in riverine N flux (Hirsch et al 2010). SON content in 1980 was estimated from a soil organic carbon assessment (Liang et al 2019) assuming a C/N ratio of 14. Inorganic N fertilizer, manure N application, biological N fixation and wastewater discharge were directly obtained from local yearbooks. Atmospheric N deposition was obtained from published reports (Liu et al 2013). Landsat images were interpreted to acquire annual land-use distributions that were subsequently used to construct watershed land-use trajectories (www.gscloud.cn/). Mean groundwater travel time was derived from previous results obtained using stable isotopes of oxygen and hydrogen (Hu et al 2020).

The model was calibrated using the 31 year record of annual riverine total N fluxes. Nash–Sutcliff coefficient (NSE) and determination coefficient (R²) were used to evaluate goodness-of-fit. The Markov Chain Monte Carlo procedure with the Differential Evolution Adaptive Metropolis (DREAM-ZS) algorithm was used to calibrate parameters and determine uncertainties (SI, Part B). Then, we replaced the constant SON mineralization rate with the dynamic exponential function in the source code. A similar procedure was repeated to calibrate the newly introduced parameters (i.e. k_min and b). Global parameter sensitivity analysis was used to determine the parameters having the most significant influence on riverine N flux for the modified model. All model coding, sensitivity analysis and statistical analyses were conducted in MATLAB (Ver. 10.0, MathWorks, USA).

Model sensitivity analysis indicated that simulated riverine N fluxes were most sensitive to denitrification rates in human waste, soil and groundwater (λ_pop, λ_s and γ, table S1). Calibrated values for basal SON mineralization rate (k_min) and temperature sensitivity (b) were 0.04 (95% CIs: 0.01–0.39) and 0.55 (95% CIs: 0.02–1.30), respectively. Based on the calibrated parameters, the SON mineralization rate displayed considerable variations, ranging from 0.09 to 0.12 yr⁻¹ (average = 0.10 yr⁻¹, figure 3(e)). This indicates that using a fixed SON mineralization rate may impose a maximum 16.7% over/underestimate in a given year (i.e. (0.12–0.10)/0.12), highlighting the benefit of using a variable SON mineralization rate based on actual temperature variability.

Zoom In Zoom Out Reset image size

Compared to previous applications of ELEMeNT-N, some of the calibrated parameters presented large differences. Specifically, we found a lower protection coefficient for cropland (h_c: 0.14, 95% CIs: 0.01–0.27), lower soil denitrification rate (λ_s: 0.07, 95% CIs: 0.01–0.23), higher soil porosity (n: 0.64, 95% CIs: 0.54–0.70) and higher water content (s: 0.63, 95% CIs: 0.53–0.70; table 1). We attribute these differences primarily to the cropping frequency and types. In Europe and America, the cropping index is ∼105%, whereas it often exceeds 200% in this region (Liu et al 2021b). The higher cropping frequency leads to more frequent tillage resulting in greater soil disturbance, which hinders the formation of the protected SON pool and increases soil porosity (Morugan-Coronado et al 2020). Soil oxygen exchange with the atmosphere increases with higher soil porosity, thereby limiting soil denitrification, whereas higher soil water content can promote denitrification. In addition, denitrification is strongly influenced by soil organic matter (SOM) availability. The SOM content in this region is low, with ∼50% of cropland having SOM content below 20 g kg⁻¹ (Liang et al 2019). Therefore, the combination of low SOM content and high soil porosity may have a strong effect on denitrification, limiting soil denitrification rate. Our previous stable isotope studies also demonstrated weak denitrification in this region (Hu et al 2019). No site-specific characteristics were observed for SON mineralization rate in the protected pool (k_p: 2.4 × 10⁻⁴, 95% CIs: 2.0–4.4 × 10⁻⁴), protection coefficient for non-cultivated land (h_nc: 0.40, 95% CIs: 0.28–0.52), denitrification rate in groundwater (γ: 0.16, 95% CIs: 0.02–0.40) and human waste (λ_pop: 0.83, 95% CIs: 0.59–0.99, table 1).

The modified ELEMeNT-N model incorporating temperature-dependent SON mineralization achieved a high accuracy, accounting for 87% (R²) of the variability in annual riverine N flux and a NSE of 0.80 (figure 3(b)). These metrics can be considered as 'very good' according to watershed model performance ratings (R² >0.75 and NSE >0.75; Liew et al 2003, Moriasi et al 2007), and the efficiency metrics were comparable or superior to previous ELEMeNT-N applications (mean NSE: 0.76, ranges: 0.56–0.95). Compared to the original ELEMeNT-N (R²: 0.72 and NSE: 0.71, figure 3(a)), the most significant improvement occurred during the period when the riverine N flux irregularly changed (i.e. no regular increasing/decreasing trend during 1980–2002; figure 3(c)). These results verify the efficacy of the model modification and document the applicability of the modified ELEMeNT-N for use in the Yongan watershed. Despite the overall good performance metrics for the modified ELEMeNT-N, it should be noted that several simulation results deviated appreciably from the measured riverine N fluxes (figures 3(b) and (c)). These discrepancies likely result from the lack of data before 1980 that prevented determination of pre-1980 nitrogen accumulation in watershed soils (Yang et al 2020, Hu et al 2023).

During the study period (1980–2010), a total of 1630.0 kg N ha⁻¹ was applied to the Yongan watershed (SI, Part C, figure S1). The model results indicate that only 20.4% (12.4%–34.5%) of the cumulative N input was removed via riverine N export (figure 4). This fraction is consistent with 20%–25% determined by previous studies in this watershed (Chen et al 2014a, Hu et al 2018) and falls in the mid-range determined for several global watersheds (mean: 20%, range: 10%–40%; Han et al 2009, Howarth et al 2012, Swaney et al 2012). The majority of N input was removed via denitrification, accounting for 43.2% (range: 38.0%–64.4%; figure 4). This leave an estimated 36.4% (range: 13.2%–62.5%) of surplus N (593 kg N ha⁻¹) retained within the watershed as legacy N pools in SON (26.2%; range: 12.0%–42.4%) and groundwater (10.2%; range: 1.2%–20.1%; figure 4). These results are supported by a 65% increase of N in the upper 20 cm layer of agricultural soils between 1980 and 2009 (Chen et al 2014b) and a 3.8-fold increase in groundwater N concentrations between 1980 and 2019 in the Yongan watershed (Wu et al 2022). When compared to a previous study within this watershed, the magnitude of legacy N accumulation was much higher than the ∼480 kg N ha⁻¹ estimated by the ReNuMa model (Hu et al 2018). An important distinction is that the ReNuMa model assumed that a considerable proportion of N was retained in forest biomass (∼134 kg N ha⁻¹; Hu et al 2018), whereas ELEMeNT-N does not explicitly distinguish N accumulation in forest biomass from that in SON/groundwater. Therefore, the legacy N accumulation calculated by the modified ELEMeNT-N is generally consistent with that of ReNuMa when considering legacy total N accumulations in soil, groundwater and forest biomass.

Zoom In Zoom Out Reset image size

Modeled SON accumulation presented a ∼2.5-fold spatial variation (219–539 kg N ha⁻¹), while a smaller variation was observed for groundwater (143–188 kg N ha⁻¹; figures 5(a) and (b)). Specifically, SON accumulation in forest lands (average: 385 kg N ha⁻¹) was significantly lower than that in other land uses (410 kg N ha⁻¹). This can be attributed to the lower N inputs and higher denitrification rate in forest soils (Wang et al 2018). Monitored groundwater N concentrations and simulated groundwater N accumulation showed a strong positive correlation (figure 5(a), R² = 0.99, p< 0.01). A recent study documented a considerable amount of nitrate accumulation in deeper soil horizons and vadose zone of subtropical red soils (Yang et al 2020). These results indicate that the deep soil/vadose zone can function as a transient N storage zone, absorbing N leaching from upper soil horizons and thereby retarding its transport to groundwater. We attribute the smaller spatial variations in groundwater N accumulation to this transient storage in the deeper soil/vadose zone.

Zoom In Zoom Out Reset image size

As expected, the spatial distribution of legacy N accumulation was closely related to the magnitude of surplus N (figure 6(c)). Specifically, the correlation between total legacy N and surplus N was strongest (figure 6(c), R² = 0.65, p< 0.01), followed by SON (figure 6(c), R² = 0.52, p< 0.01), and weakest for groundwater N (figure 6(c), R² = 0.47, p< 0.05). The direct groundwater N source is soil N leaching rather than N inputs to the watershed (Gu et al 2013). Therefore, the accumulation of groundwater N is also regulated by SON mineralization rate and soil permeability/leaching dynamics (Wu et al 2016, Liu et al 2020).

Zoom In Zoom Out Reset image size

Both SON and groundwater were important N sinks throughout the study period. The accumulation rate of legacy N storage displayed a continuous decreasing trend (figure 6(a)), with a rate of −0.32 kg N ha⁻¹ yr⁻¹. The trend for SON accumulation can be divided into two stages: a relatively constant (16.5 ± 0.6 kg N ha⁻¹) stage during 1980–1999 versus a decreasing trend from 2000 to 2010 (from 10.5 to 8.3 kg N ha⁻¹, figure 6(a)). We assign two primary reasons for these results. The primary reason is the increasing temperature in recent years (figure 3(d)), which leads to a higher SON mineralization rate. Meanwhile, the N inputs decreased during these periods (figure S1). As such, SON accumulation decreased from 2000 to 2010. Another reason may derive from an increase in the area of drained agricultural land. Correlation analysis indicated that drained agricultural land area was negatively related to both legacy N and SON accumulation (figure 6(b)). Increases in drained agricultural land result in faster routing of precipitation and irrigation water to streams, thereby reducing SON accumulation (Blann et al 2009, Basu et al 2011). Theoretically, enhanced drainage systems can also lead to faster transport of subsurface N to streams and reduced accumulation in groundwater (Blann et al 2009, Buzek et al 2009). However, no significant trend was observed in groundwater N accumulation over the study period, with a mean (±SD) accumulation rate of 5.0 ± 2.7 kg N yr⁻¹. This may be ascribed to the accumulation of nitrate in the deeper soil/vadose zone (Yang et al 2020), which may provide an additional N supply to groundwater as soil moisture conditions change (Wu et al 2016).

The ELEMeNT-N modeled lag time between watershed N inputs and riverine export in the Yongan watershed was 11.9 years (95% CIs: 8.3–21.3 years). This result was consistent with the decadal lag time determined by previous studies in this watershed (Chen et al 2014a, 2015, Hu et al 2018). The majority of anthropogenic N input to a watershed is transiently stored in as SON before subsequently being released to groundwater through mineralization and leaching (Spoelstra et al 2001, Sebilo et al 2013). In terms of ELEMeNT-N results, the residence time for SON was 9.7 years (95% CIs: 6.3–18.4) in the active pool versus 4281 years (95% CIs: 2276–35 981) in the protected pool. The simulated hydrologic lag time was 2.2 years (95% CIs: 2.0–3.0), which was lower than the mean groundwater transit time (4.5 years, 95% CIs: 4.0–5.5) determined by stable isotopes (Hu et al 2020). We posit two primary causes for this groundwater hydrology deviation. First, current estimation methods for groundwater transit time are subject to large uncertainties (Yang et al 2018, Rodriguez and Klaus 2019). Thus, the deviation may be introduced by differences in characterizing transit times between the two methods. Additionally, mean travel time of legacy N in the modified ELEMeNT-N includes the surface runoff component, which has a much lower travel time than groundwater.

ELEMeNT-N indicated that riverine N load at the watershed outlet was dominated by legacy N. Current year's N application (including point and non-point sources) contributed only 11.8% (range: 6.3%–20.2%) of the riverine N flux and legacy N contributed the remaining 88.2% (79.8%–93.7%). These results were consistent with those determined by a dual stable nitrate isotopes approach, in which the N source apportionment indicated that >80% of riverine N originated from legacy N pools (i.e. soil + groundwater; Hu et al 2019).

Due to N removal via denitrification, biomass uptake/storage and riverine output, the contribution of previous years' N applications decreased with increasing time. Anthropogenic N applied 2–5 years, 6–10 years and >10 years ago contributed 35.2% (range: 25.4%–83.7%), 30.6% (range: 8.1%–34.0%) and 22.3% (range: 0%–44.5%) of the total riverine N load (figure 7(a)). Notably, limited by the model structure and data availability, N aged >10 years consisted of two components: N inputs from more than 10 years ago and the natural background N pool. This result was similar to previous estimates (Chen et al 2014a, 2015), in which the contribution of the current year's net anthropogenic N inputs was ∼10%, whereas the legacy and natural background N contributed ∼75% and ∼15%, respectively.

Zoom In Zoom Out Reset image size

Riverine export of current year's N applications exhibited contrasting trends for different time periods, namely a continuous increase during 1980–1999 followed by no significant trend during 2000–2010 (figure 7(a)). When compared to annual surplus N inputs (SI, Part C, figure S1), an obvious nonlinearity was observed (figure 7(b)). In subtropical watersheds, 10.7 kg N ha⁻¹ yr⁻¹ was suggested as a threshold value for N saturation (Howarth et al 2012, Swaney et al 2012). During the study period, observed N inputs far exceeded this threshold value, leading to a progressive N saturation. Therefore, the capacity of soil, vegetation and groundwater to transiently store N declines as N saturation progresses. The decrease in N inputs may be offset by this reduction in N storage efficiency, thereby resulting in a relatively stable riverine export of current year's N application (figure 7(b)). In addition, N aged 2–5 and 6–10 years began to decrease after their long-term continuous increase in the earlier portion of the study period. The transition occurred 4 and 7 years after the decrease in surplus N inputs (∼2000) (figure 7(a)). These results confirm those of previous studies that found reducing N inputs is essential for reducing the impact of legacy N on mitigating N pollution.

The modified ELEMeNT-N provides quantitative information for watershed N legacy effects (figures 4, 5 and 7(a)), which are critical considerations for managers designing science-based, N pollution control strategies and assessing their efficacy. For example, the spatial distribution of legacy N provided by ELEMeNT-N (figure 5) allows decision makers to prioritize nutrient management strategies to critical regions. The model predicted that it could take at least 15–21 years following a 50%–100% reduction of N surplus to meet the water quality goal for the Yongan River (1 mg N L⁻¹) (figure S3). Therefore, it is necessary to set appropriate expectations for the public in meeting long-term pollution targets (Bouraoui and Grizzetti 2014). Due to the large legacy N contributions to riverine N flux, targeted strategies should be adopted to proactively prevent legacy N accumulation. Reducing N inputs (e.g. fertilizer, manure) is the most direct approach, as they are a direct source of the legacy N (figure 6(c)). Construction of agricultural drainage systems is a necessary and effective measure to limit legacy N accumulation as well. It can lead to shorter watershed travel time, thus effectively reducing legacy N accumulation, especially as SON (figure 6(b); Blann et al 2009, Sloan et al 2017). However, drainage systems may result in more immediate N pollution delivery to rivers. Therefore, additional strategies targeting drainage water treatment (e.g. buffer zones between drainage and rivers) are warranted to avoid short-term water quality impairments (Mander et al 2005, Hu et al 2019).