Assessment of Rangeland Degradation in New Mexico Using Time Series Segmentation and Residual Trend Analysis (TSS-RESTREND)

: Rangelands provide signiﬁcant socioeconomic and environmental beneﬁts to humans. However, climate variability and anthropogenic drivers can negatively impact rangeland productivity. The main goal of this study was to investigate structural and productivity changes in rangeland ecosystems in New Mexico (NM), in the southwestern United States of America during the 1984–2015 period. This goal was achieved by applying the time series segmented residual trend analysis (TSS-RESTREND) method, using datasets of the normalized difference vegetation index (NDVI) from the Global Inventory Modeling and Mapping Studies and precipitation from Param-eter elevation Regressions on Independent Slopes Model (PRISM), and developing an assessment framework. The results indicated that about 17.6% and 12.8% of NM experienced a decrease and an increase in productivity, respectively. More than half of the state (55.6%) had insigniﬁcant change productivity, 10.8% was classiﬁed as indeterminant, and 3.2% was considered as agriculture. A decrease in productivity was observed in 2.2%, 4.5%, and 1.7% of NM’s grassland, shrubland, and ever green forest land cover classes, respectively. Signiﬁcant decrease in productivity was observed in the northeastern and southeastern quadrants of NM while signiﬁcant increase was observed in northwestern, southwestern, and a small portion of the southeastern quadrants. The timing of detected breakpoints coincided with some of NM’s drought events as indicated by the self-calibrated Palmar Drought Severity Index as their number increased since 2000s following a similar increase in drought severity. Some breakpoints were concurrent with some ﬁre events. The combination of these two types of disturbances can partly explain the emergence of breakpoints with degradation in productivity. Using the breakpoint assessment framework developed in this study, the observed degradation based on the TSS-RESTREND showed only 55% agreement with the Rangeland Productivity Monitoring Service (RPMS) data. There was an agreement between the TSS-RESTREND and RPMS on the occurrence of signiﬁcant degradation in productivity over the grasslands and shrublands within the Arizona/NM Tablelands and in the Chihuahua Desert ecoregions, respectively. This assessment of NM’s vegetation productivity is critical to support the decision-making process for rangeland management; address challenges related to the sustainability of forage supply and livestock production; conserve the biodiversity of rangelands ecosystems; and increase their resilience. Future analysis should consider the effects of rising temperatures and drought on rangeland degradation and productivity.


Vegetation Cover
NM's grassland biomes, particularly in the Desert Grassland Association-desert plain grassland and mixed grassland or mixed prairie-are dominated by black grama grass (Bouteloua eriopoda) and tobosa grass (Hilaria mutica) in the desert plains. Bluestem (Andropogon scoparius), san blue stem (A. halli), and Indian Grass (Sorghastrum nutans) are parts of mixed grassland or mixed prairie. Woodland biome can be found within a range of 1371 m to 2286 m elevation (amsl) and consists of one-seed juniper (Juniperus monosperma) and pinon pine (Pinus edulis) sometimes with oak (Quercus spp.) with an understory of grassland, forbs, and shrubs. Coniferous forest biomes can be found within a range of 2590 to 3658 m amsl in Petran subalpine and petran montane dominated by Enngelman spruce (Picea engel-mannii) and subalpine fir (Abies lasiocarpa). Petran montane forest association (2591.8 to 2896 m a.m.s.l) covers an extensive area of the state dominated Douglas fir, and white fir species. The land use/land cover data used in this analysis included the National Land Cover Dataset of 2011 [43] as well as the state's ecoregions [41] and quadrants. These datasets were overlayed with the identified breakpoints to be related to changes in productivity over these different land cover types.

GIMMS NDVI
To study temporal changes in rangelands' structure (i.e., significant changes in NPP trend), long-term NDVI timeseries was used as a proxy to NPP. Specifically, the GIMMS NDVI product based on the third generation Global Inventory Modeling and Mapping Studies NDVI (GIMM NDVIv3.1g) data was used in this study. These data was based on  (Figure 1). Three major types of vegetation biomes (i.e., forest, shrubland and grassland) comprise the respective ecoregions [42].

Vegetation Cover
NM's grassland biomes, particularly in the Desert Grassland Association-desert plain grassland and mixed grassland or mixed prairie-are dominated by black grama grass (Bouteloua eriopoda) and tobosa grass (Hilaria mutica) in the desert plains. Bluestem (Andropogon scoparius), san blue stem (A. halli), and Indian Grass (Sorghastrum nutans) are parts of mixed grassland or mixed prairie. Woodland biome can be found within a range of 1371 m to 2286 m elevation (amsl) and consists of one-seed juniper (Juniperus monosperma) and pinon pine (Pinus edulis) sometimes with oak (Quercus spp.) with an understory of grassland, forbs, and shrubs. Coniferous forest biomes can be found within a range of 2590 to 3658 m amsl in Petran subalpine and petran montane dominated by Enngelman spruce (Picea engel-mannii) and subalpine fir (Abies lasiocarpa). Petran montane forest association (2591.8 to 2896 m a.m.s.l) covers an extensive area of the state dominated Douglas fir, and white fir species. The land use/land cover data used in this analysis included the National Land Cover Dataset of 2011 [43] as well as the state's ecoregions [41] and quadrants. These datasets were overlayed with the identified breakpoints to be related to changes in productivity over these different land cover types.

GIMMS NDVI
To study temporal changes in rangelands' structure (i.e., significant changes in NPP trend), long-term NDVI timeseries was used as a proxy to NPP. Specifically, the GIMMS NDVI product based on the third generation Global Inventory Modeling and Mapping Studies NDVI (GIMM NDVIv3.1g) data was used in this study. These data was based on Remote Sens. 2021, 13, 1618 5 of 27 the Advanced Very High Resolution Radiometer (AVHRR) sensors [5,8] and it is of the most accurate datasets currently available. This version of the data corrected for calibration errors in a previous one [22]. The GIMMS NDVI data is available for the 1981-2015 period only at spatial and temporal resolutions of 8 km and 16 days, respectively. The GIMMS NDVI dataset was obtained from ECOCAST [44].

