Monitoring Landscape Dynamics in Central U.S. Grasslands with Harmonized Landsat-8 and Sentinel-2 Time Series Data

Abstract

Remotely monitoring changes in central U.S. grasslands is challenging because these landscapes tend to respond quickly to disturbances and changes in weather. Such dynamic responses influence nutrient cycling, greenhouse gas contributions, habitat availability for wildlife, and other ecosystem processes and services. Traditionally, coarse-resolution satellite data acquired at daily intervals have been used for monitoring. Recently, the harmonized Landsat-8 and Sentinel-2 (HLS) data increased the temporal frequency of the data. Here we investigated if the increased data frequency provided adequate observations to characterize highly dynamic grassland processes. We evaluated HLS data available for 2016 to (1) determine if data from Sentinel-2 contributed to an improvement in characterizing landscape processes over Landsat-8 data alone, and (2) quantify how observation frequency impacted results. Specifically, we investigated into estimating annual vegetation phenology, detecting burn scars from fire, and modeling within-season wetland hydroperiod and growth of aquatic vegetation. We observed increased sensitivity to the start of the growing season (SOST) with the HLS data. Our estimates of the grassland SOST compared well with ground estimates collected at a phenological camera site. We used the Continuous Change Detection and Classification (CCDC) algorithm to assess if the HLS data improved our detection of burn scars following grassland fires and found that detection was considerably influenced by the seasonal timing of the fires. The grassland burned in early spring recovered too quickly to be detected as change events by CCDC; instead, the spectral characteristics following these fires were incorporated as part of the ongoing time-series models. In contrast, the spectral effects from late-season fires were detected both by Landsat-8 data and HLS data. For wetland-rich areas, we used a modified version of the CCDC algorithm to track within-season dynamics of water and aquatic vegetation. The addition of Sentinel-2 data provided the potential to build full time series models to better distinguish different wetland types, suggesting that the temporal density of data was sufficient for within-season characterization of wetland dynamics. Although the different data frequency, in both the spatial and temporal dimensions, could cause inconsistent model estimation or sensitivity sometimes; overall, the temporal frequency of the HLS data improved our ability to track within-season grassland dynamics and improved results for areas prone to cloud contamination. The results suggest a greater frequency of observations, such as from harmonizing data across all comparable Landsat and Sentinel sensors, is still needed. For our study areas, at least a 3-day revisit interval during the early growing season (weeks 14–17) is required to provide a >50% probability of obtaining weekly clear observations.

1. Introduction

The U.S. Great Plains provide habitats of global significance to migratory volant species [1,2,3]; important sources of food and fuel to cattle, wild ungulate species such as deer, and people; and serve a role in storing organic carbon [1,4,5]. Monitoring land surface seasonal dynamics and changes in the Great Plains is challenging because grasslands and wetlands are temporally dynamic and sensitive to disturbances and changes in weather [1,4,5,6,7,8]. Despite these challenges in monitoring, information on land surface seasonality, disturbance, and resilience to change across this prairie landscape is needed by decision makers who administer environmental resources and land management programs.

There is a long history of applying time series of satellite data to remotely monitor vegetation characteristics, particularly using data from the Advanced Very High Resolution Radiometer (AVHRR) and Moderate Resolution Imaging Spectroradiometer (MODIS) [9,10,11,12,13]. These sensors, aboard daily polar orbiting Earth observing satellites, provide the needed temporal resolution for observing seasonal dynamics. However, the data are spatially coarse, and the pixels often represent a mixture of perennial grasslands, prairie wetlands, and agricultural lands with annual crops, making it difficult to distinguish a variety of change agents on the landscape including the influence of seasonal weather, disturbance, or agricultural management [14]. To date, satellite sensors that collect observations with higher spatial resolution have not provided sufficient revisit frequency or global coverage to track within-season dynamics [15].

New advances in approaches for conducting time series analysis with the long-running Landsat archive, coupled with techniques to fuse Landsat data with other sensors that have similar characteristics, offer the potential to observe within-season landscape dynamics at 30-m spatial resolution. The dataset fusion of Landsat-8 Operational Land Imager (OLI) and Sentinel-2 MultiSpectral Instrument (MSI) [16] provides an opportunity to investigate seasonal landscape dynamics and disturbances that may be obscured at coarser spatial resolutions. However, a challenge with employing such fused data is that the respective acquisition schedules of the different sensors results in temporally (and spatially) unevenly distributed observations that may impact the detection of change and be difficult to implement with traditional algorithms that were developed for data with consistent time steps.

