A Comparison between the MODIS Product (MOD17A2) and a Tide-Robust Empirical GPP Model Evaluated in a Georgia Wetland

Abstract

Despite the importance of tidal ecosystems in the global carbon budget, the relationships between environmental drivers and carbon dynamics in these wetlands remain poorly understood. This limited understanding results from the challenges associated with in situ flux studies and their correlation with satellite imagery which can be affected by periodic tidal flooding. Carbon dioxide eddy covariance (EC) towers are installed in only a few wetlands worldwide, and the longest eddy-covariance record from Georgia (GA) wetlands contains only two continuous years of observations. The goals of the present study were to evaluate the performance of existing MODIS Gross Primary Production (GPP) products (MOD17A2) against EC derived GPP and develop a tide-robust Normalized Difference Moisture Index (NDMI) based model to predict GPP within a Spartina alterniflora salt marsh on Sapelo Island, GA. These EC tower-based observations represent a basis to associate CO₂ fluxes with canopy reflectance and thus provide the means to use satellite-based reflectance data for broader scale investigations. We demonstrate that Light Use Efficiency (LUE)-based MOD17A2 does not accurately reflect tidal wetland GPP compared to a simple empirical vegetation index-based model where tidal influence was accounted for. The NDMI-based GPP model was capable of predicting changes in wetland CO₂ fluxes and explained 46% of the variation in flux-estimated GPP within the training data, and a root mean square error of 6.96 g C m⁻² in the validation data. Our investigation is the first to create a MODIS-based wetland GPP estimation procedure that demonstrates the importance of filtering tidal observations from satellite surface reflectance data.

1. Introduction

Estimates of Gross Primary Productivity (GPP) are useful for understanding local and global carbon dynamics and are often scaled from sites to landscapes using satellite data [1,2,3,4]. A common data source for this effort is the Moderate Resolution Imaging Spectroradiometer (MODIS), which has a coarse spatial resolution but collects daily global observations. Information from MODIS has been packaged into the MOD17A2 GPP product, which is a cumulative composite of daily GPP estimates. The MOD17 algorithm uses the radiation-use efficiency concept first proposed by Monteith [3] and provides a method for estimating pixel-wise GPP and Net Primary Productivity (NPP) globally [5]. Satellite-derived GPP, such as from MODIS, is useful because it helps model terrestrial energy, carbon, and water cycle processes. However, the algorithm for MOD17 was calibrated for upland biomes and cannot be directly applied to wetlands with any accuracy. Consequently, there is a knowledge gap concerning the spatial and temporal variability of wetland GPP.

Tidal marshes, the wetland type that is the focus of this study, are ecologically productive ecosystems [6,7] and are at risk to loss from threats such as land-use change and sea level rise [6]. Estimations of tidal marsh GPP can improve our understanding of feedback mechanisms that regulate carbon dynamics in the context of global change and facilitate climate policy-making that might preserve coastal habitats.

Remote sensing techniques provide an efficient method for estimating GPP and MODIS is a useful data source for this kind of monitoring. While MODIS is spatially coarse (pixel spatial resolution of 250 m to 1 km), its daily return frequency is especially useful for estimating plant and ecosystem properties within tidal marshes [8,9]. Daily observations allow sufficient data density after cloudy and tidally influenced observations are removed. Further, MODIS products often are used to estimate GPP in terrestrial systems [10,11] and parameterizing a similar product for carbon-rich productive tidal marshes allows for landscape-scale comparisons among terrestrial and coastal systems through MODIS data.

MODIS can only assist with monitoring spatiotemporal variability in wetlands if the influence of tides on vegetation reflectance can be isolated. Water or saturated soil background in wetlands reduce the intensity of the near-infrared (NIR) reflectance, change the slope and position of red-edge, and alter the magnitude of vegetation indices such as Enhance and Normalized Vegetation Index (EVI and NDVI) and other red-edge-type indices [12,13,14]. Consequently, GPP models used for terrestrial ecosystems cannot be readily applied to wetlands without accounting for the effect of tides on surface reflectance [2,4,15].

