Evaluating the Impact of Grazing Cessation and Reintroduction in Mixed Prairie Using Raster Time Series Analysis of Landsat Data

: Great efforts have been made to manage and restore native prairies to protect native species, enrich biodiversity, protect ecological resilience, and maintain ecosystem services. Much of this has been focused on preventing degradation from overgrazing and crop conversion. Understanding the consequences of management polices is important to identify best practices. Previous research has compared restoration outcomes from variable intensity grazing, prescribed ﬁre, and grazing removal. However, few studies have explored the optimal durations of management practices and variation in restoration outcomes among vegetation communities. This study evaluates whether the impact of grazing cessation and reintroduction varies among native vegetation communities and measures the effective time periods of grazing cessation and reintroduction. Restoration outcomes were evaluated using four biophysical indicators (fresh biomass, soil organic matter, green cover, and litter cover) and two vegetation indices (normalized difference vegetation index (NDVI) and normalized difference water index (NDWI)) measured from Landsat images using seasonal Kalman ﬁlter and raster time series analysis. The results show that: (i) Grazing cessation increased soil organic matter and green cover while decreasing fresh biomass compared to moderate grazing management, while grazing reintroduction inﬂuences those indicators in an opposite direction; (ii) The effective time period for prairie conservation is about 11–14 years and varies among vegetation communities and biophysical indicators; (iii) The effective intensity of grazing cessation is highest in valley grassland, moderate in upland grassland, and mildest in sloped grassland; (iv) Grazing reintroduction returned the three native vegetation communities to the initial condition (i.e., the stage in 1985 before large grazers were removed), with less time than the time consumed for grazing cessation to restore the prairie ecosystem to the maximum changes; (v) Grazing reintroduction effectively inﬂuences upland and valley grasslands for 7 to 9 years, varying from different indicators, while it continuously affected sloped grassland with no clear time lag; (vi) The intensity of grazing reintroduction was strongest in sloped grassland, moderate in upland grassland, and mildest in valley grassland. The results of this study suggest expected time periods for prairie management methods to achieve results.


Introduction
Native prairies in the Great Plains of North America not only provide diverse habitats and quality forage for livestock and wildlife [1], but also maintain a variety of ecological functions including carbon storage, water conservation, soil stability, low-input grazing, They have enabled consistent, long-term spatial analysis, and have been widely used for temporal analysis in ecosystem monitoring [42].
One measure that continues to impart inconsistency and bias to the temporal analysis is seasonal variation due to plant phenology [42,44]. Studies that have performed temporal analysis using the Landsat series have often acquired images in the same season in recent years [45], or obtained consistent images of peak growing season conditions [46,47]. However, limited observation frequency and cloud effects (prolonged precipitation in prairie limits available optical images) compounds the issue of acquiring growing season Landsat images in prairie regions [43,45,48]. Therefore, it would advance the field of grassland monitoring to develop new methods that minimize seasonal effects with the capacity to account for missing data in long-term temporal studies in prairie regions.
Our study explored the optimum duration of grazing cessation and reintroduction on a native prairie with a known history of grazing by native and domestic herbivores and their variation among prairie types. The specific objectives were: (1) To develop a new method to minimize seasonal effects for the temporal analysis of native prairies based on Landsat images; (2) To analyze the optimum duration for grazing cessation and reintroduction; (3) To explore the variability of vegetation response to management practices among three native vegetation communities, including upland, sloped, and valley grasslands.

Study Area
The study area is located in Grasslands National Park (GNP; Figure 1) and the surrounding pastures, in southwestern Saskatchewan within the northern mixed grass prairie region. Most of the lands for GNP were acquired for the purpose of prairie conservation and restoration in 1985, and Parks Canada continuously purchased lands in native prairie regions until 2000 ( Figure 1). Prior to gazetting as a national park, the land was grazed by cattle with small parcels converted to tame pastures, hay, or limited cropping. In this study, only the lands acquired in 1985 that were previously cattle grazed native prairie were analyzed (Figure 1). responses of native prairie ecosystems to management practices that vary in frequency, intensity, and duration [24]. The Landsat series, with almost 50 years of images (since 1972) have demonstrated the potential to monitor prairie response to restoration practices [42,43]. The Landsat series also offers abundant historical information to benefit our understanding of historical disturbances, reducing uncertainty in prairie restoration outcomes [17,24]. They have enabled consistent, long-term spatial analysis, and have been widely used for temporal analysis in ecosystem monitoring [42].
One measure that continues to impart inconsistency and bias to the temporal analysis is seasonal variation due to plant phenology [42,44]. Studies that have performed temporal analysis using the Landsat series have often acquired images in the same season in recent years [45], or obtained consistent images of peak growing season conditions [46,47]. However, limited observation frequency and cloud effects (prolonged precipitation in prairie limits available optical images) compounds the issue of acquiring growing season Landsat images in prairie regions [43,45,48]. Therefore, it would advance the field of grassland monitoring to develop new methods that minimize seasonal effects with the capacity to account for missing data in long-term temporal studies in prairie regions.
Our study explored the optimum duration of grazing cessation and reintroduction on a native prairie with a known history of grazing by native and domestic herbivores and their variation among prairie types. The specific objectives were: 1) To develop a new method to minimize seasonal effects for the temporal analysis of native prairies based on Landsat images; 2) To analyze the optimum duration for grazing cessation and reintroduction; 3) To explore the variability of vegetation response to management practices among three native vegetation communities, including upland, sloped, and valley grasslands.

