The Potential of Satellite Remote Sensing Time Series to Uncover Wetland Phenology under Unique Challenges of Tidal Setting

: While growth history of vegetation within upland systems is well studied, plant phenology within coastal tidal systems is less understood. Landscape-scale, satellite-derived indicators of plant greenness may not adequately represent seasonality of vegetation biomass and productivity within tidal wetlands due to limitations of cloud cover, satellite temporal frequency, and attenuation of plant signals by tidal ﬂooding. However, understanding plant phenology is necessary to gain insight into aboveground biomass, photosynthetic activity, and carbon sequestration. In this study, we use a modeling approach to estimate plant greenness throughout a year in tidal wetlands located within the San Francisco Bay Area, USA. We used variables such as EVI history, temperature, and elevation to predict plant greenness on a 14-day timestep. We found this approach accurately estimated plant greenness, with larger error observed within more dynamic restored wetlands, particularly at early post-restoration stages. We also found modeled EVI can be used as an input variable into greenhouse gas models, allowing for an estimate of carbon sequestration and gross primary production. Our strategy can be further developed in future research by assessing restoration and management effects on wetland phenological dynamics and through incorporating the entire Sentinel-2 time series once it becomes available within Google Earth Engine.


Introduction
Remote sensing products are increasingly important to facilitate upscaling of critical ecosystem services such as carbon (C) sequestration from local monitoring sites to regional extents of management and planning [1][2][3]. Such efforts often utilize spatially explicit time series of spectral vegetation indices (or "greenness") that can provide useful proxies of plant biomass, coverage, and photosynthetic activity (e.g., [4][5][6]). Predicting such vegetation properties at different stages of a growing season enables a better understanding of various biotic and abiotic controls on plant and ecosystem functioning [7][8][9] and modeling of ecosystem services targeted by management, restoration, and conservation efforts. However, in contrast to well-studied "upland" ecosystems, this potential remains less certain in coastal environments such as tidal marshes, where high short-term spectral variability precludes accurate characterization of vegetation (e.g., [10][11][12][13][14]). This uncertainty contributes to gaps in upscaling and budgeting of these important "blue carbon" ecosystems and projecting their future under climatic, hydrological, and anthropogenic change drivers [15][16][17].
From a remote sensing perspective, the key driver behind this issue is variable tidal flooding, which leads to attenuation of vegetation signals [14,18] and limits the ability of spectral reflectance from snapshot remote observations to indicate relevant aspects of plant physiology and abundance. Furthermore, such attenuation non-linearly increases with the magnitude of flooding [14], limiting performance of vegetation indicators more substantially in frequently flooded environments such as low-elevation zones of tidal marshes. Yet, correcting for flooding is not a straightforward task because the net effect of attenuation depends not only on tidal stage and local ground elevation, but also on vegetation canopy properties (height, biomass, structure) that define to what extent vegetation can be obscured by flooding at a given point in a season [4,19]. The interplay among these factors affects how much vegetation is effectively "visible" to a remote sensor at the time of overpass in a given part of a wetland.
Characterizing this interplay in a tidal setting is extremely difficult due to the lack of relevant geospatial and earth observation data at sufficiently high temporal frequencies.
Among various instruments and platforms, light detection and ranging (LiDAR) systems are currently best suited for characterizing both ground elevation and vegetation height and for modeling local inundation based on these characteristics [20][21][22]. However, in the absence of continuous global LiDAR observations at high spatial resolutions, acquisitions of such data at local and regional scales remain patchy, sporadic, and too costly for matching temporal frequencies of tidal variation (e.g., from hourly to biweekly). Furthermore, in tidal wetlands the accuracy of LiDAR-derived metrics may be compromised by variation in the soil/sediment and by the obstruction of soil by dense vegetation and dead biomass, or litter [21,23]. These issues limit the potential for tracking spatial and temporal (seasonal) dynamics of wetland vegetation biomass and height which modulate tidal attenuation during the growing season. This limitation is amplified by the common lack of alternative up-to-date wetland vegetation products, e.g., plant community-or species-level maps at sufficiently high spatial resolution and temporal frequency to help approximate changes in biomass and canopy properties at a local scale [24].
Previous studies have offered several strategies to help circumvent these challenges to summarize spatio-temporal variability of wetland vegetation for tidal ecosystem mapping, modeling, and phenological assessments. For example, the first USA national-scale model of aboveground tidal marsh C stocks [24] focused on the relationships between field-measured biomass and temporally matching satellite products across a comprehensive range of sites and vegetation types, often sampled at or close to the peak biomass stage. Peak biomass conditions reduce the magnitude of tidal attenuation and provide a way to characterize the overall biomass and C stock representative of a given ecosystem, which can be useful for annual C budgeting [25,26]. In other cases, it may be necessary to explicitly consider seasonal variation in vegetation proxies to indicate shorter-term changes in biogeochemical fluxes and wetland habitat for specific seasonal phases of animal ecological cycles [5][6][7], which in tidal settings requires addressing flooding-induced attenuation specific to the time and location of satellite observations [16]. A study of tidal marsh vegetation phenology in three Pacific Coast sites [27] incorporated the offsets in water level height between the tidal gauge data from individual site locations into models of vegetation phenology, and additionally removed satellite images collected at high water levels. Alternatively, high temporal frequency remote sensing data can be used to identify and flag conditions when a given pixel is likely to be flooded relative to its overall phenological behavior, such as the Tidal Marsh Inundation Index [16] based on 500 m MODIS satellite reflectance products in green and shortwave infrared electromagnetic regions. Such methods, however, may be less feasible when satellite time series with optimal spatial resolution are sparser and compromised by clouds, which increases the uncertainty associated with exclusion of high-tide images, especially in sites with no accurate local tidal data. Furthermore, from an ecosystem modeling perspective, high-tide dates are still informative of how tidal flooding affects ecosystem functioning; for instance, the rate of carbon uptake by flooded vegetation during photosynthetically active daytime hours may be expected to differ between flooded and non-flooded conditions [3,28]. As the contribution of tidal vegetation biomass to remote sensing signals is expected to increase with the increment in biomass and height over the course of the growing season [28], the degree of attenuation reduces during the peak biomass stages. Thus, rather than predicting the total amount of vegetation, studies aiming to support biogeochemical and habitat modeling may need to focus on capturing the characteristic range of spectral indicator values as a "phenological envelope," the upper and lower bounds of which may change throughout the growing season. However, the uncertainty around the specific magnitude of tidal attenuation and its local spatial variability, amplified by sparseness of remote sensing data, reduces the ability to infer such a range based on a single-date or even a single-year phenological time series.
Alternatively, this uncertainty can be navigated by using phenological information from historical remote sensing observations encompassing a greater variety of ecosystem signals under different stages of tidal flooding and vegetation growth than single-year series. A greater range of historical image dates provides a more comprehensive idea on the "upper" and "lower" boundaries of vegetation indicators possible over a course of a tidal cycle at a given seasonal stage, as well as more general phenological variability and anomalies caused by unusual years or disturbances (e.g., [29,30]). To predict such a phenological range more realistically for a given analysis year, historical information can be further coupled with environmental variables specific to this year, such as climatic indicators affecting phenological timing [27]. However, the relative importance of such indicators may vary geographically [24,27] and therefore may need to be explicitly investigated in region-specific applications.
Regional-scale modeling of seasonal changes in greenness also requires understanding where historical phenological information might be less informative, e.g., where historical dynamics are not sufficiently representative of current conditions. Although spectral signals of tidal systems frequently exhibit high short-term variability due to flooding cycles, certain types of landscape processes and land use practices might create periods of more pronounced instability and/or directional change. Wetland restoration sites, for example, can experience phases of vegetation colonization and expansion after the initial perturbation [31], during which changes in vegetation signals, functions and canopy structure may be different from the later stages of more "stable" dynamics [5,6,32]. However, when and to what extent the "characteristic" phenological variability of wetlands may be sensitive to such perturbations, and whether such sensitivity exceeds the overarching effects of climatic controls on remotely sensed phenological signals, is not well understood.
This study developed and investigated a strategy to model annual trajectories of a greenness indicator in a tidally affected estuarine setting by combining per-pixel historical multi-year phenological variability with relevant landscape properties and current-year climatic characteristics from open-access data. This approach aims to provide a cost-effective strategy to predict wetland greenness phenology when usable cloud-free satellite data at suitable spatial resolution are sparse and spatially explicit, temporally frequent information on vegetation composition and structure is not available. We focused this analysis on the San Francisco Bay Estuary, California, USA, where there is an increasing interest in wetland restoration and regional-scale monitoring and modeling of restoration and management outcomes, including ecosystem function and productivity and C sequestration. These efforts require more comprehensive and seasonally detailed vegetation indicator information than may be possible from sparser satellite imagery. Regional wetlands include both remnant historical marshes and more recently restored and managed wetlands Remote Sens. 2021, 13, 3589 4 of 28 with a diversity of hydrological regimes and vegetation types affecting ecosystem-level phenology from a remote sensor perspective. This provides a useful opportunity to test the suitability of multi-year phenological variability to inform extrapolation of seasonal dynamics for a given year. Using seven years of open-access satellite data, we investigated the following questions: (1) Can phenological history of greenness inform prediction of single-year greenness from discrete satellite time series? (2) Which climatic and other environmental covariates should be included in such predictive modeling to account for the unique phenological variation of the prediction year and local wetland conditions? and (3) To what extent is such predictive modeling sensitive to wetland restoration and management status, including time since restoration?

