Flooding in Landsat across tidal systems (FLATS): An index for intermittent tidal filtering and frequency detection in salt marsh environments

Remote sensing can provide critical information about the health and productivity of coastal wetland ecosystems, including extent, phenology, and carbon sequestration potential. Unfortunately, periodic inundation from tides dampens the spectral signal and, in turn, causes remote sensing-based models to produce unreliable results, altering estimates of ecosystem function and services. We created the Flooding in Landsat Across Tidal Systems (FLATS) index to identify flooded pixels in Landsat 8 30-meter data and provide an inundated pixel filtering method. Novel applications of FLATS including inundation frequency and pattern detection are also demonstrated. The FLATS index was developed to identify flooding in Spartina alterniflora tidal marshes. We used ground truth inundation data from a PhenoCam and Landsat 8 pixels within the PhenoCam field of view on Sapelo Island, GA, USA to create the index. The FLATS index incorporates a normalized difference water index (NDWI) and a phenology-related variable into a generalized linear model (GLM) that predicted the presence or absence of marsh flooding. The FLATS equation for predicting flooding is 11 - NDWI4 68 Pheno3 , 4 , and we found that a cutoff 0.1 was the optimized value for separating flooded and non-flooded pixel classes. FLATS identified flooded pixels with an overall accuracy of 96% and 93% across training data and novel testing data, respectively. FLATS correctly identified true flooded pixels with a sensitivity of 97% and 81%, across training and testing data, respectively. We established the need to apply FLATS when conducting vegetation time-series analysis in coastal marshes in order to reduce the per-pixel reflectance variations attributed to tidal flooding. We found that FLATS identified 12.5% of pixels as flooded in Landsat 8 tidal marsh vegetation time-series from 2013 to 2020, after traditional quality control and preprocessing steps were conducted, which could then be filtered out or modeled separately in order to conduct remotely sensed vegetation assessments. Therefore, in tidal wetlands, we recommend incorporating FLATS into Landsat 8 preprocessing prior to vegetation analysis. We also demon- strated innovative applications for the FLATS index, particularly in detecting flooding frequency and flooding patterns relevant to the broader biophysical modeling framework, including mapping marsh vulnerability due to fluctuation in inundation frequency. The FLATS index represents advancements in the understanding and application of inundation indices for coastal marshes.


