Seasonality of inundation in geographically isolated wetlands across the United States

Remote sensing provides the opportunity to map spatio-temporal variations in wetland extent over large regions. Several efforts have been made in this regard, especially at local scales (Wu and Lane 2016, Hird et al 2017, Kaplan and Avdan 2018, Rezaee et al 2018, Wu et al 2019). Recently, the Global Surface Water (GSW; Pekel et al 2016) and Dynamic Surface Water Extent (DSWE; Jones 2015, 2019) maps have been used to study inundation dynamics (Allen and Pavelsky 2018, Luijendijk et al 2018, Mentaschi et al 2018). Here, we used the GSW inundation data from March 1985 to October 2015 to map the seasonality of inundated areas of GIWs. GSW was used here as it provides the opportunity to extend the methodology to other places beyond the CONUS. In addition, it is easy to use, and provides an estimate of inundation for each pixel within a month. The GSW v1.0 data was downloaded using the Google Earth Engine.

We studied the inundation characteristics in ten large wetlandscapes of the CONUS (figure 1). These include the California vernal pools (CVP), Prairie potholes (PP), Basin wetlands (BAS), Maine vernal pools (MVP), Playa lakes (PL), Cypress domes (CYD), Coastal plain wetlands (COP), Pocosins (POC), Delmarva bays (DEB), and Nebraska Sandhills (NES). Inundation characteristics were assessed within a selected rectangular region of 1000 km² area at each wetlandscape. The location of the rectangular region within each wetlandscape was determined based on the following heuristic: (a) wetlandscape locations were first identified based on the Cohen et al (2016) and Tiner (2003b), (b) within each wetlandscapes, locations with a high density of wetlands based on NWI data were subset, (c) from these subsets, locations with little or infrequent gaps in GSW data were identified, (d) rectangular regions with a long-term (>5 years, figure S1 and table S1 available online at stacks.iop.org/ERL/17/054005/mmedia) USGS discharge gaging station nearby for which the contribution area largely lies in the wetlandscape, were preferred. Post-processing of inundation data (figure S7) shows that the shape index or the ratio of the total square root of wetland area to total wetland perimeter in many of these landscapes have changed during the analysis period, thus indicating anthropogenic influences (Van Meter and Basu 2015, Krapu et al 2018).

Zoom In Zoom Out Reset image size

To assess the role of hydroclimatic fluxes on the seasonality of GIWs' inundation area, monthly forcings data such as rainfall, SM, evapotranspiration (ET), and potential evapotranspiration (PET) were obtained from the NLDAS-2 data set (Xia et al 2012). The majority of these fluxes have been validated in a range of settings (Luo et al 2003, Pan et al 2003, Nan et al 2010, Nearing et al 2016) and has been used extensively in hydrologic studies (Mitchell et al 2004, Yu et al 2013, Maxwell et al 2015, Wang et al 2018, Zhang et al 2019). Data were extracted for all NLDAS-2 grids that intersected with 1000 km² study region within each wetlandscape. Records corresponding to each wetlandscape were spatially averaged.

GIWs were derived within the selected rectangular regions. To this end, we used a method similar to the one presented in Lane and D'Amico (2016). The method involves using the NWI and National Hydrographic Dataset (NHD). To select GIWs from NWI, first riverine, marine, and estuarine wetland systems were excluded from the NWI, i.e. only the palustrine and lacustrine wetlands were considered for further analysis. Next, NHD data that includes hydrographic features such as rivers, lakes, reservoirs, wetlands, and oceans was used. Following Lane and D'Amico (2016), we obtained 10 m buffer polygons for rivers, NHD lakes larger than 8 ha, other water bodies (e.g. reservoirs, playa, etc) larger than 1.5 ha, and all other flowlines and area features (including bays/inlets, locks, levees, etc.). Palustrine and lacustrine wetland polygons, selected above from NWI, that intersect with the 10 m buffered NHD polygons were then removed, as they were considered to have surficial connections with other waterbodies. The remaining wetlands were termed as GIWs. The resulting extracted wetlands (figure 2) were the considered GIWs for ensuing analyses.

Zoom In Zoom Out Reset image size

