Wetland restoration yields dynamic nitrate responses across the Upper Mississippi river basin

We used the Soil and Water Assessment Tool (SWAT) model (version 659) to assess the effects of wetland restoration on changes in ${{\rm{NO}}}_{3}^{-}$ yields, loads, and concentrations across the UMRB (Figure S2). SWAT is a process-based watershed-scale hydrologic model commonly applied to simulate management impacts on hydrologic flows and aquatic nutrients (Gassman et al 2007). The model discretizes a watershed into a series of subbasins, which are composed of one or more hydrologic response units (HRUs), or areas of similar soil, slope, and land use properties, and simulates a daily water and nutrient balance for each HRU (Neitsch et al 2011). Transformation of nitrogen in the soil is controlled by partial differential equations that represent the nitrogen cycle (e.g., N₂ fixation, N plant uptake, denitrification, leaching, and volatilization; see Text S1 for further details).

Each HRU's simulated daily streamflow and nutrient loads then enter the river network, where nutrients are stored and processed (i.e., removed permanently or cycled through the river system) and are routed progressively downstream to the basin outlet (Neitsch et al 2011). In-stream transformation of nitrogen in the model is governed by partial differential equations representing growth and decay of algae, biological oxidation rates for different nitrogen species, water temperature, and settling of organic N with sediment (see Text S1).

We extended a recently developed wetland-integrated SWAT model for the URMB focused on river basin hydrology (Rajib et al 2020), updating the published version to simulate the response of ${{\rm{NO}}}_{3}^{-}$ to river basin-scale wetland restoration scenarios. Novel to this study, we substantially modified model inputs (Text S2; figure 1) and agricultural management operations (Text S3), calibrated wetland N removal rates (see below and Text S6), and calibrated nitrate loads at 19 sites across the UMRB (Text S4; Tables S2-S5). The model was constructed to include 279 subbasins with a mean area of 1,600 km² (SD = 1,100 km²), each having a single HRU depicting the subbasin's dominant land use (figure 1(A)), soil, and slope.

Zoom In Zoom Out Reset image size

