Long-Term Nitrogen and Phosphorus Outflow from an Instream Constructed Wetland under Precipitation Variability

Abstract

In many agricultural watersheds, surface runoff often causes unwanted nitrogen (N) and phosphorus (P) losses from croplands into stream networks. When this phenomenon is pronounced, it significantly changes N and P concentrations in streams affecting aquatic ecosystems. To protect stream water quality, the installation of instream-constructed wetlands (ICWs) for treating runoff water is often reported as a low-cost alternative to conventional water treatment systems. Indeed, ICWs have the capacity to collect and temporarily retain nutrients transported from agricultural landscapes and then slowly release them into downstream networks. However, the long-term hydrologic behavior of ICWs relative to N and P outflow control is still insufficiently reported. Especially in the context of climate change, it is relevant to investigate the effect of precipitation variability on ICWs N and P outflow. This study uses the soil and water assessment tool (SWAT) model to approximate the long-term hydrologic behavior of an experimental ICW installed in a small agricultural watershed. The model was set assuming a continuous corn and soybean rotation on croplands, then a multidecadal (period 2001–2020) simulation was used to evaluate the implication of precipitation variability on total nitrogen (TN), nitrate-N (NO3-N), total P (TP), and dissolved P (DP) outflows. Results show meaningful changes in the precipitation pattern with contrasting effects on N and P outflows. While analyses show significant trends in the maximum monthly precipitation, nutrient outflows during two consecutive decades, 2001–2010 and 2011–2020, show increases of 46% for TN, and 82% for TP. At the watershed scale, month-to- month TN and TP outflows range from 24 to 810 kg N and 26 to 1358 kg P during 2011–2020, compared with 42 to 398 kg N and 40 to 566 kg P during 2001–2010. The increase in nutrient outflow is particularly pronounced for TP and DP which show significant trends and high correlations (r > 0.70) with maximum monthly precipitation. An exception is nitrate-N outflow, which counts on average for less than 5% of TN outflow but appears more affected by the timing of N fertilization in the watershed.

1. Introduction

Worldwide, nitrogen (N) and phosphorus (P) are by far the most used fertilizers in agriculture. However, the use of these fertilizers is a cause of unwanted N and P losses towards neighboring streams affecting water quality [1,2]. With the growing global food demand, the use of N and P in crop production is expected to keep increasing to sustain high yields and profitable agriculture [3,4]. This tendency exacerbates water quality protection challenges in many agricultural watersheds across the United States (US). For water quality conservation, instream constructed wetlands (ICWs) are helpful for protecting surface water bodies from diffuse pollution sources. In wastewater management, ICWs are often considered low-cost alternatives to conventional treatment systems [5]. When installed in a water pathway, an ICW functions as a sink that traps nutrients transported from agricultural landscapes through surface runoff [5]. In the water, N and P nutrients exist in various compound forms (e.g., ammonium $N{H}_4^{+}$, nitrite $N{O}_2^{-}$, nitrate $N{O}_3^{-}$, urea $CN{H}_{4}O$ for N; $H_{2}P{O}_4^{-}$, $HP{O}_4^{2-}$, $P{O}_4^{3-}$ for P) [2,6]. As nutrients accumulate in the ICW, N and P are attached to dead tissues, sediments, and precipitated materials often deposit at the bottom of the ICW pond. The ICW slowly releases dissolved nutrients with the outflow water into the downstream network. Given the ICW outflow rates are low (except in the case of breaching), the amount of nutrients released is substantially lowered thereby reducing the risk of water quality impairment. In the ICW, nutrients are subjected to pH, temperature, and various biophysical activities (e.g., microbial, algal, and plant nutrient uptake) the interplay of which affects nutrient balance in water [6]. However, a significant concern related to ICWs is the potential disturbances caused by exceptional precipitation events. Indeed, heavy precipitation events generate abnormal runoff and water currents that disrupt the regular water and nutrient outflows patterns. As a result, substantial quantities of nutrients could instantaneously be released out of the ICWs exposing downstream water to abnormal nutrient loads. Even though this phenomenon is noted in previous studies [7,8,9], its magnitude in relation to precipitation variability is not sufficiently documented, particularly regarding N and P contamination. Yet, long-term ICWs experimental studies are needed to fully appraise the effect of precipitation variability on the ICW hydrologic behavior. Moreover, projections of precipitation patterns under climate change suggest major perturbations which are often translated as an increase in extreme events [10,11,12]. Unfortunately, long-term ICW experimental studies are seldom reported due to costs and labor requirements. In such situations, a process model such as the soil and water assessment tool (SWAT), can be used to approximate the long-term behavior of an ICW. Hence, the present study aims to use SWAT model simulations to examine the long-term hydrological behavior of an ICW relative to N and P outflow under precipitation variability.