Rangeland Productivity
Rangeland productivity for the USA was obtained from the Rangeland Production Monitoring Service (RPMS) dataset developed by the United State Forest Service [45]. The dataset was prepared using NDVI from the Thematic Mapper for the 1984-2020 period at 250-meter pixels. The dataset provides estimates of annual production of rangeland vegetation in pounds per acre, which is useful in understanding trends and variability of rangeland forage resources [45]. The RPMS dataset was coupled with the assessment framework (Section 3.2) to evaluate the significance of rangeland productivity changes compared to those identified by the TSS-RESTREND method.

Precipitation, Drought, and Fire
Gridded precipitation data developed by Parameter-elevation Regressions on Independent Slopes Model (PRISM) was obtained from PRISM Climate Group [46]. Monthly precipitation at~4 km pixels was used. To evaluate the accuracy of breakpoints in terms of timing, extent, and direction of change relative to potential disturbances, historical drought and fire events were compared with the breakpoints. Mean monthly and annual self-calibrated Palmar Drought Severity Index (PDSI) was acquired from [47]. Fire data were obtained from New Mexico Resources Geographic Information System (RGIS) [48].

Methods
This analysis followed three main steps-data preparation, the application of TSS-RESTREND, and the development and application of an assessment framework to evaluate the detected breakpoints and changes in productivity ( Figure 2). The first step involved data acquisition, projection, resampling, and extraction of pixel values of GIMMS NDVI, PRISM precipitation, and RPMS productivity. The second step (Section 3.1) involved the application of the TSS-RESTREND method to identify breakpoints in space and time, their significance, and their structural changes. The last step provided a framework (Section 3.2) to evaluate the detected breakpoints using the RPMS-an approach that was lacking in the work by [8,23].

NDVI and Precipitation Relationships
The Bimonthly GIMMS NDVI data were filtered using a quality control (QC) procedure to remove non-reliable values based on a quality flag [49,50]. Pixels with at least 75% reliable values were used. Non-vegetated pixels were excluded based on a threshold of a median NDVI of less than 0.1 [50,51]. Out of 4474 pixels, a total of 4454 pixel were used in the analysis. Complete mean monthly timeseries of NDVI (ctsNDVI) was assessed over each pixel. The PRISM precipitation data was resampled using a bilinear method to match the spatial resolution of the GIMMS NDVI and assessed at monthly time scale.
An ordinary Least Square regression (OLS) was used to develop the relationships between the ctsNDVI and the complete timeseries of optimal accumulated precipitation for each pixel for the 1982-2015 period (referred to as ctsVPR) [8]. The optimal accumulated precipitation was determined by identifying the pairs of ctsNDVI and accumulated precipitation that provide the highest correlation coefficient. The correlation coefficients of all possible OLS relationships between ctsNDVI and a matrix of precipitation accumulation period (1-12 months) with offset period (0-3 months) were evaluated following [8]. The relationships with the highest correlation coefficients were used to calculate the residual between the observed and predicted ctsNDVI (referred to as ctsVPR-residual). The relationships with the highest correlation coefficients were used to calculate the residual between the observed and predicted ctsNDVI (referred to as ctsVPR-residual).

Figure 2.
Depreciation of TSS-RESTRND analysis and validation of breakpoints. In step 1, Data preparation which includes (i.e., acquisition of NPP data, precipitation, and productivity data), projection, resampling, and extraction of pixel level values. In Step II (TSS-RESTREND): method of residual fit, p -vector, trend of productivity change (i.e., decreasing, increasing, non-significance change and indeterminant), and ecosystem structure change were detected. In Step III (Evaluation), include stratified random sampling for validation, Welch's t test on randomly selected pixels (i.e., to test the significant change of productivity before and after the break years).

Identification of Breakpoints
The BFAST method was applied on the ctsVPR-residual to list potential breakpoints. Briefly, the BFAST method decomposes the timeseries into season, trend, and reminder components-an approach that allows to detect changes in the season and trend components [52,53]. The list of potential breakpoints identified by the BFAST method [8] based on the ctsVPR-Residuals were further evaluated for their significance in the VPR-allowing to assess their impact on NPP as represented by the maximum NDVI timeseries. A Chow test was applied on the VPR-Residuals on all pixels with significant VPR (α = 0.05). Based on this test, all pixels that have no significant breakpoints in the VPR-Residual (α = 0.05) but have significant VPR (α = 0.05) were further assessed using the standard RE-STREND by developing a regression between the VPR-Residuals and time (Equation (1)) [8].
where is intercept, is slope, and x is year Pixels that showed significant breakpoints in VPR-Residuals (based on the Chow test above) and had significant VPR were further evaluated for the significance of these breakpoints but in the VPR also using Chow test. Pixels with significant breakpoints in VPR-Residuals but not in the VPR were further assessed using the Segmented RESTREND by developing a multivariate regression between VPR-Residual, time, and a dummy variable (Equation (2)) Figure 2. Depreciation of TSS-RESTRND analysis and validation of breakpoints. In step 1, Data preparation which includes (i.e., acquisition of NPP data, precipitation, and productivity data), projection, resampling, and extraction of pixel level values. In Step II (TSS-RESTREND): method of residual fit, p-vector, trend of productivity change (i.e., decreasing, increasing, non-significance change and indeterminant), and ecosystem structure change were detected. In Step III (Evaluation), include stratified random sampling for validation, Welch's t test on randomly selected pixels (i.e., to test the significant change of productivity before and after the break years).