The MODIS GPP product, MOD17A2, while useful for upland monitoring, cannot be used for tidal wetlands because it has not been parameterized for this habitat. The MOD17A2 product employs a Production Efficiency Model (PEM) to estimate GPP. PEMs are based on the original logic of Monteith [3], who suggested that GPP is linearly related to absorbed photosynthetically active radiation (APAR), such that

GPP ∝ ε × Σ (fAPAR × PAR_in)

(1)

where, PAR_in is the incident photosynthetically active radiation, fAPAR is the fraction of PAR_in absorbed by the canopy, and ε is light use efficiency (LUE). The MOD17A2 product estimates GPP by assuming a linear relationship between fAPAR and NDVI and substitutes LUE as a vegetation biome-specific constant obtained from the lookup tables [16]. However, tidal marshes are not one of the eleven biomes implemented within MOD17A2 [17], perhaps because, compared to terrestrial systems, the processes that regulate ecosystem-atmosphere carbon dioxide (CO₂) fluxes, including GPP and LUE in wetlands are not well understood [18]. Furthermore, models which do not require LUE maybe easier to implement across habitats and are worth exploring [16].

Vegetation indices (VIs) are frequently used in GPP models. Sims et al. [19] demonstrated that a simple model based on the Enhanced Vegetation Index (EVI) alone, could provide estimates of GPP. Their methods yielded results that were as good as or better than MOD17A2 for many sites. Sims et al. [16] developed another model, the temperature and Greenness (TG) model, which was based on MODIS EVI and the Land Surface Temperature (LST) products. Further, Wu et al. [20] proposed a new vegetation index GPP model that demonstrated the usefulness of VIs as reliable proxies of both LUE and fAPAR [19,21,22,23].

The vegetation index approach for estimating GPP might also be useful within coastal marshes, provided that we can overcome the distortion of conventional vegetation indices by periodic tidal flooding. Such an approach needs to begin with daily surface reflectance information, because composited MODIS products used for calculating MOD17A2, such as fAPAR and leaf area index (LAI), are derived from multi-day composites of MODIS NDVI or EVI products. As such, observations that make up multi-day composites may contain noise from observations collected during high tide and cannot be removed post hoc [24]. Thus, we need a procedure for filtering daily MODIS surface reflectance for tidal observations, calculating a vegetation index related to GPP, and then creating tide-free composites from that filtered time series. This procedure should be a step forward for estimating GPP from MODIS within tidal wetlands.

Others have studied productivity within coastal wetlands through remote sensing data [1,14,25,26]. However, these studies have either not accounted for tidal influence [14,25,26], focused on biomass or net ecosystem exchange (NEE) rather than GPP [14,26], or evaluated canopies that are not submerged by tides, (e.g., mangroves; Barr et al., 2013 [1]). Indeed Yan et al. (2010) [26] observed a loss of phenological pattern in pixels that were most tidally influenced, which might be caused by the use of composite MODIS products without tidal filtering [24]. Thus, additional work to solve the problem of tidal influence on the spectral reflectance of coastal marshes, and particularly to apply this effort to GPP estimation, is still needed.

In this paper, we present a method to calibrate MODIS GPP in tidal wetlands. The objective of this study was to build a tide-robust model to predict GPP within a Spartina alterniflora salt marsh on Sapelo Island, GA. The specific objectives of this research were to (1) evaluate and compare the potential of using different MODIS vegetation indices to characterize flux tower GPP, (2) to establish the statistical relationship between vegetation indices from MODIS 500 m surface reflectance products and GPP flux, and (3) calibrate a MODIS GPP algorithm within a tidal wetland.

2. Materials and Methods

2.1. Study Site

Our study area was a salt marsh on Sapelo Island, Georgia, USA dominated by a monoculture of S. alterniflora at lat 31.444°, long −81.283° (Figure 1). This marsh is monitored as part of the Georgia Coastal Ecosystems Long Term Ecological Research (GCE-LTER) program. The tides here are semi-diurnal and have large amplitudes (>2 m). At this marsh site, there is also an in situ eddy-covariance (EC) tower which measures CO₂ fluxes between the marsh surface and the atmosphere. The EC footprint varies, but typically covers an area larger than 500 m². We selected a focal MODIS pixel (500 m × 500 m) for our GPP model parameterization which is the one that is most often within the footprint of the EC tower and contains only a salt marsh. This pixel should best represent the information measured by the EC tower, whereas other marsh pixels nearby should be less representative.