Study Area
The study area is located in Grasslands National Park (GNP; Figure 1) and the surrounding pastures, in southwestern Saskatchewan within the northern mixed grass prairie region. Most of the lands for GNP were acquired for the purpose of prairie conservation and restoration in 1985, and Parks Canada continuously purchased lands in native prairie regions until 2000 ( Figure 1). Prior to gazetting as a national park, the land was grazed by cattle with small parcels converted to tame pastures, hay, or limited cropping. In this study, only the lands acquired in 1985 that were previously cattle grazed native prairie were analyzed ( Figure 1).  After the lands were acquired for the national park, all livestock grazing ceased, with only very light use by native grazers (mule deer (Odocoileus hemionus), white-tailed deer (O. virginianus), pronghorn (Antilocapra americana)) occurring [49]. This initiated the process of secondary succession by modifying the pressures and disturbances [50] and allowed the native prairie vegetation community to rest and conserve its biodiversity. After a long period of grazing cessation (i.e., the management policy of removing larger grazers), large Remote Sens. 2021, 13, 3397 4 of 21 amounts of dead materials accumulated, increasing fire risk [51]. In 2006, partially in response to this, and to maintain biodiversity, plains bison were introduced to the study area [52].
Three native vegetation communities in GNP are upland, sloped, and valley grasslands ( Figure 2). Shrub communities also occur along the Frenchman River valley, traversing the west block of GNP ( Figure 2). This block of GNP also contains eroded and disturbed (old fields of seeded tame species-smooth brome, crested wheatgrass and sweet clover) communities. Five isolated prairie dog colonies also occur within upland grasslands ( Figure 2).
After the lands were acquired for the national park, all livestock grazing ceased, with only very light use by native grazers (mule deer (Odocoileus hemionus), white-tailed deer (O. virginianus), pronghorn (Antilocapra americana)) occurring [49]. This initiated the process of secondary succession by modifying the pressures and disturbances [50] and allowed the native prairie vegetation community to rest and conserve its biodiversity. After a long period of grazing cessation (i.e., the management policy of removing larger grazers), large amounts of dead materials accumulated, increasing fire risk [51]. In 2006, partially in response to this, and to maintain biodiversity, plains bison were introduced to the study area [52].
Three native vegetation communities in GNP are upland, sloped, and valley grasslands ( Figure 2). Shrub communities also occur along the Frenchman River valley, traversing the west block of GNP ( Figure 2). This block of GNP also contains eroded and disturbed (old fields of seeded tame species-smooth brome, crested wheatgrass and sweet clover) communities. Five isolated prairie dog colonies also occur within upland grasslands ( Figure 2).