Identification of Breakpoints
The BFAST method was applied on the ctsVPR-residual to list potential breakpoints. Briefly, the BFAST method decomposes the timeseries into season, trend, and reminder components-an approach that allows to detect changes in the season and trend components [52,53]. The list of potential breakpoints identified by the BFAST method [8] based on the ctsVPR-Residuals were further evaluated for their significance in the VPR-allowing to assess their impact on NPP as represented by the maximum NDVI timeseries. A Chow test was applied on the VPR-Residuals on all pixels with significant VPR (α = 0.05). Based on this test, all pixels that have no significant breakpoints in the VPR-Residual (α = 0.05) but have significant VPR (α = 0.05) were further assessed using the standard RESTREND by developing a regression between the VPR-Residuals and time (Equation (1)) [8].
where β 0 is intercept, β 1 is slope, and x is year Pixels that showed significant breakpoints in VPR-Residuals (based on the Chow test above) and had significant VPR were further evaluated for the significance of these breakpoints but in the VPR also using Chow test. Pixels with significant breakpoints in VPR-Residuals but not in the VPR were further assessed using the Segmented RESTREND by developing a multivariate regression between VPR-Residual, time, and a dummy variable (Equation (2)) where β 0 is intercept, x is year, z is value of dummy variable (0 or 1), β 1 is slope, β 2 is offset at the breakpoint, and β 3 is the change in the slope at the breakpoint.
Pixels that showed significant breakpoints in VPR-Residuals and in VPR may indicate the presence of significant structural changes and thus the assumption of stationarity of the accumulation and offset periods used in the calculations of optimal accumulated precipitation before and after a breakpoint [5,54]. The set if the NDVImax and precipitation timeseries before and after a breakpoint was separated to recalculate new and independent VPR on either side of the breakpoints [8]. The precipitation data were standardized to account for the differences in the accumulation and offset periods among the breakpoints (Equation (3)).
where z is the standard score, x i is observed values, µ is the mean, and δ is the standard deviation. The NDVImax and standard score timeseries were further evaluated by fitting a multivariant regression (Equation (4)) [5,8].
where x is the standardized precipitation for year i, z the value of the dummy variable (0 or 1), β 0 is intercept, β 1 is slope, β 2 is the offset at the breakpoint, and β 3 the change in the slope at the breakpoint. Pixels that did not meet any of the above conditions were classified as indeterminantmet the following conditions: (1) had no significant VPR and no significant breakpoints in VPR; or (2) no significant VPR, significant breakpoints in VPR, and no significant breakpoints in segmented VPR.

Identification of Structural Changes
Structural changes of each pixel within NM ecosystems were identified based on three properties that include the significance of the breakpoints, direction of change (i.e., increasing and decreasing in productivity), and method of detection (as described in Section 3.1.2). With the regard to the significance level and direction of change, all pixels were classified into nine categories as shown in (Table 1) following [55] and similar to other dryland degradation studies in Australia [8] and China [5]. Stable decrease in productivity if the opposite of the previous condition was detected (a condition that can be considered as an initiation of degradation) [11,13]. A summary of the pixels based on these structural change categories was presented relative to NM's total area, quadrant, ecoregions, and major land use/land cover classes.

Breakpoints Assessment Framework
The TSS-RESTREND method by [8] provided only a qualitative assessment of breakpoints with other independent data in terms of significance and direction of changes. These two properties are important to properly characterize changes in productivity. A framework was developed in this study to address this gap by proposed means to quantitatively assess these properties. This framework consisted of four steps: (1) Develop random samples within the identified significant breakpoints; (2) Select and use independent productivity data; (3) Evaluate the random samples at the pixel level; and (4) Group and evaluate all random samples that fall within identified ecoregions-allowing to identify whether the changes in productivity at the individual pixels are reflective of consistent regional changes.
Random Samples: A set of random samples in terms of size and location can be developed following [56,57]. To estimate the size of random samples, a prior knowledge of image accuracy/variability is required [57]. The degree of variability in the RPMS data is unknown, thus it was assumed that the maximum variability would be about 50% with 95% confidence level and ±5% precision [56]. This criterion helps to determine representative sample pixels using a stratified random sampling approach [58]. The strata were developed based on ecoregions, land use/land cover, and direction of change. Based on the formula from [57], the sample size was calculated as (Equation (5)).
where n o is the sample size, Z is the selected critical value of desired confidence level (1.96), p is the estimated proportion of an attribute that is present in the population (0.5), and q = 1 − p (0.5) with p and e (0.5) represent the desired level of precision. The allocation of the random sample was based on the majority of the identified breakpoints in each method (Section 3.1.2). The random samples can be broadly allocated based on the direction of change and land cover within the identified significant breakpoints in the two main categories-the Segmented VPR and Segment RESTREND. Ecoregions and land use/land cover (NLCD 2011) [43] maps were overlaid to identify the locations of the pixels within these regions for further assessment.
Selection of independent data: The productivity data from the RPMS [45] was used to evaluate the accuracy of the identified breakpoints based on the TSS-RESTREND in terms of direction and significance of change in productivity. The data was resampled to match that of GIMMS NDVI [44,58]. Using the coordinates of GIMMS NDVI pixels, the corresponding RPMS productivity was extracted. The rangeland productivity was converted from pound per acre per year to kilogram per hectare per year. this process allowed to prepare a dataset that include RPMS productivity and with the corresponding TSS-RESTREND estimates of timing and direction of change, significance of breakpoints, and method of detection for each random sample.
Pixel level assessment: Using the location of the randomly selected samples (breakpoint pixels) and their identified year of breaks, the mean annual productivity for each pixel based on the RPMS data can be calculated for before and after the break years. Over each pixel, the number of years before and after the breaks were different and thus these mean values had different sample sizes. These mean values before and after the breaks over each pixel were then statistically compared using the Welch's student-test for the significance in their difference and direction of change. The Welch's t test was conducted assuming that the variances were not equal before and after the breaks [58,59]. The obtained results based on RPMS data using the Welch's test were compared with those from the TSS-RESTREND method. A summary of the agreement of this comparison was provided over these pixels as well as over the representative land cover classes.
Ecoregion level assessment: a similar approach was followed in this step of the framework except in this case the randomly selected pixels were grouped to represent ecoregions-all sampled pixels within each ecoregion represent a single entity (analysis unit). The mean values of productivity based on RPMS of each group before their individual breaks were further averaged and compared to the corresponding average after the break. The averages were assessed for their significance and direction of change using the Welch's student-test. The results were then compared with those obtained based on the TSS-RESTREND over the equivalent ecoregions.

Observed Types of Structural Changes
In general, identified breakpoints by the Segmented VPR method indicate significant ecosystem structural change-either decrease (i.e., irreversible degradation) or increase in productivity. Those identified with the Segmented RESTREND method indicate changes-either decrease (initiation of degradation) or increase (reversal from degradation)-in productivity that are not significant enough to alter the ecosystem structure.

