Using the TSS-RESTREND Methodology to Diagnose Post-Reclamation Vegetation Trends on the Western Slope of Colorado

The extraction of crude oil and natural gas on US public lands has experienced a significant boom within the past decades. This rapid growth poses challenges for natural resource management due to the unique nature of oil and gas development (Di Stéfano et al. 2021). In many large landscapes, the development of oil and gas has resulted in a vast network of discrete disturbances stemming from well pads, pipelines, and connecting access roads. This development has direct and indirect effects including habitat degradation, increased susceptibility to invasive plant species, and altered landscape hydrology (Chambers et al. 2022; Walker et al. 2020; Yu et al. 2015). These impacts often impair or inhibit other important ecosystem values and permitted uses such as livestock grazing, wildlife habitat and connectivity, and recreation (Allred et al. 2015; Walker et al. 2007; Waller et al. 2018).

On US public land, oil and gas development is primarily administered by the Bureau of Land Management (BLM), where the agency oversees about 99 million surface ha (245 million acres) of public land and 283 million subsurface ha (700 million acres) with ∼5 million subsurface ha (12.8 million acres) producing oil and gas in profitable quantities (Allred et al. 2015; Bureau of Land Management 2021). US federal land management agencies, such as the BLM, are required to manage for multiple natural resource uses without permanent impairment of the productivity of the land and the quality of the environment, commonly referred to as the multiple-use mandate (Haskell 1976).

Reclamation is the process of repairing highly damaged lands, such as often occurs with oil and gas development, where management objectives are often guided by policy and target conditions (often referred to as “reference” conditions; Di Stéfano et al. 2021; Gerwing et al. 2022; Lima et al. 2016). In the context of reclaiming public lands, particularly for agencies such as the BLM, the assessment of reclamation success predominantly revolves around permit compliance. This process begins with an initial evaluation of soil stabilization efforts, including tasks such as topsoil removal, respreading, and restoring the land’s approximate original contour. Subsequently, the evaluation places significant emphasis on vegetation cover and, to a lesser extent, species diversity. Together, soil stability and vegetation cover have the most significant weight within the legal framework governing reclamation for the BLM. In addition, enhancing the likelihood of reclamation success may involve comparing various aspects to a reference condition, including physical characteristics, soil dynamics, ecosystem services, vegetation dynamics, and predicted responses to management activities (Lupardus et al. 2023; Twidwell et al. 2013).

Reference condition can be determined by vegetation and soil conditions before the pad was established or through a space-for-time approach using a nearby site unimpacted by oil and gas development (Di Stéfano et al. 2021). However, the majority of oil or gas pads on public land in the United States lack predisturbance condition information, necessitating the selection of a comparable and undisturbed site for reference, which can be difficult. Furthermore, assessing outcomes following treatments using pre- and postdisturbance data can be confounded by weather (e.g., wet or dry cycles), whereas the inference from the space-for-time approach is completely dependent on selecting a well-matched site to the pad in question (Fick et al. 2021). These challenges make it difficult to interpret and manage for postreclamation vegetation trends, allowing for a wide and subjective array of evaluations on reclamation outcomes (Di Stéfano et al. 2021).

Time series remote sensing is a widely used method for monitoring landscape changes, encompassing phenomena such as oil and gas development and subsequent reclamation efforts (Waller et al. 2018), and there are several time series analysis methods specifically designed to characterize plant cover trends (Symeonakis 2022). In this regard, breaks for additive season and trend (BFAST) and residual trends analysis (RESTREND) have emerged as valuable tools. BFAST is a well-established method in remote sensing applications that detects both gradual and abrupt changes in vegetation cover over various time periods by identifying significant changes or “breaks” in time series data. It has proven effective in identifying alterations associated with land degradation, fire events, and other landscape transformations (Verbesselt et al. 2010; Wessels et al. 2012). However, BFAST has been found to be overly sensitive to climactic events (e.g., drought), making it difficult to determine vegetation responses specific to management actions (Burrell et al. 2017). RESTREND has been used to detect an overall upward or downward trend in vegetative cover, but it relies on precipitation to be a linear predictor of plant growth (Higginbottom & Symeonakis 2014). Although this method has its merits, it may not be the most suitable choice for analyzing degraded areas within arid landscapes, where plant growth responses to precipitation exhibit strongly nonlinear patterns (Wessels et al. 2012).