Existing, or 'baseline', model wetlands (figure 1(B)) were identified as topographic depressions that intersected non-floodplain wetlands within the National Wetlands Inventory (Lane and D'Amico 2016, Wu and Lane 2016, US Fish & Wildlife Service 2018). From hereafter we refer to these non-floodplain wetlands as 'wetlands'. Baseline wetlands covered ∼23,000 km² of the basin (or ∼5% of UMRB area). Further information regarding physical wetland representations in the model is provided in Text S5 and Table S1. The complete, new model set-up resulted in baseline N inputs and outputs depicted in figures 1C and 1D.

Wetland ${{\rm{NO}}}_{3}^{-}$ removal was simulated at a daily time-step per each subbasin (see Text S6). For each day, the subbasin's wetlands received a constant fraction (specified as the fraction of land draining to the wetlands) of the subbasin's simulated daily flow and nutrient load, comprising surface, shallow subsurface and groundwater contributions. Tile effluent could not be routed directly to wetlands in SWAT (similar to most existing large river basin models; see section 4.2). The wetland's ${{\rm{NO}}}_{3}^{-}$ removal (kg day⁻¹) was then simulated as the product of the wetland's daily simulated surface area (ha), the wetland's daily simulated ${{\rm{NO}}}_{3}^{-}$ concentration ([${{\rm{NO}}}_{3}^{-}$]; kg m⁻³), a user-specified and calibrated wetland N-removal constant (m yr⁻¹), and a constant unit conversion factor (Neitsch et al 2011). The model's N-removal constant is intended to account for denitrification and plant uptake of ${{\rm{NO}}}_{3}^{-}$ of the aggregated wetland within each subbasin (Ikenberry et al 2017). Spillage and a corresponding [${{\rm{NO}}}_{3}^{-}$] from a subbasin's wetland to the nearest stream reach occurred whenever the wetland's daily simulated storage exceeds the wetland's maximum storage capacity.

We simulated change in ${{\rm{NO}}}_{3}^{-}$ concentrations, yields, and loads at the outlets of the 279 subbasins and the UMRB in response to restoration of potentially restorable wetlands (Figure S3). Our restoration simulations considered only wetlands outside of floodplains (i.e., non-floodplain wetlands) because of their proximity to nonpoint pollution sources (Cheng et al 2020) and the capacity to classify their previous (historic) locations for restoration based on topographic analyses. Potentially restorable wetlands (figure 1(B)) were defined as topographic depressions (Wu and Lane 2016) that did not intersect current NWI non-floodplain wetlands, i.e., those not included in the baseline model. The depressions were assumed to indicate locations that may have historically functioned as wetlands but were drained or otherwise converted to support alternative land uses (primarily agricultural production within the UMRB). Upland areas that did not contain depressions were not considered as potential locations for wetland restoration. The potentially restorable wetlands totaled ∼8,000 km² (or ∼2% of UMRB area) and ranged from 0 to 12% of subbasin areas.

We first developed and simulated a single model scenario depicting complete restoration of the basin's potentially restorable wetlands. With this scenario, we assessed the maximum capacity of wetland restoration to reduce ${{\rm{NO}}}_{3}^{-}$ levels at subbasin and river basin outlets. The scenario was constructed by revising the baseline model's wetland representation to include the summed physical attributes (i.e., storage capacity, surface area, catchment area) of both the basin's existing wetlands (used in the baseline calibrated model) and potentially restorable wetlands. This scenario increased the wetland maximum surface area from 0 to 478 km² across the UMRB subbasins (mean = +28 km²). The additional maximum wetland storage capacity in the wetland restoration scenario ranged from 0 to 2.4 × 10⁸ m³ (mean = 1.6 · 10⁷ m³), and the proportion of subbasins draining to wetlands increased by 0 to 33% (mean = 4%). The model scenario also depicted a corresponding reduction in fertilizer application rates to account for the quantity of land that was assumed to be converted from agricultural production to restored wetland (Text S7).

We additionally ran a suite of 279 model scenarios, each targeting complete restoration of a single subbasin's wetlands, in isolation. We hereafter term this approach subbasin 'targeting', i.e., targeting individual subbasins one-at-a-time for wetland restoration. For each scenario, we revised the baseline model's description of wetlands in only the targeted subbasin to depict the summed physical attributes of the targeted subbasin's existing and potentially restorable wetlands (Table S1). In turn, we reduced the quantity of fertilizer applied in only the targeted subbasin to account for the subbasin's assumed conversion of land from agricultural production to restored wetland (Text S7).

The model scenarios – depicting either basin-wide or subbasin-targeted restoration – were compared against the baseline model, which incorporated representation of only existing or 'baseline' wetlands, to facilitate assessment of restoration effects in three ways. First, we simulated the daily effect of the restored wetlands on discharge (Q; mm yr⁻¹), ${{\rm{NO}}}_{3}^{-}$ yields (kg N ha⁻¹ yr⁻¹), and ${{\rm{NO}}}_{3}^{-}$ concentrations (mg N l⁻¹), which we averaged to and presented as mean annual changes at the UMRB outlet and the UMRB subbasin outlets. Second, we simulated the restored wetland ${{\rm{NO}}}_{3}^{-}$ removal rates (kg N ha⁻¹ yr⁻¹) within each subbasin, which we calculated as the quotient of the mean annual ${{\rm{NO}}}_{3}^{-}$ load reduction from each subbasin's wetlands (kg N yr⁻¹) and the area of restored wetlands (ha) in the subbasin. Lastly, we compared the magnitude of subbasin (or 'local') outlet ${{\rm{NO}}}_{3}^{-}$ load reductions (kt N yr⁻¹) from wetland restoration to ${{\rm{NO}}}_{3}^{-}$ load reductions (kt N yr⁻¹) at the UMRB outlet.

Simulated basin-wide wetland restoration (8,000 km² of potentially restorable wetlands, or 2% of the UMRB area) reduced the mean annual river basin outlet ${{\rm{NO}}}_{3}^{-}$ load by 53 kt N yr⁻¹, or by 12%. Approximately 93% of this load reduction was directly attributed to restored wetland storage and removal of ${{\rm{NO}}}_{3}^{-};$ the remainder resulted from the accompanying reduction in fertilizer applications atop converted land (see Text S8). Basin-wide wetland restoration also reduced the mean annual ${{\rm{NO}}}_{3}^{-}$ yield at the river basin outlet by 1.2 kg N ha⁻¹ yr⁻¹, with year-to-year yield reductions ranging from 0.77 to 1.66 kg N ha⁻¹ yr⁻¹ (figure 2(A)). Mean annual discharge (Q) at the basin outlet decreased by 15 mm yr⁻¹ (5%), with year-to-year reductions ranging from 10 to 18.4 mm yr⁻¹ (figure 2(B)). Year-to-year reductions in nitrate concentrations, hereafter symbolized as [${{\rm{NO}}}_{3}^{-}$], at the river basin outlet fluctuated between 0.16 and 0.38 mg N L⁻¹ (figure 2(C)), and the basin's mean annual [${{\rm{NO}}}_{3}^{-}$] declined by 0.25 mg N L⁻¹ (7%).

Zoom In Zoom Out Reset image size

The effects of subbasin-targeted wetland restoration were varied. Restoration reduced ${{\rm{NO}}}_{3}^{-}$ yields at the 279 subbasin outlets by 0 to 10.6 kg N ha⁻¹ yr⁻¹ and Q by 0 to 135.6 mm yr⁻¹ (figure 3(A)–(B)). While restoration reduced [${{\rm{NO}}}_{3}^{-}$] at a majority of subbasin outlets, it tended to increase [${{\rm{NO}}}_{3}^{-}$] in subbasins with higher baseline [${{\rm{NO}}}_{3}^{-}$] simulations (figure 3(C)), i.e., of the 25 subbasins with mean annual [${{\rm{NO}}}_{3}^{-}$] ≥ 10 mg N L⁻¹, 23 (92%) had increased [${{\rm{NO}}}_{3}^{-}$] following restoration. Changes in [${{\rm{NO}}}_{3}^{-}$] from wetland restoration ranged from −1.25 to 9.1 mg N L⁻¹ across subbasin outlets.

Zoom In Zoom Out Reset image size

Subbasin-targeted restoration revealed spatial variations in wetland ${{\rm{NO}}}_{3}^{-}$ reductive capacity. Simulated wetland ${{\rm{NO}}}_{3}^{-}$ removal rates ranged from 0.13 to 546 kg N ha⁻¹ yr⁻¹ across UMRB subbasins (mean = 50 kg N ha⁻¹ yr⁻¹). The restored wetlands removed ${{\rm{NO}}}_{3}^{-}$ at higher rates in subbasins with higher baseline simulated ${{\rm{NO}}}_{3}^{-}$ yields, i.e. the basin's critical source areas or 'hotspots' (Evenson et al 2021); figure 4. These subbasins tended to be in the basin's mid-section where agricultural activities and subsurface tile drainage are prevalent (figure 1(A), (D)). In fact, all subbasins above the 75th percentile ${{\rm{NO}}}_{3}^{-}$ yields were simulated as tile drained agriculture.

Zoom In Zoom Out Reset image size

Subsurface tile drains limited the restored wetlands' ${{\rm{NO}}}_{3}^{-}$ reductive capacity. Tile drains were an important pathway for ${{\rm{NO}}}_{3}^{-}$ loss – or export – from tile drained subbasins; model structure prevented routing this effluent directly to wetlands. The median percentage of total subbasin ${{\rm{NO}}}_{3}^{-}$ export via tiles was 67% (figure 5(A)). We accordingly observed that restored wetland ${{\rm{NO}}}_{3}^{-}$ removal rates increased as the percentage of non-tile ${{\rm{NO}}}_{3}^{-}$ (i.e., ${{\rm{NO}}}_{3}^{-}$ contributions from surface, shallow subsurface and groundwater flows) in the total ${{\rm{NO}}}_{3}^{-}$ yield became larger (figure 5(B)).

Zoom In Zoom Out Reset image size

We observed a clear linear relationship between wetland restoration-affected ${{\rm{NO}}}_{3}^{-}$ load reductions at the subbasin outlets and those of the river basin outlet (figure 6). The relationship was strong at higher magnitude load reductions at the subbasin scale (moving towards the upper-right corner of the figure 6 scatterplot). However, approximately equivalent magnitude ${{\rm{NO}}}_{3}^{-}$ load reductions at the basin outlet were less likely to occur if wetland restoration resulted in low-magnitude subbasin ${{\rm{NO}}}_{3}^{-}$ load reductions. Restoration-affected subbasin scale load reductions resulted in no change or even minor increases in ${{\rm{NO}}}_{3}^{-}$ loads at the basin outlet in 19, or 7%, of the subbasins (see figure 6 points and subbasins outlined in black). These subbasins were mostly in the northern portion of the basin, further from the basin outlet – affording long travel times through the river network and thereby elevated potential for storage and transformation processes to occur (Text S1) – and with a relatively high number of existing wetlands, as represented in the baseline model (figure 1).

Zoom In Zoom Out Reset image size

It is clear that wetlands can effectively remove nitrogen across landscapes (Johnston 1991, Cheng and Basu 2017). However, apart from recent empirical approaches that assume stationary watershed conditions (Hansen et al 2018, Cheng et al 2020), the cumulative effects of wetlands on downgradient aquatic nutrient conditions remain largely elusive (Golden et al 2019). Watersheds are nonstationary systems, requiring dynamic, process-based simulations to quantify the link between wetland restoration and watershed nutrient conditions, particularly across large river basins. Using this process-based approach, we found that wetland restoration across substantially reduces ${{\rm{NO}}}_{3}^{-}$ yields and loads (and to a lesser degree, concentrations) across the UMRB – with caveats.

Restoration of 8,000 km² of the URMB's historic wetlands would lead to a 12% reduction in mean annual ${{\rm{NO}}}_{3}^{-}$ loads at the river basin outlet. These benefits were achieved by an average 4% increase in subbasin area draining to wetlands (i.e., the fraction of land draining to wetlands following restoration) and without additional conservation practices (e.g., tile-effluent management, described more below). Our results suggest that higher magnitude load reductions may be attained via minor increases in the quantity of land draining to wetlands. We therefore underscore the capacity of wetlands to benefit water quality at river-basin scales (Golden et al 2019, Cheng et al 2020), in the greater Mississippi River Basin specifically (Mitsch et al 2001), in portions of the UMRB (Hansen et al 2018), and at smaller watershed-scales (Arheimer and Wittgren 2002, Czuba et al 2018).

Yet our results indicate that wetland restoration does not uniformly reduce ${{\rm{NO}}}_{3}^{-}$ in subbasins across this large river basin. Additional or complementary management practices will be needed to meet N load reduction targets set for the greater Mississippi River Basin – and likely other large river basins. Our simulated 53 kt yr⁻¹ (12%) reduction in mean annual ${{\rm{NO}}}_{3}^{-}$ loads following restoration of 8,000 km² of wetlands falls short of the 20% total nitrogen (TN) load reduction target set for the year 2025 (Hypoxia Task Force 2015) – noting that ${{\rm{NO}}}_{3}^{-}$ comprised 86% of the baseline model mean annual TN load. These results also highlight a need for cost-benefit analyses (à la Singh et al (2019) and Hansen et al (2021)) to inform restoration decision making processes.

Further, our simulations provide a low estimate of the extent to which wetland restoration reduces river basin nitrogen loads at the outlet, compared to other recent studies. For example, Cheng et al (2020) predicted targeted wetland restoration resulting in a 467 kt yr⁻¹ N load reduction for the greater Mississippi River Basin – meaning a 115 kt yr⁻¹ decrease in N loads for the UMRB, assuming restoration of the same extent of potentially restorable wetland area considered in our analysis (see Text S9 and Table S2). This is almost double our simulated ${{\rm{NO}}}_{3}^{-}$ load reductions in response to wetland restoration across the UMRB. Based on our findings, we suggest that wetland restoration should be implemented in concert with additional in-field conservation practices (e.g., cover crops), edge-of-field practices that treat tile-effluent (e.g., constructed wetlands or denitrifying bioreactors (Crumpton et al 2020)), or improved fertilizer management with reductions in fertilizer inputs to the basin.

Our methodological approach may, in part, explain our lower ${{\rm{NO}}}_{3}^{-}$ reductions to wetland restoration relative to previously published studies (Mitsch et al 2001, Cheng et al 2020). Cheng et al (2020), for example, used the spatial distribution of wetlands and wetland nutrient removal rate estimates to calculate cumulative wetland removal potential – and this removal was assumed to be directly realized at the basin outlet as a steady-state basin condition. In contrast, our model simulated dynamic wetland ${{\rm{NO}}}_{3}^{-}$ removal, estimated via the daily expansion and contraction of wetland surface area and daily changes in [${{\rm{NO}}}_{3}^{-}$].

Our results also suggest potentially muted downstream impacts of wetland restoration because of landscape and in-stream N processing that is represented in the model (Text S1), highlighting the importance of considering these processes when assessing river basin-scale effects of wetland restoration. For example, tile drainage reduced wetland capacity to remove ${{\rm{NO}}}_{3}^{-},$ figure 5; Text S3). Further, wetland restoration produced dissimilar subbasin and river basin outlet load reductions if subbasins had low-magnitude ${{\rm{NO}}}_{3}^{-}$ decreases (figure 6) – and these differences resulted from the model's simulation of in-stream N processes and removal of additional N (figure 6; Text S1).