To evaluate GIW dynamics, first, the maximum extent of GIWs was derived from the GSW maps. To this end, a multi-step process was followed (figure 3). By counting all the water pixels in GSW data within a wetlandscape, the inundation area for every month from March 1984 to October 2015 was evaluated. From this, the wettest month during the analysis period, characterized by the maximum inundation area within a wetlandscape, was identified. For the wettest month in the GSW map, we then distinguished wet pixels into GIWs and non-GIWs. Non-GIWs were the contiguous wet pixels that touched or intersected with water bodies that are not classified as GIWs in the NWI data (identified in section 2.4) or, in other words, water bodies that are termed as NWI → other waterbodies in figure 2. If contiguous wet pixels did not touch any GIWs in the NWI data (or NWI → GIWs in figure 2), these pixels were classified as non-GIWs as they are likely to be either non-GIW wetlands or ephemeral pondings or misclassifications in the satellite data. It is also possible that commission errors were also due to poorer data resolution and inconsistent mapping protocol of older NWI data (Tiner 1997, Wu and Lane 2017), as we assumed the NWI to be the most accurate baseline dataset for these analyses. The rest of the wet pixels in the GSW map were identified as GIWs. Each contiguous wet pixel region was considered a GIW henceforth. The derived map represented the maximum extent of selected GIWs during the analyses period. It is to be noted that there is a possibility that the identified extent may not be the absolute maximum extent of the selected GIWs in reality, as the analyses could only be performed for the dates for which GSW data exists.

Zoom In Zoom Out Reset image size

The maximum extent of GIWs, extracted above was then used to obtain the temporal dynamics. To this end, for each GIW, the inundation area was evaluated for all months, with no-data values being less than 5% of the maximum extent area. When no-data pixels exceeded 5% of the maximum extent area of a given GIW for a given timestep, the record was neglected. Next, the inundation area fraction (IAF) of a GIW was calculated as the ratio of the wet pixel area in a given timestep divided by the all-time maximum wet pixel area of the GIW. The IAF varies from 0 to 1, thus providing a standardized inundation fraction for GIWs. IAF equal to zero indicates a completely dry GIW or a GIW with no water detected within it. A value of one indicates its maximum extent. For many GIWs, the minimum is greater than zero, indicating that those GIWs do not fully dry up. During wet periods, it is expected that individual GIWs will have larger IAFs. Next, the average IAF for a given month overall GIWs within each wetlandscape, hereafter called allGIW_IAF, was obtained. Finally, mean monthly allGIW_IAF (allGIW_IAF_m hereafter) and monthly standard deviation of allGIW_IAF (allGIW_IAF_std hereafter) were evaluated. In other words, a single allGIW_IAF_m for January was obtained based on allGIW_IAF estimates from all January months in the data period. The intra-annual variation in allGIW_IAF_m indicates the seasonality of GIW dynamics within a wetlandscape, while allGIW_IAF_std provides a measure of inter-annual variations in averaged monthly inundation fractions.

For certain months and locations, GSW data are non-existent (see figures 1 and S2). Following cues from recent studies (e.g. Walker et al 2020) that have demonstrated a correspondence between wetland area and regional water availability descriptors such as streamflow, regressions between inundation characteristics, and monthly streamflows were developed. Here, the allGIW_IAF_m (or allGIW_IAF_std) for the months with data were considered the response variables. The explanatory variable was the mean monthly streamflow per unit contribution area at a nearby streamflow gaging station. It is to be noted that other meteorological variables such as precipitation (SM + rain) and (precipitation—ET ) were also used to generate these relations, but streamflow turned out to be the most explanatory variable overall (see table S2). More information regarding the contribution area of streamflow gaging stations is provided in figures 1, S1 and table S1. Using the regression equations derived in figures 4 and S3, estimates of allGIW_IAF_m (allGIW_IAF_std) for the months with missing data were obtained. The filling was performed in wetlandscapes CVP, MVP, PL, POC, and DEB. For the wetlandscapes where the regression performance was not satisfactory, i.e. with p-value > 0.05 (e.g. wetlandscapes COP and NES), the filling was not performed. As the goal of this study is restricted to assessing the seasonality of wetlands, instead of the exact magnitude of allGIW_IAF_m (or allGIW_IAF_std), we qualitatively evaluated our overall results with respect to observation data from isolated observations in the considered wetlandscapes. Notably, an alternative regression configuration was also considered. Herein, the response variable was the average inundation fraction of all GIWs for each time step or allGIW_IAF, and the explanatory variable was the corresponding streamflow for these time steps. Results from the analyses indicate that the selected regression configuration (shown in figure 4) is overall more effective than one shown in figure S4 and hence was used for the ensuing analyses.