Vegetation indices, such as the Soil Adjusted Total Vegetation Index (SATVI), have often been used within these time series analysis methods to measure trends in total vegetative productivity (Marsett et al. 2006). The focus on total vegetative productivity is an issue, however, for well pads, where invasive weeds are a persistent postreclamation issue and may inflate estimates of vegetative productivity (Nauman et al. 2017). In contrast, the greenness-to-cover index (GCI) has been used to detect recovery on arid lands after complete removal of vegetative cover and evaluate changes in plant community composition (Sesnie et al. 2018; Villarreal et al. 2016). GCI is calculated as the normalized difference between normalized difference vegetation index (NDVI) and total vegetation fractional cover (TVFC; Eq. 1; Villarreal et al. 2016).

Here, we investigate the use of the time series segmented residual trends (TSS-RESTREND) methodology for analysis of oil and gas reclamation trends. TSS-RESTREND was developed for detecting and diagnosing land cover change in landscapes with an unstable vegetation-precipitation relationship (VPR), including degraded areas (Burrell et al. 2017). TSS-RESTREND combines BFAST and RESTREND along with the statistical Chow test to evaluate land cover change in unstable landscapes that might be missed using traditional time series analysis.

First, the VPR is quantified by a regression model of the GCI and precipitation values. A times series of the residuals from this relationship is then evaluated for possible BFAST breaks, after which the identified breaks are further evaluated with the Chow test. The Chow test has not been commonly used in remote sensing, but it has been extensively used in economics to identify temporal instability or change in time series data. In our case, we used the Chow test to determine if there were significant differences between the coefficients of the VPR residuals before and after a detected change point or BFAST break (Chow 1960). These coefficients signify the average change in the VPR residuals for each year of the time series. By integrating the Chow test with BFAST and RESTREND, the TSS-RESTREND methodology can help identify changes in plant community composition over time, enhancing our understanding of plant community dynamics in these severely degraded areas.

The general objective of this article was to determine if TSS-RESTREND could be used to detect and evaluate plant community change after reclamation and if these changes corresponded with documented management actions. More specific objectives were (a) assessing if time since reclamation was related to trends in the VPR, (b) determining if breaks or trends in the VPR could be related to management changes visible in aerial imagery or noted in records, and (c) evaluating the utility of using TSS-RESTREND for assessing reclamation outcomes. We accomplished these objectives by applying TSS-RESTREND with GCI to detect compositional plant community changes postreclamation for a set of reclaimed oil and gas well pads in northwestern Colorado, USA. We then determined possible connections between remotely sensed vegetation patterns and management actions documented in the management records of each well pad. By finding and evaluating these connections, we can better understand how management actions may be promoting or hindering recovery.

We conducted this study within the western portion of the BLM’s White River and Little Snake Field Offices in northwestern Colorado, USA (40.3°N 108.3°W). This portion of the field offices covers approximately 2 million ha, with 385,000 ha having an active or pending BLM oil or gas lease (Fig. 1; BLM Colorado State Office 2021). Mean annual precipitation at Dinosaur, CO (central location in study area) between 1981 and 2010 was 289 mm, with a mean annual snowfall of 889 mm (US Climate Data 1981). Approximately 90% of the study area is public land that is predominantly managed by the BLM, and the other 10% is privately owned. Vegetation in the study area is dominated by big sagebrush (Artemisia tridentata Nutt.) and perennial bunchgrasses such as Indian ricegrass (Achnatherum hymenoides [Roem. & Schult.] Barkworth), bluebunch wheatgrass (Pseudoroegneria spicata [Pursh] Á. Löve), and Sandberg bluegrass (Poa secunda J. Presl; Lowry et al. 2005). The soils in the study area is primarily derived from residuum and/or colluvium weathered from shale or sandstone (Soil Survey Staff 2021). The most prevalent land uses of the study area are energy extraction (due to the underlying oil shale deposits; Taylor 1987) and livestock grazing. The BLM manages most oil and gas permits in the region, which includes reviewing drilling applications, monitoring compliance with extraction regulations, and evaluating the success of well pad reclamation.