Widely used in water resources and environmental assessments, SWAT allows modelers to adjust model parameters and configurations to simulate hydrological processes and hydrological behaviors of watersheds and reservoirs. In this study, data from an experimental ICW are used in the SWAT model to reproduce and comparatively analyze N and P outflow over a multi-decadal precipitation regime. The outcomes of the study are structured in five sections including this introduction, the data and method, the results, the discussion, and the conclusion.

2. Data and Method

2.1. Data

The study uses measurements of N and P outflows along with recorded discharges from an experimental watershed namely the Herrings Marsh Run (HMR). The HMR is a small watershed draining runoff water from approximately 409 ha in the Northeast Cape Fear River Basin (US Geological Survey’s Hydrologic Unit Code 03030007) in North Carolina. As shown in Figure 1, the HMR watershed, which is highly agricultural with 66% of its area as croplands, includes an experimental ICW which has been monitored by the US Department of Agriculture’s Agricultural Research Service (USDA-ARS). Descriptions of the ICW parameters and the SWAT modeling procedure are reported by Sohoulande et al. [13]. The ICW’s N and P outflows measurements during the period 1993 to 1999, and historical discharge measurements of the stream-gage USGS-0210783240 (latitude 35.10°, longitude −77.93°, elevation 28.96 m) are used in the present study to calibrate and validate the SWAT model. Description of the ICW outflow sampling and nutrient concentration measurement procedures are reported in previous manuscripts [9,14,15]. Long-term historical climate data, including daily precipitation, maximum, and minimum air temperature over the period 1991 to 2020 were input in SWAT. This allows a model warming period (1991–1992), a calibration period (1993–1996), a validation period (1997–1999), and a multi-decadal simulation period (2001–2020). The climate data were retrieved from the National Oceanic and Atmospheric Administration (NOAA) stations USC00311881 (latitude 35.02, longitude −78.28), USC00314684 (latitude 35.19, longitude −77.54), and USW00013713 (latitude 35.34 and longitude −77.96). Geospatial data used in SWAT include a 10 m digital elevation model (DEM) obtained from the US Geological Survey (USGS) database, a 10 m gridded soil map obtained from the USDA-Natural Resources Conservation Service’s Gridded Soil Survey Geographic (gSSURGO) database, and a 30 m gridded map of the year 2011 land use land cover (LULC) obtained from the Multi-Resolution Land Characteristics Consortium database [16]. In a previous study, Sohoulande et al. [13] analyzed land use changes in the watershed based on 2001, 2006, and 2011 LULC data, then showed there was no major LULC change across the watershed justifying the use of 2011 LULC as a SWAT input.

2.2. Setting SWAT Model: Crop Management, ICW Inputs, and Performance Indicators