Study Area
The study focused on tidal wetlands within the San Francisco Bay Estuary (Bay Area) as defined by the San Francisco Estuary Institute ( Figure 1C, [33]). The study area was also constrained by the elevation dataset we used, which only includes salt and brackish tidal marshes within the Bay Area [34]. The Bay Area can be broken up into three regions, South San Francisco Bay (South Bay), San Pablo Bay, and Suisun Bay ( Figure 1B). For model training, we focused on 14 areas within the Bay Area tidal wetlands ( Figure 1B). The areas were semi-randomly selected and distributed across the Bay Area within known tidal wetland locations. The Bay Area has a dry season roughly between May and October, with the majority of the precipitation occurring in December through February [35]. Most of the freshwater supply is from the Sacramento-San Joaquin River system and enters the Bay Area through the Suisun Bay; only 10% of freshwater flows from smaller streams and creeks such as the Coyote Creek in the South Bay. This makes the South Bay more saline than the Suisun Bay [36].
Due to the availability of eddy covariance flux tower data at two locations within our study area (Figure 1), we also compared flux tower-derived gross primary production with modeled plant greenness within Eden Landing Ecological Reserve (Eden Landing) and Rush Ranch Open Space Preserve (Rush Ranch). Eden Landing is located within the South Bay and is managed by the California Department of Fish and Wildlife. The site is part of the South Bay Salt Pond Restoration project, where restoration started in 2003 and aimed at restoring former industrial salt ponds to their more natural state. Restoration at Eden Landing was completed in 2008 and includes roughly 835 acres of tidal salt marsh. The dominant plant community is Sarcocornia pacifica, however, Spartina foliosa, algal mats, and other vegetation are also present. Rush Ranch is located within the Suisun Bay and is part of the National Estuarine Research Reserve system. Rush Ranch is one of the few undiked brackish tidal marshes in San Pablo Bay and the marsh remains relatively undisturbed. The plant community here is more diverse and is vegetated with plants such as Distichlis spicata, S. pacifica, Schoenoplectus americanus, and Lepidium latifolium.