To detect changes in plant community over time and differences in reclamation management, we used a stratified random sample to select 40 of the total 979 reclaimed pads in the BLM Northwest District Office of Colorado. First, we limited our study to classes of Southwest Regional Gap Analysis Project (SWReGAP) land cover where big sagebrush was the dominant vegetation, to reduce environmental variability and because of the habitat needs of sagebrush obligate species such as the greater sage-grouse (Centrocercus urophasianus). The specific SWREGAP classes were Inter-Mountain Basins Big Sagebrush Shrubland (n = 38 of the final selected sites) and Inter-Mountain Basins Montane Sagebrush Steppe (n = 2; Lowry et al. 2005). The pads were then classified into three different aridity index classes based on quantiles of the 25th, 50th, and 75th percentiles (<0.1777, “dry”; 0.1777 to 0.2248, “semi-dry”; and >0.2248, “wet”; Trabucco & Zomer 2019). The pads were then binned within four different times since the reported abandonment date (i.e., 5, 10, 15, and 20 years; Table 1). The aridity index is a scale ranging from 0 to 1 that describes moisture availability for potential growth of vegetation (Zomer et al. 2022). Data set values for the aridity index are available for the years 1970 to 2000 (Trabucco & Zomer 2019). We treated the reported abandonment dates as reclamation dates because final reclamation typically occurs within 6 months of that date. Of the total 979 reclaimed pads within the study area, 46.4% had big sagebrush SWReGAP cover classes, of which we selected 40 using the sample_frac function in the R package dplyr (Wickham et al. 2023).

To implement the TSS-RESTREND methodology on our 40 selected pads, we used two spatial data sets: (a) a GCI index derived from Landsat 5 and 7 imagery at a 16-day interval during the growing season (March-September) using Google Earth Engine (Gorelick et al. 2017) and (b) daily precipitation values acquired from the GridMET data set (Abatzoglou 2013). Both data sets were obtained for the growing seasons of each year from 1984 to 2020 and used in their native resolution (Landsat imagery at 30 m; GridMET at 4 km). Spatial data were projected to Universal Transverse Mercator zone 13 North (UTM 13N).

To calculate the relationship between the vegetation index (GCI) and the local precipitation (VPR), we first performed an ordinary least squares regression on the max GCI for each growing season and the optimal accumulated precipitation. The optimal precipitation was calculated by finding which combination of accumulated period (1–12 months) and offset period (0-3 months) had the highest correlation coefficients with the GCI max. The difference between the observed GCI and GCI predicted from the optimal VPR is then called the VPR residuals (Fig. 2). These resulting residuals from the relationship are then assumed to be the GCI values detrended from annual climatic fluctuations, thereby isolating impacts of management actions or disturbances on the GCI.

To look at the effect of time since reclamation on vegetation, we applied a linear mixed-effect model using the lmer function in the lme4 package in R version 4.2.3 (Bates et al. 2015). The evaluated response variable was the time series of annual VPR residual values after reclamation of each well pad. The fixed effects were years since the reclamation and the mean annual aridity index value for 1970 to 2000, with the well pad identifier used as a random effect.