The HMR watershed was subdivided into 13 sub-basins (Figure 2) and the crop management scenario used in SWAT assumed a continuous rotation of corn-soybean. Indeed, the USDA’s National Agricultural Statistic Service (USDA-NASS) county-level statistics of crops show corn and soybean are the major crops grown in the area during the last two decades as they both represented more than 72% of the annual acreage of raw crops harvested [17]. To balance the spatial representation of corn and soybean in the HMR watershed, the 13 sub-basins were split into two groups (see Figure 2). The management and operations scheduling used in the SWAT model reflect common extension service guidelines for the US Coastal Plain [18,19]. The management operations include a split application of N and P fertilizers on corn crops with a total rate of 200 kg N/ha (42% at planting in April and 58% in June) and 56 kg P/ha (50% at planting in April and 50% in June). For the soybean, only a split application of P is considered at a total rate of 56 kg P/ha (50% at planting in May and 50% in June). The crop management assumes that there was no tillage and that fertilizers are surface broadcasted. The corn growing season started in April and terminated in August, while the soybean growing season started in May and terminated in November. Grains are harvested while the rest of the crop biomass remains on croplands.

In SWAT, the ICW was represented as a reservoir with a surface area (RES_ESA) of 4.52 ha, and a volume (RES_EVOL) of 25,400 m³ when filled at emergency spillway; a surface area (RES_PSA) of 3.4 ha, and a volume (RES_PVOL) of 6000 m³ when filled at the principal spillway. The ICW location in subbasin 10 (Figure 1) is considered for setting the reservoir inputs. Additional water quality and quantity parameters used to set the ICW in SWAT are reported by Sohoulande et al. [13]. Examples of these parameters are the ICW hydraulic conductivity, phosphorus settling rate in ICW, and initial concentration of DP. The SWAT model was calibrated and validated for monthly watershed discharge, ICW’s nitrate-N (i.e., NO3-N), and dissolved phosphorus (DP) outflow. Precisely, the time slices 1993 to 1996 and 1997 to 1999 were respectively used for the model calibration and validation. The SWAT calibration and uncertainty programs (SWAT-CUP) [20] were first used to narrow down sensitive parameters and ranges which were later manually adjusted to enhance the model simulations. Aligning with previous studies [21,22], targeted SWAT parameters, including soil available water capacity (SOL_AWC), runoff curve number at moisture II (CN2), saturated hydraulic conductivity (SOL_K), denitrification threshold water content (SDNCO), denitrification exponential rate coefficient (CDN), nitrate percolation coefficient (NPERCO), nitrogen uptake distribution parameter (N_UPDIS), deep aquifer percolation fraction (RCHRG_DP), phosphorus uptake distribution parameter (P_UPDIS), phosphorus percolation coefficient (PPERCO), phosphorus soil partitioning coefficient (PHOSKD), phosphorus availability index (PSP), residue decomposition coefficient (RSDCO), a fraction of algal biomass that is Phosphorus (AI2), and plant uptake compensation (EPCO).

The model performance during calibration and validation were analyzed using four different indicators including the Nash-Sutcliffe’s Efficiency NSE (Equation (1)), the index of agreement d₁ (Equation (2)), the coefficient of determination R², and the Root Mean Squared Error RMSE (Equation (3)) [21,23]. Once the model was calibrated and validated, two-decade (2001 to 2020) model simulations were considered for a long-term ICW’s N and P outflow analysis in relation to precipitation variability.

$$NSE=\hspace{0.17em}1\hspace{0.17em}\hspace{0.17em}-\hspace{0.17em}\hspace{0.17em}\frac{{\displaystyle \sum _{i=1}^n{(q_i-q_i^{\prime })}^2}}{{\displaystyle \sum _{i=1}^n{(q_i-\overline{q})}^2}}$$

(1)

$$d_1=\hspace{0.17em}1\hspace{0.17em}\hspace{0.17em}-\hspace{0.17em}\hspace{0.17em}\frac{{\displaystyle \sum _{i=1}^n\left|q_i^{\prime }\hspace{0.17em}-\hspace{0.17em}{q}_i\right|}}{{\displaystyle \sum _{i=1}^n\left(\left|q_i^{\prime }\hspace{0.17em}-\hspace{0.17em}\overline{q}\right|\hspace{0.17em}+\hspace{0.17em}\left|q_i\hspace{0.17em}-\hspace{0.17em}\overline{q}\right|\right)}}$$