Zoom In Zoom Out Reset image size

Inundation dynamics of GIWs are expected to depend on a range of water exchange fluxes (figure S6). GIWs are often recharged by falling rain and snow on them or by the fluxes indirectly supported by rain and SM, such as lateral groundwater flow between wetland and the neighboring aquifer (Neff et al 2020, Park et al 2020) and runoff from the local contributing area (Shook et al 2013, Wang et al 2021). Notably, runoff from the contributing area is also a function of rain, SM, and ET . Water is removed from wetlands by open water evaporation, ET via vegetation, and bed leakage (Hayashi et al 2003, Shook et al 2013, Liu and Kumar 2016). While fine spatio-temporal resolution estimation of these fluxes for each wetland at the scale of our analysis is prohibitive due to lack of refined data and prevalence of a range of processes in play, the overall seasonality is expected to be a function of rain (R), PET, ET, and snowmelt (SM) (Sánchez-Carrillo et al 2004, Park et al 2014, Penatti et al 2015, Liu and Kumar 2016, Zhu et al 2017, Bertassello et al 2020, Lee et al 2020). Here our goal is to assess the extent to which regional data of R, PET, ET, and SM, with inherent uncertainties, can explain the seasonal dynamics of GIW inundation at wetlandscape scale.

To assess the extent to which hydroclimatic forcings influence the seasonal wetland dynamics, multiple linear regression with variable selection was performed at each wetlandscape. The change in allGIW_IAF between consecutive months or allGIW_IAF(t)-allGIW_IAF(t-1) (hereafter ΔallGIW_IAF) was the dependent variable. Here t is the index for a given month. Use of ΔallGIW_IAF, instead of allGIW_IAF, as the dependent variables were driven by the fact that similar allGIW_IAF could be observed within a year during seasonal recharge and recession periods when the meteorological forcings could be very different. The candidate independent variables were R, ET, PET, and SM. For the analysis period of N years, regressions were generated using different combinations of data of N-1 years based on the leave-one-out method. All 15 possible subsets of independent variables were considered one-by-one. Among the subsets of independent variables, the subset producing the highest R² was selected. This subset was then used to predict ΔallGIW_IAF for the year not considered in the regression. The process was repeated for other N-1 years. This was followed by evaluation of allGIW_IAF(t) for all months during a year, for every year, using the equation allGIW_IAF(t) = allGIW_IAF(t-1) + ΔallGIW_IAF(t). Predicted allGIW_IAF for all years were then averaged to obtain the seasonality in allGIW_IAF_m.

Intra-annual variations of IAFs in wetlandscapes are shown in figure 5. Overall, wetlands that are affected by SM (hereafter SM-affected), such as the PP, BAS of Minnesota, MVP, and NES, experienced peak allGIW_IAF_m in April. In these wetlandscapes, wetlands may often be frozen in certain winter months. Here the inundation area evaluation is performed only for the months with average air temperature greater than zero in at least 50% of the observation years. It is to be emphasized that the shown temporal variations do not indicate that all GIWs in these wetlandscapes were at their largest inundation extent in April every year, but rather the mean of monthly extents over all GIWs are maximum in April. After April, allGIW_IAF_m progressively decreased till August. Following August/September, the decreasing trend usually reversed. The extent of intra-annual variations in allGIW_IAF_m, however, showed differences across wetlandscapes. For example, compared to the maximum allGIW_IAF_m, the minimum allGIW_IAF_m for PP, BAS, MVP, and NES were 1.4%, 0.0%, 23.1%, and 1.6% in magnitude. The coefficient of variation or the ratio of standard deviation to the mean, indicating the magnitude of intra-annual variations of allGIW_IAF_m, for these wetlands were 1.46, 1.51, 0.56, and 1.40, respectively. Notably, intra-annual variation of allGIW_IAF_m in MVP was found to be relatively smaller than in other wetlandscapes. Recovery of the wetland area in fall was also the highest in the MVP, with October allGIW_IAF_m being around 76% of the maximum allGIW_IAF_m.

