We modeled N retention dynamics in eight wetland regions of the continental US representing different climate and geomorphology: California vernal pools, prairie potholes in North Dakota, basin wetlands in Minnesota, cypress domes in Florida, Texas playa lakes, North Carolina pocosins, coastal plain wetlands in Georgia, and Nebraska sandhills. Within each of these eight landscapes, a 1000 km² rectangular region was selected based on the presence of a high density of National Wetlands Inventory (NWI) wetlands, and infrequent gaps in remote sensing data [25]. GIWs were then derived within the selected rectangular regions using the NWI and National Hydrography Dataset (NHD), following the methodology outlined in Park et al [25]. Briefly, the method involves filtering the NWI dataset to focus only on palustrine and lacustrine wetlands. First, the NHD was used to create 10 m buffer polygons to filter out rivers, lakes >8 ha, other water bodies >1.5 ha, other flowlines or features, and all palustrine and lacustrine wetlands within the 10 m NHD buffered polygons were removed. Bathymetric relationships for each GIW were then derived to convert observed inundated areas to volume estimates. This was determined using the 1/3rd arc-second USGS DEM [26] cropped to the GIW extent, and starting from the top of the wetland incrementally determining the areas and volumes at each vertical increment using the rasterio python package [27].

We then characterized the inundation dynamics of the GIWs using the global surface water (GSW) dataset that was derived from Landsat imagery and contains monthly observations of surface water body extent [28]. By intersecting the boundaries of the identified GIWs with the GSW dataset, we developed a monthly time series of wetland inundated areas and volumes between 1985 and 2015. The method led to the characterization of 3698 GIWs across the eight wetlandscapes, with the greatest density of GIWs in the 1000 km² block in North Dakota, and much fewer GIWs in the vernal pools and the North Carolina pocosins (table 1). We treat each year of data for each wetland separately in the subsequent analysis, and refer to each as a wetland-year. Use of the remotely sensed observed inundation data allows us to integrate climate (e.g. dry and wet years) and landscape factors (e.g. downstream and upland wetlands), as well as various human controls (managed and unmanaged wetland systems) that contribute to the inundation dynamics.

A reduced-complexity model was developed to estimate temporal fluctuations in N retention dynamics in a wetland as a function of time-varying inflows and outflows (figure 1). Assuming well-mixed conditions within the wetland, the rate of change of nitrogen mass M [M] within the wetland was described as:

where Q_in is the volumetric flux of water entering the wetland from its contributing area [L³/T], C_in is the N concentration in the water entering the wetland [M/L³], Q_out is the volumetric flux of water leaving the wetland as recharge or surface outflow [L³/T], V is the ponded wetland volume [L³], and k is a first order N retention rate constant [T⁻¹] that describes the temporary (e.g. plant uptake, burial) and permanent (e.g. denitrification) processes that contribute to N retention in wetlands. While further differentiation of the surface and groundwater pathways of water and nitrogen could provide additional insight on how wetlands behave, there is insufficient data to quantify how the flows and solutes are partitioned. As our study aims to quantify the within-wetland dynamics, this model structure can still address how nitrogen retention is affected by downstream disconnectivity.

We estimated k (d⁻¹) using an empirical relationship developed by Cheng and Basu [6], that used measured N flux data, as well as water flow and surface area information (SA; m²) from 178 wetlands across the world:

The inverse relationship was attributed to the higher ratio of reactive sediment area to water volume in smaller wetlands that increases the opportunity for N in the water column to come in contact with the sediments, where N is removed by denitrification. Of course, other factors such as temperature and soil redox conditions that are specific to the different regions will impact the k values; however, Cheng and Basu’s meta-analysis highlights that wetland morphology acts as one of the most important first-order controls [6]. This simplification allows us to estimate k values at regional scales using wetland SA data that is easily accessible. Cheng et al used this relationship to estimate wetland N removal in 30 million wetlands across the continental US [10]. However, they assumed wetland area to be invariant in time. As a small wetland dries up over the course of a season, the reactive SA to volume ratio increases, and this can increase the retention rate constant even further. Here, we use time-varying SA of the wetlands from remotely sensed data to estimate time-varying k using equation (2).