Observed Types of Structural Changes
In general, identified breakpoints by the Segmented VPR method indicate significant ecosystem structural change-either decrease (i.e., irreversible degradation) or increase in productivity. Those identified with the Segmented RESTREND method indicate changeseither decrease (initiation of degradation) or increase (reversal from degradation)-in productivity that are not significant enough to alter the ecosystem structure.
The total number of pixels identified by the different methods and their direction of change were presented in Table 3. Those showed decreased (significant or insignificant) productivity were 786 or 17.6% (46,976 km 2 ) out of 4454. About 4.9% (~14,016 km 2 ) were identified using the Segmented VPR method; 4.3% (~12,160 km 2 ) were identified using the Segmented RESTREND method; and the remaining 8.5% were identified using the RESTREND method (insignificant decrease or increase in productivity). Table 3. A summary of the pixels (in % relative to 4454 pixels) and methods used to identify the direction of change in productivity in New Mexico during the 1982-2015 period. From all pixels that showed increased productivity (12.8% or 570 pixels, 34,816 km 2 ), 5.7% (16,320 km 2 ) were identified using the Segmented VPR method (i.e., significant gradual increase), 4.6% (12,992 km 2 ) were identified using the RESTREND method (i.e., insignificant increase), and the remaining 2.5% (7168 km 2 ) were identified using the Segmented RESTREND method (i.e., reversal of degradation).

Method of Change Detection
More than half of NM's area (55.6%~158,592 km 2 ) show insignificant change in productivity based on the total number of pixels that fell within the NSC category (p ≥ 0.1). The pixels identified with the RESTREND method that accounted for 25.8% (~73,472 km 2 ), and the remaining 29% were considered indeterminate.
More than half of NM's area (55.6%~158,592 km 2 ) show insignificant change in productivity based on the total number of pixels that fell within the NSC category (p ≥ 0.1). The pixels identified with the RESTREND method that accounted for 25.8% (~73,472 km 2 ), and the remaining 29% were considered indeterminate.

Dominant Land Cover Class Changes
Significant trends (increasing or decreasing) in productivity (irreversible degradation, initiation in degradation, reversal in degradation, and initiation in reversal of degradation) were observed on NM's dominant land cover classes that include shrubland, grassland, and evergreen forest (Figure 7). From all analyzed pixels (i.e., 4454), 2.2% of NM's grassland (6336 km 2 ), 4.5% of the shrubland (12,800 km 2 ), and 1.7% of the evergreen forest pixels (4800 km 2 ) showed significant decreasing trends in productivity-either with a complete or initiation of degradation. From all analyzed pixels,~1.4% (3840 km 2 ) and 0.88% (2496 km 2 ) of NM's grasslands that showed significant decrease in productivity were attributed to initiation of degradation (i.e., the Segmented RESTREND method) and non-reversal degradation (i.e., the Segmented VPR method), respectively (Figure 7).
Similarly, from all analyzed pixels, 2.2% (6336 km 2 ) and 2.3% (6464 km2) of NM's shrublands that showed significant decrease in productivity were attributed to initiation of degradation and non-reversal degradation, respectively. The degradation of shrublands was dominant in the northwestern (2048 km 2 ) and southeastern (2112 km 2 ) quadrants, while the degradation in grasslands was notably observed in northwestern (512 km 2 ) and northeastern (1664 km 2 ) quadrants.
Similarly, from all analyzed pixels, 2.2% (6336 km 2 ) and 2.3% (6464 km2) of NM's shrublands that showed significant decrease in productivity were attributed to initiation of degradation and non-reversal degradation, respectively. The degradation of shrublands was dominant in the northwestern (2048 km 2 ) and southeastern (2112 km 2 ) quadrants, while the degradation in grasslands was notably observed in northwestern (512 km 2 ) and northeastern (1664 km 2 ) quadrants. On the other hand, 5.7% (out of the 4454 pixels) of NM's shrubland pixels (16,248 km 2 ), 1.3% of the grassland pixels (3648 km 2 ), and 0.92% of the evergreen forest pixels (2624 km 2 ) experienced either significant gradual increase in productivity or reversal of degradation ( Figure 8). From all analyzed pixels, 1600 km 2 (0.56%) and 2048 km 2 (0.72%) of the grasslands that showed significant increase in productivity were attributed to On the other hand, 5.7% (out of the 4454 pixels) of NM's shrubland pixels (16,248 km 2 ), 1.3% of the grassland pixels (3648 km 2 ), and 0.92% of the evergreen forest pixels (2624 km 2 ) experienced either significant gradual increase in productivity or reversal of degradation ( Figure 8). From all analyzed pixels, 1600 km 2 (0.56%) and 2048 km 2 (0.72%) of the grasslands that showed significant increase in productivity were attributed to reversal in degradation (using the Segmented RESTREND method) and significant increase in productivity (using the Segmented VPR method), respectively.

Assessment of Breakpoints
This section provides a summary of the results obtained based on the breakpoints assessment framework described in Section 3.2.

Identified Random Samples
The total number of breakpoints based on the Segmented RESTREND and Segmented VPR methods were about 304 and 510, respectively. From which, 384 samples were randomly selected with 165 and 219 were based on the Segmented RESTREND and Segmented VPR methods, respectively (Tables A1 and A2 in Appendix A). Since the Segmented VPR had more significant breakpoints with a noticeable change in productivity, only those pixels were subjected to the random selection. Only 155 samples out of the 219 (or 71%) showed significant difference in productivity before and after the break years (either decreasing or increasing) (Figure 9). The remaining 64 samples did not meet the criteria for significance and were not considered. From the 155 samples, 24% and 76% were obtained over grasslands and shrublands, respectively. The distribution of the grassland samples (i.e., 24%) represented the Southwestern (SW) Tablelands (46%), Arizona/New Mexico (AZ/NM) Plateau (35%), and Chihuahua Desert (16%) ecoregions. Similarly, 58%, 25%, and 11% of the shrubland samples (i.e., 76%) represented the Chihuahua Desert, AZ/NM Plateau, and SW Tablelands ecoregions, respectively.