While our findings demonstrate the benefits of targeted approaches to wetland restoration, we suggest a careful focus on processes that may limit wetland capacity to receive surplus nitrogen. We found restoration was more effective within the basin's high N yield areas (figure 4), which lends credence to targeting approaches. However, restoration targeting these high N yield subbasins was less effective in subbasins with elevated ${{\rm{NO}}}_{3}^{-}$ losses via tile drains (figure 5), enhanced by tile effluent not being routed to wetlands in the model. While wetlands are sometimes restored or constructed specifically to receive tile effluent (Crumpton et al 2020), our model reflects the UMRB as a whole: an expansively tiled landscape that promotes hydrologic connections and unmitigated nutrient transport. Thus, our results support targeted restoration of wetlands in high nitrogen yield or source areas but suggest that efforts may be less successful where tile-drains continue to bypass wetlands, limiting wetland capacity to receive and remove ${{\rm{NO}}}_{3}^{-}$ from the landscape. Tile drainage should be considered in addition to other landscape characteristics and processes - e.g., climate change (Bosch et al 2014), and nutrient legacies (Van Meter et al 2018) - in prioritizing restoration efforts and projecting river-basin scale effects.

River-basin scale simulations necessitate aggregating fine-scale processes that may affect restoration efficacy. For example, we simulated wetland NO₃–N removal as a function of a N-removal constant, depicting the aggregate effect of wetland denitrification, settling, and plant-uptake processes per subbasin. Model simulations were highly sensitive to the N removal constant (p < 0.001; Table S3). While we calibrated the parameter to in-stream ${{\rm{NO}}}_{3}^{-}$ loads, its specification could be alternatively verified using subbasin scale summaries of wetland N removal rates, when available – though we caution that individual wetland N removal rates may insufficiently represent aggregate wetland N removal per subbasin. Future research on aggregating or scaling-up individual wetland N removal rates will decrease uncertainty and improve confidence in model predictions.