(2)

$$RMSE={\left[n^{-1}{\displaystyle \sum _{i=1}^n{\left(q_i-q_i^{\prime }\right)}^2}\right]}^{0.5}$$

(3)

In Equations (1)–(3), q_i and q_i′ represent respectively the observed and simulated values at the ith month with 1 ≤ i ≤ n and n representing the number of months in each period; $\overline{q}$ is the average observation of the period.

3. Results

3.1. SWAT Model Performance

The SWAT model was calibrated and validated for monthly discharge at the watershed outlet and then for nitrate-N and DP outflow from the ICW. Evaluation of the model performance during the calibration and validation is summarized in

Table 1. SWAT model performance for streamflow, nitrate-N, and dissolved phosphorus outflow from the instream constructed wetland.

Efficiency Criteria	Discharge (m3/s)		Dissolved Phosphorus Outflow (kg P)		Nitrate-N Outflow (kg N)
	Calibration	Validation	Calibration	Validation	Calibration	Validation
NSE	0.70 *	0.83	0.46	0.38	0.80	0.61
d1	0.67 *	0.71	0.66	0.76	0.69	0.65
R2	0.80 *	0.90	0.48	0.52	0.81	0.76
RMSE	13.03 *	0.05	48.80	17.35	0.50	0.57

* Discharge calibration is for the Northeast Cape Fear River basin, an explicit description is reported by Sohoulande et al. [13].

. Indeed, the efficiency criteria values for discharge (e.g., 0.70 ≤ NSE ≤ 0.83, 0.67 ≤ d₁ ≤ 0.71), DP outflow (e.g., 0.38 ≤ NSE ≤ 0.46, 0.66 ≤ d₁ ≤ 0.76), and nitrate-N outflow (e.g., 0.61 ≤ NSE ≤ 0.80, 0.65 ≤ d₁ ≤ 0.69) were acceptable SWAT simulations based on guidelines for hydrological model calibration and validation [21,24,25]. The model was henceforth used to simulate monthly DP, total P (TP), nitrate-N, and total N (TN) over two consecutive decades, 2001–2010 and 2011–2020. For these simulations, recorded daily precipitation and temperature were used in SWAT. Graphs in Figure 3 comparatively present the average fluctuations of TN, nitrate-N, TP, and DP during different months of the year. At the first glance, these graphs show two major contrasts. The first contrast is related to the discrepancies between the two decades 2001–2010 and 2011–2020. Indeed, outflow values during 2011–2020 are higher than during 2001–2010 for all nutrients. The average annual outflow values of TN, nitrate-N, TP, and DP respectively increased by 46%, 25%, 82%, and 108% during the decade 2011–2020 compared to 2001–2010. However, these overall average increases also hide critical intra-annual discrepancies. For instance, at the watershed scale, month-to-month TN and TP outflows range from 24 to 810 kg N and 26 to 1358 kg P during 2011–2020, compared with 42 to 398 kg N and 40 to 566 kg P during 2001–2010. The second contrast is nitrate-N outflow (Figure 3b) has a maximum peak in April while TN, TP, and DP outflow (Figure 3a,c,d) have maximum peaks in September. Meanwhile, average month-to-month precipitation variability also has its peak in September. Given precipitation is the major dynamic model input throughout the entire simulation period 2001–2020, its effect on the nutrient outflow can be reasonably estimated. Yet the eventuality of precipitation variability effects on ICW nutrient outflows does not rule out factors inherent to the model. As an example, running the model over a long period assumes a cumulative nutrient build-up in the reservoir which will be discussed later in the following sections.

3.2. Signals of Precipitation Variability and ICW’s N and P Outflow