We evaluated the VPR residuals for breaks during 1984 to 2020 through BFAST modeling and then analyzed if the resulting breaks indicated instability in plant community composition using the Chow test (Fig. 2). The Chow test was performed with the CHOW function in the TSS-RESTREND package in R version 4.0.2 (Burrell 2020; Chow 1960). If no instability in plant community composition was detected by the CHOW test, we performed a single RESTREND for the entire time series of the well pad’s VPR residuals (Fig. 2). When the Chow test determined a breakpoint significant (found instability in plant community composition; p < .05), we then recalculated the VPR before and after the breakpoint. Each segment of the resulting VPR residual time series were then analyzed using RESTREND (Fig. 2). This divided analysis is called segmented RESTREND. All methods were applied using the TSS-RESTREND package in R version 4.0.2 (Burrell 2020).

To better understand how management actions on the pads relate to TSS-RESTREND outcomes, we used aerial imagery available in Google Earth (GE) up to July 2016 (Google Earth 2016) and well pad documents maintained for the public by the Colorado Oil and Gas Conservation Commission (COGCC; Colorado Department of Natural Resources 2021). Specifically, we compared break points and trends identified by TSS-RESTREND with dates of well pad abandonment, dates of vegetation reseeding, documented disputes between operators and BLM, appearance of equipment on site beyond the reported abandonment date, reported management actions, and other management or activity observed in GE imagery. We also noted if the reclaimed site was adjacent to other disturbed areas. We then grouped the well pads into six categories based on our observations of the relationship between TSS-RESTREND results and GE imagery/COGCC documents: (a) high surrounding disturbance, (b) well pad not developed, (c) inaccurate abandonment date, (d) apparent reclamation failure, (e) reclamation site redisturbed, and (f) apparent successful reclamation.

We observed a weak overall model fit between the VPR residuals and the fixed effects of time since reclamation (R² = .11) and mean aridity index values (R² = .13). We did find a slightly positive increase in VPR residuals with longer time since reclamation completion, indicating that vegetation growth response becomes more positive over time since reclamation (β = .10, p = .016). We also noticed a small positive impact on VPR residuals associated with decreases in average aridity index values, suggesting that vegetation growth response to precipitation is greater on more arid well pads (β = .13, p = .001). When considering the interaction between time since reclamation and aridity index, we found that the impact of the mean aridity index value decreases as time since reclamation increases (p = .041; Fig. 3). This indicates that the aridity of a well pad has a diminishing impact on vegetation growth with increasing time since reclamation, although the effect size remains minor (β = −.11).

The overall change in the well pads’ VPR residuals appeared variable across the study area with no visible spatial trends (Fig. 1), but we did see a positive increase in VPR residual change with an increasing reclamation age of a well pad (Fig. 4). VPR residual change is minimal for well pads within 5 to 10 years since reclamation, but there is a more positive change for well pads within 15 to 20 years since reclamation (Fig. 4).

The impact of management actions was apparent in the trends of the VPR residuals (Table 2, Fig. 5). Of the 40 well pads analyzed, we found only 5 whose times series of VPR residuals qualified for segmented RESTREND analysis, meaning there was instability in plant community composition as detected by the Chow test. In the GE imagery, we observed that all five well pads occurred in areas of high surrounding disturbance, meaning there were multiple other well pads and access roads within 30 m or less of the well pad surveyed (Fig. 6). For the other 35 well pads, trends in their VPR residuals could often be categorized by events observed in GE imagery or recorded in COGCC paperwork. Of those 35 well pads, 15 had no BFAST breaks and 20 had at least one BFAST break.

We found that well pads with either one or no BFAST breaks in their time series and minimal change in VPR residuals over time either did not have the plant community removed (not developed) or had persistent bare ground with little vegetation during the study period (n = 10; Fig. 7). Where there were multiple BFAST breaks (two or more) and higher variability in VPR residuals, BFAST breaks often corresponded with dates of interim reclamation and final reclamation observed in GE imagery but not always recorded in COGCC paperwork (n = 4; Fig. 8). For example, we detected from GE imagery that vegetation established after interim reclamation was usually later removed during final reclamation, which resulted in a less predictable VPR residual trend postreclamation. Within this unpredictability, some well pads showed an overall positive trend in VPR residuals while others were negative, suggesting an inconsistent vegetation response to precipitation (Fig. 5, 8) that we linked to management actions, as evidenced in the GE imagery and pad documentation.