Changes in Productivity at Pixel Level
A summary of the comparison of the significance of the differences in mean productivity before and after the breakpoints between the Segmented VPR and the RPMS data was shown in Table 4 along with the direction (increase or decrease) of change in productivity. The results were presented as the percent of pixels that fell within each category relative to the number of the random samples (i.e., 155). reversal in degradation (using the Segmented RESTREND method) and significant increase in productivity (using the Segmented VPR method), respectively.

Assessment of Breakpoints
This section provides a summary of the results obtained based on the breakpoints assessment framework described in Section 3.2.

Identified Random Samples
The total number of breakpoints based on the Segmented RESTREND and Segmented VPR methods were about 304 and 510, respectively. From which, 384 samples were randomly selected with 165 and 219 were based on the Segmented RESTREND and Segmented VPR methods, respectively (Tables A1 and A2 in Appendix A). Since the Segmented VPR had more significant breakpoints with a noticeable change in productivity, only those pixels were subjected to the random selection. Only 155 samples out of the 219 (or 71%) showed significant difference in productivity before and after the break years (either decreasing or increasing) (Figure 9). The remaining 64 samples did not meet the criteria for significance and were not considered. From the 155 samples, 24% and 76% were obtained over grasslands and shrublands, respectively. The distribution of the grassland samples (i.e., 24%) represented the Southwestern (SW) Tablelands (46%), Arizona/New Mexico (AZ/NM) Plateau (35%), and Chihuahua Desert (16%) ecoregions. Similarly, 58%, 25%, and 11% of the shrubland samples (i.e., 76%) represented the Chihuahua Desert, AZ/NM Plateau, and SW Tablelands ecoregions, respectively. Based on the Segmented VPR method, all the 155 randomly selected samples indicated significant difference in productivity before and after the break years. However, based on the RPMS only 55% of them showed significant differences in mean productivity before and after the break years until 2019, respectively (Table 4). From the 155 samples, 37% and 18% showed persistent decreasing and increasing in mean productivity after the break years until 2019 on RPMS data, respectively. About 36% and 9% showed increasing and decreasing trends before and after the break years on the RPMS data, respectively.
From the 37% of the sampled pixels that exhibited persistent decrease in mean productivity after the break years, 13% were in the Chihuahua Desert, and 18% in the SW Tablelands ecoregions. From the 18% of the sampled pixels that exhibited significant and consistent increase in mean productivity after break years, AZ/NM Mountains and AZ/NM Plateau ecoregions accounted for 3% and 15%, respectively. From the 36% of the sampled pixels with increased but insignificant difference in mean productivity before and after the break years, 35% were over the Chihuahua Desert ecoregion, and the remaining were equally obtained over AZ/NM Mountains and AZ/NM Plateau (Table 4). From the 9% of the sampled pixels with decreasing trend but insignificant difference in productivity before and after the break years, 5% were over AZ/NM Plateau.

Changes in Productivity at Pixel Level
A summary of the comparison of the significance of the differences in mean productivity before and after the breakpoints between the Segmented VPR and the RPMS data was shown in Table 4 along with the direction (increase or decrease) of change in productivity. The results were presented as the percent of pixels that fell within each category relative to the number of the random samples (i.e., 155). Table 4. The percentages of pixels with decreasing and increasing trends in productivity as estimated from the RPMS data relative to randomly selected pixels (155) based on the Segmented VPR method categorized based on their significance of difference in mean productivity before and after the break years over major ecoregions in New Mexico. Based on the Segmented VPR method, all the 155 randomly selected samples indicated significant difference in productivity before and after the break years. However, based on the RPMS only 55% of them showed significant differences in mean productivity before and after the break years until 2019, respectively (Table 4). From the 155 samples, 37% and 18% showed persistent decreasing and increasing in mean productivity after the break years until 2019 on RPMS data, respectively. About 36% and 9% showed increasing and decreasing trends before and after the break years on the RPMS data, respectively.  Table 4. The percentages of pixels with decreasing and increasing trends in productivity as estimated from the RPMS data relative to randomly selected pixels (155) based on the Segmented VPR method categorized based on their significance of difference in mean productivity before and after the break years over major ecoregions in New Mexico. In the sampled grasslands and shrublands pixels (i.e., 24% and 76% of the 155 samples, respectively) with either decreased or increased productivity, the difference in mean productivity before and after the break years was insignificant in 5% and 40%, respectively on RPMS (Table A3 in Appendix B). The grassland (12% of the samples) and shrubland pixels (26% of the samples) with decreased productivity had significant lower mean productivity before the break years on RPMS data. Similarly, out of the pixels that showed increased trend on productivity, 7% from grasslands and 11% from shrublands had significantly higher mean production than that before the break years on RPMS (Table A3).

Changes in Productivity at Ecoregion Level
A summary of the comparison between the Segmented VPR and RPMS at the ecoregion level (Section 3.2) including grassland and shrubland cover classes based on the sampled pixels (i.e., 155) was provided in Table A5 (Appendix C)-allows to highlight whether the changes at the individual pixels are reflective of those at the regional level.
Based on the RPMS data, a continuous and significant decrease in shrublands' mean productivity after the break years in the Chihuahua (Welch's test p ≤ 0.0001) and the AZ/NM Plateau (p = 0.0397) ecoregions was observed. Similarly, significant decline in mean productivity of grasslands of the SW Tablelands (p ≤ 0.0001) and AZ/NM (p = 0.019) ecoregions were observed after break years. The shrublands within the AZ/NM Mountains and the SW Tableland ecoregions and the grassland in Chihuahua and the AZ/NM Plateau ecoregions exhibited insignificant differences in mean productivity between before and after the break years-stable ecosystem productivity during the study period.
In contrary, the sampled pixels over the shrublands in the Chihuahua Desert (p = 0.0194) and the AZ/NM (p = 0.00765). Mountain ecoregions showed significant increase in mean productivity after the break years (Table A5). There was a significant increase in mean productivity in the grasslands within AZ/NM Plateau ecoregions after the break years. However, sampled shrubland and grassland pixels in the AZ/NM Plateau and the Chihuahua Desert, respectively, showed insignificant difference in mean productivity before and after the break years, suggesting a negligible increase in mean productivity after the break years during the study period.