As mentioned earlier, actual daily precipitation data have been used in the SWAT model to simulate ICW outflows over two consecutive decades, 2001–2010 and 2011–2020, respectively. Hence, precipitation variability is a major dynamic input for the model. The setting allows an appraisal of precipitation effects on the long-term capacity of the ICW at controlling P and N outflow toward the downstream network. Figure 4 presents an overview of precipitation patterns in the watershed by reporting not only the relationship between annual maximum monthly precipitation (MaxMonPREC) and annual total precipitation (AnnPREC), but also the pattern of MaxMonPREC anomalies during the past two decades (2001–2020). The linear relationship (R² = 0.74) in Figure 4a indicates that MaxMonPREC explains 74% of the variance contained in AnnPREC. Hence, the precipitation pattern in the watershed is significantly driven by MaxMonPREC. In Figure 4b, the MaxMonPREC anomaly of individual year y [$MaxMonPREC\hspace{0.17em}\hspace{0.17em}anomaly\hspace{0.17em}(y)\hspace{0.17em}$] is estimated relatively to the average MaxMonPREC ($\overline{MaxMonPREC}\hspace{0.17em}$) of the entire period 2001–2020 as given by Equations (4) and (5).

$$MaxMonPREC\hspace{0.17em}\hspace{0.17em}anomaly\hspace{0.17em}(y)\hspace{0.17em}=\hspace{0.17em}MaxMonPRE{C}_y\hspace{0.17em}-\hspace{0.17em}\overline{MaxMonPREC}$$

(4)

$$\overline{MaxMonPREC}\hspace{0.17em}\hspace{0.17em}=\hspace{0.17em}\hspace{0.17em}\frac{1}{20}{\displaystyle \sum _{y=2001}^{2020}MaxMonPRE{C}_y}$$

(5)

An analysis of Figure 4b shows the frequency of positive MaxMonPREC anomalies shifted from 20% in 2001–2010 to 60% in 2011–2020. Explicitly, only two years (i.e., 2003 and 2006) show positive MaxMonPREC anomalies during the decade 2001–2010 while six years show positive MaxMonPREC anomalies during the decade 2011–2020. The positive anomalies indicate the MaxMonPREC values are above the long-term average suggesting changes in the regular precipitation distribution during the year and an occurrence of extreme precipitation events. In alignment with the study objective, this sub-section examines how the observed precipitation signals translate into ICW N and P outflow. Hence, Figure 5 presents the annual variation of TN (AnnTN), and TP (AnnTP) outflows in respective comparison with monthly maximum TN (MaxMonTN), and TP (MaxMonTP). Both graphs in Figure 5a,b show high TN outflow years have high TP outflows, and low TN outflow years have low TP outflows. This suggests strong similarities in annual nutrient outflow patterns. However, these annual patterns may hide meaningful intra-annual signals as the analysis of Figure 3 shows the maximum peak of nitrate-N outflow occurs in April while the maximum peaks of TN, TP, and DP occur in September. Interestingly, the peak of precipitation in the watershed also occurs in September as shown by the boxplots in Figure 6. Precipitation distribution during the year is uneven and the large interquartile range of September’s boxplot in Figure 6 indicates high precipitation variability.

3.3. Effect of Precipitation Variability on ICW’s N and P Outflow

Over the period 2001–2020, the observed shift of positive MaxMonPREC anomaly frequency (Figure 4a) is likely a signal of an ongoing change in precipitation patterns across the study region. Indeed, during the decade 2011–2020, the Southeast Coastal region of the US recorded unusual extreme precipitation events which could explain changes in the precipitation pattern [26,27]. However, Kendall’s Tau statistics in

Table 2. Trends in annual total and annual maximum monthly precipitation, TN, NO3-N, TP, and DP outflow over the period 2001 to 2020.

Variable	Annual Total		Maximum Monthly
	Kendall’s Tau	p-Value	Kendall’s Tau	p-Value
Precipitation	0.19	0.26	0.33 *	0.04
TN outflow	0.21	0.21	0.25	0.13
NO3-N outflow	0.18	0.29	0.17	0.32
TP outflow	0.39 *	0.02	0.39 **	0.02
DP outflow	0.47 **	0.00	0.51 **	0.00

* Indicates a significant statistic at p-value = 0.05; ** indicates a significant statistic at p-value = 0.01.