We found that consistently positive trends in well pad VPR residuals were indicative of an increased vegetation response to precipitation through time, which may imply ongoing plant community recovery postreclamation. The situations resulting in consistently positive trends in VPR residuals (>0) were where (a) an inaccurate abandonment date was recorded in the COGCC database (n = 2) and (b) final reclamation was performed and vegetation appeared to be successfully reestablished (n = 7; Figs. 8, 9). For those well pads with inaccurate abandonment dates, the correct date was identified by BFAST and corroborated by the GE imagery (Fig. 9). For one pad, this included equipment being visible 4 years after the reported abandonment date (Fig. 9c) but not visible after the BFAST identified break (Fig. 9d). In instances of final reclamation being performed, BFAST breaks corresponded with well pad establishment and abandonment, and the reclamation plant community appeared to be reestablished and blend in with the surrounding area in the GE imagery (Fig. 10).

The diminishing effect of aridity on reclamation outcomes over time suggests that water availability has a reduced impact on long-term reclamation success. This result highlights the potential for reclamation efforts to overcome initial arid conditions and implies that, with time, plant communities can become increasingly self-sustaining, even under water-limited conditions (Shackelford et al. 2021). In addition, the gradual positive increase in plant community response over time indicates the possible long-term recovery of plant communities in response to reclamation and that these outcomes become more pronounced as time passes. It emphasizes the importance of long-term monitoring and continued investment in reclamation practices to achieve successful and sustainable outcomes in disturbed landscapes (Shackelford et al. 2018).

Even though individual oil and gas well pads are relatively small in the study region (∼0.45 ha; Martinez & Preston 2018), the combined footprint of oil and gas operations is substantial, necessitating the use of tools that can be applied to many discrete disturbances across a landscape, such as TSS-RESTREND. By assessing the time series based on time since reclamation completion and comparing it with events identified from GE imagery, we uncovered specific management actions that yield distinct reclamation outcomes. However, further studies are required to establish whether differences in management actions consistently contribute to variability in reclamation outcomes across larger management units such as BLM field offices or districts.

The effectiveness of the TSS-RESTREND methodology for detecting changes in VPR residuals, particularly in relation to significant changes in the postreclamation plant community, warrants further discussion. The methodology frequently faced challenges in attributing these changes to substantial shifts in the plant community because the Chow test did not classify them as significant breaks unless there was a high amount of surrounding disturbance unrelated to the specific well pad. This observation suggests that the Chow test may set too stringent a statistical threshold for identifying breaks in the VPR residuals indicative of reclamation management actions, which could have been used in a subsequent segmented RESTREND analysis. Segmented RESTREND is generally more suitable for areas experiencing high levels of degradation, such as those during and following energy development (Burrell et al. 2017; Waller et al. 2018). Furthermore, our RESTREND results for the remaining pads exhibited weak or no relationship between time and VPR residuals, confirming that linear analysis is often inadequate for capturing frequent and subtle vegetation cover changes on disturbed and sparsely vegetated arid landscapes (Jiang et al. 2017; Lawley et al. 2013; Liu et al. 2019). Further work is needed to determine how the Chow test in TSS-RESTREND can be best used to detect important changes in the plant community postreclamation such as the establishment of invasive plant species on a reclaimed well pad.