Landsat Images for Temporal Analysis
A total of 126 cloud-free Landsat images (Table 1), including 2 Thematic Mapper (TM) images from Landsat 4, 87 TM images from Landsat 5, 12 Enhanced Thematic Mapper Plus (ETM+) images from Landsat 7, and 25 Operational Land Imager (OLI) images from Landsat 8, were collected from 1986 to 2020 (i.e., only images acquired from April to October were used for temporal analysis to avoid ice and snow cover). All Landsat images were preprocessed (i.e., with both geometric and atmospheric correction; Level 2 data from Landsat Collection 2) when downloaded from the United States Geological Survey (USGS) website (https://earthexplorer.usgs.gov/, accessed on 14 January 2020).

Landsat Images for Temporal Analysis
A total of 126 cloud-free Landsat images (Table 1), including 2 Thematic Mapper (TM) images from Landsat 4, 87 TM images from Landsat 5, 12 Enhanced Thematic Mapper Plus (ETM+) images from Landsat 7, and 25 Operational Land Imager (OLI) images from Landsat 8, were collected from 1986 to 2020 (i.e., only images acquired from April to October were used for temporal analysis to avoid ice and snow cover). All Landsat images were preprocessed (i.e., with both geometric and atmospheric correction; Level 2 data from Landsat Collection 2) when downloaded from the United States Geological Survey (USGS) website (https://earthexplorer.usgs.gov/, accessed on 14 January 2020).
According to the Level 2 data information from Landsat Collection 2 (https://www. usgs.gov/core-science-systems/nli/landsat/landsat-collection-2-level-2-science-products, accessed on 20 April 2021), the algorithm for generating surface reflectance data from Collection 2, Level 1 products of Landsat 4-5 TM and Landsat 7 ETM+ was the Landsat Ecosystem Distribution Adaptive Processing System (LEDAPS) algorithm (Version 3.4.0), and that of Landsat 8 OLI was the Land Surface Reflectance Code (LaSRC) algorithm (Version 1.5.0). The spatial resolution of all Landsat images was 30 × 30 m, which were projected as UTM Zone 13N, WGS1984. The pixel values of multispectral bands from the Level 2 products were converted to surface reflectance with a scaling factor of "0.0000275 × pixel value − 0.2" (e.g., if the pixel value is 32,000, the reflectance is 0.68). Pixel values of thermal bands from Landsat Collection 2, Level 2, were converted to surface temperature with a scaling factor of "0.00341802 × pixel value + 149.0 − 273.15" (e.g., if the Remote Sens. 2021, 13, 3397 5 of 21 pixel value is 45,000, the surface temperature is 29.6609 • C; the thermal band of Landsat 8 used in this study is band 10).

Classification of Vegetation Communities
We conducted Principal Component Analysis (PCA) to the band composition of the DEM image and all eleven Sentinel-2 images (i.e., Normalized Difference Vegetation Index (NDVI) calculated for each Sentinel-2 image before PCA). The first band from PCA results was excluded before unsupervised classification because this band largely consists of DEM information instead of vegetation. An ISO unsupervised classification based on ArcGIS 10.5 software was applied to the band composition of PCA results excluding the first component. Post-classification techniques, including reclassification, majority filter, boundary clean, region group and nibble were also applied to the classification results to improve the results for vegetation communities. The information from prairie dog towns, annual crops, and some disturbed communities were incorporated into classification results based on 1983 and 1995 field surveys [53]. The overall accuracy for the final classification map of vegetation communities (Figure 2) was 89.7%, based on 100 randomly distributed sample points using visual interpretation.

Links between Surface Reflectance of Landsat Images and Grassland Biophysical Parameters
Our previous research on prairie responses to grazing cessation in GNP [50] indicated that the difference of biophysical parameters of conserved and grazed sites had a linear relationship with the difference of surface reflectance (or land surface temperature) in Landsat bands. This research showed that the difference in near-infrared (NIR) reflectance, first shortwave infrared (SWIR1) reflectance and second shortwave infrared (SWIR2) reflectance, and land surface temperature (LST) from the thermal band, are all significantly correlated with the differences of fresh biomass, soil organic matter, green cover and litter cover (Table 2), respectively [50]. Therefore, we used these relationships as indicators of the prairie response to grazing management, considering the potentially overwhelming influence of climate change.

Formation of Time Series Dataset by "Zoo Objects" and Missing Value Interpolation
Time series datasets were prepared using "zoo objects" in the R software package (https://www.r-project.org/, accessed on 16 March 2020) for each pixel of NIR, SWIR1, and SIWR2 reflectance and LST. First, a raster time series object was generated from a raster stack object in the 126 Landsat images and the corresponding acquisition date for each image in R software. Then, a time series dataset was formed for each pixel using a "for loop" in the zoo package. To do this, the temporal dataset for each pixel was aggregated Remote Sens. 2021, 13, 3397 7 of 21 into monthly data using the "aggregate" function (i.e., if more than one image was acquired in the same month, their data were aggregated by their mean value). This aggregated temporal dataset had missing values in the months for which cloud-free Landsat images were unavailable. To form a regular "zoo object", each pixel was assigned a specific frequency of 7 (the number of acquisition months from April to October). The missing values for each pixel were then interpolated using a seasonal Kalman filter (the na.StructTS() function in the zoo package). The seasonal Kalman filter has been shown to be a functional method to filter missing values in a dataset with repeated seasonal variability [54]. As a final step, the regular zoo object for each pixel of NIR, SWIR1, SWIR2, and LST was generated without any missing values. NDVI is the most popular vegetation index for representing green vegetation information [55,56], and Normalized Difference Water Index (NDWI) is a good indicator for plant water content [57][58][59]. Therefore, in addition to these four bands, the regular zoo objects for NDVI [60] Equation (1) and NDWI [61] Equation (2) were also generated.
where NDVI refers to the normalized difference vegetation index, NIR refers to the nearinfrared band (i.e., band 4 for TM and ETM+ images and band 5 for OLI images); R refers to the red band (i.e., band 3 for TM and ETM+ images and band 4 for OLI images).
where NDWI refers to normalized difference water index, SWIR1 refers to the first shortwave infrared band (i.e., band 5 for TM and ETM+ images and band 6 for OLI images).

Pixel Based Time Series Analysis
Based on the filtered regular zoo object for each pixel, a time series object was formed for each pixel (i.e., start in 1986 with a frequency of 7) using the "ts" function in the raster time series (rts) package from the R software. Time series analysis was then applied to decompose the temporal object of each pixel into trend component, seasonal component, and random noise component using the "stl" function in the "rts" package. With time series analysis, seasonal variation is minimized for maintaining temporal consistency of longterm studies [47]. The median value into the trend component of each pixel was extracted to calculate a new raster for each year. The NIR, SWIR1, SWIR2, LST, NDVI, and NDWI, with limited seasonal effects for these rasters, were used for further temporal analysis.

Segmented Linear Regression
Zonal statistics in ArcGIS 10.5 were applied to calculate the mean values of newly generated rasters for each year from 1986 to 2020. We did this using polygons formed from the intersection of different vegetation communities ( Figure 2) and the land acquisition year ( Figure 1). Then, temporal trends (for NIR, SWIR1, SWIR2, LST, NDVI, and NDWI independently) were calculated based on the mean value of zonal statistic results (see the selection of three vegetation types in the surrounding pastures in Supplementary, Figure S1). The temporal trends of the difference in each multispectral band and index for a certain vegetation community were calculated by using the temporal trend of the GNP sites (acquired in 1985) minus that of the surrounding pastures. Then, the temporal trends of the difference in NIR, SWIR1, SWIR2, LST, NDVI, and NDWI between 1985 and 2020 was smoothed by the loess function in R software. Finally, the segmented linear regression in the "segmented" package in R software was applied to each smoothed temporal trend to evaluate the significant turning points. The segmented linear relationship (i.e., broken-line relationship) has a variety of applications in ecological research [51]. In this study, the segmented linear relationship was used to analyze thresholds for the optimum duration of prairie restoration practices including grazing cessation and reintroduction. After the turning points were measured for each smoothed trend, the slope was calculated for NIR reflectance of upland, sloped, and valley grasslands in GNP was higher than that of surrounding pastures in 1985 (the year GNP was first acquired as a national park), based on the temporal trend of the difference in NIR reflectance between GNP sites and the surrounding pastures (Figure 3a-c). The difference in NIR reflectance gradually decreased until it reached a relatively stable stage when GNP was under grazing cessation from 1985 to 2006 (Figure 3a-c). However, the effective time periods for the impact of grazing cessation on three vegetation communities were different. The time periods when grazing cessation caused a significant change in the difference in NIR for upland, sloped, and valley grasslands were 1985-1999, 1985-1997, and 1985-1996, respectively (Table 3; Figure 3a-c). During the effective time periods, valley grassland had the strongest response to grazing cessation (the absolute slope is 0.0029 of the first linear segment; Table 3), while upland grassland (the absolute slope is 0.0014 of the first linear segment; Table 3) and sloped grassland (the absolute slope is 0.0021 of the first linear segment; Table 3) had a milder response to grazing cessation based on the temporal dynamics of the difference in NIR reflectance. After these time periods, grazing cessation had no statistically significant influences on the three vegetation communities anymore (Figure 3a-c). After bison were reintroduced in 2006, the temporal trend of the difference in NIR reflectance for upland grassland changed in reverse of the initial changes when it first began to be conserved in 1985, while it became stable after 2014, 8 years after grazing reintroduction was implemented ( Figure 3a). However, the difference in NIR reflectance remained stable for the first four and six years for sloped and valley grasslands, respectively, and then increased afterwards when they were under grazing reintroduction management (Figure 3b,c). Grazing reintroduction had a stronger impact on sloped and valley grasslands (absolute slope of second linear segment is 0.002 for sloped grassland and 0.0019 for valley grassland; Table 3) than upland grassland (absolute slope of first linear segment is 0.0013; Table 3).
Except in the valley grassland, the difference in SWIR1 reflectance for upland and sloped grasslands in GNP in 1985 was lower than the surrounding pastures (Figure 3d-f). Grazing cessation caused a significant decrease in the difference in SWIR1 reflectance between GNP and the surrounding pastures in the first 14 years for upland grassland (1985-1999; Figure 3d) and in the first 13 years for sloped grassland (1985-1998; Figure 3e). According to the temporal changes on the difference in SWIR1 reflectance, the intensity of the impact of grazing cessation was slightly stronger in upland grassland than in sloped grassland (i.e., the absolute slope of the first linear segment for upland and sloped grassland were 0.0016 and 0.0014, respectively; Table 3). After these time periods, grazing cessation caused the difference in SWIR1 reflectance increasing with a much milder intensity (Figure 3d,e; the absolute slope of the second linear segment for upland and sloped grassland were 0.0004 and 0.0003, respectively). The impact of grazing cessation on the difference in SWIR1 reflectance did not show a time lag effect for valley grassland; instead, it caused a continuous decrease in the temporal trend of SWIR1 difference from 1985 to 2006 until grazing reintroduction broke the temporal trend (Figure 3f). After bison were reintroduced in GNP, the difference in SWIR1 increased for upland grassland during 2006-2013, for sloped grassland during 2006-2020, and for valley grassland during 2006-2014 (Figure 3d-f). The impact of grazing reintroduction was stronger for upland and sloped grasslands than it was for valley grassland (i.e., the absolute slope of the first linear segment was 0.0014, 0.0017, and 0.0009 for upland, sloped and, valley grasslands, respectively; Table 3). However, grazing reintroduction was no longer effective for upland grassland since 2013, and was much less effective for valley grassland after 2014 (Figure 3d-f).
Remote Sens. 2021, 13, x 10 of 21 the difference in LST to increase until 2014, with a linear slope of 0.1, 0.15, and 0.25 for upland, sloped, and valley grasslands, respectively (Table 3).    Note: The effective periods were analyzed by the turning points from segmented linear regression: "Slope" refers to the slope of a certain line segment from segmented linear regression, which indicates the intensity of the influences of grazing cessation and reintroduction on three native vegetation communities. "NA" means there was significant turning points when the vegetation community was under grazing cessation, and reintroduction or "NA" means the slope of the linear segment is not statistically significant.

of 21
In 1985, the difference in SWIR2 reflectance for upland and valley grasslands in GNP was higher than the surrounding pastures, and sloped grassland was almost the same for both GNP and surrounding pastures (Figure 3g-i). Grazing cessation caused a significant decrease in the difference in SWIR2 reflectance between GNP and the surrounding pastures in the first 13 years for both upland and sloped grasslands (1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998); Figure 3g,h). However, the intensity of the impact of grazing cessation on the difference in SWIR2 reflectance was greater for upland grassland than for sloped grassland (i.e., the absolute slope of the first linear segment was 0.0012 and 0.0008 for upland and sloped grasslands, respectively; Table 3). The temporal trend of the difference in SWIR2 during grazing cessation implementation period did not show a statistically significant turning point for valley grassland, instead, grazing cessation caused a continuous decrease in the difference of SWIR2 reflectance from 1985 to 2006 (Figure 3i). Since 2006, when the grazing reintroduction was implemented, the difference in SWIR2 increased for upland grassland during 2006-2013 and valley grassland during 2006-2015 (Figure 3g,i), while it continued to increase for sloped grassland from 2008 to 2020 (Figure 3h). Grazing reintroduction had the strongest effects on the difference in SWIR2 reflectance for sloped grassland (the absolute slope of the linear regression was 0.0017; Table 3) than upland and valley grasslands (the absolute slope of the first linear segment for upland and valley grasslands was 0.0014 and 0.0009, respectively; Table 3).
Grazing cessation had no significant impact on the difference in LST for upland grassland in the first 6 years, for sloped grassland in the first 9 years, and for valley grassland in the first 10 years (Figure 3j-l; Table 3). However, after these periods, the difference in LST between GNP and surrounding pastures slightly dropped while GNP was continuously under grazing cessation (Figure 3j-l). Valley grassland had the greatest decline in the difference in LST between 1995 and 2005 (i.e., the absolute slope of the second linear segment is 0.11; Table 3), followed by upland and sloped grasslands (i.e., the absolute slopes of the second linear segment are 0.04 and 0.01 for upland and sloped grasslands, respectively; Table 3). Grazing reintroduction had much stronger influences on the temporal trend of the difference in LST (Figure 3j-l). Since 2006, grazing reintroduction caused the difference in LST to increase until 2014, with a linear slope of 0.1, 0.15, and 0.25 for upland, sloped, and valley grasslands, respectively (Table 3).

Temporal Trends of the Differences in NDVI and NDWI between GNP and Surrounding Pastures
The difference in NDVI gradually decreased in upland during 1985-2000, in sloped grassland during 1985-1997, and in valley grassland during 1985-1997 (Figure 4a-c). The decreasing pace for valley grassland is larger than that for upland and sloped grasslands (i.e., the absolute slope of the first linear segment is 0.0007, 0.0012, and 0.0015 for upland, sloped, and valley grasslands, respectively; Table 4). After these periods, the difference in NDVI of upland and sloped grasslands became stable, while grazing cessation influenced the difference in NDVI of valley grassland in an opposite direction of the initial changes after 1997 (Figure 4a-c; Table 4).
After bison were reintroduced in GNP in 2006, the difference in the NDVI increased in upland grassland (Figure 4a Table 4).
The temporal trends of the difference in NDWI show similar trends with NDVI, though different turning points and intensity between sloped and valley grasslands when both communities were under grazing cessation and reintroduction managements (Figure 4e,f; Table 4). In upland grassland, the difference in NDWI gradually decreased during the whole period of grazing cessation and increased grazing reintroduction from 2006 to 2020 (Figure 4d).  Note: The effective periods were analyzed by the turning points from segmented linear regression: "Slope" refers to the slope of linear regression for the certain line segment from segmented linear regression, which indicates the intensity of the influences of grazing cessation and reintroduction on three native vegetation communities. "NA" means there was significant turning points when the vegetation community was under grazing cessation, and reintroduction or "NA" means the slope of the linear segment is not statistically significant.  Note: The effective periods were analyzed by the turning points from segmented linear regression: "Slope" refers to the slope of linear regression for the certain line segment from segmented linear regression, which indicates the intensity of the influences of grazing cessation and reintroduction on three native vegetation communities. "NA" means there was significant turning points when the vegetation community was under grazing cessation, and reintroduction or "NA" means the slope of the linear segment is not statistically significant.

Monitoring Prairie Management Effects Using a Landsat Imagery Raster Time Series Analysis
Long-term temporal analysis has been shown to be essential to monitor prairie responses to management change. However, temporal inconsistency in image series challenges long-term analysis when using remotely sensed imagery in prairie research. Previous studies have used images acquired in the maximum growing season to analyze temporal dynamics of ecosystems [42]. However, the low temporal resolution of Landsat images may miss the reflectance of the peak growth due to the short northern prairie growing season [43]. The maximum growing season corresponds to rainfall peaks, resulting in large intervals of cloud cover and missing data during key growing periods. The short lags between rainfall and grassland community growth confounds the consistent use of remote sensing in prairie ecosystems.
Seasonal variation also affects the temporal monitoring of prairie ecosystems because of vegetation phenology. We applied a seasonal Kalman filter to fill missing values caused by cloud effects based on the seasonal cycles in long-term Landsat images. These vary across vegetation communities with different penology [54]. The pixel-based seasonal Kalman filter was applied to generate the regular time series raster object and minimize the effects of seasonal variation among vegetation communities. The filter generated better time series data than spline interpolation in particular during a year with many missing images. Temporal trends were separated by time series analysis to eliminate the seasonal effects [47], which greatly enhanced the temporal consistency of the reflectance of Landsat images.

The Relationship between the Difference in Biophysical Parameters and Multispectral Reflectance
Green vegetation cover [62][63][64] and biomass [65][66][67] are two of the most common biophysical indicators used to evaluate management responses in prairie ecosystems. These have been estimated using statistical models based on field-measured data and spectral indices based on red, NIR, or SWIR1 bands [68]. In GNP, a large amount of non-photosynthetic materials challenges estimates of fresh biomass and green cover from Landsat images [69]. We were able to take advantage of our previous research relating to the difference in NIR and SWIR2 reflectance and the difference in fresh biomass and green cover between GNP sites and the surrounding pastures, respectively [50]. Soil organic matter and litter cover are two indicators that are not easily measured using Landsat reflectance [58], yet play a strong role in prairie management. Previous research shows that SWIR1 reflectance is well correlated to soil organic matter [50,70,71]. Litter accumulation decreases LST [50], which is consistent with the ecological function of litter that enhances soil moisture and mitigates soil surface temperature [51]. Therefore, we take advantage of the results from our previous research [50] to use the difference in NIR, SWIR1, SWIR2, and LST to predict the difference in fresh biomass, soil organic matter, green cover, and litter cover between GNP and the surrounding pastures, respectively (Table 2). More importantly, it highly reduced the effects of climate change by using the temporal trends of the difference in NIR, SWIR1, SWIR2, and LST for the vegetation communities in GNP and the surrounding grazing pastures, because the same vegetation type in GNP and the surrounding pastures was under the same intensity of climate change for a given time, albeit under different management policies.
Our previous research [50] used the five paired ecological comparison sites (i.e., all the sites are upland grasslands) in GNP and the surrounding pastures to form the linear relationship between the difference in NIR, SWIR1, SWIR2, and LST and the difference in fresh biomass, soil organic matter, green cover, and litter cover, respectively (Table 2). Each paired site has similar vegetation composition and soil properties, therefore, the difference in NIR, SWIR1, SWIR2, and LST were mainly caused by difference management policies between GNP and the surrounding pastures. However, there are still some biases and uncertainties when using the difference in NIR, SWIR1, SWIR2, and LST between GNP and the surrounding pastures. First, the linear models of the difference in Landsat reflectance and the difference in biophysical parameters did not exactly fit the assumption of homoscedasticity, which would reduce the statistic power. Second, it would bring uncertainties by directly applying those linear models on Landsat TM, ETM+ and OLI images, because the initial linear models were built based on Landsat TM images and field data ( Table 2). Third, the limited samples of field data also make the linear models less reliable (i.e., 8 paired samples for fresh biomass and soil organic matter collected in 2004 and 2005, and 13 paired samples for green cover and litter cover collection in 2003, 2004, and 2005). Fourth, when this study was extended into the whole study area, the heterogeneity of the mixed prairie makes it impossible to find the pixels in the surrounding pastures, which have similar vegetation composition and soil properties to pair each of the pixels of GNP sites, including upland, sloped, and valley grasslands. The difference data used in this study was calculated by the mean value of upland, sloped, and valley grasslands in GNP minus the mean value of typical upland, sloped, and valley grassland in the surrounding pastures, respectively. However, the inconsistency of species composition, soil properties and canopy information would bring biases to the difference in the temporal trends between GNP and surrounding pastures. Finally, the accuracy of vegetation classification (Figure 2) also brought bias to the results because the calculation of the mean value of difference vegetation communities relies on the classification results. The accuracy of our classification map (89.7%) might have some biases because it was conducted by visual interpretation of 100 random points instead of field validation data.

The Impact of Grazing Cessation on Three Native Vegetation Communities
According to the positive linear relationship between the difference in fresh biomass and NIR reflectance (Table 2), the difference of NIR reflectance between GNP sites and the surrounding pastures indicates a higher fresh biomass in the park than that of the surrounding pastures in 1985, when GNP was first acquired for all three vegetation types (Figure 3a-c). The impact of grazing cessation on the temporal dynamics of the difference in NIR reflectance indicates that grazing cessation, compared to moderate grazing management in the surrounding pastures, decreased fresh biomass in upland, sloped, and valley grassland within a certain time period of 11-14 years, varying from different vegetation communities (Figure 3a-c; Table 3). The intensity of grazing cessation on fresh biomass varies from three vegetation communities as well. Under the implementation of grazing cessation, the difference in fresh biomass between GNP and surrounding pastures decreased 30.24 g for upland grassland during 1985-1999, 38.88 g for sloped grassland during 1985-1997, and 49.22 g for valley grassland during 1985-1996 (i.e., calculated based on slope of the linear segment and positive linear relationship between the difference in fresh biomass and the difference in NIR reflectance; Tables 2 and 3). The stability of the temporal trend of the difference in NIR reflectance after those time periods suggests no effectiveness of grazing cessation on fresh biomass (Figure 3a-c; Table 3).
The decrease in the difference of SWIR1 reflectance reveals the increase in the difference of soil organic matter between GNP and the surrounding pastures, caused as a result of the negative linear relationship between the difference in soil organic matter and the difference in SWIR1 reflectance ( Table 2). Compared to the moderate grazing management in the surrounding pastures, grazing cessation increased soil organic matter in upland grassland for the first 14 years and in sloped grassland for the first 13 years, while it continuously increased soil organic matter in valley grassland for all 21 years under grazing cessation. The intensity of the impact of grazing cessation on GNP varies from different vegetation types as well. Grazing cessation caused an increase in soil organic matter by 0.51% for upland grassland, 0.55% for sloped grassland, and 5.73% for valley grassland (Tables 2 and 3).
Our previous research shows that the difference in green cover between GNP sites has a negative linear relationship with the difference in SWIR2 reflectance (Table 2). Thus, the decrease in the difference of SWIR2 reflectance indicates an increase in the difference of green cover between GNP and the surrounding pastures for upland and sloped grasslands during 1985-1998, as well as for valley grassland during all 21 years under grazing cessation (Figure 3g-i), compared to the moderate grazing management in the surrounding pastures. However, the intensity of the influences of grazing cessation on green cover varies for three vegetation communities. Grazing cessation caused an increase in green cover by 5.43% for upland grassland, 3.62% for sloped grassland, and 9.50% for valley grassland (Tables 2 and 3).
The differences in NIR, SWIR1, and SWIR2 between GNP and the surrounding pastures have shown similar temporal trends as those seen in previous research [50]. However, the difference in LST does not show similar temporal trends with our previous research, which indicates that other factors besides litter cover, including species composition, vegetation canopy, water content, and soil properties, might override the effects of litter cover difference on LST.
When using fresh biomass, soil organic matter, and green cover as indicators, grazing cessation had the strongest impact on valley grassland and the mildest impact on sloped grassland, mostly within a time period of 11 to 14 years instead of the whole management period. This may be because the vegetation in valley grasslands have a higher aboveground biomass and height, with more accumulation of soil organic matter. It is also possible that greater moisture in the valley grassland would increase the decomposition of dead material, resulting in more green cover due to less competition for growing space with standing dead materials [51]. The increase in the difference of soil organic matter and green cover suggests that grazing cessation is more effective than long-term moderate grazing management (i.e., a management policy for the surrounding pastures). The decrease in the difference of fresh biomass when GNP was under grazing cessation management is due to a large amount of accumulation of non-photosynthetic materials.
The decrease in the difference of NDVI and NDWI when GNP was under grazing cessation is also due to dead material accumulation, as it is generally a restoration objective of grazing cessation in prairies [8]. Based on the temporal trends of the difference in NDVI and the difference in NDWI, grazing cessation influenced upland grassland with the longest period but the mildest intensity, while it impacted valley grassland with the shortest time period but the largest intensity among all three vegetation communities ( Figure 4 and Table 4). However, the difference between GNP and the surrounding pastures was not visibly obvious in NDVI or NDWI images, while the difference was clear in NIR, SWIR1, and SWIR2 bands, especially the NIR band (see Supplementary, Figure S1 for the images acquired in 2006).

The Impact of Grazing Reintroduction on Three Native Vegetation Communities
Based on the temporal trend in the difference of NIR reflectance (Figure 3a-c) and the relationship between fresh biomass and NIR reflectance (Table 2), grazing reintroduction increased fresh biomass by 16.05 g during 2006-2014 for upland grassland and then became stable (Figure 3a), while it was not effective for fresh biomass in sloped grassland for the first 4 years and valley grassland for the first 6 years (Figure 3b,c). After the first stable stage, grazing reintroduction increased fresh biomass by 30.86 g in sloped grassland from 2010-2020 and by 23.45 g in valley grassland from 2012-2020 (Tables 2 and 3). The impact of grazing reintroduction was slightly better at enhancing fresh biomass in sloped and valley grasslands than in upland grassland, compared to moderate grazing management.
Grazing reintroduction decreased soil organic matter by 1.78% in upland grassland from 2006 to 2013 and by 1.31% in valley grassland between 2006 and 2014, while it continuously decreased soil organic matter by 3.71% in sloped grassland until 2020 (Tables 2 and 3). Compared to grazing cessation, soil organic matter restored in at the slowest pace in sloped grassland (absolute slope of the first linear segment is 0.0014; Table 3). Among three vegetation types, however, it degraded at the highest speed while it was under grazing reintroduction management (absolute slope of the first segment is 0.0017; Table 3). The impact of grazing reintroduction on soil organic matter for valley grassland was the mildest among the three vegetation types.
Since 2006, the year grazing reintroduction was implemented in GNP, green cover decreased by 2.44% in upland until 2013, by 5.01% in valley grassland until 2015, and by 7.52% in sloped grassland until 2020 (i.e., no clear time lag effects showed for sloped grassland). For sloped grassland, green cover recovered in the slowest pace under grazing cessation, but it decreased in the highest speed under grazing reintroduction among the three vegetation types (Table 3).
According to the changes of the fresh biomass, soil organic matter, and green cover, grazing reintroduction had a shorter effective time period and milder intensity than grazing cessation for upland and valley grasslands (Figure 3; Tables 2 and 3). Among all three types, grazing reintroduction had the greatest intensity of impact on sloped grassland and the smallest intensity of influences on valley grassland. For sloped grassland, grazing reintroduction had stronger effects than grazing cessation in a reversed direction.
Based on the temporal trends of the difference in NDVI and NDWI between GNP sites and surrounding pastures, grazing reintroduction influenced upland grassland for the longest period but with the mildest intensity, while it was effective on sloped and valley grassland, with a time lag of 7 to 8 years ( Figure 4 and Table 4). After 7 to 8 years under grazing reintroduction management, the differences between NDVI and NDWI in sloped and valley grasslands changed in an opposite direction, just as GNP did under grazing cessation (Figure 4). This phenomenon was shown in the temporal trend of the difference in NIR reflectance, as well in sloped and valley grasslands (Figure 3b,c). This may be due to the selective grazing by bison in the first few years when they were first introduced into GNP [50].
The dynamics of the most temporal trends of the differences in Landsat bands or indices indicates that grazing reintroduction returned upland, sloped, and valley grasslands to their initial stage in 1985 in less time than the time of the maximum changes caused by grazing cessation (Figures 3 and 4).

The Influences of Climate Change on Three Native Vegetation Communities
It is generally a challenge to separate the influences of climate change and prairie management with remote sensing applications, using biophysical parameters as indicators to evaluate the response to climate change or prairie management. With long-term Landsat images, the temporal trends of different vegetation indices were often used to analyze the impact of prairie management or climate change. However, vegetation communities have similar temporal trends of NDVI and NDWI ( Figure 5), suggesting that climate effects have a high potential to override the impact of prairie management (i.e., the grazing cessation and reintroduction in this study).
When analyzing the temporal trends of NIR, SWIR1, SWIR2, and LST for three native vegetation communities in GNP, the increase in NIR, SWIR1, and SWIR2 reflectance and LST reveals an increase in fresh biomass, and a decrease in soil organic matter, green cover, and litter cover, respectively, to some extent. However, other factors, including water content and global warming, would influence the variation of the temporal trends of Landsat reflectance and LST. The vegetation communities in GNP have similar temporal trends of NIR, SWIR1, SWIR2, and LST, as those in the surrounding pastures (Supplementary, Figures S2-S5); however, the different change rates still reveal the impact of grazing cessation and reintroduction. Therefore, the differences in Landsat reflectance and vegetation indices are good options with which to analyze the influences of prairie management only. LST temporal trends were similar for the three native vegetation communities in GNP and the surrounding pastures (Supplementary, Figure S5), which suggests that their dynamics may be dominated by climate change. The air temperature recorded at the closest meteorological station did not show an annual mean temperature increase from 1985 to 2020, however, the monthly mean temperature for prairie growing season (i.e., June to August) increased by 0.46°C per year from 1985 to 1999 and stabilized from 2000 to 2020. The continuous increase in LST was consistent with air temperature. closest meteorological station did not show an annual mean temperature increase from 1985 to 2020, however, the monthly mean temperature for prairie growing season (i.e., June to August) increased by 0.46 ℃ per year from 1985 to 1999 and stabilized from 2000 to 2020. The continuous increase in LST was consistent with air temperature.

Conclusions
Grazing cessation has been suggested as a superior prairie management policy to long-term moderate grazing for enhancing soil organic matter and increasing green cover. Among all three native vegetation communities, grazing cessation had a stronger impact on valley grasslands and the mildest effects on sloped grasslands. The significantly effective time period of grazing cessation to occur differed among biophysical indicators and vegetation communities in GNP, which were 11 to 14 years. Grazing cessation increased soil organic matter for 14 and 13 years for upland and sloped grasslands, respectively; increased green cover for 13 years for both upland and sloped grasslands; decreased fresh biomass for 14, 12, and 11 years for upland, sloped, and valley grasslands, respectively. Figure 5. The temporal trends of normalized difference vegetation index (NDVI) and normalized difference water index (NDWI): (a-c) are the temporal trends of NDVI for upland, sloped, and valley grasslands in Grasslands National Park (GNP), respectively; (d-f) are the temporal trends of NDVI for upland, sloped and valley grasslands in the surrounding pastures, respectively; (g-i) are the temporal trends of NDWI for upland, sloped, and valley grasslands in GNP, respectively; (j-l) is the temporal trend of NDWI for upland, sloped and valley grasslands in the surrounding pastures, respectively.

Conclusions
Grazing cessation has been suggested as a superior prairie management policy to long-term moderate grazing for enhancing soil organic matter and increasing green cover. Among all three native vegetation communities, grazing cessation had a stronger impact on valley grasslands and the mildest effects on sloped grasslands. The significantly effective time period of grazing cessation to occur differed among biophysical indicators and vegetation communities in GNP, which were 11 to 14 years. Grazing cessation increased soil organic matter for 14 and 13 years for upland and sloped grasslands, respectively; increased green cover for 13 years for both upland and sloped grasslands; decreased fresh biomass for 14, 12, and 11 years for upland, sloped, and valley grasslands, respectively.
Grazing reintroduction affected the prairie ecosystem with a shorter effective time period than grazing cessation for both upland and valley grasslands. The intensity of grazing reintroduction was strongest in sloped grasslands and mildest in valley grasslands. Grazing reintroduction had weaker influences on upland and valley grasslands, but stronger effects on sloped grassland than on grazing cessation. The effective time period of grazing reintroduction was 7 to 9 years for upland and valley grasslands, while there was no clear time lag effect in sloped grassland.
Generally, the increase in NDVI and NDWI is an indicator of restoration outcome for grazing cessation. However, the opposite temporal trends of NDVI and NDWI were observed in this study when GNP was under grazing cessation. This is because the continuously accumulated non-photosynthetic materials highly reduced the reflectance in NIR, an important band for calculating NDVI and NDWI.
The results of our study provide fundamental information for prairie mangers on the effects of grazing cessation and reintroduction compared to long-term moderate grazing management. More importantly, the results of this study suggest the time periods for management to take effect among grassland communities.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/rs13173397/s1, Figure S1: Selection of three vegetation types in the surrounding pastures; Figure S2: NIR reflectance of upland, sloped, and valley grasslands from both GNP and the surrounding pastures; Figure S3: SWIR1 reflectance of upland, sloped, and valley grasslands from both GNP and the surrounding pastures; Figure S4: SWIR2 reflectance of upland, sloped, and valley grasslands from both GNP and the surrounding pastures; Figure S5: LST of upland, sloped, and valley grasslands from both GNP and the surrounding pastures.