show trend signals are significant with MaxMonPREC but not with AnnPREC. This suggests the complexity of the changes affecting precipitation in the region given the sole consideration of annual total precipitation is not sufficient to evidence changes. Yet the shift of MaxMonPREC anomalies outlined in Figure 4a and the positive trend of MaxMonPREC reported in Table 2, sustain an overall increase of extreme precipitation events during the decade 2011–2020. Table 2 also reports trend analyses of ICW nutrient outflows including annual and maximum monthly TN, nitrate-N, TP, and DP over the period 2001 to 2020. Yet, significant trends are only observed with AnnTP, AnnDP, MaxMonTP, and MaxMonDP. Figure 7 further reports correlation gradients between precipitation and ICW nutrient outflow variables. The precipitation variable MaxMonPREC, is highly correlated with AnnTP (r = 0.76), AnnDP (0.76), MaxMonTP (r = 0.69), and MaxMonDP (r = 0.71). However, the correlations between MaxMonPREC and AnnNO3 or MaxMonNO3 are insignificant. These results substantiate the tendencies outlined in Table 2. However, high correlations are also noted between AnnTN and the ICW’s P outflow variables (0.92 ≤ r ≤ 0.99), but the correlations with AnnNO3 and MaxMonNO3 remain insignificant. Note, the Pearson correlation thresholds for p-value < 0.05 and p-value < 0.01 with df = 19 are 0.43 and 0.55 respectively. On average the annual nitrate-N outflow counts for less than 5% of the TN outflow. Hence, monthly nitrate-N outflow variations are likely outweighed by the overall TN fluctuations. This could explain in part the contrast between month-to-month TN and nitrate-N outflow observed in Figure 3a,b.

4. Synthesis and Discussion

This study uses SWAT simulations to analyze an ICW’s hydrological behavior under a changing precipitation pattern. The modeling approach assumes that the essential parts of N and P originate from croplands, considered the primary non-point sources of nutrients in most agricultural landscapes [28]. ICWs are generally valued as alternative solutions for mitigating the impact of diffuse sources of nutrients on stream water quality. Indeed, ICWs serve as buffer ponds that accumulate and temporarily retain nutrients washed-off from croplands. As nutrients are washed-off from croplands, N and P constituents are transported into the ICW in dissolved and undissolved forms. Substantial fractions of the undissolved N and P are attached to organic tissues, sediments, and precipitated materials which often accumulate at the bottom of the ICW pond. The dissolved nutrients are slowly released from the ICW into the downstream network. However, the behavior of N and P in the ICW is different [29,30]. For instance, P attached to organic tissues and inorganic particles frequently sink and accumulate at the bottom of the pond. Eventually, the accumulated P is flushed by water currents caused by extreme precipitation events. In steady flow conditions, however, DP is slowly released from the ICW outlet at a rate that generally does not cause a significant rise in P concentrations downstream. Unlike P, N is unstable and its fate in the ICW is quite complex. In water, N can be present in various organic and inorganic compound forms such as urea, amino acids, ammonium ($N{H}_4^{+}$), nitrite ($N{O}_2^{-}$), nitrate ($N{O}_3^{-}$), nitrous oxide ($N_{2}O$), ammonia ($N{H}_3$), nitrogen gas ($N_2$), and nitrogen dioxide ($N{O}_2$). Gaseous forms of N (e.g., $N_2$, $N_{2}O$, $N{O}_2$, $N{H}_3$) are emitted from the water surface into the atmosphere. The magnitude of this gaseous emission is highly dependent on the interplay between N constituents in the water and biophysical factors such as microbial activities (i.e., denitrification), pH, temperature, and wind. Of the N constituents in the water, nitrate-N is the essential plant growth nutrient [2,6], which in excess can trigger algal blooms causing eutrophication. In addition, nitrate in drinking water is subjected to regulatory criteria (e.g., in the U.S. the allowed nitrate concentration for drinking water supplies is 10 mg/L) due to subsequent nitrite toxicity. As shown in the results, nitrate-N only counts for a fraction (5%) of the ICW’s TN outflow. Yet, a separate analysis of TN and nitrate-N shows the peak nitrate-N outflow occurs independently of the regular TN and P outflow pattern. Indeed, the average month-to-month nitrate-N outflow, particularly the peak observed in April (Figure 3b) has a relationship with the crop management scenario which specified N fertilization at the beginning of the corn season in April. Though the early application of N is still recommended for corn production, a recent study has highlighted discrepancies between available soil N following N application and actual plant N need causing substantial nitrate loss from soils [31]. Hence, the average nitrate pattern portrayed in Figure 3b is likely driven by the crop management scenario used.

