Study Area

The Great Plains region is composed mainly of grasslands (≈36% of the area), intermixed with cropland (about 40% of the area) in the central United States between eastern and western forests (see Fig. 1) (see table 1, Augustine et al. this issue for breakdown of land cover categories). Since the focus of our work was an assessment of grassland production trends, we only focused on grassland vegetation and necessarily removed all other land cover types (see Fig. 1). The major grassland types in the study area include northern mixed-grass prairie, southern mixed-grass prairie, tallgrass prairie, and shortgrass steppe identified using the Landfire Existing Vegetation Type (EVT) version 1.3 (USDA and US Department of Interior, Washington, DC); Comer et al. (2003). The extent of these vegetation types used as an analysis mask generally reflects the major grassland divisions shown in Lauenroth et al. (1999) as part of the coordinated effort in this special issue.

Table 1

Validation characteristics of the annual net primary production (ANPP) dataset against 3 areas where observations have been made within the Great Plains.

The climate of the Great Plains exhibits wide ranges in temperature and precipitation. A temperature gradient occurs from a 4°C mean in the north to 21°C mean in the south, and an increasing precipitation gradient occurs from 20 cm/yr in the west to 110 cm in the east (Shafer et al. 2014; Ojima et al. this issue). Large annual, seasonal, and even intraseasonal variation in both temperature and precipitation is common in the Great Plains. Historical vegetation consisted of shortgrass prairie in the semiarid west, mixed-grass prairie in the central region (both north and south), and tallgrass prairie in the mesic eastern region (Küchler 1964).

Remotely Sensed Data

The remote sensing data are from the Thematic Mapper (TM) archive found on the Google Earth Engine. To obtain data from 1984 to 2017, we accessed TM (Landsat 5), Enhanced Thematic Mapper (ETM) (Landsat 7), and the Operational Land Imager (OLI) from Landsat 8. The TM data were preferentially used from 1984 to 2011, while ETM data were used only when no cloud free data could be obtained from the TM data. The OLI data were used from 2013 to 2017. Henceforth, we refer to data from these three satellite sensors as “TM data,” recognizing that they represent three different sensors. The TM data are offered at a nominal spatial resolution of 30 m in the spectral channels used in this study (the red and near infrared bands). Before 1999, the repeat frequency, or the time in between successive images of the same area, was 16 d. However, since 1999, data from at least one of the Landsat platforms have been available every 8 d.

Table 2

Validation study site data characteristics. Methods for collection of data used for validating annual net primary production estimates from three sites in various major grassland types.

These data were converted to the commonly used NDVI. The NDVI is formulated as:

where red and NIR stand for the spectral responses acquired in the red (visible) and near-infrared regions (band passes), respectively. This ratio is probably the most widely used “vegetation index” since satellite remote sensing began. The annual maximum NDVI for each year from 1984 to 2017 formed the basis for estimating rangeland ANPP of the conterminous United States. The data used to derive NDVI had 3 preprocessing steps applied. The flow for the preprocessing generally follows that of Kennedy et al. (2018). First, the data used are considered Tier 1 (a description can be found here:  https://www.usgs.gov/land-resources/nli/landsat/landsat-collection-1?qt-science_support_page_related_con=1), representing the highest quality data available in the Landsat data series. These data are processed to terrain-corrected surface reflectance. Second, these Tier 1 data are further screened for clouds, snow, and shadows using the CFMask algorithm (Zhu and Woodcock 2012) version 2.0 (US Geological Survey, Washington, DC). The CFmask process is often used by algorithms or programs seeking to produce consistent time-series data such as the Analysis Ready Data (Dwyer et al. 2018) and LandTrendr (Kennedy et al. 2018). Data for 2012 are missing from our Landsat time series because Landsat 5 become inoperable in 2011, leaving only Landsat 7 data, which are full of missing data (stripes) due to the Scan Line Corrector problem for Landsat 7 that is widely known throughout the remote sensing community (Chen et al. 2011). The 2012 data were filled using the Version 6 Mod13Q1 NDVI from the Moderate Resolution Imaging Spectroradiometer (MODIS), which represents a 16-d maximum value composite product at 250-m resolution. For these data we allowed only those pixels with quality control flags representing the highest quality annual maximum NDVI values. Next, the MODIS NDVI data were calibrated to those of Landsat 8 (OLI) by comparing NDVI values from 2013 and 2014 to the MODIS data from the same time in a similar manner as Ke et al. (2015). The valid range of the MODIS-based NDVI product is –2 000 to 10 000 (Didan et al. 2015), while the range for NDVI derived from Landsat 8, based on surface reflectance, is –1 to +1. To calibrate the MODIS NDVI values in 2012 to be very like those from Landsat 8 (OLI), we evaluated the response across 110 vegetation types for each sensor for 2013 and 2014. In addition, since the MODIS NDVI has a larger pixel size (here 250 m on a side) compared with the TM suite (here 30 m on a side), the OLI pixels were resampled to match the MODIS pixel size using a cubic convolution (e.g., OLI pixels were resampled to 250 m on a side). In addition to these preprocessing steps, intersensor calibration was performed to ensure consistency between the 3 Landsat sensors (each representing different eras) including Landsat 5, 7, and 8. The sensor characteristics of Landsat 5 and 7 are similar enough to not need intercalibration, but calibration between Landsat 7 and 8 is required given the differences in the spectral channels across which data are recorded. The OLI has spectral channels within the red and near-infrared wavelengths that are distinct from previous TM sensors, and oftentimes a crosswalk between the two sensors' data is performed. Correspondingly, we converted the OLI NDVI data to those in the ETM+ using the coefficients developed by Roy et al. (2016). This transformation is given as follows:

where OLI represents the converted NDVI from the OLI sensor based on NDVI from the ETM+ sensor derived using surface reflectance data.

Figure 1.

The location of the three spatial scales (major grassland divisions, Ecological Subsections, and allotments within the National Grasslands managed by the US Department of Agriculture, Forest Service) used to conduct the assessment of production trends on the Great Plains study area. The grassland subdivisions come from the Ecological Systems represented in the Landfire EVT version 1.3 (Comer et al. 2003). In addition, the locations of the Central Plains Experimental Range (CPER), High Plains Grasslands Research Station (HPGRS), and Konza field study areas are shown.

Once the sensor intercalibration was complete, an analysis mask representing the major grassland types (see Fig. 1) was applied. From this spatial subset the annual maximum NDVI value was chosen at each pixel to be retained for further analysis. Choosing the maximum NDVI value enables some advantages, but also disadvantages, as other types of processing, such as integrating the NDVI values of the growing season. The advantages are threefold. First, this is a simpler approach that requires no a priori assumptions of growing season length over which to integrate the NDVI, which is often done (e.g., Gaffney et al. 2018). Second, when an integration is applied, this usually requires a smoothing algorithm to decrease noise and enable a curve to be fit between missing data points, which creates synthetic data that may or may not be realistic to retrieve the maximum value from. Third, these curves, depending on how they are they fit, may create conditions that don't exist. Although using the annual maximum values has benefits, there are also tradeoffs with this approach. The biggest tradeoff can exist when two distinct growing seasons are found (so-called bimodal systems). In most cases, for our study area, we did not find significant production occurring at two times in the year, especially at more northern latitudes. Instead, most of the production occurred during one of the two growing season events. Another factor potentially confounding our NDVI ANPP relationships is that shrubs can have different phenological patterns than herbaceous species. Other shortcomings of relying on maximum NDVI values are that cloud shadows can inflate NDVI values and masking of cloud shadows is challenging. Here we used the method of Housman et al. (2018) to identify and remove cloud shadows. In addition, no correction for the viewing and illumination geometries using a Bidirectional Reflectance Distribution Function (BRDF) was made, which may also affect results presented here (Roy et al. 2016).

Calibrating NDVI to Annual Production

To calibrate maximum NDVI values to ANPP, four steps were employed. First, all NDVI data were spatially subset to the extent of conterminous US rangelands representing about 662 million acres using the spatial data from Reeves and Mitchell (2011). Second, Ecological Sites were spatially represented, where they exist, using the Soil Survey Geographic Database (SSURGO 2018). The annual production data associated with Ecological Sites and stored in the SSURGO database were collected by hundreds of people by destructive harvest or double sampling (USDA, NRCS 2017; Ken Spaeth, personal communication, 2018). The general rule for collecting annual production on each site was to collect 30–50 quadrats of varying size from 0.18 to 0.89 m² (depending on vegetation type) for each site being described (Ken Spaeth, personal communication, 2018). These site locations were uniform in vegetation, soils, and landform and large enough to include the complete plot (USDA, NRCS 2017). Because of the long history and large geographical extent over which the ANPP estimates have been obtained, exact dates for all acquisitions cannot be known (Ken Spaeth, personal communication, 2018). Ecological sites generally have ecological site descriptions, which, among other things, contain information on average, above-average, and below-average values of ANPP. These data represent an unprecedented dataset for calibrating remotely sensed data in rangeland environments because of their comprehensive coverage, and the wide range of ANPP is useful for calibrating the full range of TM NDVI values found across the extent of conterminous US rangelands. The ANPP estimates from ecological sites were spatially aggregated to biophysical settings (BpS) from the Landfire Project (Rollins 2009; Landfire 2018) such that for each BpS evaluated, there were three data points representing the below-average, average, and above-average ANPP estimates. Third, the TM maximum NDVI values were also aggregated to these same BpS classes across the extent of US rangelands enabling direct comparisons with the ANPP data. The mean, minimum, and maximum annual NDVI values were compared with the average, below-average, and above-average ANPP observations for each associated ecological site yielding 3 points for each of 110 vegetation types evaluated across the extent of the study area. From these overlays, nonlinear regression models were developed. Fourth, these models were applied to NDVI data in each yr from 1984 to 2017 across the extent of US rangelands and clipped to the boundaries of the Great Plains (see Fig. 1).

Validating Annual Production Estimates

The validation process comparing the predicted ANPP with observations from three expressions of North American grasslands (shortgrass steppe, northern mixed-grass prairie, and tallgrass prairie) was based on four attributes including r² (coefficient of determination), bias, mean absolute error (MAE), and the years of minimum and maximum observed and predicted ANPP values (Table 1). The MAE was computed as the absolute value of the residuals from the observed minus the predicted production values, while the bias was computed as the summation of residuals/sample size (n).

Three different datasets representing ANPP observations in the study area were used to quantify error rates and identify where ANPP estimates may be less reliable (Table 2). The first dataset was from a shortgrass steppe ecosystem at the Central Plains Experimental Range (CPER) near Nunn, Colorado administered by the USDA–Agricultural Research Service (ARS). The data used from the CPER represent the yr between 2000 and 2017, and a description of methods and study sites for these data can be found in Irisarri et al. (2016). The second validation dataset came from the USDA-ARS High Plains Grassland Research Station (HPGRS) in the northern mixed-grass prairie near Cheyenne, Wyoming. The data used from the HPGRS represent the yr between 2001 and 2016, and a description of the study sites and methods can be found in Derner and Hart (2007). The third dataset came from the Konza Prairie Long-Term Ecological Research site, tallgrass prairie region, representing yr from 1984 to 2015. These data were collected across a series of watersheds representing different land treatments, such as burning and grazing. A description of these methods and study sites can be found at Konza (2018) and in Table 2. For each year in the assessment, the ANPP values were averaged across treatments since this assessment focuses on ANPP, not the differences in ANPP due to treatments. The net result of the methods for all three datasets is that each one produced direct observations of ANPP for varying temporal periods. For validation we compared these observations of ANPP to our predictions enabling a MAE and bias while identifying the highs and lows and interannual variation in both the observed and predicted time series.

In addition to field validation, we also conducted a validation by holding out 93 points from the calibration dataset of production values derived via ecological sites using the SSURGO soils database. From the comparison of ANPP observations from ecological sites with ANPP estimates from the RPMS, the mean absolute error (MAE), bias, and coefficient of determination (r²) were calculated from this holdout dataset.

Conducting the Assessment

The main goal of the assessment was to characterize the ANPP response of rangelands throughout the Great Plains from 1984 to 2017 through analysis of the RPMS across three spatial scales. The three spatial scales include major grassland types, ecological subsections, and management units of the National Grasslands (see Fig. 1). The major grassland types evaluated in the assessment included northern mixed-grass prairie, southern mixed-grass prairie, tallgrass prairie, and shortgrass steppe. We also portrayed trends at the subregional level within each grassland type with ecological subsections, which are part of the USDA Forest Service's National Hierarchical Framework of Ecological Units (Cleland et al. 2007). The allotments are units within National Grasslands managed by the USDA Forest Service, and the latitudinal gradient separating these grassland units enables comparison of more site-specific climatic events, such as droughts that occurred in 1984, 2002, 2011, and 2015. In addition, the allotments were used to demonstrate the use of the RPMS at local scales and to provide assessment results that are relevant to managers and related stakeholders in the region.

ANPP response was quantified at each spatial scale as correlation with respect to time, interannual variation, year of minimum production, and spatial variation. The correlation was computed as Pearson's r. The interannual variation was calculated as the coefficient of variation about the mean for each of the three spatial scales. To illustrate, consider the shortgrass steppe major grassland type. To compute the coefficient of variation for the shortgrass steppe, the average and standard deviations of ANPP were computed across all years within that spatial domain from 1984 to 2017, thereby permitting an assessment of variation though time. Similarly, the spatial variation was estimated as the coefficient of variation about the mean across the spatial domain of each spatial scale. For example, for each allotment, for each year, the mean and standard deviation of ANPP was quantified, enabling the calculation of coefficient of variation about the mean. This represents the spatial variation or heterogeneity across each analysis unit. However, the estimate of spatial variation may be muted due to positive spatial autocorrelation in the residuals. Similarly, the year of minimum production was quantified from the time series of ANPP for each of the three spatial scales.

In addition to these metrics, linear regression was used to estimate the general rate of change and the significance of trends for each of the spatial scales. The Yule-Walker approach (Gallant and Goebel 1976), which accounts for temporal autocorrelation, was used to regress ANPP and seasonal precipitation against time, which also allowed determination of significance for each trend. The Yule-Walker method augments the standard least-squares regression model with a first-order autoregressive term for the error. The resulting coefficient estimates are statistically efficient and correct standard errors and significance tests are generated. Analyses were conducted using SAS PROC AUTOREG in version 9.4 of the SAS System for Windows (© 2018 SAS Institute Inc., Cary, NC).

Precipitation Data

We compared seasonal precipitation trends from 1984 to 2017 with ANPP trends as a means of explaining potential causality of ANPP trends. To prepare and analyze precipitation data, monthly records were obtained from the Parameter-elevation Regressions on Independent Slopes Model (PRISM) project (Daly et al. 2001). Gridded monthly precipitation records were aggregated to represent the growing season (April, May, June, July, August, and September). Trends in these seasonal data were evaluated across the major grassland divisions. The Yule-Walker regression approach was used to detect significance in seasonal trends of precipitation.

Evaluating ANPP Trends on the Great Plains Across Three Spatial Scales

Major Grassland Divisions

Using the RPMS, we found significant (P ≤ 0.05) positive trends in ANPP across all the major grassland types since 1984 (Fig. 2, Table 3). The greatest annual rate of increase was observed in the northern mixed-grass prairie, estimated at about 20 kg ha^–1 with a correlation of 0.62 (P < 0.0001; Pearson's r) with respect to time. This is followed closely by southern mixed-grass prairie with a rate of increase of about 14 kg ha^–1 and a correlation of 0.59 (P < 0.0001). The shortgrass steppe and tallgrass prairie types exhibited increases of about 8 and 11 kg ha^–1 and correlations with respect to time of 0.54 (P = 0.0234) and 0.47 (P = 0.0035), respectively. Note the significant drop in production in the shortgrass steppe during the 2012 season (see Fig. 2). 2012 also exhibited uncharacteristically early phenology, or a “false spring,” which left vegetation susceptible to late frost (Ault et al. 2013) and may have aided the desiccation of grasslands. Excluding this year of significant reduction as an outlier, the correlation of production with time is estimated at 0.7.

In addition to the trend of increasing ANPP with time, all grassland types seem to exhibit stable to increasing heterogeneity, over the study period, except for the tallgrass prairie, which seems to be increasing in homogeneity (decreasing spatial variation) through time (see Fig. 2). The interannual variation was greatest on the northern mixed-grass prairie at about 21%. The shortgrass steppe exhibited average variation of 18%, followed by southern mixed-grass prairie and tallgrass prairie with 14% and 8%, respectively.

Ecological Subsections

Results for Ecological Subsections reveal subregional differences and indicate that the northern-most areas have experienced considerable increases in production through time, especially in Montana and North Dakota (Fig. 3). This northern part of the region exhibited positive trends with respect to time (correlation ranging from 0.4 to 0.7). In the southeastern part of the Great Plains, production exhibited a lower but still positive correlation with respect to time (correlation ranging from 0.45 to 0.65). Areas with lowest correlations include parts of New Mexico and Texas and the far eastern extent of the study area (correlations ranging from 0.01 to 0.35). Average coefficients of variation in ANPP (amount of interannual variation) range from 6% to 37%, see Fig. 3, B).

The year of the lowest estimated ANPP was likely due to drought (see Fig. 3, C). Roughly 40% of the study area experienced minima during the droughts of 1988 and 1989. Another notably low period of production occurred in 1984 and 1985 across central Texas and the Black Hills region of South Dakota. In 2002 much of eastern Colorado and surrounding lands exhibit a local minimum. The droughts of 2011 and 2012 produced the lowest recorded ANPP since 1984 across much of the southern plains, especially northern Texas and western Oklahoma.

The amount of spatial variation across each ecological subsection generally decreases from west to east (see Fig. 3, D). Regions such as the Montana, North Dakota, and the Sand Hills of Nebraska exhibit the greatest spatial variation, while the easternmost extents tend to be more homogeneous across the landscape relative to more xeric sites to the west.

National Grasslands

Within the grasslands, 1 774 grazing allotments were evaluated but, because many were of small size (i.e., < 5 pixels), only 1 572 allotments were retained. We found three noteworthy issues regarding ANPP on grazing allotments administered by the US Forest Service in the National Grasslands of the Great Plains. First, like the results for the major grassland types, the grazing allotments have generally exhibited increasing trends of ANPP since 1984. Most of the allotments with the highest correlation of production with respect to time (r ≥ 0.7) are in western North Dakota on the Little Missouri National Grasslands (Fig. 4, A and C). About 8% of the allotments exhibited a correlation of 0.6 or more of production with respect to time across the extent of the Great Plains. Only 10 allotments exhibited decreasing trends, while most were slightly increasing to flat trajectories.

Figure 2.

Time series and correlation coefficients (Pearson's r) with annual production for each of the major grassland divisions within the Great Plains study area including tallgrass prairie, shortgrass steppe, northern mixed-grass prairie, and southern mixed-grass prairie. The error bars represent the standard deviation about the mean derived across the associated grassland type in each year. In this way, the error bars represent the spatial variability in annual net primary productivity. Also shown are the trends in growing season precipitation across the domain of the major grassland types.

Our second notable finding was that the allotments with the highest interannual variation are located in western South Dakota, south of the Black Hills area (Buffalo Gap National Grasslands), and northwestern Nebraska (Oglala National Grasslands (see Fig. 4, B) and the grasslands of north Central Colorado (Pawnee National Grasslands). Allotments in these areas exhibit variation of > 40% on an interannual basis and, across all the allotments in the National Grassland system, about 38% exhibit variation > 30%.

Table 3

Results of Yule-Walker regression to evaluate trends in growing season precipitation and annual net primary production (ANPP) for all 4 major grassland types.

Our third main finding was related to allotments in the southern Great Plains, particularly on the Kiowa Rita-Blanca National Grasslands in Texas and Oklahoma, which exhibited extraordinary changes in production after the 2011 and 2012 drought period. The trajectory (see Fig. 4, D) indicates that production increased about 3.5-fold within just 5 yr after the drought. To add context to the trends revealed through assessment of ANPP across these three spatial scales, we quantified growing season precipitation from 1984 to 2017.

Precipitation

Across all major grassland types, except the shortgrass steppe, growing season precipitation has significantly (P ≤ 0.1) increased since 1984 (Table 3). The tallgrass prairie exhibited the greatest increase in growing season precipitation through time of 3.1 mm per season with a correlation with respect to time of 0.32 (P = 0.081). The southern and northern mixed-grass prairies exhibited similar significant increases of 2.2 and 2.1 mm per season with correlations of 0.33 (P = 0.098) and 0.41 (P = 0.012), respectively. The trend in growing season precipitation in the shortgrass steppe since 1984 was 1 mm per season, but this was not significant (P > 0.1). The interannual variation, expressed as the coefficient of variation of growing season precipitation since 1984, was greatest on the northern mixed-grass prairie and shortgrass steppe, both exhibiting a coefficient of variation of about 20%. The southern mixed-grass and tallgrass prairies both exhibited a coefficient of variation of 17%.

Figure 3.

An analysis of trend with respect to time (A, quantified with Pearson's r), interannual variability (B, variability with respect to the mean), year of minimum ANPP (C), and spatial variability of annual net primary productivity (D) in the ecological subsections (Cleland et al. 2007).

Figure 4.

An analysis of trend (quantified with Pearson's r) and variability (coefficient of variability, with respect to the mean) in the National Grasslands within study area. A, Trend in production with respect to time. B, Interannual variability. C, The temporal trajectory of the area in the Little Missouri National Grassland within the box. D, Variability and 3.5-fold increase in annual net primary productivity after the drought of 2011 and 2012 is shown for the Kiowa-Rita Blanca National Grassland.

Figure 5.

Relationship between annual maximum normalized difference vegetation index (NDVI) and annual production from ecological site descriptions used in this study from across the western United States. There are two ranges or relationships defined here, one for a lower end of production and one for the higher end. The upper end range begins at an NDVI value of 0.46.

Calibrating NDVI to Annual Production

The relationships between ANPP estimates and maximum NDVI were divided into two groups to enable different models to be fit to the lower and upper ends of production given as:

where y is the estimated ANPP in kg ha^–1 and X is the NDVI for the upper range (X ≥ 0.46) and

where y is the estimated ANPP in kg ha^–1 and X is the NDVI for the lower range (X < 0.46). The division into 2 sections was done, in part, because of the asymptotic nature or “saturation” feature (Santin-Janin et al. 2009) of NDVI with respect to ANPP (Fig. 5).

The lower range is calibrated such that negligible production is estimated below NDVI values of 0.13, which is consistent with the observations in Huete (1988) where the upper range of NDVI for bare ground, depending on local factors, should be about 0.1.

The comparison of MODIS-based NDVI to NDVI from the OLI across 110 vegetation types over 2 yr yielded the following transformation:

where NDVI_new is the MODIS-based NDVI in 2012 calibrated to OLI NDVI and NDVI_MODIS is the MODIS-based NDVI value in 2012. This linear formulation resulted in an r² of 0.93, which is near that reported by Ke et al. (2015) where r² values averaged about 0.94 when comparing OLI and MODIS NDVI. This conversion for data in 2012 completed an intercalibrated NDVI dataset from 1984 to 2018 that was ready to be calibrated to annual production.

Validation

Overall validation results indicate that the amount of variation in ANPP explained (r²) by NDVI was 89% (Fig. 6) with the following fitted equation:

where ANPP_pred is the predicted ANPP and ANPP_obs is the observed ANPP.

Within each individual validation comparison, results were considerably better for the shortgrass steppe and northern mixed-grass prairie (CPER, HPGRS) than the tallgrass prairie site (Konza) (see Fig. 6). Of these, the HPGRS exhibited the greatest amount of variation explained at 90%. For the CPER, the amount of variation explained was 79% while at Konza, the variation explained was the lowest at 30%. At the CPER and HPGRS, the bias or amount of overprediction or underprediction was 278 and 308 kg ha^–1, respectively, while at Konza bias was higher at 404 kg ha^–1. Likewise, the MAE was also highest at Konza at 660 kg ha^–1, while at CPER and HPGRS, MAE values were 280 and 308 kg ha^–1, respectively. The final metric used for direct validation was to evaluate the estimated year of minimum and maximum ANPP that corresponds to the time frame of the validation data. Across the CPER and HPGRS, the estimated ANPP matched the years of minima and maxima exactly, but for Konza, the maximum yr in the observed data was 2007 while the estimated yr was 2015. In addition to directly comparing the RPMS estimates with observations of ANPP obtained on shortgrass steppe, northern mixed-grass prairie, and tallgrass prairie, validation using a holdout dataset was conducted. For the lowest values of production (< 3 200 kg ha^–1), the mean absolute error and bias are 236 and –139 kg ha^–1, respectively (Fig. 7) and for the highest values of production (≥ 3 200 kg ha^–1), the MAE and bias are 1 268 and –1 128 kg ha^–1, respectively. The coefficient of determination across the full range of production values examined in the holdout process was 0.76. The fitted equation corresponding to Figure 7 was as follows:

where CrossVal_pred is the predicted ANPP and CrossVal_obs is the observed ANPP.

Figure 6.

Validation of the annual net primary productivity estimates using the Konza, High Plains Grasslands Research Station (HPGRS) Long-Term Ecological Research (LTER) data, and the Central Plains Experimental Range (CPER) (Derner et al. 2007) as observations.

Quantifying Trends and Variation in ANPP Across the Great Plains

ANPP exhibited an increasing trend through time across each of the four major grassland types evaluated. To our knowledge, this has not been demonstrated for the entire Great Plains but has recently been documented in grasslands of central Asia. Yang et al. (2018) reported that 44% of grassland area in west central China exhibited increasing biomass between 2000 and 2016. Although increases in ANPP were found in the Great Plains, those overall averages belie subregional patterns within each type. Large areas within the Great Plains that exhibited no increasing trend in production can be adjacent to regions with strongly increasing ANPP (see Fig. 3, A).

Interannual variation can also be similarly asymmetrical across the landscape. For example, shortgrass steppe exhibited extreme interannual variation in some areas, such as northeastern Colorado, which is equal to or greater than variation seen anywhere in the grassland divisions of the Great Plains (see Fig. 3, B). Further, variation at localized areas may exceed 50% on an interannual basis. In contrast, most of the subsections within the tallgrass prairie domain tend be considerably less variable than the other grassland types (see Fig. 3, C and Table 1). Interannual variation in the Konza validation data is about 20%, which matches our assessment of interannual variability at that site using the RPMS. Overall, interannual variation decreases from west to east following the increasing precipitation gradient. Similarly, the spatial variation in production (see Fig. 3, D) indicates that the tallgrass prairie is considerably less variable from a spatial perspective as well, again reflecting more stability owed to greater moisture balance relative to drier areas to the west. A few subsections exhibit extremely high spatial variation including the Sand Hills of Nebraska, the northwestern edge of the northern mixed-grass prairie, and along the eastern flank of the southern Rocky Mountains in southern Colorado (see Fig. 3, D).

Average variation from either a spatial or temporal perspective only tells part of the story, especially from a manager's standpoint. Another way to evaluate variation is to estimate the distance between local minima (low production) and maxima (high production) across a set of years. Evaluating the rate of increase in ANPP following a drought may reveal recovery (and possibly resiliency) in each grassland. To illustrate, across the shortgrass steppe, in 2011 and 2012, the growing season precipitation was the lowest since 1984 at 174 and 182 mm, respectively, and the estimated ANPP was commensurately low (see Fig. 2). Subsequently, precipitation increased to a maximum in 2015 of 443 mm and a maximum in ANPP of 1 587 kg ha^–1. Percent changes in precipitation and ANPP were 142% and 81%, respectively. In a similar manner, the southern mixed-grass prairie also reached a minimum of precipitation in 2012 of 272 mm while reaching a maximum in 2015 of 537 mm. During that same time, estimated ANPP also reached a maximum in 2015 representing an increase of 28% from 2012. In 2011, the northern mixed-grass prairie experienced the highest amount of growing season precipitation (414 mm) since 1984, followed by a 50% decrease in 2012 producing the fifth lowest precipitation amount (273 mm). By 2014, the precipitation had increased from this minimum by 76% and ANPP rebounded commensurately by about 65%. Finally, on the tallgrass prairie, the period from 2006 and 2007 represents the single greatest interannual change in growing season precipitation since 1984. In 2006, the precipitation was the fourth lowest on record at 464 mm while 2007 exhibited 768 mm representing a 65% increase. The increase in ANPP from 2006 to 2007 was estimated at 22%. These recoveries from drought do not necessarily reflect resilience, especially since growing season precipitation is not the only determinant of yield and the relationships between precipitation and ANPP vary considerably (Petrie et al. 2018). However, if these major grassland types were not functioning reasonably well, they probably would have smaller recoveries after drawdowns of these magnitudes due to drought or other perturbation (Reeves and Baggett 2014).

Although novel, this assessment of recovery after drought, or following large reductions in ANPP, is more relevant to managers when applied to individual allotments or sites. A specific example of this capability (see Fig. 4, D) is detection of increased ANPP by 3.5-fold within just 5 yr after the notable drought that enveloped the region in 2011 and 2012. Conversely, in a nearby unit (about 2 km) on the Kiowa Rita Blanca National Grasslands, the recovery over the same period was only 40%. Not all allotments in this region experienced such dramatic increases following the 2011/2012 drought, but 62% experienced postdrought production (2013 to 2017) of at least 1.5 times the values seen during the drought. This type of resilience or rebound in central US grasslands (see Figs 2 and 3), following extreme experimental drought and heat waves, has been documented by Hoover et al. (2014).

Overall, postdrought increases in ANPP have been observed across the Great Plains, but this does not indicate that there are no lasting impacts of drought. Species change is a significant issue that should be evaluated after drought. Hoover et al. (2014) found that losses in forb abundance and production were compensated for by an increase in grass production in an extreme drought study. Following the 2011/2012 drought on the southern Great Plains in 2013 (see Fig. 3, C and Fig. 4, D), Texas experienced the lowest estimated calf crop for the period of record (USDA-NASS 2018) at 3.8 million, followed by the smallest inventory of cattle and calves since 1958 in 2014. Likewise, public land managers were also challenged during these periods of drought. For example, it can be difficult to authorize fewer head months than a permit allows because this can strain relations between managers and permittees (Safranek, personal communication, 2018). The situation can be further complicated because control of native ungulates is administered by the states, and thus the amount of available forage during a drought can be reduced and out of the control of the manager of a grassland unit (Safranek, personal communication, 2018). In other circumstances, managers may be attempting a new technique or grazing strategy with a permittee and severe reductions in forage from drought or fire can reduce the efficacy of these experiments (Safranek, personal communication, 2018).

Potential Drivers and Impacts of Increasing ANPP

All four of the major grassland types exhibited significant (P ≤ 0.05) increases in ANPP since 1984. About 35% of the grazing allotments in the National Grasslands have exhibited correlations of ≥ 0.5 in ANPP with respect to time (see Fig. 4, A). The greatest increases generally appear in the northernmost region of the Great Plains and particularly within the Little Missouri National Grasslands. We suggest that increasing growing season precipitation and woody encroachment are two potential factors influencing the ANPP patterns across the Great Plains.

Growing season precipitation has increased significantly (P ≤ 0.1) on three of four major grasslands, offering compelling evidence as a potential driver of ANPP. Despite the increasing trends in growing season precipitation, there is a considerable amount of unexplained variability in the trends of ANPP. The linkage between growing season precipitation and ANPP varies, and the amount of precipitation in the prior year can also influence ANPP in the current year (Petrie et al. 2018). Thus, it is conceivable that woody expansion, or other factors like CO₂ enrichment, are also acting to increase ANPP across the extent of Great Plains grasslands (Morgan et al. 2008).

Woody expansion in grasslands is noted globally (Knapp et al. 2008; Van Auken 2009), in Canada (Bai et al. 2009; Madden et al. 2000) and the United States (Houghton et al. 1999; Gaskin et al. this issue). Especially in the southern extent of grasslands, encroachment has been noted for a long time, in some cases since the early 1930s (Browning and Archer 2011), but the species and situations of woody components on the landscape vary considerably. In the northern mixed-grass prairie, western snowberry (Symphoricarpos occidentalis) is among those species that have been documented as increasing in density, but the scope and scale are presently unknown (Bai et al. 2009). Where western snowberry increases have been documented, concomitant increases in ANPP and phytomass are often observed, in some cases sixfold compared with non-encroached areas (Bai et al. 2009). In the southern mixed-grass prairie, especially toward the southern extent, mesquite and juniper (Prosopis and Juniperus spp.; Ansley et al. 1995; Ansley and Rasmussen 2005) are increasing in density. This is significant because ANPP increased substantially where encroachment has occurred and sites with high woody vegetation densities had ANPP values up to 4× greater than sites with only herbaceous vegetation (Hughes et al. 2006). In the tallgrass prairie, eastern redcedar (Juniperus virginiana) is a primary species responsible for encroachment. Other species, such as roughleaf dogwood (Cornus drummondii) and smooth sumac (Rhus glabra), are also having significant impacts on the tallgrass prairie region (Ratajczak et al. 2009). Stands dominated by roughleaf dogwood have ANPP values about 3× those of uninvaded adjacent tallgrass prairie sites (Lett et al. 2004). The situation for encroachment in the shortgrass steppe is less clear, and the evidence for widespread shrub encroachment on shortgrass steppe is lacking.

While this assessment does not directly link increased ANPP to increasing shrub density, the results suggest potential causality. An assessment of changes in shrub cover across rangelands (e.g., Xian et al. 2015) in tandem with ANPP data developed here will provide answers at a national scale regarding the role that shrubs are playing in the increasing ANPP signal. Reasons for increased abundance of shrubs and trees include land management history, such as reduced fire frequency and changing herbivory, CO₂ enrichment, and increasing interannual variability of precipitation (Morgan et al. 2008; Van Auken 2009). Gherardi and Sala (2015) have posited a newer theory whereby increasing interannual variation in precipitation can lead to a greater proportion of ANPP being produced by shrubs, at least under experimental conditions even if overall productivity shows little or declining response.

Regardless of the causes of encroachment, woody expansion is considered “one of the greatest contemporary threats to mesic grasslands of the central United States” (Briggs et al. 2005). Increasing shrub density has significant implications for ecosystem processes, such as a shift in biomass from belowground to aboveground. From a producer's standpoint, woody cover usually displaces herbaceous cover, thereby reducing the amount of forage. Reduced herbaceous production has direct consequences for sustainability of herbivory, which is a significant economic, social, and ecological consideration throughout the Great Plains.

Validation Considerations

The consistent ANPP time series developed here shares many of the characteristics with the observations of production, especially in the temporal aspects and in the way the ANPP estimates match the highs and lows on the landscape. Our temporal results are consistent with others such as Yang et al. (2018), who reported the coefficients of determination (r²) in grasslands of the Three Rivers Headwaters Region on the Qinghai-Tibet Plateau as ranging from 0.79 to 0.82, which are like the ranges found in the present study.

Although the RPMS is a reliable monitoring system for most US rangelands, there is clear “saturation” evident in the relationship between NDVI and ANPP. As a result, increasing NDVI does not necessarily produce commensurate increases in ANPP and reliable calculation of ANPP may not be possible when ground cover approaches 100% (Paruelo et al. 1997; Schmidt et al. 2016). For US rangelands, saturation did not become significant until ≈3 258 kg ha^–1, and most US rangelands (about 85% as revealed in the RPMS) have mean ANPP < the 3 258 kg ha^–1 threshold. In the Great Plains, about 25% of the rangelands exhibit average ANPP exceeding the threshold, particularly in tallgrass prairie. The difficulty in calibrating vegetation indices to ANPP or biomass in tallgrass is widely noted and well documented (e.g., fig. 1 in Sharma et al. 2018). The implication for the present work is that the ANPP estimates across the tallgrass prairie are less accurate than the other vegetation types over which the RPMS offers data.

Use of other indices with greater sensitivity on sites capable of producing high biomass, such as the EVI, Wide Range Dynamic Vegetation Index (Viña and Gitelson 2005), or the Generalized Difference Vegetation Index (GDVI) (Wu 2014), could be used, especially in tallgrass prairie. It has not been definitively shown, however, which of these indices work best throughout the full range of the tallgrass prairie region and this is an area of research that could be pursued in the future.