Datasets
We used a variety of spatial datasets within this study (Table 1) to capture the potential precursors of wetland phenological dynamics and their regional variability. To select and classify wetland types, we used San Francisco Bay Area EcoAtlas Modern Baylands (Baylands) data ( Figure 1C, [33]). The Baylands data are a shapefile containing wetlands within the Bay Area. The wetlands are further classified into different categories based on characteristics such as tidal influence, management, and elevation. This dataset allowed us to model the EVI within known tidal wetland areas and analyze if wetland type played a role in model error. To further test for model error, we used a habitat restoration dataset (Habitat Projects). The Habitat Projects dataset is a spatial dataset and includes restoration projects within California, with a focus on wetland and coastal ecosystems [37]. To estimate vegetation greenness, we used Landsat 8 (Landsat) surface reflectance time series data accessible via Google Earth Engine (GEE), the open-source, cloud-based, web-enabled applied programming interface for multi-scale geospatial analyses [40]. This 30-meter spatial resolution product includes data from 2013 to the present and is available as atmospherically corrected using Land Surface Reflectance Code (LaSRC; [41]). With the image Quality Assessment (QA) band, we further reduced image noise by masking out clouds (bit attribute 5) and shadows (bit attribute 3). We selected all Landsat scenes ( Figure A1) within the study area between 2013 and 2019. To estimate plant greenness, we calculated the Enhanced Vegetation Index (EVI). The EVI is less sensitive to soil and oversaturation than the commonly used Normalized Difference Vegetation Index, and several researchers have successfully used the EVI in wetland settings [5,[42][43][44]. We calculated the Enhanced Vegetation Index using the following equation: where NIR is the near infrared band, red is the red band, and blue is the blue band. To calculate EVI history, we subset the Landsat data between 2013 and 2018 and calculated the two-week EVI mean regardless of year. To do this within GEE, we forced all the years to be the same placeholder value, while keeping the month and day values as is. We then summarized the data by taking the mean EVI value within 14-day periods. This gave us 26 spatially explicit EVI mean values summarized from the historical dataset.
To help account for interannual variations in plant growth, we included temperature and drought variables in our model, also using datasets accessible via GEE (Table 1).
Previous studies found that aboveground biomass and plant greenness linearly correspond to temperature [45]. Similarly, researchers found that drought, or Palmer Drought Severity Index (PDSI), and aboveground biomass within wetlands are linearly correlated [46]. For weather data, we used PRISM Climate Group's Daily Spatial Climate Dataset data (which includes temperature) and used the PDSI provided through the GRIDMET Drought dataset [39,47]. PDSI values range from −10 to 10 where values less than 0 indicate drier than normal conditions and values more than 0 indicate wetter conditions. Larger magnitudes indicate higher dry or wet signals.
We used two elevation datasets for this study. The first digital elevation model (DEM) is a 10-meter spatial resolution dataset resampled to 30 meters to match Landsat resolution. The dataset is a composite of DEMs acquired by the United States Geological Survey and the Department of Water Resources and is referenced to elevations above NAVD88 [38]. The second elevation dataset is also a composite of elevation datasets derived mostly from LiDAR but is referenced to elevations above mean higher high water (MHHW) and has a spatial resolution of 5 meters resampled to the 30-meter resolution of Landsat [34]. Tidal range is not constant across the Bay Area; using elevations above MHHW is a way to account for this difference. It is also a better determinant of flooding duration and the tidal influence on vegetation. Furthermore, wetlands within the Suisun Bay region have experienced a large amount of land subsidence. This gives the wetlands a lower elevation above NAVD88 than other wetlands within the Bay Area. However, due to the influence of dikes and levees, these wetlands do not experience the same tidal flooding as another wetland at a similar elevation above NAVD88.

Model Building and Predicting
In Google Earth Engine, we created a point at each Landsat pixel within the 14 model development sites, giving us a total of 239 points. Within this point dataset, we attached the elevation data, EVI history, and the 2019 Landsat 8-derived EVI time series. We then found the PDSI and temperature data with dates that most closely matched the 2019 Landsat acquisitions dates and attached them to the point dataset. The data were then exported into the R statistical program where we built and tested a series of predictive linear models to estimate the EVI at each Landsat pixel within the study area. Given the lack of previous studies and novelty of this analysis, we chose a linear regression approach as the main framework due to its relative ease of interpretation and comparison of different predictors (their coefficients and contributions to the model), and more transparent model error structure and performance assessment than with more complex and less transparent algorithms. Please note that even though we used linear multi-variate regression, not all predictors were included in a linear form, some were transformed to quadratic functions based on their statistical distribution, model diagnostics, and/or physical interpretation. For example, day of year (DOY) was used as a non-linear term because solar radiation follows a non-linear trajectory during a growing season. We then compared the models using Akaike information criterion (AIC) analysis [48]. We built the model using the following variables, a subset of the variables, and/or a variable transformation: where DOY is day of year, EVI h is historically observed EVI, PDSI is Palmer Drought Severity Index, Temperature is minimum, maximum, or mean temperature ( • C), and Elevation is elevation (cm) above NAVD88 or mean higher high water (cm). β 0 is the intercept while β 1 − β 5 are the fitted model coefficients corresponding to the predictor variables.
Using the best fit model and climate variables from 2019, we created a modeled 2019 phenological envelope by using the upper and lower bounds of the 95% confidence interval around predictor variable coefficients. This gave us maximum and minimum expected EVI values for 2019 ( Figure 2). We did this for both the 14 study areas and for the tidal wetlands within the entire Bay Area. To assess model accuracy, we calculated the root mean square error (RMSE) modified Willmott's index of agreement (dmod) for the 2019 mean modeled EVI valu sus Landsat observed EVI values [49]. The modified Willmott's index of agreem way to calculate model performance. Values range from 0 to 1, and 1 indicates a model fit. This index improves upon the original Willmott's index. We used 95 r validation points throughout the Bay Area to calculate RMSE and dmod. We generat dom points in ArcGIS, extracted the Landsat observed EVI and the mean model then exported the data into R. We calculated the RMSE using the following formu where n is the sample size, EVIobs is the Landsat observed EVI, and EVImod is the m mean EVI. We calculated dmod using methods outlined by [50]. We also estimated model error within the entire Bay Area based on the pheno envelopes. We found the mean difference between the predicted EVI values and ob EVI values obtained from Landsat in 2019. To calculate the difference, we took the L observed EVI minus the predicted values and set the difference to zero if the ob EVI values fell within the 95% confidence interval, or the phenological envelope predicted values. We calculated the difference for each available 2019 Landsat sce found the mean difference by taking the average of the difference maps, giving spatial map of mean difference. This gave us one spatially explicit mean differenc

Additional Satellites in EVI Histories
To test if increasing the number of satellites used for creating EVI history im our results, we expanded our satellite imagery to include Landsat 7 and the Harm Sentinel-2 imagery. We used the same calibrated model described above but use ferent EVI history input based on what satellite(s) were used to calculate the mean To assess model accuracy, we calculated the root mean square error (RMSE) and the modified Willmott's index of agreement (d mod ) for the 2019 mean modeled EVI values versus Landsat observed EVI values [49]. The modified Willmott's index of agreement is a way to calculate model performance. Values range from 0 to 1, and 1 indicates a perfect model fit. This index improves upon the original Willmott's index. We used 95 random validation points throughout the Bay Area to calculate RMSE and d mod . We generated random points in ArcGIS, extracted the Landsat observed EVI and the mean modeled EVI, then exported the data into R. We calculated the RMSE using the following formula: where n is the sample size, EVI obs is the Landsat observed EVI, and EVI mod is the modeled mean EVI. We calculated d mod using methods outlined by [50]. We also estimated model error within the entire Bay Area based on the phenological envelopes. We found the mean difference between the predicted EVI values and observed EVI values obtained from Landsat in 2019. To calculate the difference, we took the Landsat observed EVI minus the predicted values and set the difference to zero if the observed EVI values fell within the 95% confidence interval, or the phenological envelope, of the predicted values. We calculated the difference for each available 2019 Landsat scene and found the mean difference by taking the average of the difference maps, giving us one spatial map of mean difference. This gave us one spatially explicit mean difference map.

Additional Satellites in EVI Histories
To test if increasing the number of satellites used for creating EVI history improved our results, we expanded our satellite imagery to include Landsat 7 and the Harmonized Sentinel-2 imagery. We used the same calibrated model described above but used a different EVI history input based on what satellite(s) were used to calculate the mean historical fortnight EVI.
Landsat 7 is available within GEE, while the Harmonized Sentinel-2 imagery currently is not. The Harmonized Sentinel-2 data are available on NASA's website, atmospheri-Remote Sens. 2021, 13, 3589 9 of 28 cally corrected and calibrated to match Landsat 8 in pixel resolution (30 m) and spectral response [41]. Downloading and uploading the Harmonized Sentinel-2 data must be done manually, therefore we did not include all available images. We selected 49 Harmonized Sentinel-2 images that had no or low cloud cover and focused on dates that had fewer Landsat scenes. We then calculated the mean historical EVI using the combination of Landsat 7 and Landsat 8 then an additional mean EVI history using Landsat 7, Landsat 8, and the Harmonized Sentinel-2 imagery. This gave us two additional EVI history datasets, allowing us to analyze if adding additional satellites to the historical datasets improved our original model using EVI history based on Landsat 8 only. Using the two additional EVI history datasets, we predicted the EVI for 2019 throughout the Bay Area and calculated the difference similarly as described above.

Eddy Covariance Flux Tower Data
Eddy covariance data were collected at Rush Ranch and Eden Landing following similar approaches which allowed for calculation of 30-min average fluxes of CO 2 between the marsh and atmosphere at the ecosystem scale [51]. Each tower was equipped with an open path CO 2 and H 2 O infrared gas analyzer (LI-7500A/RS, LI-COR, Lincoln, NE, USA) and a sonic anemometer (Windmaster Pro, Gill Instruments Ltd., Lymington, Hampshire, UK and CSAT3, Campbell Scientific, Logan, UT, USA for Eden Landing and Rush Ranch, respectively). The instruments were mounted at 4.3 m and 3.4 m above the marsh surface for Eden Landing and Rush Ranch, respectively. Instruments were mounted on the northwestern side of each tower as both sites have predominant wind from the west (~270 • ). Raw data were recorded at 20 Hz frequency and the high-frequency data were processed to half hourly fluxes using EddyPro, with a missing samples allowance of 10%. Data and metadata are publicly accessible on the AmeriFlux website (US-Edn and US-Srr; https://ameriflux.lbl.gov/). Tower footprints were calculated using the Kljun model [52,53]. Currently, Rush Ranch eddy covariance data are available for between 2013 and 2017 and Eden Landing eddy covariance data are available for between 2018 and 2020.
Partitioning of net ecosystem exchange (NEE) into gross primary productivity (GPP) and ecosystem respiration at Rush Ranch followed Knox et al. [54]. Partitioning of NEE at Eden Landing was performed using an artificial neural network with input variables of sensible and latent heat flux, water table height, date, and tidal difference (the change in water table height between each half hour).

PEPRMT GPP Model
We used EVI predictions from Rush Ranch and Eden Landing to run a simple light use efficiency model previously applied to restored wetlands in the Sacramento-San Joaquin River Delta (PEPRMT model; [5,7]). The model was parameterized using Rush Ranch GPP measurements in 2016 and then applied to Eden Landing 2019 in a validation exercise [7].

Satellite Image Count
The number of satellite images used for assessing the EVI history varied based on location within the Bay Area. Overall, the San Pablo Bay region had the highest image availability, and the South Bay had the fewest ( Figure 3A). When factoring in the fortnight number, we further see that the South Bay consistently has the lowest number of available Landsat scenes ( Figure 3B). There is also a seasonal pattern in available images, with the beginning and end of the year having the lowest image availability.

Model
Based on the model with the lowest AIC value (Table A1, dmod = 0.74), we found the following regression model to best predict the data: DOY is the day of year, EVIh, is the EVI history, PDSI is the Palmer Drought Severity Index, Tmin is the minimum temperature (°C), and MHHW is elevation above mean higher high water (cm). The DOY and MHHW variables have a polynomial model fit. β0 represents the model intercept. Each subsequent β coefficient corresponds to the predictor variable. The β coefficients for each predictor vary based on if we are calculating the lower bound, mean, or upper bound EVI values ( Table 2). Using the upper and lower bounds of the 95% confidence interval around the mean β coefficients, we modeled the 2019 EVI throughout the entire Bay Area tidal wetland landscape. This gave us 26 images of spatially explicit EVI envelopes, each representing a specific date range. Mean EVI values for the specific date range were also mapped, allowing for a spatial visual representation of plant greenness throughout the Bay Area for a given point in time (Figure 4). The summertime snapshot shows distinct areas of lower and higher plant greenness. The eastern portion of the San Pablo Bay has lower greenness than the rest of San Pablo Bay, the Suisun Bay has the overall highest greenness, while the South Bay has the lowest greenness.

Model
Based on the model with the lowest AIC value (Table A1, d mod = 0.74), we found the following regression model to best predict the data: DOY is the day of year, EVI h , is the EVI history, PDSI is the Palmer Drought Severity Index, Tmin is the minimum temperature ( • C), and MHHW is elevation above mean higher high water (cm). The DOY and MHHW variables have a polynomial model fit. β 0 represents the model intercept. Each subsequent β coefficient corresponds to the predictor variable. The β coefficients for each predictor vary based on if we are calculating the lower bound, mean, or upper bound EVI values ( Table 2). Using the upper and lower bounds of the 95% confidence interval around the mean β coefficients, we modeled the 2019 EVI throughout the entire Bay Area tidal wetland landscape. This gave us 26 images of spatially explicit EVI envelopes, each representing a specific date range. Mean EVI values for the specific date range were also mapped, allowing for a spatial visual representation of plant greenness throughout the Bay Area for a given point in time (Figure 4). The summertime snapshot shows distinct areas of lower and higher plant greenness. The eastern portion of the San Pablo Bay has lower greenness than the rest of San Pablo Bay, the Suisun Bay has the overall highest greenness, while the South Bay has the lowest greenness.

Phenological Envelope
Extracting the EVI envelopes for each of the 14 model development sites, we mined the phenological curves capturing seasonal EVI variability for each site (Figu

Phenological Envelope
Extracting the EVI envelopes for each of the 14 model development sites, we determined the phenological curves capturing seasonal EVI variability for each site ( Figure 5). This gave us the expected minimum and maximum EVI values throughout 2019 for each site. Phenological envelopes varied for each site, but overall sites had the lowest EVI during the winter with the EVI increasing in the spring, leveling off in the summer, then decreasing in the fall. Spatially, summer EVI values were highest in the Suisun Bay region ( Figure 5A), lowest in the South Bay region (Figure 5C), and summer EVI values in the San Pablo Bay region ( Figure 5B) were in between. Across these sites, the minimum modeled EVI was −0.02 and occurred in the winter while the maximum EVI was 0.63, occurring in summer.
Using the mean predicted EVI for Rush Ranch and Eden Landing, we ran a simple light use efficiency model (PEPRMT) [5,7] to test its ability to help improve predictions of GPP in tidal wetlands. We parameterized the model using 2016 data from Rush Ranch and validated the model using 2019 data from Eden Landing. Eddy covariance data was not available for 2019 within Rush Ranch, therefore we used 2016 data for GPP. Model performance was strong at Rush Ranch (d mod = 0.72, RMSE = 1.74; Figure 6B) and Eden Landing (d mod = 0.63, RMSE = 0.96; Figure 6A). This gave us the expected minimum and maximum EVI values throughout 2019 for each site. Phenological envelopes varied for each site, but overall sites had the lowest EVI during the winter with the EVI increasing in the spring, leveling off in the summer, then decreasing in the fall. Spatially, summer EVI values were highest in the Suisun Bay region ( Figure 5A), lowest in the South Bay region ( Figure 5C), and summer EVI values in the San Pablo Bay region ( Figure 5B) were in between. Across these sites, the minimum modeled EVI was −0.02 and occurred in the winter while the maximum EVI was 0.63, occurring in summer.
Using the mean predicted EVI for Rush Ranch and Eden Landing, we ran a simple light use efficiency model (PEPRMT) [5,7] to test its ability to help improve predictions of GPP in tidal wetlands. We parameterized the model using 2016 data from Rush Ranch and validated the model using 2019 data from Eden Landing. Eddy covariance data was not available for 2019 within Rush Ranch, therefore we used 2016 data for GPP. Model performance was strong at Rush Ranch (dmod = 0.72, RMSE = 1.74; Figure 6B) and Eden Landing (dmod = 0.63, RMSE = 0.96; Figure 6A).

Model Accuracy
We assessed model accuracy both for validation points throughout the Bay Area and for tidal wetlands within the entire Bay Area. Model accuracy was high for the validation points; d mod = 0.67, RMSE = 0.11 (Figure 7). Throughout the Bay Area, the mean difference, or model error, varied but overall mean difference was low ( Figure 8A). In general, the model tended to overestimate the EVI, with distinct areas of higher error within the San Pablo Bay region. When using EVI history that included Landsat 7 and Landsat 8 (EVI L7L8 ) or using the EVI history with Landsat 7, Landsat 8, and Harmonized Sentinel-2 (EVI L7L8Sent ), we see minor differences within model error (Figures 6 and 8). The EVI mean difference using EVI L7L8Sent history appears to be closer to zero (i.e., more yellow on the map, Figure 8A), indicating lower model error, with a few regions in the South Bay having higher model error ( Figure 8A). Furthermore, the EVI envelopes have minor changes in shape when using multiple satellites to calculate EVI history versus using only Landsat 8 ( Figure 6).

Model Accuracy
We assessed model accuracy both for validation points throughout the Bay Area and for tidal wetlands within the entire Bay Area. Model accuracy was high for the validation points; dmod = 0.67, RMSE = 0.11 (Figure 7). Throughout the Bay Area, the mean difference, or model error, varied but overall mean difference was low ( Figure 8A). In general, the model tended to overestimate the EVI, with distinct areas of higher error within the San Pablo Bay region. When using EVI history that included Landsat 7 and Landsat 8 (EVIL7L8) or using the EVI history with Landsat 7, Landsat 8, and Harmonized Sentinel-2 (EVIL7L8Sent), we see minor differences within model error (Figures 6 and 8). The EVI mean difference using EVIL7L8Sent history appears to be closer to zero (i.e., more yellow on the map, Figure 8A), indicating lower model error, with a few regions in the South Bay having higher model error ( Figure 8A). Furthermore, the EVI envelopes have minor changes in Figure 6. Top graphs are the 2019 predicted EVI values for Eden Landing (A) and Rush Ranch (B) using three different EVI history datasets. L8 indicates modeled EVI using Landsat 8 for EVI history, L8S indicates modeled EVI using Landsat 8 and Sentinel-2 for EVI history. L7L8S indicates modeled EVI using Landsat 7, Landsat 8, and Sentinel-2 for EVI history. The bottom graphs are eddy covariance measured gross primary production (GPP) at Eden Landing in 2019 (A) and Rush Ranch in 2016 (B). The green lines represent the predicted EVI values, the black dots and lines indicate the eddy covariance measured GPP, and the blue dots represent the PEPRMT modeled GPP. Eddy covariance data were not available for 2019 within Rush Ranch, therefore we used 2016 data. Root mean square error (RMSE) and modified Willmott's index of agreement (d mod ) for the observed versus predicted GPP are indicated in the figures.
shape when using multiple satellites to calculate EVI history versus using only La ( Figure 6). When summarizing model error within restoration project sites, we further shape when using multiple satellites to calculate EVI history versus using only Landsat 8 ( Figure 6).  When summarizing model error within restoration project sites, we further see the highest model error within locations in the northwestern parts of the regions, i.e., San Pablo Bay ( Figure A2B). For most restoration project sites, EVIs were overpredicted for 2019, while the averaged model error was −0.06 across the sites, the highest model over When summarizing model error within restoration project sites, we further see the highest model error within locations in the northwestern parts of the regions, i.e., San Pablo Bay ( Figure A2B). For most restoration project sites, EVIs were overpredicted for 2019, while the averaged model error was −0.06 across the sites, the highest model over prediction was −0.213, and the highest model underprediction was 0.054. The areas with the highest model over-and underprediction were both located within the San Pablo Bay. A few small areas of underprediction also occurred within the Suisun Bay, while all the areas within the South Bay were overpredictions.
Comparing mean model error within bayland types and restoration projects, we further see that model error differs between types ( Figure A4A Model accuracy can also be influenced by tides and time since the end date of restoration. We analyzed model accuracy within the 14 model development locations and found the water level above mean higher high water observed during the Landsat image acquisition. This allowed us to determine that during the highest water levels, there was also higher calculated model error ( Figure A5). We also examined how time since the end date of a restoration project was related to model error. The model restoration end date and model error were significantly correlated (Spearman correlation coefficient = −0.34), and model error was larger in areas with more recent restoration (Figure 9).
Remote Sens. 2021, 13, x FOR PEER REVIEW prediction was −0.213, and the highest model underprediction was 0.054. The ar the highest model over-and underprediction were both located within the San Pa A few small areas of underprediction also occurred within the Suisun Bay, whi areas within the South Bay were overpredictions.
Comparing mean model error within bayland types and restoration projects ther see that model error differs between types ( Figure A4A Model accuracy can also be influenced by tides and time since the end date ration. We analyzed model accuracy within the 14 model development locati found the water level above mean higher high water observed during the Lands acquisition. This allowed us to determine that during the highest water levels, th also higher calculated model error ( Figure A5). We also examined how time since date of a restoration project was related to model error. The model restoration e and model error were significantly correlated (Spearman correlation coefficient and model error was larger in areas with more recent restoration (Figure 9). . Restoration project end date versus mean difference. Mean difference is the model error calculated for Lan observed minus predicted EVI. The project end dates for 128 restoration projects were obtained from EcoAtlas from t Habitat Projects dataset. A smoothing curve was fitted to the data using loess to estimate the trend within the data.

The Potential of Multi-Temporal Phenological Information to Facilitate Greenness Interpolation in Tidal Systems
Using historical EVI at given time steps in conjunction with elevation, temp and the Palmer Drought Severity Index, we were able to model minimum and m expected plant greenness, important inputs into greenhouse gas flux modeling, t out an entire year. Remote sensing and particularly multi-spectral satellite data a tools in measuring plant greenness on a landscape scale and for establishing tim data capturing phenological dynamics of plant coverage and biomass. Howeve spectral satellites cannot "see" through clouds, and within tidal wetlands, vegeta Figure 9. Restoration project end date versus mean difference. Mean difference is the model error calculated for Landsat observed minus predicted EVI. The project end dates for 128 restoration projects were obtained from EcoAtlas from their Habitat Projects dataset. A smoothing curve was fitted to the data using loess to estimate the trend within the data.

The Potential of Multi-Temporal Phenological Information to Facilitate Greenness Interpolation in Tidal Systems
Using historical EVI at given time steps in conjunction with elevation, temperature, and the Palmer Drought Severity Index, we were able to model minimum and maximum expected plant greenness, important inputs into greenhouse gas flux modeling, throughout an entire year. Remote sensing and particularly multi-spectral satellite data are useful tools in measuring plant greenness on a landscape scale and for establishing time series data capturing phenological dynamics of plant coverage and biomass. However, multispectral satellites cannot "see" through clouds, and within tidal wetlands, vegetation signals are influenced by tidal flooding as a function of elevation, tidal height, and vegetation height [14,25,26]. This makes it difficult to assess both the degree of flooding throughout the entire landscape at a given point in time, and the degree to which remotely sensed plant greenness is representative of the vegetation biomass at that time. Our modeling approach circumvented the need to consider tidal data as well as cloud-related limitations and the lack of repeated 3D assessments of vegetation at the landscape scale, allowing us to predict likely plant greenness. We also were able to obtain a full-season approximation of phenological variability (i.e., phenological envelopes), by estimating plant greenness throughout the entire year.
The historical two-week mean EVI values have the greatest importance in the model, indicated by the largest standardized beta coefficient. The Palmer Drought Severity Index has the second highest standardized beta coefficient and has the second largest influence on the model. This indicates that plant growth is dependent on drought conditions and the positive PDSI beta coefficient indicates that plant greenness is higher under wetter conditions. Other studies have also shown that tidal wetland plant growth decreases under years of drought [31,55]. With decreased rain, plants can experience higher stress through altered hydroperiod, reduced soil water availability and increased salt burden. The day of year beta coefficient has a negative squared polynomial term, meaning that EVI follows a concave parabolic relationship with DOY; the max EVI is in the middle of the year and the min is at the beginning and tail ends of the year. For minimum temperature, there is a negative beta coefficient, indicating that with warmer temperatures, the EVI is higher, reflecting the dynamics of plant biomass, photosynthetic activity, and metabolism following seasonal changes in solar radiation and atmospheric conditions. When looking at elevation above mean higher high water, there is a positive squared polynomial beta coefficient, indicating that the EVI follows a convex parabolic relationship with elevation.

Model Performance
We observed the highest model error within diked wetlands and restoration projects (Figures 8 and 9). Diked wetlands are often highly managed for waterfowl, salt production, or agriculture, and thus their plant communities may not predictably respond to seasonal environmental conditions in their human-regulated state [56]. For example, the majority of the wetland areas within the Suisun Bay are managed by private Duck Clubs that control the tidal influence within the marsh and aim to grow food for waterfowl [57]. This manipulation of the landscape may lead to the larger model error. Within restoration project sites, we also see higher model error, which likely resulted from inter-annual transitions in vegetation coverage, biomass, and composition over the course of post-restoration plant colonization and re-establishment, making historical dynamics of greenness less representative of current phenology [31,32]. We also see high model error in four restoration projects within the San Pablo Bay region ( Figure A2). Two of these projects were completed in 2009 and the other two were completed in 2015 ( Figure A2). For the projects completed in 2015, the EVI history included values from before, during, and after the restoration project. The projects completed in 2009 and 2015 are likely still changing and revegetating. One of the restoration sites, the Napa-Sonoma Ponds 3, 4, 5, was formerly a salt pond with hypersaline conditions ( Figure A3). It is possible that the system is still recovering and still approaching natural functioning observed in other similar wetlands. Furthermore, the elevation data from this site were collected in 2003 and have inevitably changed. Since elevation is a model parameter, the inaccurate elevation for this site may also be a factor in model error. To create a more accurate elevation map for restoration project sites, a more current LiDAR elevation could be used and converted to elevation above mean higher high water [21,27]. Since our study region is large, has different tidal regimes, and added complexity with diking, converting the more recent LiDAR dataset to elevation above mean higher high water across different types of wetlands and coastal structures was beyond the scope of our study. On a smaller scale of individual sites or management units, this could more feasible, and it would be useful for future research to dedicate efforts to create an accurate elevation and tidal inundation map within the Bay Area tidal wetlands.
Model performance slightly improved when using Landsat 7, Landsat 8, and Harmonized Sentinel-2 data to calculate the mean historical EVI dataset. This is evident by the overall lower model error when using these three datasets. We did not include all available Harmonized Sentinel-2 data since the dataset is not currently available on GEE and the process to add all available images is cumbersome, however, in the future the harmonized product may become available on GEE which would streamline the process. Another important thing to note about the Harmonized Sentinel-2 data is that the cloud mask is not optimal. We found that the cloud mask was often overly sensitive and removed sections within wetland areas that did not contain clouds. This is a noted issue within the Harmonized Landsat 8 Sentinel-2; and future work could be done to improve automation of cloud removal for the Harmonized Sentinel product [41].
Lastly, we estimated model performance by comparing the predicted EVI to the satellite observed EVI. As mentioned previously, remotely sensed plant greenness might not accurately represent actual greenness. This is especially concerning in tidal wetlands since flooding may completely or partially submerge vegetation. Since it is difficult to ensure that the satellites acquired imagery only during low tide (this may not always be possible at lower revisit frequencies), tidal height varies both throughout the different images and throughout the landscape. After accounting for tidal height in the 14 model development sites, we observed higher model error during the highest tides. In a portion of the satellite images, the vegetation signal is likely reduced by tidal influence. This indicates that actual model error might be lower than estimated.

Vegetation Greenness throughout the Bay Area
There is also notable regional variation in greenness phenology across the study region, which becomes evident in the modeling outcomes and has important implications for modeling biogeochemical fluxes. During the summer, tidal wetlands within the Suisun Bay in general have higher plant greenness, or EVI, than the South Bay or the San Pablo Bay ( Figure 5). Suisun Bay is fresher than San Pablo Bay and South Bay due to freshwater inputs from the Sacramento and San Joaquin Rivers. Salinity stress can decrease plant growth rates while not impacting the photosynthetic rates, and this would lead to lower plant greenness of the same species and minimal difference in observed CO 2 flux [58]. Salinity also helps dictate the vegetation assemblage. Byrd et al. also found that the Suisun Bay region has higher aboveground biomass within tidal wetlands than San Pablo Bay [24]. Though we did not estimate aboveground biomass, plant greenness and aboveground biomass are often related [4,24,59]. The eastern portion of the San Pablo Bay region has lower plant greenness than the rest of the San Pablo Bay. This area was historically tidal marshes, converted into salt ponds and farms, and subsequently restored into tidal marshes again. Restoration projects are still ongoing within this area, but several projects were completed in the early 2000s. Though many restoration projects are complete within this area, the sites are likely still recovering. Similarly, a large portion of the tidal wetlands within the South Bay were converted into salt ponds. Restoration projects have restored part of the salt ponds back into tidal wetlands, however, industrial salt ponds remain within the area. The high salinity and management practices within the salt ponds make it difficult for plant growth, and typically only a narrow range of algae can grow in their hypersaline water [60].
Although there currently is no single consistent vegetation map product for the whole study region, information from existing mapping efforts and field studies (e.g., [31,[61][62][63][64]) strongly suggests the role of varying composition of plant communities in the observed magnitude of greenness and its phenological patterns. The dominant vegetation type within many tidally controlled bay area wetlands is Sarcocornia pacifica, or pickleweed; and this is especially true within the South Bay. S. pacifica is a succulent and does not completely turn brown like other Bay Area tidal wetland plants such as Spartina folisosa, Pacific cordgrass [65]. In a study in Spain, researchers found that Sarcocornia fruticosa had peak aboveground biomass within the winter months, while belowground biomass peaked in the spring and summer [66]. This may be one reason why the greenness curve for Eden Landing does not appear to decrease at the end of the year. However, another study looked at aboveground biomass of S. pacifica (also known as Salicornia virginica) and separated the biomass between the woody and the green succulent part. They found the succulent fraction was highest in the summer and lowest in the winter [67]. The succulent portion of S. pacifica is the green fraction while the woody part is the older brown portion; and this would indicate that there is higher photosynthesis during the summer months than the winter months and thus a summer growth signal should be observed within satellite data. However, the study also mentions the high spatial and temporal variability in S. pacifica aboveground biomass observed in other studies [67]. Furthermore, the greenness signal within Eden Landing, and other wetlands within the South Bay, is lower than wetlands in San Pablo Bay and Suisun Bay, making it more difficult to decipher a greenness curve. This is further compounded by algal films that are sometimes present within the area, further skewing the signal from the plants [31,65].

Restoration Areas
Given the varying and often uncertain time scales of wetland restoration projects before equilibrium stages [68], it is possible that the restored wetland areas are not yet fully restored, and the extent and density of vegetation are increasing. There are unaltered tidal wetland areas located adjacent to these restored sites, and the unaltered wetlands have higher plant greenness values than the restored sites ( Figure A3). Model error is variable between each of the different restoration project sites, but we do see a trend that older, more established restored sites have less model error than more recently restored sites. Restoration projects ending in roughly 2000 or earlier all have similar model error and projects that ended after 2000 have higher model error ( Figure A2). Wetland areas within restoration projects should be further investigated to determine the time frame it takes for the system to recover, and how wetland type plays a role in recovery time, which remains an important question in restoration efforts globally [31,32,68].

GPP Modeling
Wetland plant greenness is an important input into carbon sequestration models such as PEPRMT. Pairing our method of modeling the EVI with PEPRMT will allow for landscape-scale analysis of carbon sequestration within tidal wetlands over a given time frame. We have demonstrated that using our modeled EVI within the PEPRMT model provides an accurate estimate of gross primary production. This approach allows researchers to upscale ecosystem function across a large landscape. Future work could enhance these predictions further by incorporating in situ phenological observations from cost-effective phenocams and integrating them with satellite data [5,6], which would require expanding the scope of phenocam monitoring in this region. Engaging phenocam data to improve broader-scale predictions of phenology as satellite-based greenness envelopes would be especially beneficial for wetland sites showing high uncertainty in EVI modeling, such as actively managed or recently restored wetlands.

Limitations and Future Research Directions
Several limitations encountered in this study provide useful insights on the future research directions for advancing the use of remotely sensed phenological information in tidal marsh ecosystem analyses. One notable challenge concerned the variable availability of cloud-free Landsat time series across the study area. Although most of our focal region fits within one Landsat scene (WRS-2 path 44, row 34, Figure A1), the number of "usable" Landsat observation dates highly varies due to regular influxes of fog from the Pacific Ocean and its uneven distribution across the San Francisco Bay region. The South Bay experiences more cloud cover than the Suisun Bay and San Pablo Bay, which reduced the number of usable images to quantify historical variability in the EVI for wetlands in this part of the region. In addition, two Landsat tiles covered most of the Suisun Bay, three tiles covered the western part of the San Pablo Bay, and one Landsat tile covered the South Bay ( Figure A1). The low image availability, especially within the South Bay, makes it particularly difficult to estimate tidal wetland plant growth cycles without using modeling approaches such as the one outlined in this paper. When data are limited, including the Sentinel Landsat harmonized product can be beneficial. Currently, the harmonized product is not available on GEE, however, data can be manually uploaded and incorporated into the EVI history calculation. Since this product is not currently available on GEE, we did not use this approach for the main model, however, the ability to find and process large amounts of Landsat data within GEE makes it ideal for time series analysis. In the future, Sentinel Landsat harmonized data may become available within the GEE database, which will likely improve our modeling approach.
Another notable challenge concerns lower applicability of the historical EVI to currentyear seasonality in sites that have undergone a recent perturbation, such as restoration. Restoration projects can drastically alter the tidal regime, plant assemblage, and plant density. If the historical EVI dataset includes data from pre-restoration, the EVI signal can greatly differ from that after completion of the restoration project. Furthermore, the landscape continues to change even after completion of the restoration and there is not a universal timeline for when the ecosystem reaches an equilibrium state. For areas within known restoration projects, there is little to be done to improve our modeling approach, although this issue underscores an important and widely recognized research need to better understand the time frame of project dynamics, recoveries, and equilibria [68,69]. However, a shorter EVI history can be used which would only include data after completion of restoration. For historical data, emerging time series analysis tools such as Landsat-based Detection of Trends in Disturbance and Recovery (LandTrendr) can also help to identify characteristic time frames of post-restoration perturbation and recovery [70,71].
Lastly, using our modeling approach, we did not consider less common rapid unpredictable environmental changes. California is fire prone, so it is possible that fires will rapidly alter vegetation within wetland areas. For example, within Southern California, a fire resulted in decreased plant greenness within S. pacifica dominated wetland; vegetation eventually re-established but on the order of years [72]. It is more obvious which wetland areas experience fires, however, another natural phenomenon that is less spatially predictable and would influence vegetation coverage is salt marsh dieback. Salt marsh dieback is most noted in Spartina alterniflora dominated wetlands, has been recorded in California marshes, and causes regions of vegetation to thin or completely die [73]. The region eventually revegetates but the cycle takes years. Factors such as salinity, elevation, and drought may play a role in the vulnerability of a marsh to experiencing salt marsh dieback, however, the phenomenon is unpredictable [74]. Salt marsh dieback and burned areas would cause a lower observed plant greenness signal than predicted within our model.
Satellite data can be further validated or gap-filled through ground-based information such as through phenocams. Phenocams are low-cost cameras set up in the field that continuously take multi-spectral images. The cameras are set up horizontally across the landscape, which differs from satellite images which capture data from above. A relationship between the phenocam data and the satellite data can be established, allowing for gaps to be filled within the satellite time series. Phenocams also do not cover as wide of an area as satellites, which makes it more difficult to monitor vegetation across a large landscape while using phenocam data alone [6,54,75]. However, more widely collected phenocam data could improve our modeling approach in the future or provide for better validation of the growth curve envelopes that we model with satellite-based data.

Conclusions
Using EVI history and site characteristics such as elevation, the Palmer Drought Severity Index, and temperature, we can estimate the two-week Enhanced Vegetation Index, a popular spectral indicator of plant greenness, across the landscape for a given year. This allows us to estimate an EVI times series without having to consider tidal flooding, vegetation type, or unavailable satellite data due to poor atmospheric conditions. It also provides a practical strategy to predict proxies of vegetation phenological changes for modeling of plant-driven ecosystem processes in complex tidally influenced estuarine environments in the absence of repeated LiDAR observations. At the same time, our predictive modeling approach is sensitive to factors impacting the representativeness of phenological history of current ecosystem dynamics, particularly time since restoration, where more recently restored wetland sites tend to have larger model error. However, even among such projects, the degree of error varied, indicating that both site and time since restoration are important parameters to consider. Within restoration projects sites, the type of restoration, plant assemblage, and location within the Bay Area may all be factors in restoration success and system establishment. Overall, our findings indicate that this modeling approach can streamline estimating a vegetation time series across a large spatial extent. The predicted EVI output can then be incorporated into additional models such as PEPRMT which estimates carbon sequestration within wetlands. Future research could develop strategies to alleviate some of the key detected uncertainties. This includes more explicitly assessing restoration and management effects on wetland phenological dynamics within the regional context. This would allow for incorporating the status of these human interventions into phenological modeling.    Appendix E Figure A5. Model error (mean difference) versus water level (m) above MHHW. Mean difference is the predicted EVI minus 2019 observed Landsat EVI. We estimated the water level by finding the water level at the tide gauge closest to the study sites used for model development. We only calculated the mean difference versus water level for the 14 study areas used for model development. We found the mean difference versus water level for all the available 2019 Landsat images. We fit a smoothing trendline with standard error to better understand how mean difference is influenced by water level.  Appendix E Figure A5. Model error (mean difference) versus water level (m) above MHHW. Mean difference is the predicted EVI minus 2019 observed Landsat EVI. We estimated the water level by finding the water level at the tide gauge closest to the study sites used for model development. We only calculated the mean difference versus water level for the 14 study areas used for model development. We found the mean difference versus water level for all the available 2019 Landsat images. We fit a smoothing trendline with standard error to better understand how mean difference is influenced by water level.
(A) (B) (C) Figure A5. Model error (mean difference) versus water level (m) above MHHW. Mean difference is the predicted EVI minus 2019 observed Landsat EVI. We estimated the water level by finding the water level at the tide gauge closest to the study sites used for model development. We only calculated the mean difference versus water level for the 14 study areas used for model development.
We found the mean difference versus water level for all the available 2019 Landsat images. We fit a smoothing trendline with standard error to better understand how mean difference is influenced by water level.