Future analyses would benefit from a refined model representing interactive processes between tile-drainage and wetland restoration. A substantial portion of UMRB subbasins (49%) were simulated as tile-drained (figure 1(A)), and these drains were an important pathway for ${{\rm{NO}}}_{3}^{-}$ losses (the mean percentage of total subbasin ${{\rm{NO}}}_{3}^{-}$ losses via tiles was 60% [STD = 33%]). UMRB wetlands are often restored by 'breaking' existing tile lines to partially restore their hydrology. In addition, wetlands may be purposefully positioned downslope from fields to receive effluent from working tiles and more efficiently capture and filter nutrients (Zedler 2003, Crumpton et al 2020). SWAT (SWAT2012 ver. 659) does not currently facilitate representation of tile-effluent entering wetlands. Thus, our analysis will have underestimated the efficacy of restoration without tile-effluent routing to restored NFWs or at downslope treatment wetlands. We recommend structural modifications to the SWAT model to facilitate representation of these dynamics.

River-basin scale analyses of wetland restoration – and other water quality mitigation strategies, for that matter – relate a land management action (e.g., restoration) to nutrient loads and concentrations at an outlet-point separated by considerable distance. The storage and transformation of nutrient across these distances represents a gulf of uncertainty. Because a host of interacting nutrient storage and transformation processes are embedded in our model (see Neitsch et al (2011) and Text S1), our results are suggestive that nutrient dynamics along these pathways − i.e. between subbasins and the UMRB outlet − have important implications in determining the relative efficacy of restoration. We propose that there is a critical need to validate model simulation of these dynamics, particularly as applied at river-basin scales where they have greater capacity to impact simulation outputs. Subbasin scale summaries of river-lengths (e.g., initial nutrient storage, mineralization rates) could inform the parameterization and/or verification of the model's simulation of in-stream nutrient dynamics.