Characteristics of Change
The significance of the breakpoints as identified by the TSS-RESTREND methods can be interpreted relative to observed ecosystem structural changes [8,13]. Out of 67.9% of the pixels that met the criteria of RESTREND [13], 8.5% and 4.6% showed decreased and increased productivity, respectively. These pixels exhibited gradual change as their VPR remined consistent over time with no major ecosystem structural changes [13,60]. The behavior of the pixels that were identified by the Segmented VPR method (11.6%) experiencing irreversible degradation (4.9%) or increased productivity (5.7%) can partially be attributed to abrupt land use changes (decreeing or increasing) induced by human activities and climate variability [13,61] such as overgrazing or easing of drought conditions [62]. More than half of NM that did not experience significant change in productivity (NSC = 55.6%- Table 2) was dominantly in the western part of the state. This can partially be explained by the fact that western NM is the driest region in the state. Thus, it experiences weaker interactions related to climate variability and human activities-thus minimal effects on productivity. The pixels that were identified as indeterminant (10.8%) were those that none of the methods was able to fit their observed behavior [8,55] and there was no clear explanation for such behavior.
There were relatively higher human activities in northeastern, northwestern and southeastern NM represented by crude oil and natural gas production and livestock grazing practices. The prevalence of irreversible degradation and initiation of degradation were apparent in these regions [63]. On the other hand, significant changes in productivity (i.e., increase) was observed mostly in southwestern and northwestern NM, where forest were minimally influenced by human activities. On landscapes where human activities were dominant, i.e., southeastern NM, significant increase in productivity can be associated with parallel restoration activities and land use management practices [55,64].
Overall, there was a consistent increase in the number of detected significant breakpoints during the 1980s, 1990s, and 2000s, with 9.5%, 19.2%, and 71.4%, respectively ( Figure 10). Again, the combined effects of climate change and anthropogenic factors can explain this increase [65]. For example, NM's precipitation remained close to the long-term average (only recently showed increased variability) since 1920 while air temperature showed an increasing trend since 1970s [63,66,67] with a simultaneous increase in fossil fuel production activities since 1980s.
ticed that from all the detected breakpoints during the study period , 67.9% were observed only in 2000s from which 38.2% exhibited decreased trend in productivity (Figure 10a). Coincidently, these breakpoints overlapped with frequent and extended periods of drought events in NM that were also observed regionally in the southwest USA during 2000s. The impact of this regional drought had resulted in a decrease in productivity in more than 30% of the coterminous USA, out of which ~15% (equivalent to more than 41 million ha) was rangelands [37].  It is challenging to identify a single factor that can be considered as the direct cause of these breakpoints and ecosystem structural changes [8,68]. However, it was possible to compare these breakpoints with some observed ecosystem structural changes whether gradual (e.g., land use management and climate change) or abrupt (e.g., wildfire). Such a comparison can help in evaluating the accuracy of these breakpoints in terms of timing, distribution, and direction of change [8,53].

Land Cover Changes Relative to Drought and Wildfire
Some of the major drivers of change in NM's dryland ecosystems were identified as climate (e.g., drought) which is influenced by increased concentration of atmospheric greenhouse gases (GHG) (e.g., CO 2 ); wildfire; grazing practices (i.e., livestock density); and land use conversion [21,35,68] were used here to highlight their effects on the breakpoints.

Detected Changes Compared to Previous Studies and Current Restoration Activities
The significant decrease in grasslands and shrublands was mostly observed in northwestern, northeastern, and southeastern NM-consistent with [4] that indicated that degradation was mostly over grassland-savanna. Degradation of grasslands, in some cases, was directly linked to increased productivity of shrublands. In NM, increased productivity in shrublands was noticed southwestern deserts and plains. This increase was attributed to land cover conversion from black grama and other valuable grasses-dominated areas to bushes due to over grazing and drought [69]. The evident climate warming, expected increase interannual precipitation variability, and projected increase in aridity can further enhance the growth of shrubs (i.e., encroachments of woody plants) over grasses since the later are heavily depended on transient surface moisture [70,71]. It was also argued that increased CO 2 concentration in the atmosphere can improve plant CO 2 uptake and reduce water loss through plant's stomata thus increasing the photosynthesis process and this mechanism actually favor woody vegetation over non-woody ones [72,73].
The general notion that was suggested in some studies was that encroachments of shrublands to grasslands can be considered as one stage of degradation that can threaten the integrity of the rangelands [74]. However, it was argued that, shrublands are not necessarily degraded, nor do they necessarily represent "degradation" due to their ability to support valued ecosystem services [71,75], and they also have a long-term mean annual above ground NPP-equivalent to that of grasslands [76]. In NM, there have been a number of restoration and brush management efforts undertaken to control invasive species, and retain favored shrub species [64]. This suggested that the attribution of drivers to changes in productivity in environmental studies remains challenging owing to the diversity of driving forces and limited sources of ground truth data for validation [77].
The increased productivity in grasslands that was observed in western NM can partially be attributed to local scale successful restoration efforts by the Bureau of Land Management [64]. These efforts targeted the replacement of Creosote and mesquite by healthy grasslands, and reclamation of surfaces resulted from oil and gas extraction operations. According to Powell [78], gradual improvements in range conditions (increased productivity) in southwestern NM was not only related to the results of better moisture, but also cumulative efforts on rangeland management, such as proper stocking, vastly improved grazing distribution, and brush management.