The use of the long-term imagery archive and spatial scale of Landsat was critical for the application of TSS-RESTREND to oil and gas development in the study area. All of the study wells were developed before 2000, and the Landsat archive going back to 1984 allowed us to collect predevelopment data on the locations of the well pads, which is not possible with most sensors. In addition, the oil and gas pads in this area are typically far too small for the use of the 1.1-km AVHRR sensor used in previous TSS-RESTREND studies (e.g., Burrell et al. 2017). The long temporal resolution and a finer spatial resolution of the Landsat can detect subtle reflectance changes over time and provide a predisturbance baseline for evaluating the impacts of many land uses, including reclamation outcomes (Maynard et al. 2016; Zhu et al. 2019). Newer satellites with finer spatial resolution and more frequent imagery will certainly be important tools for land managers in evaluating the progress of reclamation going forward and new energy development, especially when field monitoring data are not available (Masek et al. 2020; Maynard et al. 2016). Similarly, purpose collected unmanned aerial imagery could supplement field or satellite data but cannot replace the wealth of information provided by the longer retrospective view of the Landsat or other longer-term satellite sensors (Ren et al. 2019; Taddeo 2022).

Although the Landsat sensor resolution is finer than other long-term satellites, the typical size of pads is only 4 to 14 Landsat pixels (0.36 ha to 1.26 ha), depending on the drilling technology and number of wells (Martinez & Preston 2018). The moderate spatial resolution of Landsat, the discrete nature of a pad’s footprint, partial reclamation of the pad footprint during production (interim reclamation), and influence of adjacent vegetation on target pixels may have prevented us from capturing the finer vegetative characteristics typical of early successional plant communities and decreasing our ability to detect compositional changes in the postreclamation plant community (Di Stéfano et al. 2020; Madonsela et al. 2017; Walker et al. 2014). In addition, there is a higher risk of an inaccurate classification when a smaller pad covers 1-4 Landsat pixels, limiting the use of Landsat for management decisions.

We observed in the GE imagery multiple redisturbances not recorded in the paperwork available in the COGCC database. Sometimes these redisturbances appeared to be from the performance of interim reclamation, where the well pad footprint was reduced from what was needed to establish the pad to the footprint needed for active extraction (e.g., 0.12 ha to 0.09 m²). In these types of redisturbances, interim reclamation appeared to establish an intermediary plant community that provided site stabilization, which was then completely removed at the time of final reclamation (USDI & USDA 2007). The interim process often includes multiple instances of topsoil disturbances, which can result in the loss of soil organic matter and hinder soil recovery and plant community stabilization after final reclamation (Costantini et al. 2016; Rottler et al. 2019).

For the oil and gas operators, recovery of the disturbed surface is often a secondary consideration in reclamation, with the bulk of reclamation cost going toward plugging the well to prevent seepage and pollution of underground soil and water (Alboiu and Walker 2019; Di Stéfano et al. 2021). Notices of management actions (sundry notices) over the life of the well and postreclamation do not often include important details on reclamation actions (seed mix and seed rates used, how seed was applied, etc.) and progress (monitoring data), which limits the ability of operators and land managers to pass down reclamation knowledge to future generations and requires more reliance on individual experience (Ladouceur et al. 2022; Samuel et al. 2022). In these circumstances, the TSS-RESTREND could be used to pinpoint events over a pad’s life cycle and to evaluate the long-term effects of management decisions, augmenting existing reclamation knowledge.

When applying TSS-RESTREND to postproduction oil or gas wells, we were able to identify changes in the plant community resulting from documented management actions, but it did not necessarily capture long-term recovery of the plant community. The TSS-RESTREND methodology was helpful in determining the types of management actions, recorded in documentation, that likely contributed to limited plant community establishment and successful well pad reclamation. Our analysis on the utility of TSS-RESTREND was constrained though by limited information available on specific management practices such as seeding and topsoil storage. Accurately assessing recovery trends may be improved, however, by adding a focus on consistent oversight and improved records management. Both of these aspects of reclamation were found to confound our analysis. The types of analysis conducted here coupled with better documentation of reclamation practices and outcomes (e.g., seeding dates and methods, seed mixes used, conditions of approval and liability release, etc.) could greatly improve our understanding of what management approaches are working, when, and where. This centralized approach to reclamation information can both improve management and facilitate development of new knowledge.

The USDA-ARS, Plains Area, is an equal opportunity/affirmative action employer, and all agency services are available without discrimination. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the US government. This work was supported by the US Geological Survey Ecosystems Mission Area.