Traditional algorithms such as logistic regression [17,18,19], polynomial regression [20,21], and harmonic analysis [15,22] are widely used to model seasonal dynamics. The coefficients of such mathematical models may be sensitive to the timing and frequency of observations in highly dynamic land cover types, but few studies have investigated how observation frequency influences model output from such algorithms. Given the timing and irregular frequency of Landsat-8 and Sentinel-2 (HLS) time series data, we assessed the performance of previously established monitoring methodologies by using HLS time series data. The selected study sites represented grassland processes that typically pose challenges when monitoring with coarser-resolution data. Site selection was based on the availability of HLS data, phenological camera data, a known history of grassland fires, and an intensively studied wetland area with in-situ data. We asked: (1) Does the observation frequency obtained with the addition of Sentinel-2 data with Landsat 8 data enhance our ability to detect the start of the vegetation growing season, occurrence of grassland fires, and wetland dynamics in Great Plains grassland landscapes? (2) How does uneven temporal and spatial coverage of observations influence detections of these processes?

To address the questions, two previously established methodologies were applied to the HLS time series data and the Landsat-8 time series data. One method used weighted least-squares regression to develop smoothed time series models from observations and applied a delayed moving average (DMA) to model vegetation phenology [9,13]. The other method used harmonic analysis to detect land cover changes [22]. The DMA model has been employed operationally by the U.S. Geological Survey’s Earth Resources Observation and Science (EROS) Center to produce a land surface phenology database for the conterminous United States based on a consistent weekly time series of composited MODIS normalized difference vegetation index (NDVI) data [9,21]. We used that algorithm with HLS data to model vegetation phenology for an entire growing season (2016) and to estimate the timing of the start of the season (SOST) at a 30-m scale. We applied a harmonic analysis called Continuous Change Detection and Classification (CCDC) [22] to test the detection of fire disturbance and wetland dynamics. We evaluated the results by comparing the differences in using Landsat-8 time series data alone versus HLS time series data. We then performed comparisons with standard products and/or ground-based observations.

2. Study Area

2.1. Vegetation Phenology

We estimated the start of the 2016 growing season over one HLS tile located along the eastern border of North Dakota and western Minnesota, representing 12,000 km² (Figure 1). The monthly average temperature for this study site ranges from −12 °C to 20 °C, and monthly precipitation ranged from 0.5 cm to 16 cm [23] in 2016. Precipitation in the form of snow is common during winter months (November-April), and spring green-up generally occurs in April when the average low temperature rises above 0°C [23]. The study area is dominated by cropland and pasture (78%) and grassland covers about 2% according to the 2016 Cropland Data Layer (CDL) [24]. The site contains cool and warm season grasses and, importantly, includes the location of the Oakville phenological camera (“PhenoCam”) [25]. The Oakville PhenoCam provided near-ground collected phenological data for a native tallgrass prairie site.

2.2. Grassland Fire

Eastern Kansas (Figure 1) contains some of the largest areas of unplowed native tallgrass prairie in North America [26]. The selected study site is a subset of one HLS tile and covers approximately 8.08 km² (361 × 746 pixels of 30-m spatial resolution). The monthly average temperature ranges between −3°C and 26°C and average annual precipitation was about 95 cm [23] from 2013 to 2016. The growing season is from April to November [9]. This region also undergoes frequent prescribed burning (Figure 1), and most of the fires are initiated during March and April for the purpose of managing grassland species [27]. In 2016, multiple early growing season fires (15 March, 8 April, and 24 April) and two late growing season fires (17 October and 25 October) were indicated by the 30-m product Burned Area Essential Climate Variable (BAECV) data [28].

2.3. Wetland Dynamics

The Cottonwood Lake Study Area (CLSA) in east-central North Dakota (Figure 1) contains various classes of wetlands, ranging from very ephemeral to nearly permanent, that have been intensively researched [29]. The CLSA covers 1.2 km², or 200 × 200 pixels (30-m pixel resolution) area within one HLS tile. The average annual precipitation for this location is 50 cm, and the monthly average temperatures ranged between −16°C and 21°C [23] during the years 2013 to 2016. Ephemeral wetlands at the site range in size from less than a single pixel (30 m × 30 m) to several pixels. Permanent wetlands can cover dozens of pixels or more, especially during wet years. Water depths of the most permanent wetlands range from 1 m to 3 m, whereas depths in most ephemeral wetlands are less than 0.5 meters when they are inundated [30].

3. Source Data

3.1. Harmonized Landsat-8 and Sentinel-2 (HLS)

We acquired HLS surface reflectance and data quality bands for the North Dakota–Minnesota grassland phenology site (HLS tile 14UPU), the Kansas grassland fire site (HLS tile 14SQH), and the North Dakota wetland site (HLS tile 14TMT). The HLS tiles are based on the Sentinel-2 tiling system that measures 109.8 km by 109.8 km in the UTM/WGS84 projection [16].

Table 1. The temporal range, number of images, and extracted spectral bands (based on Landsat-8) for each study area.

	Temporal Range		Number of Images		Input Bands (µm)
Study Topic	Landsat-8	Sentinel-2	Landsat-8	Sentinel-2
Phenology	7 January 2015–25 April 2017	4 September 2015–25 July 2017	157	60	Red 0.64–0.67 Near-infrared (NIR) 0.85–0.88
Grassland fire	16 April 2013–27 April 2017	6 August 2015–30 August 2017	185	70	Blue 0.45–0.51 Green 0.53–0.59 Red 0.64–0.67 NIR 0.85–0.88 Shortwave Infrared-1 1.57–1.65 Shortwave Infrared-2 2.11–2.29
Wetland seasonal dynamics	12 April 2013–23 April 2017	4 October 2015–15 June 2017	182	55	same as the grassland fire study