Breakpoints and Drought
The sc-PSDI for NM and the detected significant breakpoints with decreasing and increasing trends were shown in Figure 10. A weighted average of the sc-PDSI was calculated over the breakpoints with decreasing and decreasing trends separately to evaluate their timing against drought events. Some previous findings suggested that extended periods of drought (or relatively dry conditions) can introduce lasting negative impacts (degradation) on rangelands and other ecosystems [8,79]. Based on Figure 10, it was noticed that from all the detected breakpoints during the study period (1982-2015), 67.9% were observed only in 2000s from which 38.2% exhibited decreased trend in productivity (Figure 10a). Coincidently, these breakpoints overlapped with frequent and extended periods of drought events in NM that were also observed regionally in the southwest USA during 2000s. The impact of this regional drought had resulted in a decrease in productivity in more than 30% of the coterminous USA, out of which~15% (equivalent to more than 41 million ha) was rangelands [37].
The increased number of breakpoints after 2000 suggested a permanent degradation or damage to ecosystems which mechanistically can occur when these ecosystems have little to no time to recover from a previous consecutive drought event. A recovery period from a drought event as defined by [79] is the return of an ecosystem to pre-drought values of GPP and can vary from immediate to multiple years depending on vegetation, climate, disturbance, and drought. It was indicated that dryland ecosystems such as those in NM, had experienced increased recovery periods that were highly sensitive to precipitation and temperature [79]. With the expected future projections of rising temperature and drier conditions, the drought recovery of dry ecosystems would even be longer. When droughts have shorter return period (or more frequent) this makes ecosystems more susceptible to drought or ecosystem degradation continues and builds up until a threshold is reached and the degradation become irreversible (or permanent), a condition that was referred to in [79] as a tipping point. This study suggested that these thresholds may have been reached as represented by the detected breakpoints. This explains, to some extent, the appearance of breakpoints compared to drought events. The timing of these breakpoints suggested the accumulation of dry conditions that impact ecosystem functions until reaching tipping points. Further analysis is needed to understand the timing or emergence of the breakpoints and drought accumulation periods and the response of ecosystems relative to vegetation type and climate.
Moreover, these recent exceptionally drought conditions in NM (during 2000-2015) also varied in spatial extent, duration, and intensity [80]. These variations may have contributed to the occurrence of 38% breakpoints with decreasing trends (i.e., initiation and irreversible degradation) in southeastern (Eddy, Lea, Chaves, and Otero counties); northeastern (Colfax, De Bacca, and San Miguel counties); and northwestern NM (San Juan county). From all pixels that experienced significant decrease in productivity, 35% were observed in De bacca, San Miguel, San Juan, and Otero counties as drought might have a profound influence in reducing vegetation productivity (e.g., plant mortality) [81,82].
On the other hand, reversal from degradation and significant increase in productivity was dominantly observed in the northwestern, southeastern, and southwestern NM (Figure 8). In NM, wetter years were observed from 1986 to the end of 1990s, followed by drier years from 1999-2003 and 2005 ( Figure 10). Significant breakpoints with increase in productivity were identified during 2000s dry years in Otero, Dona Anna, Luna, San Juan and Socorro counties. About 13.5% of these breakpoints were observed in the first four of counties. This can partially be attributed to the fact that precipitation remained near the long-term mean after the drier years over these counties that resulted in reversal from degradation or significant increase in productivity, respectively ( Figure 10).

Breakpoints and Wildfire
Wildfire generally reduces plant cover, alters habitat structures, decrease rangeland conditions, and requires much longer recovery period [83]. The significant breakpoints with decreased productivity (i.e., irreversible degradation based on Segmented VPR) that coincided with fire events are shown in Figures 11 and 12. Major wildfire incidents occurred in 1989 and 1994 followed by a consistent increase in frequency each year since 2001. This increase was noticeably concurrent with drought events. The number and time of some breakpoints mostly followed these of the fire incidents. The alignment of fire incidents with breakpoints can indicate the ability of TSS-RESTREND to detect the timing of ecosystem changes as suggested by [8]. However, based on Figures 11 and 12, it appeared that only a small number of breakpoints coincided with fire incidents during the study period-suggested that fire may have a limited contribution to the development of breakpoints. While, the timing of these fire incidents can partly explain the occurrence of the breakpoints, it was not rationale to state that these fire incidents were the only direct cause of the breakpoints. Other factors need to be considered to provide a rational explanation of the remaining breakpoints such as vegetation types and their response to environmental and climate disturbance (drought). The accumulation of successive dry conditions can push an ecosystem to pass a threshold of permanent damage of vegetationideal fire-prone conditions. Further analysis is needed to identify the causes of breakpoints relative the fire and was out of the scope of this analysis. conditions can push an ecosystem to pass a threshold of permanent damage of vegetation-ideal fire-prone conditions. Further analysis is needed to identify the causes of breakpoints relative the fire and was out of the scope of this analysis.

Breakpoints and the RPMS
Based on the RPMS data, a decrease in the mean annual productivity was detected over the pixels with significant breakpoints (i.e., identified by Segmented VPR method) during the 1984-2019 period ( Figure 13). The mean productivity of these pixels was lower than that of the long-term of the 36 years which was about 613.2 kg/ha. Moreover, the variability in the annual productivity of these pixels ranged from 588 kg/ha in 1995 to 641 kg/ha in 2006.

Breakpoints and the RPMS
Based on the RPMS data, a decrease in the mean annual productivity was detected over the pixels with significant breakpoints (i.e., identified by Segmented VPR method) during the 1984-2019 period ( Figure 13). The mean productivity of these pixels was lower than that of the long-term of the 36 years which was about 613.2 kg/ha. Moreover, the variability in the annual productivity of these pixels ranged from 588 kg/ha in 1995 to 641 kg/ha in 2006. Of the sampled pixels that were evaluated for changes in productivity before and after the break years, 38% with decreasing trend showed a significant difference in mean productivity after the break years with 12% in grasslands (Table A3) and 26% in shrublands (Table A4). Similar to this finding, [84,85] indicated that there was an apparent decline in NPP in the southwestern USA that was attributed to the response of these ecosystems to the combination of a warmer temperature and a decline in precipitation. The pixels that showed significant decrease in productivity in NM's rangeland (i.e., in Chihuahua Desert and AZ/NM plateau ecoregions) coincided with those from [45] that showed a similar behavior (Table A5). The difference in the productivity before and after the break years might be substantiated by the xeric nature of the southwestern USA and southeastern Great Plains in response to interannual variability of the precipitation [86]. On the other hand, 8% and 36% of sampled pixels with decreased and increased productivity (i.e., based on the Segmented VPR), respectively showed insignificant difference in mean productivity before and after the break years based on the RPMS data [45]. These findings were in agreement with those by [83] which suggested that the structure and composition of semi-arid and arid regions of the southwestern USA have undergone noticeable changes over the last two decades. Of the sampled pixels that were evaluated for changes in productivity before and after the break years, 38% with decreasing trend showed a significant difference in mean productivity after the break years with 12% in grasslands (Table A3) and 26% in shrublands (Table A4). Similar to this finding, [84,85] indicated that there was an apparent decline in NPP in the southwestern USA that was attributed to the response of these ecosystems to the combination of a warmer temperature and a decline in precipitation. The pixels that showed significant decrease in productivity in NM's rangeland (i.e., in Chihuahua Desert and AZ/NM plateau ecoregions) coincided with those from [45] that showed a similar behavior (Table A5). The difference in the productivity before and after the break years might be substantiated by the xeric nature of the southwestern USA and southeastern Great Plains in response to interannual variability of the precipitation [86]. On the other hand, 8% and 36% of sampled pixels with decreased and increased productivity (i.e., based on the Segmented VPR), respectively showed insignificant difference in mean productivity before and after the break years based on the RPMS data [45]. These findings were in agreement with those by [83] which suggested that the structure and composition of semi-arid and arid regions of the southwestern USA have undergone noticeable changes over the last two decades.