Remotely sensed wetland inundation data was also used to estimate the time-varying wetland volumes, as well as inflows and outflows Q_in and Q_out from the wetland. Specifically, we used estimated monthly ponded volumes of the wetland (V), and available data for rainfall, evapotranspiration and runoff to estimate the inflows and outflows with the following mass balance [29]:

where P is the direct precipitation rate on the wetland [L/T], ET is the evapotranspiration rate from the wetland surface [L/T], and SA_max is the maximum inundated SA of the wetland over the 30 years of available remote sensing data [L²]. Although the ponded area of the wetland changes over time, we assumed direct precipitation to be intercepted by the maximum area. Additionally, though direct ET occurs only from the inundated area, we used the maximum SA to estimate the overall ET flux, given ET in the moist margins of the wetland can induce losses from the ponded area through local groundwater exchange [30].

Monthly precipitation (P), evapotranspiration (ET) and runoff (Q) from the TerraClimate dataset [31] was used to estimate Q_in and Q_out (figure 1). We assumed a catchment area to wetland area ratio of 8.2 based on the geometric mean of 33 000 wetland catchment delineations conducted by Wu and Lane [32], and estimated Q_in’ (m³/month) as the product of Q [L/T] and the catchment area of the wetland. We then estimated changes in the monthly wetland inundated volume (dV/dt) from the satellite data. For each monthly timestep, when the total water inputs to the wetland was greater than the observed change in volume (Q_in’ + (P–ET)*SA_max ⩾ dV/dt), outflow was estimated as the difference between the water inputs and the change in volume (Q_out = (Q_in’ + (P–ET)*SA_max−dV/dt) and Q_in was assumed to be equal to Q_in’. Conversely, when the total water inputs were less than the observed change in water volume (Q_in’ + (P–ET)*SA_max < dV/dt), Q_out was assumed to be zero while Q_in’ was increased to meet the observed change in volume (figure 1). It is important to note here that the estimated inflows to and outflows from the wetland include flows through both surface and subsurface pathways, given that they are estimated from the measured wetland inundation dynamics that intercept flows from all pathways. These estimates are of course uncertain, and can be refined using better site specific information. However, our goal was not to develop an exact model for a particular landscape, but rather to explore how consideration of inundation dynamics can alter N retention patterns in these small, ephemeral wetlands across the landscape.

We assumed the concentration C_in to be temporally invariant to better isolate the effects of hydrologic variability on N retention within the wetland. We did consider C_in to vary as a function of wetland size, given that small wetlands are often located in upland areas adjacent to nitrogen sources (e.g. agricultural land) and intercepts flow with the highest N concentrations [15], whereas large wetlands generally have larger catchments, and thus lower concentrations due to dilution. Here, we used a concentration of 30 mg l⁻¹ for C_in of the smallest wetlands and the largest wetlands receiving 0.3 mg l⁻¹, and a linear interpolation on a logarithmic scale to determine the input concentrations across wetland sizes. It should be noted that the percent N retention, which is the focus of this study, is independent of the input concentration, when C_in is temporally invariant. Assumption of temporal invariance of C_in is justified given observations of chemostatic response for N in human-dominated landscapes with large N inputs [33].

The hydrologic and N retention model (equations (1) and (3)) were numerically solved with the forward Euler method using a daily time step to ensure model stability for each wetland and year. To accommodate the smaller time step, linear interpolation was used to convert the wetland volume time series to a daily time scale.

Model results were used to estimate the N retention under transient scenarios for each of the 3698 wetlands for all available years (16 362 wetland-years). The transient retention efficiency (R_TS) was calculated by normalizing the N mass retained (M; equation (1)) by the mass entering the wetland for each wetland-year. For analyses relating to wetland size, each wetland-year was assigned to a size bin based on the maximum annual wetland area (i.e. 10^2.5–10³, 10³–10^3.5 m², etc). The steady state N retention efficiency R_ss was estimated using a commonly used model for lentic water bodies under well-mixed conditions [34]:

where k (d⁻¹) and τ_ss (d) are estimated from the wetland SA using empirical relationships developed by Cheng and Basu (equation (2)); τ_ss = 1.51SA^0.23 (r² = 0.4, SA is the size of wetland in m²) [6]. We used the annual median SA from the GSW dataset for our SA estimate in the steady-state scenario.

To compare residence times between the two flow models, an effective residence time τ_TS was estimated for the transient state models using equation (4), given R_TS and the median of the time-varying k estimated from remotely sensed SA using equation (2). Finally, the Damkohler number Da (=kτ) was estimated for the steady and transient models. The Da is a dimensionless number that captures the chemical reaction and water residence timescales, where Da> 1 indicates that the system has sufficiently large residence times to allow for more nutrient retention, while Da < 1 indicates that the water and nutrients are flowing through the system at faster timescales than the reaction timescales and thus limiting nutrient retention [35].

We first explored the behavior of prairie pothole wetlands in the North Dakota region to quantify how flow transience controls N retention. This was done by running the N retention model (equations (1)–(3)) for each of the 2395 wetlands over the 30 year timeframe (1985–2015). There were some years within this timeframe when data was not available for some wetlands, possibly due to drier conditions and detection limits, and this led to 12 570 wetland-years of simulation (table 1). A median N retention efficiency of 84% (interquartile range IQR: 79%–91%) was observed across these 12 570 wetland-years, in contrast to the steady state scenario where only 40% (IQR: 39%–42%) of the N entering the wetland was retained (figure 2(a)). For the steady state scenario, a single N retention efficiency was calculated for each of the 2395 wetlands, as a function of their size (equation 4). The range of retention efficiencies for the transient scenario was greater than the steady scenario, given it explicitly considered that retention efficiencies of a single wetland can vary across years, as a function of transient hydrology.

When disaggregated by size, we found N retention to increase with wetland size under the steady state scenario, from 38% for small wetlands (10^2.5–10³ m²) to 50% for large wetlands (10^5.5–10⁶ m²). Interestingly, the relationship between N retention and size is reversed under the transient scenario, with greater N retention (88%) in the smaller wetlands compared to the larger ones (71%) (figure 2(b)). We hypothesize that this dependence of N retention on size and flow transience can be attributed to differences in flow dynamics between the small and large wetlands.

To explore this hypothesis further, we first studied the hydrologic dynamics of the wetlands. Small wetlands (<10^3.5 m²) in the prairie pothole region (PPR) have a median hydroperiod of four months (figure 2(c)), and typically fill up during the spring freshet and dry up over the summer months [30]. Such seasonal wetting-drying dynamics imply that a large proportion of water entering the wetland leaves as evapotranspiration. Our flow analysis highlights that in these small wetlands, only 28% of the water that comes in as snowmelt or runoff actually leaves as groundwater or surface water outflow, while the remaining 72% of water leaves as evapotranspiration (figure 2(d)). This is consistent with field studies in the PPR that have found that, in small wetlands, over 80% of the water can leave the system as ET, when considering losses in both the wetland and the surrounding uplands [30, 36]. Large wetlands (>10^4.5 m²), on the other hand, often have water for the entire year, and a greater proportion of incoming water leaves as lateral fluxes (43%) compared to ET based on our model results (figure 2(d)). The greater proportion of ET in small wetlands most likely contributes to a more terminal behavior for N retention in these small ephemeral systems.

Next, we compared the residence times and retention rate constants in small and large wetlands as a function of flow transience. For the steady state assumption, residence times increase as a function of wetland size [6], and the combined effect of a negative k–SA and a positive τ–SA relationship contributes to an increase in the N retention efficiency with increasing wetland size (figure 2(e)).

Transient systems are more complicated since a large proportion of the water entering the system might only leave the system through evapotranspiration, while the N is retained within the system—leading to different residence times for water and nitrogen. We estimated an effective nitrogen residence time (figure 2(e)) for the transient scenario by using the transient N retention efficiency (R_TS) (figure 2(b)) and median k to solve for τ_TS using equation (4). We found that the effective residence time in the transient scenario to be greater than the steady state scenario, and the τ–SA relationship shows a unique U-shaped behavior, with the smallest and largest wetlands having the largest residence times (figure 2(e)). While residence times in larger wetlands are larger because of their larger volume to lateral flow ratio, residence times in smaller wetlands are larger under the transient assumptions due to the terminal nature of their flow dynamics. Effective residence times in these small, ephemeral systems can be 10–20 times greater than under steady state assumption, while for the larger wetlands residence times are similar under transient and steady state conditions (figure 2(f)).

To further explore the impact of hydroclimatic variability on wetland N retention, we expanded the analysis from North Dakota to seven other wetlandscapes across the US. We found wetland inundation dynamics to be strongly driven by climate and wetland size. In the semi-arid PPR of North Dakota (aridity index AI = PET/P = 2.31), wetlands are frozen during the winter months and have peak inundation levels in April (figures 3(a) and (b)). This behavior is typical of snowmelt driven systems like the basin wetlands of Minnesota, prairie potholes in North Dakota and Nebraska Sandhills [25]. The climate drivers are further modified by wetland size, where small wetlands in the PPR typically fill up during the spring snowmelt and dry up over the summer months (figure 3(a)) while large wetlands often lose only a smaller proportion of their total volume during summer and rarely completely dry up (figure 3(b)). This contrasts wetlands in the more rainfed, humid systems like the pocosins in North Carolina (AI = 1.23), which have water for most of the year in both small and large wetlands, albeit the smaller wetlands have a greater drawdown magnitude compared to its volume (figures 3(c) and (d)).

This pattern is consistent across the hydroperiods in all eight wetlandscapes (black lines in figure 4), with smaller wetlands in semi-arid landscapes (regions with AI >2, i.e. TX, NE, CA, ND) being more ephemeral (shorter hydroperiods), while humid landscapes have a greater proportion of more permanent wetlands. For example, the median hydroperiod for the largest size class (10⁵–10^5.5 m²) in the Texas playas (AI = 4.71) is 6.8 months, while for the more humid North Carolina wetlands (AI = 1.23), the median hydroperiods for all wetlands >10^3.5 m² is one year, indicating that a greater proportion of the wetlandscape is composed of more permanent wetlands.

Evapotranspiration dominates the water partitioning in the semi-arid wetlandscapes (regions with AI >2, i.e. TX, NE, CA, ND), and typically accounts for greater than 50% of the incoming water (green bars in figure 4), with smaller wetlands having a greater proportion of ET fluxes. For example, ET loss accounted for 68% of the inflows in the small playas in Texas, but only 48% of the inflows for the larger playa wetland. More humid regions with aridity indices closer to 1 (FL, GA, NC, MN) have a greater proportion of lateral fluxes compared to ET, and no discernible flow partitioning differences with wetland size (figure 4). These wetlands may have extended periods when inflows to the wetland exceed evapotranspiration and allow for more continuous surface outflows regardless of wetland size.

Finally, these hydrologic differences translate to different patterns of N retention. While for most wetlands, the transient assumption leads to greater N retention compared to the steady state assumption, the difference is greater for the semi-arid wetlandscapes, like the prairie potholes in North Dakota compared to the more humid wetlands in Georgia (figure 4). Overall, N retention percentages for the semi-arid systems are 1.63–1.82 times greater in the transient compared to the steady state scenario, while for the humid systems N retention is only 1.43–1.65 times greater in the transient scenario (figure 5(a)). Indeed, while the Damkohler number Da is greater than 1 for all regions in the transient cases (figure 5(b)), Da is consistently higher in the semi-arid systems. This suggests that the loss of N in the semi-arid wetlands are more driven by the biogeochemical processes rather than transport out of the system.

Of course, differences in the N retention patterns between wetlandscapes cannot be completely explained by aridity alone. For example, North Carolina and Florida wetlands have similar AI values, but different N retention dynamics (figure 4), and this can be attributed to differences in temporal patterns of rainfall, as well as differences in landscape attributes that drive N retention. While the patterns of N retention dynamics across humid and semi-arid systems is interesting, the exact magnitudes should be interpreted with caution, especially for regions with relatively few wetlands, making the size dependent patterns less reliable. Furthermore, it is important to note that the residence time estimation in the steady state scenario is based on an empirical relationship, and is thus invariant across the different wetlandscapes, while the transient scenarios are driven by the local inundation dynamics from remote sensing datasets. Despite these limitations, our analyses highlight that small, ephemeral wetlands have a greater potential to sequester nutrients and protect downstream waters, and this difference is more significant in semi-arid systems.