Introduction
Remote sensing has been widely used as an effective means to expand our understanding of vegetation and ecosystem dynamics in landscapes, such as wetlands, that are difficult to monitor on the ground using traditional field monitoring techniques .
Coastal marshes present a distinct challenge because of their periodic inundation due to the ebb and flow of tides, which dampens spectral reflectance. Thus, there is a need for an accurate flagging mechanism to identify flooded pixels (Han and Rundquist, 2003;Kearney et al., 2009;Mishra et al., , 2012O'Connell et al., 2017). For example, intermittent variability in the spectral signal caused by tides can have major effects on the interpretation of data or parametrization of biophysical models, particularly in the case of vegetation time-series analysis (Cho et al., 2008;O'Connell et al., 2017;O'Connell and Alber, 2016). Vegetation indices that are commonly used to examine changes in vegetation over time, such as the Normalized Difference Vegetation Index (NDVI) or the Enhanced Vegetation Index (EVI) created for use in terrestrial systems and do not account for tidal interference in the spectral signal (Cho et al., 2008;Han and Rundquist, 2003;Liu and Huete, 1995;Rouse et al., 1974). An index that is designed specifically to flag tidally flooded pixels is necessary in order to accurately capture changes in vegetation and flooding dynamics in coastal marshes over time.
Our previous work created the Tidal Marsh Inundation Index (TMII) to identify and subsequently remove flooded observations in tidal marshes (O'Connell et al., 2017). However, we calibrated the TMII for the Moderate Resolution Imaging Spectroradiometer (MODIS). MODIS has high temporal resolution but has a much coarser spatial resolution (250 -1,000 m) which is less suited to irregular and fragmented vegetation patches that are often found in coastal marsh landscapes (Belluco et al., 2006;O'Connell et al., 2017;Rogers and Kearney, 2004). MODIS, managed by the National Aeronautics and Space Administration (NASA), provides a complete view of Earth every 1-2 days (Justice et al., 1998;Running et al., 1994). Coarse pixels typically contain more than one landcover type. These mixed pixels alter spectral dynamics so that MODIS optimized indices do not necessarily directly translate to sensors with finer spatial resolution and more homogeneous pixels (Belluco et al., 2006;O'Connell et al., 2017). Therefore, new flooding detection techniques should be developed for use with moderate-high spatial resolution satellite data.
The ability to utilize higher resolution data, such as that from Landsat 8, will allow for a more detailed understanding of tidal marsh vegetation dynamics. Landsat 8 is a collaborative mission by NASA and the United States Geological Survey (USGS) that provides multispectral remote sensing images at a 30-m spatial resolution, which is considerably higher than MODIS. Landsat 8 was launched in 2013 and has a 16day temporal resolution (Lulla et al., 2013). In this paper, we present a new index called the FLATS (Flooding in Landsat Across Tidal Systems) index, which is optimized for the use of these higher spatial resolution imagery. This expands the analytical possibilities for more accurate vegetation time-series analysis and allows us to detect finer-scale changes in flooding patterns that could impact the ecology of a given area.
In this paper, we describe the method used to create the FLATS index using Landsat 8 data in order to detect coastal inundation in a Spartina alterniflora marsh on Sapelo Island, GA, USA. We also establish best practices for applying FLATS for vegetation time-series analysis in coastal marshes, including cutoff thresholds. Finally, we provide a demonstration of innovative applications for FLATS, particularly in the area of detecting flooding frequency and flooding patterns. FLATS transforms the future of remote sensing applications in coastal wetlands due to the high-resolution inundation flagging it provides. FLATS represents advancements in the understanding and application of inundation indices for current and future coastal marshes.
Similar to TMII, the FLATS index incorporates a normalized difference water index (NDWI) that differentiates water from other land cover types (Ji et al., 2009;O'Connell et al., 2017). FLATS also incorporates an index that takes seasonal phenology into consideration, as the relative differences between flooding and vegetation differ between the growing season when vegetation is green and senescence when vegetation is brown. This phenology variable allows us to discern seasonal changes in the ratio between water, vegetation, and other land cover types (O'Connell et al., 2017). Although the FLATS index utilizes the architecture of TMII, it differs because it uses different band combinations that are more suited to finer spatial resolution data that have a higher frequency of homogeneous pixels. It, therefore, should have greater scalability to other remote sensing platforms such as near-identical Landsat 9, and even higher resolution sensors such as Sentinel 2A and 2B. FLATS can be cross-calibrated with a number of other satellites such as Landsat 5 and Sentinel-2 to examine long term trends in flooding frequency and intensity in coastal wetlands threatened by sea level rise.

Study site
We developed the FLATS index from training data that came from eighteen Landsat 8 pixels located on Sapelo Island, GA, USA (henceforth denoted as "Sapelo Island") within the Georgia Coastal Ecosystems (GCE) Long Term Ecological Research (LTER) site and the Sapelo Island National Estuarine Research Reserve (NERR) (Fig. 1A). This salt marsh is mainly composed of homogenous covers of Spartina alterniflora in three height forms (short, medium, and tall) at varying densities ranging from dense to sparse (Hladik et al., 2013). Tidal channels infiltrate the marsh platform and are flooded by semi-diurnal tides that often rise more than 2 m (tidesandcurrents.noaa.gov). An eddy covariance carbon flux tower affixed with a PhenoCam (see section 2.2.1 below for more detail) maintained by the GCE-LTER is located within the study area marsh (31.441 • N, 81.284 • W) and measures carbon dioxide fluxes between the marsh and the atmosphere (Fig. 1B). The flux tower also provides ground truth environmental data, including water level and temperature sensors. Upon optimization of FLATS, the index was applied to a 4.39 km 2 area of Sapelo Island (as seen in Fig. 1A) that contains a large monoculture of S. alterniflora marsh, the Duplin River, forested upland, and hammock areas.

Analytical approach
To develop the FLATS index, we utilized techniques similar to O'Connell et al. (2017) which are outlined in Fig. 2. First, we acquired ground truth flooding data within the focal area (Fig. 2, step 1a). Next, we obtained Landsat 8 data via Google Earth Engine (GEE), which was post-processed to filter out pixels with cloud shadow and/or cloud cover (Fig. 2, step 1b). To develop FLATS, we tested numerous band combinations in order to identify those combinations which relate best to flooding status (Fig. 2, step 1c). Once appropriate band groupings were selected, the FLATS index was created, and an optimum cutoff value was identified (Fig. 2, step 1d). The adjusted cutoff value was utilized when applying the index to new pixels within the study area (Fig. 2, step 2a). Depending on the desired application, either vegetation analysis or flooding frequency analysis, a different series of steps is recommended (Fig. 2, step 2b). If vegetation analysis is the intended application, pixels flagged by FLATS as containing values above the flooding classifier cutoff are identified as flooded and should be removed (Fig. 2, step 2bV.1). You then are left with only vegetated pixels and can proceed with your desired vegetation time-series analysis (Fig. 2, step 2bV.2). If flood frequency is the objective of your analysis, then you should keep all available data (Fig. 2, step 2bF.1) and plot the FLATS flooding classified data by pixel through time to examine trends in flooding frequency (Fig. 2, step 2bF.2).

Ground truth tidal inundation data
First, accurate tidal inundation data were obtained for the focal area in order to ground truth the FLATS index. Data was acquired from a PhenoCam (StarDot NetCam SC 5MP IR, StarDot Technologies, Buena Park, CA, USA) which is located on the GCE-LTER flux tower on Sapelo Island and is identified as the "GCESapelo" site within the PhenoCam Network (phenocam.unh.edu). The PhenoCam is a digital camera which has been programmed to automatically take oblique angle photos of the marsh surface every 30 min during daylight hours (Fig. 1C). This frequency of image collection allows for flooding and drainage patterns to be observed across the marsh platform during a range of tidal sequences. In order to discern flooding on the marsh platform, the PhenoCam field of view (black outline; Fig. 1B) was approximated based on the height and angle of the PhenoCam as well as field surveys and natural reference points (i.e., the creek, upland area, etc.). Landsat 8 pixel outlines and a habitat map were superimposed within the PhenoCam field of view ( Fig. 1B) in order to assist in image assessment (Hladik et al., 2013). The PhenoCam images were then visually examined to detect flooding in coincident Landsat 8 passovers between September 2013 and December 2020, when the PhenoCam was operational. Pixels were identified as "flooded" when water was present (i.e., when marsh vegetation was observed to be either partially or fully inundated).
Environmental data were also available from various sensors associated with the GCE-LTER flux tower to assist in validating tidal flooding, including two pressure transducers (Campbell Scientific Instruments Model CS455, Logan, UT, USA) which measured water levels within a creek adjacent to the flux tower as well as on the marsh platform surface north of the flux tower. These water level data, collected in 5-minute increments, were aggregated to match the Phe-noCam data by averaging over the preceding half-hour of data. The PhenoCam images and water level data were filtered to between 15:00 and 16:30 Coordinated Universal Time (UTC) in order to match with Landsat acquisition times. This water level data served as a second form of ground truth data and was used in conjunction with the PhenoCam images to indicate the presence of flooding. The combined PhenoCam and water level data were pivotal in the development of FLATS because tidal flooding can be spatially variable through time with respect to wind patterns and tidal creek migrations, and is not always straightforward to estimate (O'Connell et al., 2017;Wu et al., 2021). Together, these datasets provided on the ground information of tidal marsh flooding that we coordinated with remote sensing data and knowledge of the field site (including elevation and phenological characteristics) to conservatively estimate individual pixel flooding status.

Obtain and process Landsat 8 data
Remote sensing data was obtained in the form of band information from the Landsat 8 satellite from March 2013 through December 2020   for the entire study area (Fig. 1A). Multispectral Level-1 Landsat 8 surface reflectance data were acquired from GEE. The data were subset into a focal area dataset that contained only the 34 pixels found within the PhenoCam field of view (Fig. 1B). These pixels were categorized based on proximity to the flux tower resulting in near (1-3), mid (4-9), far (10-18), and very far (>18) pixels, with the very far pixels being excluded from index creation. During model development, the decision was made to exclude the pixels located furthest away from the Pheno-Cam (pixels 19-34; Fig. 1B) due to the oblique angle of the camera's field of view. These excluded pixels were located greater than 150 m away from the PhenoCam, which prevented a clear view and distorted any visual indications of flooding status.
Landsat 8 has an overpass time of approximately 16:00 UTC at this site. Prior to acquisition, the Landsat 8 data were atmospherically corrected using the Land Surface Reflectance Code (LaSRC) algorithm (USGS, 2020). The quality assessment (QA) band, inherent in the Level-1 Landsat Surface Reflectance products, flags pixels using the C Function of Mask (CFMASK) algorithm that contains conditions that often have adverse effects on image processing, such as clouds (Foga et al., 2017). Post-processing of the Landsat 8 data included multiplying band data by a scale factor as recommended by the USGS and filtering to remove clouds and cloud shadows, as well as "no data present" pixels, as indicated by the pixel QA band. We also extracted water status from the QA band, though this is known to be inaccurate for shallow inundation such as tidal flooding. Multispectral outliers, e.g., wavelengths > 0.4 µm, were eliminated from the data pool as well as < 0 which were likely scene edges.
We established that using environmental data (e.g., water level data from the pressure transducer) in conjunction with the PhenoCam images provided a credible source for flooding information on the marsh surface. We obtained 295 observations per pixel within the focal area (i.e., PhenoCam's field of view) after excluding "very far" pixels (pixels 19-34 in Fig. 1B) between March 28 th 2013, and December 27 th 2020. After post-processing this Landsat 8 data by filtering out poor quality pixels as indicated by the QA band (those with cloud cover, cloud shadow, and filled), an average of 148 occurrences per pixel remained (2,663 pixels total), indicating that approximately 50% of the available data had to be filtered out from the 7-year multispectral surface reflectance dataset.

Create training, testing, and validation sets and identify predictors of flooding through band combinations
Prior to model fitting, the focal area dataset was divided into three separate datasets. The first dataset was considered training data and accounted for a random 70% of the data between March 2013 and December 2017. The second dataset consisted of testing data and contained the remaining 30% of the data between March 2013 and December 2017. The independent validation dataset consisted of two full years of data (January 2018 through December 2020). The testing, training, and validation data included 18 pixels per Landsat scene over 66 scenes and 7 years, utilizing 984 data points in total after filtering for clouds and quality. The model was trained on the training set, and then verified against the novel testing and validation sets, which were not used in model training. Each of these datasets contained a binary "wet" variable (true/false) from the Landsat QA band that indicated the presence of water (true) or not (false). This wet variable was used to compare against our final FLATS index to help evaluate its efficacy. A second binary variable called "flooded" was created based on the ground truthed flood status (e.g., the visual interpretation of the PhenoCam images and water level data) (0 = no flooding observed, 1 = flooded, NA = flood status was not clear). A spectral plot was created to examine temporal determinants of spectral reflectance, in order to help identify band predictors for our model (Fig. 3). For this, the mean and standard deviation of each spectral band were plotted based on flooding status (dry or flooded) and by season (growing: mid-March to September, dormant: October to mid-March).

Create FLATS index and optimize cutoff value
A binomial Generalized Linear Model (GLM) was used to create FLATS, initially through the utilization of the training dataset, where Landsat 8 bands were used as the predictors of the response variable (binary ground truthed flood status). Many binomial GLMs relating band combinations to flood status were explored prior to the selection of the final FLATS model (Appendix A). In order to determine which band combinations to use in our FLATS index, we referenced the spectral plot ( Fig. 3) as well as tested a variety of established indices such as the NDVI, the Green Normalized Difference Vegetation Index (GNDVI), EVI, the Soil-Adjusted Vegetation Index (SAVI), variations of the Normalized Difference Water Index (NDWI) including the modified NDWI (mNDWI), Wide Dynamic Range Vegetation Index (WDRVI), and various versions of phenology indices (see Appendix A for full list of indices) (Gitelson, 2004;Gitelson and Merzlyak, 1998;Huete, 1988;McFeeters, 1996;O'Connell et al., 2017;Xu, 2006). We also compared bands most directly translatable from the original TMII equation to ensure FLATS was an improvement.
Seasonal differences played an important role in the variation of the spectral signal for the original TMII and continue to affect the FLATS model. A phenology index was employed in order to account for these naturally occurring seasonal changes. After consulting the spectral plot and trying many band combinations, the best option as indicated by our model fit criteria (see next paragraph below) was a phenology index developed from two bands in the visible spectrum (band 3 -Green & band 4 -Red) (see results).
In order to select the best model to predict flooding, we used the following model fit criteria: model explained deviance (ED), flood accuracy as identified by contingency table output, and classification Receiver Operating Characteristic (ROC) curves (Hanley and McNeil, 1982). ED is used to describe the strength of the relationship between a model and actual data where ED is equal to the null devianceresidual deviance / null deviance multiplied by 100 (Zuur, 2009). ED has also been described as a pseudo coefficient of determination (R 2 ) (Zuur, 2009). The contingency table was the second form of model assessment and contains accuracy, sensitivity, and specificity values. Accuracy is defined as the number of correctly classified flood status pixels divided by the total number of pixels multiplied by 100. Sensitivity indicates the percentage of correctly identified flooded pixels (true positives), while specificity indicates the percentage of correctly identified non-flooded pixels (true negatives). The third and final model assessment utilized a ROC curve and sought to identify the maximized Area Under the Curve (AUC) by using the R package "pROC" (Robin et al., 2011). The ROC curve helps the user compare true positives (sensitivity) to false positives, aiding in the distinction between whether a model is poorly optimized or lacks the power to explain the data (Hanley and McNeil, 1982;O'Connell et al., 2017). Flooded and dry pixels were distinguished based on an optimized cutoff boundary which was selected so that the model's specificity (ability to correctly identify flooding) was maximized (>90) while maintaining a balance with accuracy and specificity. Any FLATS value below the determined cutoff value is considered dry while any FLATS value above that of the cutoff is considered flooded.
Following the confirmation of model performance, the model equation was then simplified to a spectral reflectance index equation by rounding the intercept and coefficient values to the nearest tenths place. A logit transformation was also applied to the model in order to convert the model back into the predictor units. This created the final FLATS equation. The final FLATS equation was then reapplied to both the training and testing datasets to ensure that the same cutoff value could be applied and render the same result as the more complex model in both cases. The model was also verified against the novel validation data that was withheld during model development.

Application of FLATS
Next, we applied the index to the entire study area (Fig. 1A) in order to examine its proficiency at detecting flooded pixels at a larger scale. We also explored the various applications of FLATS, including timeseries analysis of marsh vegetation biophysical characteristics and flooding dynamics.

Apply FLATS to new pixels within the study site
The model was initially applied to a small subset of pixels found within the PhenoCam field of view (in order to have coincident ground truth data to verify model accuracy) but was later applied to an expanded area at the study site (Sapelo Island, GA; Fig. 1A) for application exploration. This study area covers over 4,800 pixels (approximately 4.39 square kilometers) and contains many different land cover types, including open water (the Duplin River), upland vegetation (pine trees, etc.), mudflat areas, and the three different height forms of S. alterniflora, which were mapped by Hladik et al. (2013). Prior to the application of the index to the entire study site, we used the Hladik et al. (2013) habitat map to first remove non-marsh cover types from our analysis. The FLATS index was then applied to marshes that spanned the study area from 2013 to 2020.

Selection of application type
The ability to identify flooded pixels allows the FLATS index to be used in a variety of different applications related to tidal marsh systems. We have established two workflows showing the applicability of FLATS for tidal marsh vegetation analysis and identifying flooding frequency and patterns through time.

Vegetation biophysical analysis application.
Producing a reliable vegetation analysis of coastal wetlands involves flagging tidally inundated pixels as an initial step in preprocessing (O'Connell et al., 2017). Prior to FLATS, there were no quick approaches to accurately flag tidal flooding in Landsat 8 data. To show the utility of FLATS in a vegetation analysis, we used NDVI, one of the most widely used vegetation indices. An assessment was conducted on both the focal area and the entire study area to determine how many pixels were removed by each filter type. Our goal was to showcase enhancements FLATS provides over nonfiltered data and other water filtering methods such as the Landsat QA band. For simplicity we used a single pixel (pixel 12; Fig. 1B) from within the PhenoCam's field of view to illustrate the utility of FLATS through the comparison of three NDVI time-series. This analysis compared the raw, unfiltered NDVI values, the Landsat QA band filtered NDVI values, and the FLATS filtered NDVI values for one pixel through time (2013-2020). We used ground truth data to highlight flooded and dry observations in each NDVI time-series and fitted a smoothed spline to each of the NDVI time-series to visualize the phenological signal produced by individual filter types. We also calculated the monthly mean NDVI value of all marsh vegetation pixels within the study area and then compared the NDVI time-series based on filter type (i.e., unfiltered, Landsat QA band filtered, and FLATS filtered).

Flooding analysis application.
Mapping the frequency of flooding is an alternative way to utilize FLATS. Unlike the vegetation timeseries, analysis of flooding does not involve removing the flooded pixels indicated by the FLATS index. Instead, one uses the FLATS information to visualize patterns in flooding and highlight relationships between verified flooding frequency and more commonly available variables such as water levels from local tide gauges. To demonstrate this, we analyzed tidal flooding frequency derived from FLATS to discern changing patterns through time. This flooding frequency analysis was conducted by normalizing the percent of the time a pixel was flooded, as determined by the FLATS index, by the number of times that pixel occurred through the time-series. The per-pixel normalization of flooding emphasizes marsh areas that flood more frequently as compared to other areas within the marsh that have stayed dryer over the 7 years of this study. We then conducted sub-analysis examining yearly flooding frequency normalized per-pixel on a subset of 9 pixels containing short and medium form S. alterniflora located north of the flux tower outside of the PhenoCam's field of view. This location was chosen due to its lack of migrating creeks, which are known to influence flooding patterns, and its interior marsh structure, i.e., higher elevation with medium and short form plants. The analysis serves as a snapshot into how flooding patterns can change yearly and potentially result in ponding with future sea level rise in an area not directly impacted by creek migration or low elevations. We also conducted an analysis of the flood frequency trends on a monthly basis for the entire study area. For this analysis, the flooding frequency was normalized per pixel for each month throughout the study period, where the monthly flooding frequency was determined by subsetting all dates falling within the specified month (i.e., all pixels dated January 1st -January 31st within the 2013-2020 dataset).
We also feature two ways to explore relationships between flooding and other variables: water level above the marsh surface through time as a proxy for sea level rise and monthly Mean Higher High Water (MHHW) as a function of mean monthly flood proportion. The water level and MHHW data were obtained from the National Oceanic and Atmospheric Administration's (NOAA) Fort Pulaski, GA station (ID: 8670870) that was established in 1935. The water level was adjusted to reflect the local elevation via a spectrally corrected light detection and ranging (LiDAR) derived digital elevation model (DEM) at our site on Sapelo Island, GA (vertical datum: NAVD 88) (Hladik et al., 2013;Hladik and Alber, 2012). We extracted elevation data for each marsh pixel within our study area from Hladik et al. (2013)'s DEM and then subtracted the observed water level (relative to NAVD88) at Fort Pulaski that matched closest to the Landsat 8 scene acquisition time. The Fort Pulaski water level data is available in 6-minute intervals with times that closely match Landsat 8 acquisition (15:54 and 16:00 UTC). The average water level above the surface of the marsh was calculated for each date (scene) and was examined over time. To better understand how a variable such as MHHW can relate to flooding at our study site, we plotted the normalized mean flood proportion for each Landsat 8 scene (what proportion of the entire study area was flooded on a given date) then averaged the flood proportion per month. We then plotted these mean monthly flood proportions (derived from FLATS) against the monthly mean MHHW values at Fort Pulaski between 2013 and 2020. We did not alter MHHW levels to reflect marsh elevations.

Identify predictors of flooding through band combinations
The Landsat 8 spectral plot was created using training reflectance data (2013)(2014)(2015)(2016)(2017) and demonstrated clear distinctions between seasons as well as inundation status across portions of the electromagnetic spectrum (Fig. 3). The mean reflectance was greater in the growing season than in the dormant season. This difference was most prevalent in the visible and near infrared (NIR) portion of the spectrum for dry observations. The flooded observations in each season depressed the mean reflectance when compared to their dry counterpart, particularly in the NIR and shortwave infrared portion (SWIR) of the spectrum (band 5 -865 nm and band 6 -1,608.5 nm). The spectral plot indicated that the normalized difference between a visible band (e.g., band 3 -561.5 nm or band 4 -654.5 nm) and a NIR or SWIR band (bands 5 or 6, respectively) could help discern pixel inundation status as these had the greatest change in reflectance between flooded and dry observations. The plot also emphasized the need to include a variable that accounts for phenology since there is an overlap between dry dormant season observations and flooded growing season observations in the visible portion of the spectrum.

Create FLATS model and optimize cutoff value
Although many band combinations were tried, ultimately, we selected the grouping of an NDWI-type index with a phenology index to create the final FLATS model due to the presence of bands previously identified as important in Fig. 3 and the model's high ED (see Appendix A for full list of indices). The NDWI selected for the FLATS model consisted of band 4 (the red band centered at 654.5 nm) and band 6 (a SWIR band centered around 1608.5 nm), where the equation is: NDWI 4,6 = (band 4 -band 6) / (band 4 + band 6) (1) The selected phenology index is the normalized difference of band 3 (the green band centered at 561.5 nm) and band 4, expressed as:

Pheno 3,4 = (band 3 -band 4) / (band 3 + band 4)
( The binomial GLM model selected for best classifying the flood status of tidal marsh pixels in Landsat 8 data was: Flood Status = -1.57 + 19.97*NDWI 4,6 + 68.55*Pheno 3,4 (3) This model explained 82% of the difference between the ground truth flooding data and the model's predictions in the training dataset (Table 1). This model (Eq. 3) was then transformed into the units of the predictors to create the final FLATS equation: The FLATS equation (Eq. 4) outperformed NDWI 4,6 (ED -73%) alone, which emphasizes the need for incorporation of a phenology variable. The FLATS model also significantly outperformed a model (ED -52%) which used the TMII equation while substituting for similar Landsat 8 bands (bands 3 and 6 for flooding, bands 5 and 6 for phenology) (O'Connell et al., 2017). The optimized cutoff value for FLATS was determined to be 0.1. Once the FLATS index was applied with the optimized cutoff value to the training dataset, the sensitivity of the model on the training data was 97%, while the specificity and accuracy were both 96% (Table 1). The AUC for the training data was 0.96 (Table 1).
We also used contingency table parameters to examine the testing and validation data for model efficacy. The sensitivity ranged from 81% to 86%, while specificity was 97% for both datasets (Table 1). The overall accuracy remained high, fluctuating between 93% and 95% (Table 1).

Vegetation biophysical analysis application
The FLATS index was applied to a focal dataset containing pixels within the PhenoCam's field of view (pixels 1-18) from 2013 to 2020, which originally consisted of 5,310 observations (full dataset). After filtering for clouds, cloud shadow, fill, and water using the Landsat QA band as well as filtering out scene edges, 2,663 observations remained. Finally, after applying and filtering with FLATS index using the optimized cutoff value, 1,970 dry observations were left, indicating a 63% reduction in data from the full dataset. A similar pattern in pixel reduction was detected in the study area (Fig. 1A), where there were originally 902,405 marsh pixels between 2013 and 2020. After Landsat QA filtering, 456,061 pixels remained, representing a 49.5% reduction in pixels. Filtering with FLATS further reduced the number of dry pixels remaining to 346,170, an additional 12.5% reduction bringing the total percentage of data removed from the analysis to 62%.
Examples to showcase the benefits of filtering data with FLATS vs other filtering methods can be seen in Fig. 4. Fig. 4A shows the raw NDVI phenological signal plotted for one pixel through time with ground truthed flooded pixels highlighted in teal. Fig. 4B is the pixel's NDVI values plotted after the Landsat QA band filtering was applied, still leaving 6 flooded pixels in the time series. Finally, Fig. 4C has no flooded pixels remaining after FLATS was applied to the dataset. Fig. 4D is a comparison of the smoothed splines of the 3 pixel filtering methods showcased in Fig. 4A-C where the raw NDVI phenological signal (black dotted line) was depressed compared to either of the filtered signals (orange Landsat QA and teal FLATS), particularly in 2015, 2016, and 2020 when substantial flooding occurred in the pixel. The Landsat QAfiltered and FLATS-filtered phenological signals had similar curves at the single pixel level, although the Landsat QA-filtered data still retained some pixels with minor flooding, which caused a dampening of signal when minor flooding events occurred as compared to the FLATS filtered signal. When results such as these are extrapolated over an entire study area (~3,000 or more pixels per date) significant dampening of the spectral signal will occur. The last panel in Fig. 4 (E) shows the smoothed splines of the three filter types that have been fitted to the monthly mean Table 1 Goodness of fit measures for the Flooding in Landsat Across Tidal Systems (FLATS), compared with the Landsat 8 Quality Assessment (QA) band. ED is explained deviance and AUC is Area Under the Curve. Sensitivity and specificity are the accuracy of classifying flooded and dry observations, respectively; accuracy is overall classification accuracy; N is total sample size. Flooded ID'd is the number of true flooded instances that were flagged and flooded missed is the number of true flooded observations that were missed (classified as dry rather than flooded). NDVI values for the entire study area through time (NDVI values were averaged per scene, if multiple scenes occurred in a single month those were also averaged to produce one NDVI value per month throughout the time series). Any scene that was determined to be 50% or more flooded using the Landsat QA and/or FLATS filters was removed from the analysis. Although the mean monthly NDVI splines (Fig. 4E) showed a similar overall pattern to single pixel splines (Fig. 4D), the major difference can be seen in the raw NDVI's diminished signal (particularly the dragging effect a flooded pixel has at the end of 2020) while the FLATS filtered values performed as expected.

Flooding analysis application
The majority (57%) of the pixels were flooded less than 20% of the time between 2013 and 2020, while a small percentage (9%) were flooded more than 50% of the time (Fig. 5A). These frequently flooded pixels were examined and were determined to be major creeks that cut through the marsh platform. The pixels that were flooded greater than 20% of the time but less than 50% are considered at higher risk for prolonged flooding and accounted for 34% of the area on the marsh surface. When a subset of 9 pixels containing short and medium form S. alterniflora located north of the flux tower were considered, the average flooding frequency of those pixels ranged from 10.35% (2014) to 27.78% (2020) and showed a positive trend when plotted through time. Flooding frequency was also considered on a monthly basis. The driest months were March, April, and December, with 31%, 42%, and 37% of marsh pixels being flooded, respectively (Fig. 5B). September, October, and November showed the most site-wide flooding, where an average of 99%, 93%, and 99% of marsh pixels experienced flooding in those months across all years (Fig. 5B). Flood maps were also created for each scene (with 3,000 or more pixels; i.e. mostly complete scenes without major cloud cover) to visually demonstrate FLATS' usefulness over using the Landsat QA band filtering method on a site wide basis but were not included here as they only represent brief snapshots in time and do not showcase the larger applications of FLATS over a time series. A subset of these maps can be found in the supplement (Fig. S1 A-C).
We also used Fort Pulaski water level data, which is relative to the North American Vertical Datum of 1988 (NAVD88), to study the depth of flooding through time relative to the LIDAR-derived DEM of the Sapelo Island marsh surface (Hladik et al., 2013). Fig. 6A shows the mean water level above the marsh surface averaged by scene from 2013 to 2020. There is a significant positive trend with a rise in the average water level over the 7-year study period (Fig. 6A). The FLATS flood frequency application can also be used to investigate relationships between flood proportion and other easily accessible variables. In our case, we used mean flood proportion by scene to explain MHHW values also obtained from Fort Pulaski (Fig. 6B). The relationship produced an R 2 of 0.305 and had a significantly positive relationship (p < 0.0001) while not showing distinct clustering by year. With every millimeter increase in the MHHW level we expect 0.125% more of the study area to be flooded. Our positive relationship mirrors the long-term increasing trend in sea level rise of 3.44 mm/year seen at the NOAA Fort Pulaski, GA station.

Discussion
This paper describes the creation of the Flooding in Landsat Across Tidal Systems (FLATS) model, which detects intermittent flooding within coastal marshes using Landsat 8 data. FLATS detected flooding within a S. alterniflora tidal marsh on Sapelo Island, GA, with an accuracy ranging between 93 and 96%. These high accuracies, along with high sensitivity values (81-97%), suggest that our model is effective at identifying flooded pixels. Once an optimized cutoff value was established, we then used FLATS to flag flooded pixels for removal to improve vegetation time-series analysis and to examine patterns of flood frequency.

Developing and optimizing FLATS
FLATS is an index based on spectral reflectance modeling used to predict binary flood status. FLATS was ground truthed through the use of marsh images provided by the PhenoCam and water level data from a pressure transducer located on the marsh surface. Water level data, while important in determining flood status, could not be solely relied upon to capture inundation levels over large areas of the marsh platform due to other factors that influence the movement and depth of water, such as wind speed and direction, vegetation type, density, and structure, as well as drainage patterns, and disturbance (e.g., crab burrows) (Hladik et al., 2013;Hughes et al., 2009;Nyman et al., 2009;Valentine and Mariotti, 2019). The flooding index was intended to accentuate two drastically different bands; a band that provides the greatest separation between the flooded and dry spectra (i.e., band 6) and a band that shares many of the same spectral characteristics (i.e., band 4). In addition to flood status, a phenology variable was deemed important when identifying model predictors, as emphasized by the spectral plot (Fig. 3). The inclusion of a phenology variable assists in distinguishing the slight differences in the visible spectral properties between the dormant season dry observations and the growing season flooded observations which, otherwise, would be difficult to discern. FLATS detects flooding in tidal marshes without the need for external water level data and is based on a combination of spectral bands that indicate the presence of inundation while accounting for seasonal variations in the vegetation.
We found that the normalized difference between Landsat 8 bands 4 and 6 as well as the addition of a phenology variable composed of the normalized difference between bands 3 and 4 constituted the best model for detecting flooding while accounting for seasonal variation (i.e., the FLATS model). We advocate for the use of the FLATS model instead of adjusting bands in the TMII because not only are there combinations of bands that more accurately detect intermittent tidal flooding (band 4 and band 6) but also because, unlike TMII, the use of a rolling mean function on a phenology variable is not optimal for Landsat 8 data due to its low temporal resolution (i.e., approximately 1 image every 8 days with combined revisit time frequency) as opposed to the MODIS revisit time (once per day). The NDWI index was chosen because it appropriately distinguished between flooded and non-flooded pixels, while the phenology index accurately captured the seasonal trends of the marsh vegetation. The selected NDWI band combination was similar to both McFeeters (1996), who used the green and near-infrared (NIR) bands to develop an index aimed at separating open water from vegetation and  soil in terrestrial environments, and Xu (2006), who modified McFeeters' index to include a SWIR band for build-up land areas instead of the NIR. Like O'Connell et al. (2017)'s NDWI, FLATS' chosen NDWI 4,6 consisted of a visible band and a shortwave infrared band because these wavelengths had the greatest separation between wet and dry plots.
Other attempts have been made to create and use indices to detect flooding in tidal marshes using satellite imagery, but none have successfully created a stand-alone index for moderate resolution Landsat 8 prior to FLATS. While TMII was created to utilize MODIS imagery which shares some bands centered at similar wavelengths with Landsat 8 (e.g., Landsat 8 band 5 with a center at 865 nm and MODIS band 2 with a center at 858.5 nm), there is not an equivalent for every band (O'Connell et al., 2017). Other studies such as Campbell and Wang (2020) have utilized the technique of selecting Landsat 8 equivalent bands to be used in the MODIS based TMII model in order to eliminate flooded pixels from their Landsat-based analysis. We, however, found that the simple substitution of Landsat 8 bands creates a poorly explained model (ED -52%) when compared to our FLATS model (ED -82%). The transferability of TMII from MODIS to Landsat 8 is likely inhibited by factors such as the sensors' field of view, signal-to-noise ratio, and, subsequently, its spatial resolution. As the field of view increases the proportions of vegetation to water change, which alters the spectral reflectance, and, ultimately, changes the model's performance. The type of scanning conducted by each satellite (whisk-broom for MODIS and push-broom for Landsat 8) alters the amount of time each pixel is observed (Ren et al., 2014;Wolfe et al., 2002). The push-broom scanner on Landsat 8 ′ s Operational Land Imager (OLI) sensor leads to a higher signal-to-noise ratio than that of MODIS' whisk-broom scanner Morfitt et al., 2015). The average signal-to-noise ratio of the bands used in TMII is 194.5 as compared to FLATS' 256.25, indicating that lower noise is inherent in Landsat 8 data (Morfitt et al., 2015;Xiong et al., 2008). FLATS outperforms TMII when using moderate resolution imagery, suggesting that these indices have spatial resolution cutoffs in addition to limitations brought about by the satellites they were initially created under.
Evaluation of the FLATS model involved the examination of goodness of fit measures (Table 1). The model did well in training, testing, and validation datasets, particularly when compared to data that were only Landsat QA-filtered. The sensitivity of the Landsat QA mask was less responsive to flooding than that of the FLATS model. One indication that the FLATS model performed best was the number of pixels correctly identified as flooded which was higher for FLATS (ranging 30-59) than for flooding identified from Landsat QA information (ranging 23-45). The FLATS optimized cutoff value of 0.1 was selected due to its high sensitivity while maintaining a balance with accuracy and specificity.

Applications of FLATS
A traditional application of FLATS involves filtering out flooded pixels from a vegetation time-series for future analysis. We used NDVI to compare the relative difference between unfiltered, Landsat QA-filtered, and FLATS-filtered Landsat 8 data through time (Fig. 4). This analysis emphasizes the impact various types of filtering can have on the data used when conducting remote sensing vegetation biophysical analysis. We do not claim that NDVI is the ideal index for marsh vegetation analysis, but instead denote it as a widely used and easily recognizable index, particularly for comparative analysis. Failure to filter out inundated pixels, as seen in Fig. 4A, leads to a distorted phenological signal demonstrated clearly by the flooded pixel present at the end of 2020 but also by those in 2015 and 2016. Accurately flagging flooded pixels is key to generating an accurate satellite-based vegetation time-series for subsequent analysis as partially flooded pixels dampen the spectral signature of vegetation. Filtering with FLATS provides a phenological curve at a magnitude expected in an S. alterniflora marsh. The smoothed spline of the Landsat QA filtered data generally mimics that of the FLATS filtered data for the single-pixel, but with dampened peak magnitude (Fig. 4D). When comparing the splines of the single-pixel analysis ( Fig. 4D) with that of the mean of all pixels per scene per month (Fig. 4E), we can see more variation between each of the filter types at a greater magnitude. This variation can be explained in part by the larger overall area represented in Fig. 4E (~4 km 2 ) than that in Fig. 4D (~900 m 2 ). These two scales also experience different ranges in elevation, which can lead to variations in the intensity of the NDVI values due to the changes in vegetation type and structure (e.g., tall form vs. short form S. alterniflora). Overall, Fig. 4D and Fig. 4E indicate that flood filtering via FLATS arrives at the most accurate representation of the vegetation's signal, and the three filtering types perform comparably across scales, i.e., FLATS is the best filtering method at the landscape scale and the pixel scale, suggesting FLATS provides an ideal flood filtering method prior to vegetation biophysical analysis such as predictive modeling of biomass and gross primary production (GPP).
Mapping the frequency of flooding is another potential application for FLATS. Though FLATS does not capture every flood event due to the intermittent observations of Landsat 8, a sufficiently long time series of the kind we show here results in a random tidal sampling frequency which sufficiently accounts for sample variations and provides meaningful patterns of flooding frequency on an inter-annual scale (see Fig. S2 in the supplement). With this application in mind, we distinguished creeks by their noticeably linear features and high frequency of flooding. We consider the correct classification of these linear features (creeks, streams, etc.) as water in the flood frequency maps (Fig. 5) to be further validation that our model readily identifies flooded instances. These creek areas on the marsh surface flooded greater than 50% of the time during the 7-year study period. Creeks often appear as mixed pixels in Landsat imagery due to their relatively smaller size and the density of tall form S. alterniflora that usually borders them. The majority of the pixels found on the marsh platform are only flooded 20% of the time or less. The remaining pixels (whose flooding frequency ranges between 20% and 40%) are not located in established creeks but yet are still impacted by relatively high rates of flooding, potentially leaving these areas more vulnerable to the impacts of sea level rise (Vitousek et al., 2017). An advantage of FLATS is that it can highlight the cumulative frequency of flooding, which helps distinguish areas vulnerable to ponding when areas traversed heavily by streams are excluded. Wetland fragmentation via migrating creeks and ponding presents a major challenge to continued marsh resilience (Duran Vinent et al., 2021;Mariotti, 2016). FLATS can assist in the vital effort to monitor the size and distribution of ponding. This water connectivity and distribution monitoring are important for decerning the severity of marsh fragmentation and establishing patterns to help determine if marsh recovery or drowning will occur (Day et al., 2011;Duran Vinent et al., 2021). The average water level relative to the marsh surface appears to be increasing over time which could, when coupled with the flood frequency data, indicate potential sea level rise related complications for any marsh (Fig. 6A).
Although our model does not produce water level data, MHHW helps predict the mean flood proportion derived from FLATS (Fig. 6B). This analysis approach is particularly useful if a high-resolution digital elevation model is not available, as it was in our study, because the mean flooding proportion obtained through the use of FLATS can be used as a proxy to estimate MHHW levels. MHHW can also be used as a low signal covariant for mean flood proportion in other types of statistical models (e.g., machine learning). However, flooding as estimated from FLATS itself would be better in such models, as it can provide direct observations of spatially explicit flooding impacts.
As a tidal flooding indicator, FLATS is a useful tool for a variety of applications and disciplines. In addition to those discussed here (vegetation biophysical time-series analysis and flood frequency detection), FLATS could also be used by restoration scientists, ecologists, site planners, and others for a variety of coastal flooding applications. This index could also be used to examine the long-term effects of disturbances such as wrack on flooding in tidal wetlands. The next steps include applying FLATS to both Landsat 8 and upcoming Landsat 9 to build a denser time-series and re-parametrizing FLATS to other satellites, including higher resolution Sentinel-2 as well as those with long-term datasets such as Landsat 5. One of the limitations is that we developed FLATS using data from one site and one marsh species due to the lack of established coastal wetland PhenoCam data in other marshes. However, there is great potential for this index to be scalable at other locations because of its vegetation index centric formulation, including those with varying flood regimes, and where other tidal marsh species are dominant (e.g., Juncus roemerianus or S. patens), emphasizing another direction in need of empirical exploration. Our method is easily reproducible for researchers with access to repeat photography in other marsh environments. Future improvements to this index could include model formation via machine learning and the addition of other biophysical parameters, which could both enhance the scalability and applicability of FLATS. An implementation of the index on high-resolution Sentinel-2 data could improve our understanding of ponding and creek migration in coastal marshes, while use in Landsat 5 could provide extensive flood frequency trend data for analysis over the past 35 years.

Conclusions
Filtering Landsat 8 data via FLATS can improve remote sensing models of coastal marsh vegetation by reducing the presence of spectral noise from inundation. These tide-filtered models can show phenological signals that more accurately match the dynamics at play in these indispensable coastal environments. Consequently, researchers should be able to generate improved estimates of vegetation features, such as habitat characteristics or productivity and carbon storage potential. To this end, we recommend the incorporation of FLATS into Landsat 8 preprocessing routines for tidal systems. FLATS also provides insight into shifting inundation patterns, including the migration of creeks and ponding on the marsh surface. The presence, frequency, and intensity of these flooding configurations can be an indicator of marsh resiliency or drowning. FLATS improves our capability to accurately study tidal marshes, which will have lasting impacts on our ability to better understand our changing environment.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
Data will be made available on request.