Limitation and Future Work
As the study aimed to detect degradation of vegetation, rangelands (i.e., shrubland and grasslands), the authors acknowledged some limitations that need to be considered when using these findings. The TSS-RESTREND was able to identify the breakpoints, their location, time, and type and direction of change in productivity at the state level. Due to the limited availability of ground biomass data, pixels that exhibited significant changes were validated using derived productivity from the RPMS [45] that showed only 55% agreement. Ground-based data can be more accurate means in assessing the identified breakpoints. Another important factor that can contribute to degradation is the temperature which can significantly affect NM's dryland ecosystems.
For future, the authors would address some of these limitations using ground observations (when available) for example from the long-term rangeland monitoring observatories at the Jornada Experimental Range. The obtained breakpoints seemed consistent with some of the disturbances (drought and fire). However, the authors would consider the combined use of temperature and precipitation to improve the detection of breakpoints such as the pixels that have been classified as indeterminant. Additional quantitative analysis would be conducted to determine how drought may affect the timing of breakpoints, rangeland condition, and productivity, which can be useful in developing rangeland management practices to adapt and mitigate drought impacts.

Conclusions
This study evaluated the degradation of New Mexico's rangelands during the 1984 2015 period with respect to climate using an NDVI timeseries-as a surrogate of NPPand precipitation to represent climate variability. These datasets were evaluated using the TSS-RESTREN method to detect breakpoints in the NDVI timeseries, direction and significance of change in productivity at each pixel. The study developed a breakpoint assessment framework that allowed to quantitively evaluate the identified changes that used an independent productivity data (i.e., RPMS). A qualitative assessment of the changes against land management activities, drought, and wildfire was also conducted.
The study indicated that about 17.6% of New Mexico experienced a decrease in productivity while 12.8% of the state experienced an increase in productivity. More than half of the state (55.6%) had insignificant change productivity, 10.8% was classified as indeterminant, and the remaining 3.2% was considered as agriculture was not evaluated.
The degradation in productivity was observed in about 2.2%, 4.5%, and 1.7% of NM's grassland, shrubland, and evergreen forest land cover classes, respectively. Simultaneously, about 5.7%, 1.3%, and 0.92% of NM's shrublands, grasslands, and evergreen forests were characterized by an increase in productivity, respectively. Regionally, significant decrease in productivity was observed in the northeaster and southeastern quadrants of the state while significant increase was observed in northwestern, southwestern, and a small portion in southeastern quadrants. The timing and number of detected breakpoints coincided with NM's drought frequency and severity and some fire incidents.
The TSS-RESTREN showed 55% agreement with the RPMS data over areas with significant decrease in productivity based on randomly selected pixels. Regionally, there was an agreement between the TSS-RESTREND and RPMS on the occurrence of significant degradation in productivity over the grasslands and shrublands within the Arizona/New Mexico Tablelands and in the Chihuahua Desert ecoregions, respectively.
A long timeseries assessment of rangeland productivity in New Mexico is critical to support decisions related to ecosystems management and conservation. The findings of this study can be used to address some of the rangeland degradation challenges that directly impact rangeland conditions, forage supply in the region, and New Mexican's livelihood as these systems support livestock production. In areas where degradation was prevalent (i.e., northeastern, and southeastern NM), special intervention would be needed to conserve the biodiversity of rangeland and increase the resilience of these ecosystems. As this study considered only precipitation, it was noticed that the rising temperature can play a significant role in vegetation. Thus, future analysis should consider its effects on rangeland productivity. Future research should also pay more attention to the association of degradation and productivity with recurring droughts.

Data Availability Statement:
The data used in this study are openly available in data [43][44][45][46][47][48]. Table A1. Allocation of the randomly selected pixels as identified with Segmented VPR method over the different land cover classes based on the National Land Cover Dataset of 2011 (NLCD 2011) along the corresponding ecoregions.

Appendix B
The percentages of the randomly selected pixels (i.e., 155 identified by Segmented-VPR method) over grasslands and shrublands that were evaluated for their significance of the difference in mean productivity between before and after the break years at the pixel level using the Rangeland Production Monitoring Service (RPMS) data. Table A3. The percentages of the randomly selected grassland pixels with respect to their significance in the difference in mean productivity before and after the break years based on the RPMS data along with their identified direction of change over the major ecoregions in New Mexico.  Table A4. The percentages of the randomly selected shrublands pixels with respect to their significance in the difference in mean productivity before and after the break years based on the RPMS data along with their identified direction of change over the major ecoregions in New Mexico.

Appendix C
A summary of the statistical significance test result of the difference in mean in productivity before and after the break years using the randomly selected pixels (i.e., 155 identified by Segmented-VPR method) based on the Rangeland Production Monitoring Service (RPMS) data at the ecoregion level. Table A5. Significance of the difference in mean productivity before and after the break years based on the on increasing and decreasing shrublands and grassland randomly selected pixels using RPMS data at the ecoregion level.