Zoom In Zoom Out Reset image size

Seasonal variations obtained above, for the most part, align with the groundwater table or water depth variations observed at isolated gaging stations within the concerned wetlandscapes. For example, in the PP region, Johnson et al (2004) observed peak water table elevation during April/May and a decrease in the following months. In some years, such as the deluge years of 1994 and 1995 in their observation period, they also reported a second peak in mid-summer or early fall. In most years, the decreasing trend reversed in September/October. It is to be noted that while our results indicate an increase in variability in mid-summer, the data does not show a significant increase in allGIW_IAF_m during this period. Observation sites in BAS and NES also corroborate the overall reported intra-annual variations of allGIW_IAF_m at these wetlandscapes, with peak timing in April followed by a reduction in summer and subsequent recovery in fall. For example, BEC (2010) and MBWSR (2013) reported the highest water table during April/May, a decreasing trend in summer, and recovery in fall in wetlands within Saint Louis and Dodge counties of Minnesota.

Similarly, Gosselin et al (1999) and Gosselin et al (2000) observed higher water elevations during March to May within the sandhills of Nebraska. While studies detailing water table dynamics for vernal pools in Maine are scant, its description as 'small, ephemeral wetlands that typically fill in spring with SM and precipitation, or fall with rising water tables, and are dry by summer's end' aligns with the reported temporal variation of inundation area in figure 5 (Calhoun et al 2014, p 11 002). Kifner et al (2018) also described MVP as 'typically at their highest water level in the spring, dry down by mid-summer, and refill in the autumn, while some dry on cycles longer than a year.'

Among the wetlands that are minimally affected by SM or, in other words, those mostly driven by rain (hereafter rain-fed), the intra-annual variations of inundation were generally similar to each other, but some did show stark differences (figure 5). For example, the allGIW_IAF_m peaked in February in CVP, Florida CYD, COP, POC of North Carolina, and Delmarva bays in Maryland (DEB). In contrast, the peak of allGIW_IAF_m in PL of Texas was in June. Except for PL, all rain-fed wetlandscapes experienced a decrease in allGIW_IAF_m in the months following February, with minimum allGIW_IAF_m usually occurring during May to August. COP and POC showed a second smaller peak during fall, with allGIW_IAF_m magnitude being around 93% and 68% of its maximum value, respectively. In PL, the minimum allGIW_IAF_m usually occurred between February to April. In terms of the extent of variations, the ratio of minimum to maximum allGIW_IAF_m for CVP, CYD, COP, POC, DEB, and PL were 18.7%, 3.1%, 60.6%, 40.4%, 6.5%, and 2.5%, respectively. The coefficient of variation of allGIW_IAF_m for these wetlands were 0.61, 1.07, 0.15, 0.33, 0.99, and 0.82, respectively. The result indicates that the intra-annual variation of the allGIW_IAF_m was relatively moderate in POC and DEB.

Overall, the intra-annual variations in rain-fed wetlands reported in figure 5 align with observations at isolated sites in these wetlandscapes, albeit with a few exceptions. For example, Rains et al (2006) observed significant variations in stages in vernal pools following precipitation events, with an overall decreasing trend from February to June. Similarly, Brown (1981) and Kasischke et al (2003) reported a water level peak in January in CYD of Florida. Notably, these studies also noted additional isolated peaks between August and October in response to intense precipitation that oftentimes originated from hurricanes or tropical storms. Along the lines of reported results for COP shown in figure 5, Williams et al (2002) and Amatya et al (2020) observed water table depths in wetlands near the southeastern Atlantic coast peak in spring and autumn and be relatively dry in summer. The distinct seasonality of inundation in PL (see figure 5) is in line with the water depth observations reported by Ganesan (2010) and Weinberg et al (2015), where water table peaks were reported in July/August with smaller values in fall and winter.