2.2. Data Sources

The data for this study included MODIS version 5 products, including surface reflectance data and MODIS GPP products. The MODIS tile covering Sapelo Island is h11v05. Additionally, we used the EC tower data of CO₂ flux to calculate GPP data in conjunction with tidal data obtained from the Georgia Coastal Ecosystems Long-Term Ecological Research (GCE-LTER) project. For each dataset, we acquired information beginning 1 January 2014 and ending 31 December, 2015. All data processing and analysis steps were conducted through the use of the R language for statistical analyses.

2.2.1. MOD17A2

MOD17A2 is an 8 day composite GPP product at a 1 km spatial resolution. It provides a measure of the growth of the terrestrial vegetation using daily MODIS land cover, fAPAR, LAI, and surface meteorology for the global vegetated land surface [17]. MOD17A2 was acquired from the National Aeronautics and Space Administration (NASA) website (http://modis–land.gsfc.nasa.gov) and covered 2014–2015, representing 90 composited data points. While this pixel was at a coarser scale than the focal 500 m surface reflectance pixel we used to build our salt marsh GPP model, the 1 km scale was the only resolution available for the MOD17 product at the time of this study and still samples S. alterniflora marsh close to our EC tower. In addition, the marshes surrounding the tower is coved my multiple 500 m pixels and depending on the wind speed, the flux footprint varies in size and direction, and may very well contain fluxes from adjacent areas beyond the pixel.

2.2.2. MOD09GA

The Terra MODIS surface reflectance atmospheric correction algorithm product (MOD09GA, 500 m) is daily data which includes bands 1, 2, 3, 4, 5, 6, and 7 centered at 0.648 μm, 0.858 μm, 0.470 μm, 0.555 μm, 1.24 μm, 1.64 μm, and 2.13 μm respectively. MOD09GA also provides important pixel quality and viewing geometry information for MODIS products, including an accurate cloud mask that can be used to flag cloudy pixels. MOD09GA was acquired from the National Aeronautics and Space Administration (NASA) website (http://modis–land.gsfc.nasa.gov) and covered 2014–2015.

2.2.3. Flux GPP Data

GPP was calculated using an eddy covariance flux system installed in the study site and maintained by the GCE-LTER. Two eddy-covariance systems were installed on a flux tower established at the head of a small tidal creek on the Duplin River (lat: 31.441°, long: −81.283°) in July 2013, Figure 1. The systems were mounted at 5 m above the ground facing South (180°) and North (0°), respectively. Two EC systems were utilized to minimize the gaps in the data due to maintenance and calibration, instrument malfunction, and changing wind direction. The south system covers the angle from 90° to 270° and the rest of the area is covered by the north system. Each EC system consists of a 3D sonic anemometer (CSAT3, Campbell Scientific Inc., Logan, UT, USA) to measure the three-wind components and virtual air temperature, and a closed-path CO₂/H₂O gas analyzer (Li-7200, Li-Cor Biosciences, USA) to measure the concentration of CO₂ in the air. These measurements were recorded at 10 Hz with the control of data logger (CR3000, Campbell Scientific Inc., USA). Flux data collected by these sensors were posted to the website of the GCE-LTER.

Flux Data Processing

Raw 10 Hz data from the infrared gas analyzer and sonic anemometer were processed for 30-minute runs using the EddyPro 6.1 software (Li-Cor Inc., Lincoln, NE, USA). The data was checked for data spikes and they were removed, following the methods of Vickers and Mahrt [27]. A coordinate rotation using the planar fit method [28] was applied to the sonic anemometer data to remove tilt errors as well as linear trends. The Webb-Pearman-Leuning (WPL) correction for density effects due to heat and water vapor transfer [29] was applied to correct the calculated CO₂ fluxes.

Quality control was applied to each 30-min eddy-covariance data block and data from north and south systems were combined. The processed NEE data for each sensor was filtered and combined into one data set based on the prevailing wind direction and the sensor’s quality control parameters. Data from the north facing tower was used for wind directions between 315 and 45 degrees, the south facing tower datasets were used for wind directions between 135 and 225 degrees. For the remaining wind directions, the average of both sensors was used.

The footprint model was then applied to the map with estimated fetch values and angle to screen out measurements that did not originate from the study area. Finally, only measurements which had footprints where at least 70% of the CO₂ flux originated from within the marsh area surrounding the flux tower were used in this study, which mostly occurred during daytime conditions. This combined NEE was further filtered based on marsh inundation values. The filtered data was then used to calculate the GPP.

The friction velocity threshold (u*) for nighttime fluxes at our study site was 0.10 ms⁻¹ and determined per Papale et al. [30]. We used the median of the absolute deviation method by Papale et al. [30] to detect outliers within each 30 min block of flux data. A total of 37% of the observations were removed, with a greater percentage of removed observations from nighttime conditions.

Flux Partitioning

Daytime NEE is usually modeled through nonlinear light response curves (i.e., Falge et al. [31]). Nighttime fluxes which represent respiration are modeled by temperature dependent models such as that proposed by Lloyd and Taylor [32]. Before applying these models, we first partitioned the flux data into the daytime and nighttime NEE based on light level. These partitions equate to ecosystem respiration (R_eco) and GPP per Reichstein et al. [33]. The NEE for the system is derived as follows:

$$EE=\{\begin{array}{c}{R}_{eco} Fornighttimeperiods\\ R_{eco}-GPPFordaytimeperiods\end{array}$$

(2)

The temperature dependence of R_eco is described by Lloyd and Taylor’s (1994) [32] model as follows:

$$R_{eco}=R_{ref}.e^{E_0\left(\frac{1}{T_{ref}-T_0}- \frac{1}{T_{air}-T_0}\right)}$$

(3)

where, R_eco (µmol m⁻² s⁻¹) is the ecosystem respiration, R_ref is the rate of respiration at a reference temperature, T_ref set at 10 °C, E₀ (K) is the activation energy of R_eco when T₀, the base temperature is set to −46.02 °C [32], T_air is the air temperature.

2.2.4. Marsh Inundation Data

Associated with the flux tower were additional instruments, including the Campbell Scientific CS455 pressure transducer sensor, deployed in a creek adjacent to the study marsh. The pressure transducer measures the mean water level relative to NVD88 datum every 5 min. The measurements were calculated from water pressure in mH₂O and sensor elevation of −0.23 m from RTK GPS measurements. An automated weather station with 5-min average data output was also installed to measure air temperature, relative humidity, solar radiation, precipitation, and wind speed and direction.

Additionally, there was a PhenoCam mounted on the flux tower, a digital camera which auto-collects oblique-angle images of the marsh every half-hour during daylight. We used a smart classifier to extract the percent cover of water over the marsh from each PhenoCam image [34]. Then, to estimate percent cover of marsh flooding for each 5 min period of the day, including nighttime, we used a generalized linear model with binomial error to predict percent cover flooding from creek water height, wind speed, and wind direction. This provided a better estimate of marsh flooding than creek water height alone [24].

2.2.5. MODIS Data Preprocessing

The cloud state information from the “QC_500 m: 500 m Reflectance Band Quality” layer within the MOD09GA product was used to filter cloudy pixels from the spectral reflectance dataset. Tidal inundation information was generated from our new Tidal Marsh Inundation Index (TMII) [24]. This index creates a flag that can identify tidal flooding within vegetated marshes by accounting for mixed pixels that reflect both vegetation and inundation. TMII also accounts for seasonal vegetation development and maximizes the separation of vegetation and tides across the annual cycle.

2.3. Accounting for Tides in Flux and MODIS Data

2.3.1. Creating a Tide-Free GPP Time Series from Flux Data

The data preprocessing steps for the NEE collected at the flux tower included the following steps: generating tide free NEE data, calculating tide free GPP, gap-filling, and finally, smoothing GPP over five-day windows (Figure 2). The data flow produced NEE and respiration values (R_e) measured during tidal periods and a separate dataset where tidal influences were below a determined threshold value based on the marsh inundation data. Both data sets allowed for the calculation of GPP using tidally separated NEE and nighttime R_e. We used the method developed by Reichstein et al. [33] to gap fill GPP data, where the gap filling was based on radiation and temperature inputs gathered on site.

Half hour values of tide-free and gap-filled flux tower GPP (μmol C m⁻² s⁻¹) were summed using Equation (4) to provide daily GPP in units of μmol C m⁻² per day. The coefficient of 0.0216 in Equation (4) converts units of μmol CO₂ m⁻² s⁻¹ to g C m⁻² per each 30 min flux averaging interval [1].

$$GPP=0.0216{\displaystyle \sum }_{i=1}^{48}GP{P}_{30\min }$$

(4)

As a final step, we summed the daily values to create 8 day sums of GPP to match the temporal scale in the MOD17A2 GPP product, allowing a comparison of these two datasets. We assumed that all the carbon flux from the atmosphere into the marsh ecosystem was positive.

2.3.2. Creating Cloud-Free and Tide-Free MODIS VI Composites

The steps for creating high-quality MODIS data included filtering cloud and tide-influenced observations, as well as filtering observations with view zenith angles >50° (Figure 3).

The cloud mask from the MOD09GA QC layer was used to filter MODIS data products to produce cloud-free images. Similarly, we used the TMII to identify and filter tidal inundation from the MODIS dataset, following O’Connell et al. [20]. We indexed observations with TMII > 0.2 as flooded and filtered these from the dataset.

The next step was to calculate VIs and composite those in order to match the temporal scale of the MOD17A2 product, as well as the 8 day sums created from the flux GPP. To create these 8 day VI composites, we followed the workflow we previously developed [24], which is based on the guidance outlined in Huete et al. [35]. To summarize this procedure, we composited Vis by first ordering VI observations by view zenith angle within each 8 day window. Up to five observations, in order of lowest view zenith, were used to create the average VI value for the window, provided that these observations all had view zenith angles <35°. If low view zenith observations were not available, the function selected the single VI value corresponding to the lowest view zenith angle <50°. If no low view zenith values were available <50° within the window, then the window the VI value was set as missing data. The entire processing tree is depicted in Figure 4.

2.4. Select the Best Vegetation Index as a Proxy for GPP

Despite the usefulness of VIs in other habitats as proxies of both LUE, fAPAR [19,21,22,23], and GPP [16,36,37], these methods need to be adapted for tidal marshes. To create the best MODIS GPP model, it was important to test a variety of vegetation indexes, including the Normalized Difference Vegetation Index (NDVI) [38], the Soil Adjusted Vegetation Index SAVI [39], the Normalized Difference Moisture Index (NDMI) [40], the Visible Atmospherically Resistant Index (VARI) [41], the Enhanced Vegetation Index (EVI) [35], and the Wide Dynamic Range Vegetation Index (WDRVI) [42]. Equations for these indices are

$$NDVI=\frac{{\rho }_{NIR}-{\rho }_R}{{\rho }_{NIR}+{\rho }_R}$$

(5)

where $\rho$ is spectral reflectance at the wavelength indicated by the subscript, NIR is near-infrared as indicated previously and R is the red wavelength.

$$SAVI=\frac{{\rho }_{NIR}-{\rho }_R}{{\rho }_{NIR}+{\rho }_R+L}\left(1+L\right)$$

(6)

where L is a constant and L = 0.5 is the default value and works well in most situations.

$$DMI=\frac{{\rho }_{0.86\mathsf{\mu }\mathrm{m}}-{\rho }_{1.24\mathsf{\mu }\mathrm{m}}}{{\rho }_{0.86\mathsf{\mu }\mathrm{m}}+{\rho }_{1.24\mathsf{\mu }\mathrm{m}}}$$

(7)

where the MODIS bands involved are 2 (near infrared at 0.858 µm) and 5 (1.24 µm).

$$VARI=\frac{{\rho }_G-{\rho }_R}{{\rho }_G+{\rho }_R}$$

(8)

$$EVI=G\frac{{\rho }_{NIR}-{\rho }_R}{{\rho }_{NIR}+C_1\times {\rho }_R-C_2\times {\rho }_B+L}$$

(9)

where L = 1, C₁ = 6, C₂ = 7.5, and G (gain factor) = 2.5.

$$WDRVI=\frac{\alpha {\rho }_{NIR}-{\rho }_R}{\alpha {\rho }_{NIR}+{\rho }_R}$$

(10)

where α is between 0 and 1 and is a weighting parameter for NIR reflectance.

To select the best VI for estimating flux GPP, we divided the data into training (2015) and validation (2014) data in order to build the model and then estimate its performance on the most novel data available for our target pixel, e.g., that from the new year. For training data, we developed simple linear regression models where each candidate VI (NDVI, SAVI, NDMI, VARI, EVI, WDRVI) was used, in turn, to predict flux GPP data. We also estimated Root Mean Square Error (RMSE) of the model prediction as

$$RMSE=\sqrt{\frac{{\displaystyle \sum }{\left(\widehat{y_i}-y_i\right)}^2}{n}}$$

(11)

where $\widehat{y_i}$ was the predicted GPP for y_i, the observed GPP value from the flux 8 day sum, and n was the sample size. For validation, we used the best model from the training regressions as indicated by R² to predict GPP. We then compared that to the “true” 8 day flux GPP sums to get RMSE for the de novo data without refitting the model. This RMSE for validation data provides an estimate of how well our training model works to predict new observations at our study site. Finally, for the validation GPP prediction, we compared this to the RMSE of a simple linear regression where the MODIS product MOD17A2 was the predictor of the flux GPP.

3. Results

3.1. Creating a Tide-Free GPP Time Series from Flux Data

Flux GPP measurements during some periods were missing due to instrument malfunctions. Windows with missing data were excluded from the final 8 day sums, resulting in 304 (2015 training) and 359 (2014 validation) daily GPP observations. These daily GPP sums ranged from −0.08 to 6.22 g C m⁻² across both years. After processing to 8 day sums, the GPP ranged from −0.16 to 39.76 g C m⁻², while the maximum 8 day GPP was similar among both years, the minimum 8 day GPP was much higher in 2014 than 2015 (4.35 vs. −0.16 g C m⁻²).

3.2. Creating Cloud-Free and Tide-Free MODIS VI Composites

From our focal pixel, we acquired daily surface reflectance. Some dates had missing data resulting in 358–360 surface reflectance observations annually. Of these, 44–61 were estimated as tide and cloud-free observations with view zenith <50°. These were composited into regularly spaced 8 day estimates (45 obs each year), but some composite windows contained no valid data, resulting in 20–27 observations with data (

Table 1. The MODIS samples used for GPP modeling and validation.

	Modeling (2015)	Validation (2014)
Daily surface reflectance	360	358
Cloud-free, low view zenith	52	67
Cloud and tide free	44	61
8 day composite	20	27

). The compositing process resulted in higher quality data averages from low view zenith observations which could be used for GPP modeling.

3.3. Selection of the Best Vegetation Index as a Proxy for GPP from the Training Data

Composited and smoothed vegetation indices varied in terms of their range and signal to noise ratios (Figure 5). Of these, NDMI had the least noise across the phenological cycle. When related to GPP, NDMI also a higher R² than other potential VI proxies and, thus, was selected as the best proxy for modeling variation in GPP (

Table 2. The goodness of fit for the relationship between the Vegetation Indices and GPP for the Training Data.

Index	R2
NDVI	0.04
SAVI	0.38
NDMI	0.46
VARI	0.05
EVI	0.38
WDRVI	0.03

).

Therefore, the final equation for predicting GPP was based on an 8 day composited cloud and tide-free NDMI (Equation (12)).

$$GPP=207.5\times NDMI+32.3$$

(12)

3.4. Modeling GPP from the Validation Data

When applied to the 2014 validation dataset, NDMI had a stronger relationship with GPP than MOD17A2, such that RMSE decreased from 7.46 g C m⁻² (training) to 6.96 g C m⁻². This 6.96 g C m⁻² validation RMSE represented 20% of the observed range in GPP 8 day sums (e.g., 5.29 g m⁻²–39.75 g m⁻²) (Figure 6).

During 2014, the GPP modeled from the NDMI relationship had a well-resolved phenological pattern and was more similar to flux GPP than the MOD17A2 product (Figure 7). The MOD17A2 GPP product also lacked an obvious summer maximum. Our new GPP model exhibited a more reasonable phenology than MOD17A2 in the observed year, especially during the growing season (from May to October). MOD17A2 also overestimated GPP in the non-growing season. This demonstrates the ability of the proposed model to improve the data quality of MODIS GPP products for wetlands.

4. Discussion

Combining satellite data and EC tower data is a common approach to estimate ecosystem-scale or regional-scale GPP. Remote sensing data are valuable when used with EC data because they allow us to model and map GPP beyond the tower footprint. Even though the LUE-based GPP model is considered a robust approach, simpler empirical approaches such as the correlation between GPP and vegetation indices are often valuable [16], in part because they do not require many parameters and often are geographically scalable [23]. O’Connell et al. [24] demonstrated the problems (i.e., dampened phenology) with standard vegetation indices derived from MODIS surface reflectance data for wetlands and proposed a tide flagging method as a solution. Based on that method, we developed several cloud-free and tide-filtered vegetation composites to be correlated with tower derived GPP in a simple empirical model. Our best model was able to explain 46% of the variation in tidal wetland GPP in the validation data and had a lower RMSE than MOD17. Therefore, this represents an improvement over the existing method for tidal wetland GPP calculation.

Accounting for tidal effects in wetland remote sensing data is important. We found that 20–40% of the MODIS daily surface reflectance data were affected by either clouds or tides or both in GA wetlands (Table 1). Remote sensing derived vegetation indices in wetlands have several other problems, namely reduced NIR reflectance due to saturated soil background and lack of a pronounced red-edge due to scattering at the red-absorption band [14]. These problems often introduce significant errors in remote sensing based biophysical or GPP models [1,14,17]. Moreover, the errors will increase sharply if we do not filter out tide and cloud affected pixels [20]. The procedure presented in this paper overcomes these problems by demonstrating how to tidally-filter daily surface reflectance information from MODIS, and subsequently generate a higher quality MODIS GPP model within tidal wetlands than the existing MODIS product.

The best vegetation index for estimating GPP in our study was the Normalized Difference Moisture Index (NDMI). Compared to other indices, NDMI served as a reasonable spatial proxy for wetland GPP but only after accounting for tidal conditions. We observed that NDMI, a NIR, and shortwave-IR (SWIR) based index performed better than red-NIR based indices (Table 2). As noted in O’Connell et al. (2017) [24], even though the tide affected pixels are removed, the red-NIR contrast is lower compared to NIR-SWIR contrast. This is because of the perpetually moist soil background in wetlands. Longer wavelength bands are attenuated by water and soil moisture more rapidly than visible wavelengths. Therefore, indices composed of both visible and longer wavelengths are very noisy within wetlands because the soil moisture and water depth vary over time. Further, as noted above, NDMI is composed of bands related LAI (e.g., the NIR band) as well as bands related to plant moisture content [43]. The amount of leaf area in the canopy and the water status of these leaves is positively related to plant photosynthesis rates in S. alterniflora [44]. Perhaps consequently, NDMI has a strong phenological signal and was high during summer when plant biomass was high and low during the winter when plants were dormant (Figure 5). Therefore, NDMI showed a stronger correlation with GPP compared to other indices, likely through its numerical robustness to tidal flooding and its physical relationship to plant physiology.

We applied the NDMI based model to data from Sapelo Island, GA and demonstrated that the GPP modeled from tide and cloud-free NDMI composites results in a more reliable product than the coarser scaled existing MODIS MOD17A2. We can derive two important conclusions from this observation. First, the MOD17A2 GPP product available at the time of this study is not suitable for wetlands. This may be in part because of its coarser scale though we think this contributed only to minor differences because the entire marsh area was an S. alterniflora monoculture. However, more importantly, the MOD17A2 GPP product uses an LUE-based approach by incorporating ground-based meteorological measurements, a grassland biome-properties look-up table (BPLUT), and fPAR and LAI (from MOD 15) [17,45,46,47,48,49,50]. These variables need to be parametrized and recalibrated for coastal marshes. The main improvement our modeling approach provided was to show that a tidally filtered vegetation index can be more accurate than a complex modeling approach that is not parametrized for local conditions or species. Second, although significant, NDMI only explained 46% of the variability in observed GPP which is much lower than the LUE-based PEM models developed for mangroves wetlands (e.g., Barr et al., 2013) [1]. These mangrove canopies stick above high tide and therefore do not require spectral tidal filtering that salt marshes require. However, unlike LUE, vegetation indices are often used as proxies for greenness and not carbon uptake.

In this study, we have shown that a simple vegetation index based GPP model can be more accurate than a poorly calibrated complex LUE based model (MOD17A2). Ge et al. [25] also observed that LUE-based models without tidal filtering can give a poor performance in tidal salt marshes. Therefore, we conclude that existing MODIS-based PEM or VPM models for wetlands can be improved significantly by filtering and re-compositing the vegetation index used in model parametrization. While a well-parameterized LUE-based model may ultimately give better results, this study is a step forward because this simple VI MODIS models of GPP explains a large fraction of GPP variation within the pixel best representing the flux tower footprint and is easy to calculate. The high temporal resolution of MODIS coupled with the moderate spatial resolution of a 500 m algorithm can result in an immense potential for studying and monitoring both long and short-term tidal wetland health and physiological statuses, and for allowing GPP comparisons with the upland areas that MOD17A2 was designed to estimate.

5. Conclusions

The present study represents a first attempt to design and calibrate simple MODIS GPP models for salt marsh ecosystems through the integration of remotely sensed information, flux tower observations, and tidal data. The study describes how to composite MODIS daily surface reflectance data or derived products within tidal marshes, allowing us to model tidal wetland productivity. We certified that LUE-based MOD17A2 does not accurately reflect tidal wetland GPP compared to a simple empirical vegetation index-based model where tidal influence was accounted for. The results of this study also point to the need for more robust calculations of the tidal effects on ecosystem gross primary productivity.

Mapping such productivity data gives insight into ecosystem function and can establish sites as sinks or sources of carbon in the environment. The MODIS-based models developed in this study can serve as a baseline for developing satellite-based models of GPP of the tidal wetlands. This is particularly important because these productive coastal ecosystems are vulnerable to climate change induced sea level rise and perpetual developmental pressures.

One of the limitations of this study is that we used a single MODIS pixel to calibrate and validate our model. Although a single pixel does not represent the Georgia wetlands, any study aimed at developing a coast-wide GPP product needs to start at a pixel-level. The initial pixels used to parameterize such a model need to be those in the footprints of EC towers so that the tower data can serve as ground-truth data for parameterizing models and estimating model error. We used the best such pixel in this study. The biggest factor that affects MODIS surface reflectance and, hence, the GPP product is tide. Unless we have a robust tide flagging mechanism to filter the tide affected pixels, the supposedly robust LUE models are always going to be noisy, no matter how much re-parametrization is done. Again, the tide-flagging model needs to start at a pixel level where we have tower-based and phenocam observations to serve as ground-truth data. Fortunately, MODIS produces a dense temporal dataset and even if the tower pixel is affected by cloud and tide 20–40% of the time, we will still be left with a good size dataset to calibrate and validate the GPP model. With that in mind, our goal in this study was to examine the effect of tide on the performance of a simple GPP model before investing time on a complex but more robust LUE based GPP model. Our goal was not to develop a coast-wide tide-filtered GPP map with a simple empirical model which may have other problems (site- or species-specific etc.). For that, we need more pixel-based marsh inundation observations across multiple sites. Only then can we start re-parametrizing LUE based model for the coast-wide application. This will allow us to produce a coast-wide GPP map for salt marshes and revisit the vast MODIS archive to generate an 18+ year time-series composite. It also would be valuable to extend these time series to investigate environmental drivers and patterns of change in marsh GPP and the phenology of GPP at this site and across the broader GA coast.