The strong similarities between TN, TP, and DP outflows (Figure 3a,c,d) are likely related to precipitation as the peaks of TN, TP, and DP occur in September which also corresponds to the precipitation peak (Figure 6). In fact, unusual precipitation events have been recorded in the study region during the Atlantic hurricane seasons encompassing September [26,27]. Even though the present study did not focus on extreme events such as hurricanes, results have outlined critical changes in precipitation patterns over the period 2001–2020. These changes are statistically evidenced through the MaxMonPREC time series while AnnPREC data did not show significant trends. Trends and correlation analyses show a significant connection between the MaxMonPREC and P outflow. Under normal conditions, the ICW contributes to P outflow abatement thereby protecting downstream water from excess nutrient loads. This P outflow abatement is accomplished through sedimentation, microbial interactions, and uptake by vegetation [5]. In the ICW, P compounds are mainly found as dissolved P, particulate mineral P, or particulate organic P [6]. However, infiltration rates at the bottom of wetlands are generally poor causing a cumulative P build-up in the reservoir over the years [29,30]. Indeed, precipitation events generate runoff which transports substantial amounts of nutrients, organic materials (e.g., dead tissues), and sediments into the ICW. Exceptionally, extreme precipitation events generate water currents which flush out material accumulated at the bottom of the ICW causing unusual releases of P into downstream networks [7,8,9]. This phenomenon is well represented in the long-term TN, TP, and DP outflows simulations.

As the SWAT model was run over the continuous period 2001–2020, the amounts of organic materials, precipitated compounds, and sediments built in the ICW during the decade 2011–2020 are expected to be higher than the amount built-up during the preceding decade 2001–2010. This somewhat explains the discrepancy in nutrient outflows between the two consecutive decades. However, nutrients in the ICW are also subjected to complex interplays of pH, temperature, microbiota, algae, and plants which affect P and N balances [6,29]. Overall, this study only provides an approximate picture of the long-term behavior of an ICW. Yet, the hydrologic behavior of watersheds can vary substantially in time and space. Thus, additional case studies are needed, and future studies are encouraged to use long-term and finer time-step ICW water quality measurements.

5. Conclusions

This study highlighted critical contrasts between ICW nitrate-N, TN, DP, and TP outflows in relation with precipitation variability. Despite the capacity of the ICW to abate N and P outflow, long-term SWAT simulations show sporadic disturbances of nutrient outflow patterns. These disturbances are characterized by unusual TN, DP, and TP releases into the downstream network. Changes in precipitation distribution, translated as an increase in extreme precipitation events, are linked to abnormal signals of nutrient releases. These signals highlighted in the long-term SWAT simulation of TN, DP, and TP outflows, substantiate the role of extreme precipitation events, as they generate strong water currents which flush nutrients accumulated at the bottom of the ICW [7,8,9]. On the one hand, this sporadic nutrient flushing lowers the ICW nutrient load and could be seen as a self-dredging process. On the second hand, the nutrient flushing causes an excess nutrient release into the downstream network, thereby exposing it to water quality impairment. Under the spectrum of climate change and an increasing trend of extreme precipitation events [15], disturbances of ICW nutrients outflow could be aggravated by ICW breaching due to structural damages caused by a frequent overflowing dike, abundant runoff, and strong water current.