Overall, during the analyses period, the long-term seasonality or the intra-annual variation of allGIW_IAF_m showed changes in terms of the inundation extent (figure 6). Following the strategy used in Borja et al (2020), we divided the analyses period into two halves to track the long-term change. Seasonality was then obtained for both these periods. Eastern CONUS wetlandscapes, including MVP, COP, POC, and DEB, showed an overall increase in wetland extent in the later period. In contrast, western wetlandscapes such as CVP and NES showed an overall decrease in inundation extent. These results are consistent with previously reported open surface water trends (Zou et al 2018, Borja et al 2020), although it is to be noted that our analyses capture the changes in the inundation extent of the same GIWs over time while the previous studies evaluated the total open surface area over a landscape. This could possibly be the reason why our results show a decrease in inundation, especially during the spring melt, in BAS and PP wetlandscapes, which is in contrast to the results reported in Borja et al (2020). Notably, even though the IAFs for these wetlandscapes have changed over time, the dry and wet periods during the year have remained almost the same. For PL, the month with the largest and smallest allGIW_IAF_m shifted by a month.

Zoom In Zoom Out Reset image size

The previous section discussed the seasonality or the intra-annual variation of allGIW_IAF_m for each wetlandscape. However, as IAF shows significant variations from one year to next (see the monthly mean standard deviation shown in figure 5), it is expected that the maximum (minimum) allGIW_IAF in a given year may not occur in the month with maximum allGIW_IAF_m. To assess inter-annual variations in inundation seasonality, its consistency, quantified by the fraction of years the maximum (minimum) allGIW_IAF occurs in the top three (bottom three) months of allGIW_IAF_m, was evaluated. For calculating the consistency, only the years with at least one record among the top three or the bottom three months were used.

In most wetlandscapes, the consistency of month with the maximum inundation area was more than 50% (figure 7), i.e. in more than half of the years, the maximum inundation area lies within the top three months in terms of the magnitude of allGIW_IAF_m. For SM-affected wetlandscapes, i.e. PP, BAS, MVP, and NES, corresponding consistencies were all greater than 75%. High consistency in SM-affected wetlands is unsurprising to some extent, as the SM is dominantly driven by annually recurring higher temperatures and radiation in late spring and early summer. Arguably, an even higher consistency should perhaps be expected for SM-affected wetlandscapes. However, this is not so due to a few large rainfall events in late summer and fall in some of these wetlandscapes, as noted in Neill (1993) and Johnson et al (2004), which may cause widespread inundation. The consistency of month with the maximum inundation area was relatively smaller for rain-fed wetlandscapes CYD and COP, with consistency magnitude being smaller than 50%. Notably, both CYD and COP experience occasional flooding due to hurricanes and/or tropical storms during June to November, thus leading to seasonal peaks in months that do not coincide with the three months with mean maximum allGIW_IAF_m. The lower consistency of the month with minimum allGIW_IAF_m for COP can again be partially attributed to these storms. Overall, SM-affected wetlandscapes experience high consistency for the month with minimum allGIW_IAF_m, in part due to distinct seasonality of net soil water input that is dominated by SM in late spring/early summer and relatively higher ET rates in summer.

Zoom In Zoom Out Reset image size

Qualitatively, the monthly variation of both long-term average R + SM and R + SM-ET (see figure S5) shows some correspondence with the intra-annual variation of allGIW_IAF_m (figure 5 and table S2). Notably, at least one variable between long-term average R + SM and R + SM-ET with lag of 0 and 1 month shows a coefficient of determination or r² value of more than 0.5 in seven out of ten sites. The other three sites where r² < 0.5 are CYD, COP, and POC. All of these three sites are affected by rapid inundation due to hurricanes and tropical storms. Overall, R + SM showed better r² to allGIW_IAF_m at more sites than R + SM-ET. The aforementioned results indicate that long-term average R + SM and/or R + SM-ET show similar intra-annual variations as allGIW_IAF_m for most sites.

Next, we quantitatively assess the extent to which inundation seasonality can be predicted by the seasonality of controlling fluxes. Here the assessment was performed for the two halves of the analysis period, viz. 1984–1999 and 2000–2015, using the methodology outlined in section 2.7. The reason to do the separate assessment for the two periods is that smaller period constraints the extent of anthropogenic change, if any, thus allowing a better understanding of the influence of controlling fluxes on inundation seasonality. Results (figure S8 and table 1) indicate that regressions based on R(t), SM(t), PET(t), and ET(t) can capture more than 52% (57%) of the intra-annual variations in allGIW_IAF_m in 8 out of 10 wetlandscapes during 1984–1999 (2000–2014).