shows acquired imagery and spectral bands for each site. The only year with complete temporal coverage of both HLS data for our study sites was 2016. The HLS data were preprocessed for atmospheric correction and cloud masking, geometric resampling, geographic registration, bidirectional reflectance distribution function (BRDF) normalization, and band-pass adjustments to harmonize the images from the two sensors [16]. A known issue with the acquired version (1.3) of HLS products was that the cloud masking performed poorly under various surface and atmospheric states such as bright targets, thin clouds, cloud edges, and hazy conditions [16]. The preprocessing for HLS data differs from standard Landsat-8 and Sentinel-2 products; thus, Standard Landsat-8, Sentinel-2 and HLS products are not directly comparable. From here on, we will use the terms “Landsat-8” and “Sentinel-2” to refer to their respective subsets of the HLS data.

The spatial and temporal distribution of observations varied for Landsat-8 and HLS time series data (Figure 2). Because of different acquisitioning schedules and path side-lap versus non-side-lap areas, the frequency of Landsat-8 and Sentinel-2 observations during 2016 varied from 16 to 80 at different pixel locations within a tile. In addition, observations were not distributed at equal time steps. The Landsat-8 high frequency region (yellow swaths in Figure 2a) was influenced by the observation frequency from two adjacent orbital paths. Incorporating observations from Sentinel-2 enriches the data frequency (Figure 2b) over Landsat-8 alone, but the different revisit intervals between Landsat-8 and Sentinel-2 result in a more complex spatiotemporal distribution of data.

3.2. Reference Data for the Vegetation Phenology Analyses

3.2.1. Oakville PhenoCam

We utilized ground-collected daily imagery from a PhenoCam to support the site-specific comparison between ground-observed and the HLS and USGS eMODIS phenology metrics [31] (remote sensing phenology data are available from the USGS EROS Center at doi.org//10.5066/F7PC30G1). The Department of Geography at the University of North Dakota (UND) operates a PhenoCam on the Oakville Prairie Biology Station at 47.898° N latitude and 97.316° W longitude. Oakville is a 960-acre native tallgrass prairie site managed by the UND Department of Biology. In recent years, a controlled burn regime has been introduced to reduce detritus, brush, and invasive plant species.

The PhenoCam is a StarDot NetCam SC 5.0 IR with 1296 × 960 pixel (5 megapixel) resolution. The footprint of PhenoCam is about 11 Landsat-8 pixels that form a sector shape within 5 by 3 pixels and covers one MODIS 250 m pixel. It is sensitive to both color and near-infrared (NIR) radiation, allowing for the derivation of many spectral vegetation indices. Four-channel multispectral image data (blue-, green-, and red-visible and NIR wavebands) were collected and sent through a cellular modem every 30 min between 04:00 and 23:30 Central Standard Time. NDVI was calculated for the time series imagery for 2016, the only year for which the PhenoCam data were available. The time-series algorithm for determining SOST, however, required at least 2 years of data for calculating the metric. To address this issue, the 2016 data were duplicated and used as a surrogate for the missing 2015 data, as the current year’s SOST is minimally impacted by the previous year’s trend [13].

3.2.2. eMODIS NDVI and Phenology Metrics

We used expedited MODIS (eMODIS), 250-m weekly composited NDVI as an additional source of information to generate phenology data across the entire HLS tile [20,32]. The eMODIS weekly composited data were generated from daily MODIS data. The NDVI time series product was based on a selection of high-quality observations (i.e., highest NDVI) and was temporally filtered and smoothed over the growing season to reduce the influence of cloud contamination and noise. The full description of the eMODIS processing can be found at [20]. The start of the growing season (SOST) calculated from time series eMODIS NDVI at 250 m resolution identifies the estimated day of year when photosynthetic activity in the canopy (at the scale of a pixel) increases above a winter baseline. We implemented a common processing flow such that SOST based on three different source NDVI time series data sets (eMODIS NDVI, HLS NDVI, and Oakville PhenoCam NDVI) could be compared to illuminate the influence of data frequency on the calculation. This allowed for comparison over a greater number of grassland locations.

3.2.3. CDL and NLCD

To compare the phenology results from HLS data to eMODIS for grasslands, we employed a spatial mask at 250-m resolution. The Cropland Data Layer (CDL) and 2011 National Land Cover Database (NLCD) are both 30-m products that have thematic classes from which we identified dominant land cover types in eMODIS pixels. We acquired the CDL for 2016 from the U.S. Department of Agriculture’s National Agricultural Statistics Service [24] to identify grassland-dominant pixels. We acquired the 2011 NLCD [33] to mask water pixels after the SOST estimation.

3.3. Reference Data for the Grassland Fire Analyses

3.3.1. Monitoring Trends in Burn Severity (MTBS)

The MTBS product mapped fires greater than 1000 acres in size (western United States) using multi-temporal image differences and image ratios with Landsat data [27]. Currently, MTBS products map fires from 1984 to 2015 and contain (1) polygon fire perimeters, (2) a Normalized Burn Ratio (NBR) map, (3) a differenced NBR (dNBR) map, and (4) a thematic classification of burn severity [27]. We acquired the fire perimeter product from 2013 to 2015 to identify the historical fire burns in the study area.

3.3.2. Burned Area Essential Climate Variable (BAECV)

BAECV is a Landsat scene-based product developed from machine-learning and MTBS products [28]. For each pixel, the annual product contained the (1) maximum burn probability, (2) number of Landsat scenes classified as burned within the year, (3) Julian date of the first fire burn, and (4) number of scenes with unobstructed observations [28]. The products for 2016 reported omission and commission errors for fire detection at 42% and 33%, respectively, but 62% and 57%, respectively, for grassland fires in the Great Plains of the central United States. [34]. We acquired the Julian date of the first fire burn product for 2016 to compare with the fire detection results from our study.

3.4. Reference Data for the Wetland Dynamics Analyses

3.4.1. Climate Data

We compared climate data with the modeled wetland seasonal dynamics. Monthly precipitation data for the county where the study site is located were acquired from National Oceanic and Atmospheric Administration (NOAA) climate data online servers [23] for comparison with the modeled wetland seasonal dynamics.

3.4.2. Water Surface Elevation Data

We used water surface elevation data to support the within season differences between wetland types (https://www.sciencebase.gov/catalog/item/579bc8ede4b0589fa1c982c8) [29]. The water surface elevation data were measured by field gauges (field sites in Figure 1) installed in each wetland. The gauges were read once per week during the ice-free portion of each year from 1979 to 2016. For easier interpretation, we converted the surface elevation data to water depth following the procedure from LaBaugh et al. [29].

4. Methods

To investigate whether the HLS data enhanced the models and captured short-term events and seasonal dynamics for grasslands, we conducted the three analyses described below with HLS and Landsat-8 datasets. The results were then cross-compared as well as compared with the reference data.

4.1. Phenology

We incorporated a multi-year remote sensing phenology product at 250-m resolution based on eMODIS NDVI data that currently cover the period 2001 to 2017 [9]. In this study, we applied the same methodology, including the smoothing algorithm [21] and phenology analytic techniques [13], to the HLS and Landsat-8 data to estimate the 2016 SOST at 30-m resolution. The methodology requires NDVI input data at consistent time intervals, so the HLS and Landsat-8 datasets were temporally normalized using a compositing strategy described in the steps below. For the comparison, we first compared the spatial pattern of Landsat-8 SOST, HLS SOST, and eMODIS SOST at the tile scale to evaluate the consistency of the algorithm. Then we focused on the grassland pixels and compared the variation of Landsat-8 SOST and HLS SOST within the corresponding eMODIS pixels. At last, we compared the PhenoCam SOST with single pixels from Landsat-8 SOST, HLS SOST, and eMODIS SOST. The steps below outline the data processing used to perform the SOST analysis.

Step 1. NDVI was calculated using the equation (NIR – Red-visible)/(NIR + Red-visible). The Landsat NIR band represents 0.85 to 0.88 μm. Sentinel-2 also provides an NIR band covering this same range. The red-visible band represents 0.64 to 0.67 μm for both Landsat-8 and Sentinel-2. The compositing strategy selected maximum NDVI values from either Landsat-8 or Sentinel-2 for each pixel for each week across the entire time period.

Step 2. Anomalous pixels, where the red band equaled −1 or less, were designated as “no data.”

Step 3. Data were stacked into a multiband time series covering 2015, 2016, and 2017 and processed with a smoothing algorithm [9] to remove missing NDVI observations resulting from cloud contamination and other noise. The result is filtered time series trajectories that better characterize the underlying seasonal characteristics of the vegetation at the land surface.

Step 4. The smoothed NDVI time series data were then processed using the DMA algorithm [13] to identify the SOST for each pixel.

Step 5. Resulting SOST dates were rescaled to −150 to 365 to represent all 365 days in the year of analysis (2016) and the prior year back to August 4 (equivalent to −150 or day of year (DOY) 216) to accommodate situations where the SOST was identified in the previous year. Values outside that range were set to −1000 (“no data”).

Step 6. The two sets of SOST results (calculated from time series HLS and Landsat-8) were compared with the eMODIS and PhenoCam start of season results.

Step 7. Grassland cells (total 126) for investigation were identified using the CDL (at a 250-m cell size to correspond with the eMODIS cell size). Cells that were at least 80% occupied by class 176 (Grass/Pasture) were selected for further analysis. For comparison between the eMODIS SOST, the Landsat-8 SOST, and the HLS SOST, the moderate (30-m) resolution SOST results were summarized by calculating the mean of the grassland (CDL class 176) pixels within the 250-m cells. We also calculated standard deviation to understand the variability in the moderate-resolution results. The PhenoCam SOST was also compared with Landsat-8 SOST and HLS SOST using the 30-m pixels that corresponded to the PhenoCam footprint.

4.2. Detecting Burn Scars from Grassland Fires

We used CCDC to detect changes in the spectral characteristics of vegetation following grassland fires. CCDC is a land cover change detection method based on Landsat time series [22]. The method models the seasonality and linear trend of Landsat spectral time series and detects changes depending on the model’s root mean square error. We conducted the CCDC analysis using observations that were not flagged as cloud, cloud shadow, or snow/ice in the HLS quality band. We ran CCDC first with Landsat-8 data only, then with both Landsat-8 and Sentinel-2 data. The change dates from CCDC results were then compared with each other as well as with the BAECV product.

4.3. Characterizing Seasonal Wetland Dynamics

We used tasseled cap (TC) indices [35] generated from the reflectance data to help us interpret how wetland seasonal dynamics models might be improved with more frequent data. To estimate the TC seasonal dynamics, we first derived the seasonal dynamics model for each spectral band through harmonic analysis and the CCDC method. CCDC detects different seasonal patterns from year to year and builds independent models, if needed. Moreover, CCDC automatically adjusts its highest harmonic frequency from annual to intra-annual bimodal or intra-annual trimodal change based on the number of available observations (see Figure A1) [22]. Given the uneven distribution of observations, the frequency adjustment information allowed us to investigate the model performance in a dynamic environment. Figure A1 exhibits how the CCDC model performs when different frequencies are used, given the same set of observations as input.

CCDC assumes the multi-year time series trend is linear, whereas depending on water inundation and vegetation occurrence, wetlands can have substantial variability at both multi-year and within-season scales. To capture potential nonlinear multi-year trends, we first extracted the multi-year spectral trends from the time series reflectance data using harmonic analysis (Equation (1)) and delineated the within-year spectral variation from the residuals using CCDC. The sum of the harmonic model and CCDC model represented the seasonal wetland dynamics from which the daily spectral time series were calculated and then converted to TC indices time series. The detailed procedure for each pixel was as follows:

Step 1. We discarded all observations that were flagged as cloud, cloud shadow, or snow/ice in the HLS quality band to reduce noise in the analysis.

Step 2. The multi-year dynamic model was built for each spectral band using Equation (1):

$$F\left(i, x\right)={\displaystyle \sum }_{n=2}^4\left(a_{n, i}cos\frac{2\pi x}{nT}+b_{n,i}sin\frac{2\pi x}{nT}\right)+{\epsilon }_{i,x}$$

(1)

where F is the fitted model for band i and time x (Julian date), T is the number of days per year (T = 365.25), a_n,i and b_n,i are coefficients for the n year cycle at band i, and ε_i,x is the model residual for the observation at time x and band i.

Step 3. The within-year dynamic model was built for each spectral band by applying CCDC to the residuals of the multi-year dynamic model.

Step 4. The multi-year and within-year models were added for each spectral band to represent the seasonal dynamic model.

Step 5. The seasonal dynamic models were converted to daily spectral time series and transformed to TC indices time series [35].

5. Results

5.1. Phenology

The frequency of observations in the input time series data influenced phenology detection. Figure 3 illustrates the geographic patterns of the SOST results for tile 14UPU. eMODIS showed that the majority of SOST for the study site occurred during April (DOY 92–121; Figure 3c). Although results from the Landsat-8 and HLS data also agreed with this finding, their spatial patterns illustrate how the number and distribution through time of clear observations influenced the calculation of the SOST (Figure 3g,h). In low data frequency locations, more than 80% of pixels had one to two clear Landsat-8 observations from March to May, compared with three to five clear observations in HLS. In the high data frequency zone, Landsat-8 had two to six clear observations for 88% of pixels, whereas HLS data had three to eight clear observations. Generally, the algorithm estimated an earlier SOST in low data frequency zones than in high data frequency zones. In addition, less spatial variation was apparent in the low data frequency zone than in the same area in the eMODIS output. The results derived from Landsat-8 time series data showed more pixels with an SOST after DOY 160 in the high data frequency zones (shown in dark green in Figure 3a). The SOST results from both the HLS and eMODIS time series did not show this pattern (Figure 3b,c). Moreover, the subpanels (Figure 3e,f) suggest that spatial patterns at the local scale are better captured by the 30-m data in the high data frequency zone while the eMODIS is more generalized.

The scatter plot in Figure 4 shows the relationship between eMODIS grassland SOST vs mean SOST from corresponding cells in HLS data. The HLS SOST and eMODIS SOST are similar but do not have a strong linear relationship in most of the grassland pixels and some scattered outliers in the lower right corner come from low data frequency regions. The root mean square difference was 11.6, indicating a moderate uncertainty for this relationship. The linear relationship suggested that HLS SOST are generally earlier than eMODIS SOST. It is also worth noting that both eMODIS SOST and HLS SOST range about 40 days, as the results of this study only represents 2016 within a limited study area.

The box plots in Figure 5 show the distribution of SOST values calculated for grassland 250-m cells within the study site. The average and median values for all three SOST results compared well and were all within two weeks of one another. The distribution of SOST from Landsat-8 showed a greater range from the lower to upper quartile than either the SOST derived from HLS or from eMODIS. The Landsat-8 SOST outliers (shown in black circles) were much earlier than the HLS and eMODIS outliers and occurred in December and January, which we assume to be false. Compared with the Landsat-8 SOST results, the grassland SOST results derived from the HLS time series showed a much tighter range, one being comparable with, though earlier than, the SOST range derived from the eMODIS data (Figure 5). We also calculated the standard deviation within the Landsat-8 and HLS grassland cells to illustrate the differences in spatial variability. Figure 6 shows that the standard deviation was much lower for the HLS data (in orange), showing greater spatial consistency of the phenology calculation at the 30-m resolution. In this case, we assume that the lower standard deviation in SOST means there are fewer anomalous SOST pixels in the HLS results compared with the Landsat-8 results.

We utilized ground-collected daily imagery from the PhenoCam to support a site-specific comparison, similar to that performed by Vrieling et al. [37]. The PhenoCam imagery was geographically limited, as there was only one grassland site within this study area, however, it was useful for illustrating the challenges encountered in estimating SOST using time series observations. Figure 3e shows that the PhenoCam location is within a high-frequency zone of the HLS tile. Figure 7 shows the smoothed NDVI from Landsat-8, HLS, eMODIS, and the camera at the PhenoCam site, as well as the derived SOST. The Landsat-8 NDVI time series exhibited relatively more variability during the growing season. These fluctuations in NDVI were apparently caused by limited high quality data frequency in the Landsat-8 time series even though the site is in a location where data frequency for Landsat-8 was relatively high. The HLS time series NDVI trajectory corresponded well with the PhenoCam NDVI, especially in the first half of the year. The eMODIS NDVI showed a lag of about one week in the early growing season compared with the HLS NDVI time series, a result that has been observed in a previous study [38].

5.2. Grassland Fires

The CCDC model based on both Landsat-8 and HLS time series data was able to detect burned areas when fires occurred late in the growing season but did not detect most early season burned areas (Figure 8a,b). Results from Landsat-8 and HLS data suggested that the more northerly of the two late-season fires occurred on December 1, in contrast with the 25 October date estimated by the BAECV (Figure 8c). The results for the late-season fire at the south end of the study area showed heterogeneous timing in the detection of change that ranged from 30 August to 13 December for the Landsat-8 CCDC model and from 24 September to 13 December for the HLS CCDC model. The BAECV identified 17 October as the time of change for that fire. Most areas burned in 2016 were also historically burned based on the MTBS data. The Landsat-8 CCDC also detected square-shaped changes in the northern part of the study area that were symptomatic of a cloud detection issue with the HLS quality band (Figure 8a). The HLS CCDC appeared more robust to this issue and did not show these spurious patterns (Figure 8b). Figure 9 showed a pixel trajectory of an early season fire in each spectral band. Observations following the dashed line suggested the fire event detected by BAECV caused some spectral differences that were distinct from the previous observations although the spectral differences were not lasting.

5.3. Wetland Dynamics

The Landsat-8 only time series data were insufficient to establish models in some low data frequency areas and in highly dynamic areas such as land-water boundaries typically associated with seasonal changes in inundation and vegetation (Figure 10). The eastern portion of the study area, especially the southeast corner, had the lowest data frequency in Landsat-8, although the same area had improved frequency with the HLS data. Overall, the HLS data had enough observations to support CCDC modeling with a full range of frequency in most of the study area. HLS data also generated slightly more annual models than did the Landsat-8 data in the high data frequency areas (Figure 10).

To illustrate how seasonal variation differs for wetland types, time series for three example pixels are shown: a permanent lake (referred to as P4); an open basin, flow-through wetland (referred to as T3); and a recharge wetland (referred to as T8). According to water surface elevation data, P4 contained water throughout the entire year and had water depths ranging from 1.67 m to 1.89 m [30]. T3 was inundated most of the year, but water depths were usually less than 0.42 m. T8 had water for only two weeks, during late April and early May of 2016, and had the shallowest water depths (≤0.13 m) [30]. In Figure 11, orthophotography of the three wetland types indicated spectral response differences between the wetland types. The TC seasonal variation from models based on HLS or Landsat-8 also reflected these differences. With Landsat-8 data, CCDC built models for annual and intra-annual bimodal and trimodal (full model) frequencies for P4 and T8, but no models were generated for T3. Using the HLS data, CCDC built full models in P4, T3 and T8. Models built from HLS data also appear to reflect more of the expected seasonal dynamics in T3 and T8 when compared with seasonal precipitation [23,39] than the models based on Landsat-8 only data (Figure 11). The wetland P4 contained water throughout the season, and the TC indices indicated less seasonal variation as the TC-wetness values remained higher than those in the more dynamic systems that had less inundation. The three observations in Figure 11 with remarkably high TC-brightness were due to thin cloud contamination that was missed by the HLS cloud mask. Both T3 and T8 show pronounced seasonal variation in TC indices, and T8 had slightly higher TC-greenness than T3, reflecting less inundation and more visible growth of aquatic vegetation (Figure 11).

6. Discussion

In this study, we investigated two methods for monitoring important environmental processes in Great Plains grasslands and whether improvements were gained by increased observation frequency in HLS time series data in highly dynamic grassland systems. Although this study was limited to a single year (i.e., 2016), our results showed model improvements with the addition of Sentinel-2 data to Landsat-8 time series. However, temporal frequency and distribution of observations resulting from different operational schedules for the two sensors and from cloud cover still posed challenges for our applications, especially for sites located near the centers of orbital paths where there is no side-lap. The specifics and magnitudes of these challenges varied by application and location.

HLS data frequency usually enhanced our ability to model the seasonal variation of grasslands and wetlands, for which Landsat-8 time series frequency alone was insufficient. SOST calculated from Landsat-8 time series had much greater variability and more extreme outliers than SOST calculated from the HLS time-series source. SOST derived from HLS compared more closely with eMODIS-derived SOST and had fewer outliers. Closer comparison with the NDVI time series and SOST at the PhenoCam site also demonstrated that HLS data, especially in areas having a higher frequency of observations due to orbital side-lap, could better estimate grassland phenology than Landsat-8 time series data alone. Although the PhenoCam data showed lower magnitudes of NDVI than satellite-measured NDVI, this was expected and can be explained by differences in the viewing angles of the sensors [41]. The PhenoCam’s oblique viewing angle is more likely to see grass litter and emerging grass blades than the nadir view from satellites [41]; but PhenoCam-derived NDVI still can be strongly correlated with satellite-derived NDVI in consistently indicating phenological changes over homogeneous grasslands [41].

The fire detection map derived from Landsat-8 was influenced by an artificial effect related to problems with HLS cloud masking [42]. The Landsat-8 cloud detection algorithm (HLS version 1.3) sometimes confused bright targets with clouds, then dilated the target in the Landsat-8 cloud mask, resulting in anomalous square-shaped patterns of change in the Landsat-8 data. This artificial effect, however, was not seen in the HLS results, which indicated that the CCDC model was more robust to data gaps when provided with denser time series data. Although not yet available for our study, the cloud detection issue will be improved in the following HLS version [16].

The within-season wetland models in Figure 10 and Figure 11 showed that although Landsat-8 observations revealed seasonal dynamics, the low data frequency was only enough to estimate coefficients of the annual harmonic model. The HLS TC time series captured more detailed seasonal variation in wetlands and provided the potential to build full time series models to better distinguish different wetland types [43]. HLS data, however, resulted in more simple models in the high data frequency area (Figure 10) than when using only Landsat-8 data. This was because the relatively high frequency of Landsat-8 observations was enough to initialize a CCDC model before Sentinel-2 observations were available at the end of 2015. When the Sentinel-2 observations provided data showing more seasonal variation in 2016, the variation exceeded the model-predicted root mean square error; therefore, the algorithm assumed the land surface had changed, initializing a new model using the remaining observations. These remaining observations were so few that only a simple model was fit. A similar model-fitting pattern was also observed by Roy et al. [36], who found that the linear harmonic model sometimes resulted in lower root mean square differences (RMSD) between the model and observations in low data frequency areas. This spatial RMSD inconsistency also exhibited the temporal model sensitivity issue. When linear harmonic models were built using low data frequency, the models were more likely to detect change when the data frequency and observation variation increased. This issue raises caution for harmonic analysis in model estimation or change detection related to inconsistent model sensitivity in both the spatial and temporal dimensions.

We observed that HLS data observation frequency was still limiting in some cases due to the relatively sparse observations in the early part of the growing season, which may have been linked to early season cloudiness. The CCDC model failed to find early growing season fires in both Landsat-8 and HLS time series. This is due to the CCDC model parameter that requires six consecutive anomalous observations within the time series before a change is defined. Based on visual examination of the imagery in the Kansas study site, fire scars from early growing season fires were usually visible for less than one month in this region. In the HLS time series, we typically found fewer than six observations in one month, even in locations with higher temporal frequency due to side-lap. To investigate if CCDC sensitivity could be adjusted to this shorter temporal window, we altered the parameter setting to require that only four consecutive observations be anomalous before declaring a change. This increased the detection of pixels exhibiting spectral change associated with early season fires, but also introduced considerable noise (see Figure A2). The CCDC method aims to find changes that have different seasonal variation from previous years. Thus, if historical fires occurred regularly at a similar time of the year, the algorithm would determine the fire pattern in spring to be a normal part of seasonal variation. BAECV or MTBS, on the other hand, were specifically designed to detect fire scars by comparing spectral data before and after a fire on a scene-by-scene basis. Therefore, it is not surprising that these fire-focused products provided a more complete record of fire scars than resulted from CCDC’s approach to detect a wide variety of types of change.

The SOST maps and wetland models resulting from the SOST data were less affected by the insufficient frequency of data than those resulting from the Landsat-8 data, although the impact was still visible. In particular, the results indicated the necessity of increasing observation frequency during the early growing season, when critical phenological and Earth-surface changes are underway.

A common finding in all three grassland applications was that the combination of frequency and timing of observations was critical to successfully identify key processes at a 30-m resolution. We considered the plausibility of meeting this data need by using information from the quality bands from all available 500-m MODIS Surface Reflectance Daily data (MOD09GA) [44] to estimate the probability of cloud cover in eastern North Dakota and eastern Kansas. Figure 12 shows the probability associated with different revisit intervals for acquiring at least one clear observation per week, which is needed to generate weekly composites with eMODIS data. From this, we estimated that a 3-day revisit interval would provide adequate time series data to study grassland dynamics, particularly during the early part of the growing season (weeks 14–17 in our study areas). This data frequency would lead to a greater than 50% probability of providing a weekly clear observation. Cloudier areas would require a shorter revisit period to maintain or increase the probability of a clear observation. In the future, with the combination of Landsat-8, Sentinel-2A, and Sentinel-2B, and the planned launch of Landsat-9, the global median average revisit interval could reach less than 3 days [15,45]. This speaks to the importance of maintaining the continuity of future Landsat missions and interoperability between Landsat and Sentinel-2 missions to support local to regional monitoring of landscape dynamics for informing decisions on land resources.

Our study provides an initial exploration of applying moderate resolution data harmonized from two satellite sensors to characterize highly dynamic environmental processes in grassland landscapes in the U.S. Great Plains. The results show that enhanced time series combining Landsat observations with data from similar sensors such as Sentinel-2 shows promise for enabling tracking of within-season environmental processes at higher spatial resolution than previously possible. The data available for our study was limited to only one full year of Sentinel-2 observations (2016). This restricted our analysis to a single growing season and reduced our ability to more fully understand the effects of data frequency on estimating seasonal dynamics under a broader range of annual environmental conditions. We, therefore, focused our analyses on grasslands because their within-season dynamics have been found to be challenging to model [15]. Our study was not intended to validate different approaches or models for characterizing seasonal dynamics, rather we investigated the potential improvements and challenges of modeling seasonal dynamics using HLS data. We recognize that models developed in the future specifically designed for HLS data frequency would likely provide further improvements.

7. Conclusions

In this study, we evaluated whether time series data harmonized from Landsat-8 and Sentinel-2 data could provide sufficient frequency of observations to enable monitoring seasonal dynamics and disturbances in central U.S. grasslands at a spatial resolution of 30 m. We applied two previously established algorithms to identify the timing of the start of the growing season for vegetation and to delineate burn scars from grassland fires, and we used harmonic analysis to model wetland seasonal dynamics. Our study was limited by having only one complete year (2016) of harmonized data; however, this was sufficient to illustrate that: (1) results derived from harmonized time series nearly always outperformed results derived from Landsat-8 data alone; (2) the varying observation frequency in the temporal domain could cause inconsistent accuracy in model estimates; (3) the current observation frequencies from the two satellite sensors caused spatial inconsistency in data density even after applying a weekly maximum composite process; and (4) the harmonized time series still imposed limitations in data frequency during the early part of the growing season, when landscape conditions are prone to change quickly. To have more than a 50% probability of weekly clear observations throughout the growing season for our study areas, at least a 3-day revisit interval would be required. Our results using Landsat time series enriched with data from comparable sensors showed promise for tracking within-season dynamics in grassland landscapes at a 30-m resolution. In the future, the potential to harmonize data across Sentinel-2A, Sentinel-2B, Landsat-8, and the planned Landsat-9 mission means that the global median average revisit interval for 30-m data will be less than 3 days [15,45], improving the probability for weekly clear observations to study a variety of complex landscape applications, such as the environmental dynamics of grasslands. Therefore, maintaining the continuity for future Landsat missions and interoperability between Landsat and similar sensors, such as those on Sentinel-2, will be critical.

We see many prospects for building on this research. Future work on improving the monitoring of grassland phenology and disturbance should be carried out as HLS data are improved and the period of record and data frequency are increased. A study over multiple years of HLS recently available will contribute to understanding and monitoring the range of conditions over multiple seasons with variability in climate and vegetation condition (i.e., a range of wet to dry conditions). Opportunities to learn more about how the spatial resolution of HLS enhances grassland and wetland monitoring will be supported by more dense networks of ground and near-ground ecological observations such as those being collected by the National Ecological Observatory Network (NEON) Terrestrial Observation System [46]. The integrated analyses of terrestrial observations of many different ecosystems with sensor-based, imagery, and remote-sensing data collected by NEON subsystems can facilitate the investigation of measured grassland and wetland parameters over multiple spatial scales.