Anthropogenic Climate Change Negatively Impacts Vegetation and Forage Conditions in the Greater Four Corners Region

The Southwestern U.S. has been in a megadrought for the past two decades, stressing vegetation across the region and leading to poor forage conditions on rangelands. How has human‐induced climate change affected air temperatures and atmospheric evaporative demand in the greater Four Corners region of the Southwest, and what have the corresponding impacts been to vegetative productivity? We address this question via a multi‐step climate change detection and attribution analysis. We first relate evaporative demand (quantified by Vapor Pressure Deficit, or VPD) and precipitation to the Normalized Difference Vegetation Index (NDVI) in six unique plant‐precipitation zones. VPD is strongly negatively related to NDVI, and the strength of this relation is similar in magnitude to precipitation. We then use a large ensemble of climate change simulations to estimate the magnitude of VPD increases attributable to human‐induced climate change. We find that modest increases in air temperatures have led to large increases in VPD, with the influence of climate change about twice the inter‐annual standard deviation by 2020. Finally, we estimate the reductions in NDVI due to increased VPD and present the reductions in terms of total annual Aboveground Net Primary Productivity. We conclude that human‐induced warming has roughly doubled the magnitude of recent vegetation deficits, especially in mid‐precipitation areas (receiving 250–500 mm of precipitation annually), concentrated in New Mexico. These results indicate that increased temperatures from human‐caused climate change are having persistent and damaging impacts on vegetation productivity, with significant implications for ranchers and other land users in the region.

Rangelands cover large portions of the Southwest, including 81% of federal lands in Arizona and New Mexico (USDA, n.d.). They host biodiverse ecosystems, supporting diverse plant and animal communities, as well as vegetation suitable for grazing and browsing for livestock and other animals, including grasses, forbs, and shrubs (Redsteer et al., 2013;USDA, 2009). Many inhabitants in this semi-arid region, including multiple tribal communities, engage in ranching-hence, vegetation productivity on rangelands has important implications for local livelihoods (Ferguson et al., 2016;Garfin et al., 2013;Redsteer et al., 2013). On the Colorado Plateau cattle and sheep are the most common livestock (Copeland et al., 2017;USDA, 2019). Most agricultural sales in New Mexico and Colorado counties are from cattle ranching, while sheep grazing has been a major cultural and subsistence activity for several communities in New Mexico and Arizona since their introduction in the 1700s (Frisvold et al., 2013;Milchunas, 2006).
Recent droughts, however, have degraded rangeland vegetation productivity across the region: in 2018 and 2020, over half of the region's pastureland was rated as in poor or very poor condition (NOAA, 2018(NOAA, , 2020. In desert rangelands, grazing should stay below 35% of annual vegetation production in order to maintain forage production levels, as overgrazing can harm future forage production as well as overall vegetation composition and ecosystem health (Holechek et al., 1999). Therefore, if herd sizes stay constant through drought conditions, the combination of drought-induced low productivity and constant forage demand may lead to overgrazing, resulting in further reduced forage and sometimes stressing sensitive or protected vegetation species (Frisvold et al., 2013, pp. 231-232;Grand Canyon Trust, 2022;Holechek et al., 1999). Hence, understanding the effect of drought conditions on rangeland vegetation productivity has important implications for livestock forage availability and overall rangeland management.
Here, we examine two related measures for rangeland vegetation productivity: Net Primary Productivity (NPP) and the Normalized Difference Vegetation Index (NDVI). NPP is the net rate of CO 2 gain by plants after growth and maintenance respiratory needs are met, thus representing the carbon available to increase plant biomass (Chapin et al., 2012). As such, the portion of NPP allocated to aboveground tissue (Aboveground Net Primary Productivity [ANPP]) is used to estimate forage production conditions for livestock grazing (Evans et al., 2010, p. 140;Hartman et al., 2020;Jones et al., 2021;Reeves et al., 2020). Total ANPP can be estimated for a growing season or annually by physically measuring biomass (Evans et al., 2010, p. 140). However, ANPP is difficult to measure at scale, so algorithms have been developed to estimate time-integrated ANPP on rangelands using indices such as NDVI (e.g., Hartman et al., 2020;Poděbradská et al., 2019;Reeves et al., 2020;W. Yuan et al., 2019), a satellite-based measure of vegetative greenness (Tucker, 1979). For vegetation with a distinct growing season (where leaf area is closely coupled with carbon gain), NDVI and ANPP covary (Evans et al., 2010;Gaffney et al., 2018). As NDVI is globally available and is relatively easy to derive, it has been used as a proxy for vegetation productivity for places and applications that require measurements of vegetation productivity but lack the ability to measure or calculate ANPP (e.g., Benami et al., 2021;Funk & Brown, 2006;Poděbradská et al., 2019;Thoma et al., 2002).
As photosynthesis requires the input of solar radiation and involves the exchange of carbon dioxide for water vapor, this photochemical process is generally constrained either by available energy or water. In warm, dry environments, such as the greater Four Corners region, ANPP is generally water-limited rather than energy-limited (Chapin et al., 2012, p. 171;Famiglietti et al., 2021;Funk & Brown, 2006;Jiao et al., 2021). This water limitation occurs primarily via two related yet distinct pathways-through decreasing soil moisture and increasing atmospheric evaporative demand (Novick et al., 2016). Plants can "sense" environmental conditions, including at the roots and at the leaf (Grossiord et al., 2020;Novick et al., 2016). Therefore, if soil moisture depletes or atmospheric evaporative demand at leaf level rises, vegetation may close stomata and limit photosynthesis-and therefore productivity-in order to conserve water (Dannenberg et al., 2022;Grossiord et al., 2020;Novick et al., 2016). While most vegetation in semi-arid regions is adapted to hot and dry conditions (e.g., plants have higher water use efficiencies and physiological traits to resist cavitation in low soil moisture conditions), drought conditions-or high evaporative demand and low soil moisture-can constrain photosynthetic activity (Dannenberg et al., 2022;Jiao et al., 2021;Taiz & Zieger, 2004, pp. 96-101).
The water limitation on ANPP is particularly pronounced for this study region (Dannenberg et al., 2022;Famiglietti et al., 2021). Rangelands are generally non-irrigated, so available soil water content is constrained to water inputs from precipitation and runoff (Frisvold et al., 2013, pp. 231-232). As this semi-arid region experiences relatively low climatological precipitation, precipitation acts as a primary determinant of vegetation  (Evans et al., 2010, p. 148;Jiao et al., 2021;Liu et al., 2019). Moreover, high temperatures can lead to increased atmospheric evaporative demand, stressing the already limited water resources (Dannenberg et al., 2022;Frisvold et al., 2013, pp. 231-232). Increased air temperatures lead to increased Saturation Vapor Pressure (SVP), or the amount of water vapor required to fully saturate the air. Without an increase in actual water vapor, this increases Vapor Pressure Deficits (VPD), yielding higher evaporative demand. In semi-arid regions and warmer seasons, VPD and therefore temperature become important determinants of vegetative productivity (Crimmins et al., 2017;Dannenberg et al., 2022;El-Vilaly et al., 2018;Weiss et al., 2009;A. P. Williams et al., 2020). As such, rangeland vegetation in the region is sensitive to variations in both VPD and precipitation.
The Four Corners region also has strong land-atmospheric coupling. Reduced soil moisture can increase near-surface temperatures by reducing latent heating via evapotranspiration and increasing sensible heating, thereby driving further increases in VPD (J. Zhang et al., 2008). This coupling is strongest in the summer months preceding the arrival of the monsoon, which corresponds with peak climatological temperatures (Weiss et al., 2009). This is also the season in which rangeland forage for livestock is in highest demand (Crimmins et al., 2017, p. 92;.
Previous studies have examined changing climatic conditions in the Southwest and their anthropogenic drivers. While significant reductions in precipitation have occurred since 2000, these are largely due to natural variability (Lehner et al., 2018;A. P. Williams et al., 2020). Temperatures and VPD have also been rising for the region, which studies have attributed to anthropogenic forcing (Bonfils et al., 2008;Mankin et al., 2021;A. P. Williams et al., 2020;E. L. Williams et al., 2020;Zhuang et al., 2021). It is this combination of below-average precipitation and above-average temperatures that has reduced soil moisture, increased VPD, and plunged the region into drought (Mankin et al., 2021;A. P. Williams et al., 2020. However, there has been relatively little work examining the effect of warming-exacerbated drought conditions on vegetation productivity and resulting forage conditions. In a previous analysis, we quantified the effect of human-induced temperature increases on peak annual NDVI for the 2018 water year by developing a statistical model relating VPD with NDVI and using modeled temperature fields from two large climate model ensembles. We found that anthropogenically forced increases in temperatures resulted in NDVI reductions of 18%-30% (E. L. . That study, however, was based on a spatial aggregation of the full region for a single year (2018), and hence did not account for variation in the climate-vegetation relationship across the range of elevations, land cover types, and precipitation regimes that characterize this diverse region. In this study, we expand on that work, exploring the links between temperature, precipitation, VPD, and both peak NDVI and total annual ANPP across the variable topography and climates of the Four Corners region for the past two decades, since the onset of the megadrought. By accounting for spatial and temporal variation, we can more accurately account for local conditions which may increase or attenuate local sensitivity to increasing VPD. In this study, we quantify the increase in VPD, and decrease in peak NDVI, attributable to anthropogenic increases in temperature/SVP, and, furthermore, estimate the corresponding reduction in ANPP. Thus, we link estimates of human-induced warming to a metric (ANPP) strongly related to forage availability, relevant for both ranchers and more broadly ecosystem health.

Study Region and Data
The study region was defined as the greater Four Corners (34-39°N, 112-105°W), covering the extent of exceptional drought in 2018 (as defined as the study region for E. L. , and much of the extent of exceptional drought in late 2020 ( Figure 1). Gridded (0.04° latitude × 0.04° longitude resolution) monthly precipitation, maximum and minimum temperature (T max , T min ), and maximum and minimum VPD (VPD max , VPD min ) were retrieved for 1895-2020 from the PRISM Climate Group at Oregon State University. Monthly gridded (0.05° latitude × 0.05° longitude resolution) NDVI data were retrieved from the Moderate Resolution Imaging Spectroradiometer (MODIS) (MYD13C2.006) on the USGS Earth Data platform for 2003-2020, the period of full record years (Didan & Huete, 2015). Land cover data (30 m resolution) were retrieved from the Multi-Resolution Land Characteristics (MRLC) Consortium, including the National Land Cover Database classification scheme used for primary land cover types (e.g., shrubland) for 2001 and 2016, and rangeland fractional components (for 2016) to determine percent cover of different land cover types on rangelands (e.g., percent bare) (Dewitz and USGS, 2021;Rigge et al., 2020). Elevation (∼250 m resolution) was retrieved using the R package, elevator, from AWS Terrain Tiles. Soil color and characteristics (∼200 m resolution) were retrieved from the Soil Survey from the USDA Natural Resources Conservation Service. All data were reprojected to WGS84, and furthermore resampled to the same resolution as MODIS (0.05°) using bilinear interpolation for the continuous variables and nearest neighbor interpolation for categorical variables.
Additionally, maximum near-surface air temperature was accessed from the Coupled Model Intercomparison Project Phase 6 (CMIP6) Pre-Industrial (PI), historical, and Shared Socioeconomic Pathway 245 (SSP2-4.5) experiments. SSP2-4.5 represents a "middle-of-the-road" scenario. For the historical (1850-2014) and SSP2-4.5 (2015-2100) experiments, the same 152 simulations were retrieved across 25 models, while one simulation per model was retrieved for the PI experiment (see Table S1 in Supporting Information S1 for list of models and simulations). For each simulation, the historical and SSP2-4.5 data were stacked and subset to 1895-2020, creating 152 spatio-temporal cubes, and then the cubes were resampled to 100 km resolution using bilinear interpolation. Similarly, 25 PI spatio-temporal cubes were created and resampled.

Spatially Aggregated Regression
To determine regional climatic drivers of NDVI change, relationships were explored between the spatial averages over the full study region of NDVI and precipitation (P), maximum temperature (T max ), and maximum VPD (VPD max ). Peak NDVI values are strongly associated with annually integrated NDVI and ANPP (Reeves et al., 2020). In this region, vegetation greenness typically reaches a maximum in tandem with air temperatures, during boreal summer. August was identified as the predominant peak NDVI month (Figure 2d), and hence was chosen as our predictand.
Next, a "best fit" model was identified to explain the variance in NDVI for the study region for 2003-2020. First, for each variable (P, T max , and VPD max ), 24 unique seasonal averages were created (e.g., Mar-Aug T max ). These seasonal averages were then standardized and converted to study-region-wide time series. Next, linear models were fit for all combinations of two seasonal averages (e.g., Mar-Aug T max + Mar-Aug P). Models were retained if they had a strong fit (r 2 > 0.6), a significant p-value for the temperature variable (T max or VPD max ) (<0.05), and low correlation (<0.5) and low variable inflation factor (VIF) between predictor variables. Further, multiple regression models were retained only if they included a temperature-related variable (e.g., JFM precip + Jul-Aug precip would be rejected). See Table S2 in Supporting Information S1 for the top 10 models retained.
Of the remaining models, the model with the highest r 2 value was selected as the "best fit" model. This model was further tested to ensure any covariation between independent variables was minimal in the following manner: two versions of the model were run, with and without an interaction term between the two independent variables, Figure 2. The study region, including (a) mean annual maximum Moderate Resolution Imaging Spectroradiometer Normalized Difference Vegetation Index (NDVI) value, (b) land cover type, (c) elevation, (d) the month of peak NDVI, (e) PRISM August VPD max , and (f) PRISM total annual precipitation. For NDVI, VPD max , and precipitation, data is for the climatological (2003-2020) average. Major rivers (Colorado River, San Juan River, and Rio Grande) are shown in blue. and regression coefficients and corresponding p-values were compared. Moreover, to ensure that the results of the counterfactual VPD max , NDVI, and ANPP analysis do not substantially vary based on which model was chosen, the full counterfactual experiment (Sections 2.3-2.6) was also run using the top three models, which returned similar results (as defined by overlapping confidence intervals (CIs) for model coefficients and for attributable NDVI reductions).

Pixel-Wise Regression
After the best fit model was identified for the spatially-averaged data, we identified homogeneous climate and vegetation zones within which to fit models using the following procedure. First, the best fit model was used to calculate a spatial, pixel-wise ordinary least squares (OLS) regression, in which the model was fit for each pixel using the standardized local time series of data. While the pixel-wise OLS does not account for spatial relationships between pixels (compared to approaches such as Geographically Weighted Regression), the approach still allows for a crude estimation of spatial variation that affects model performance.
Next, to ascertain the drivers of spatial variation in the pixel-wise OLS, we extracted the r 2 for each pixel along with 13 potential explanatory variables related to land cover (land cover classification, % herbaceous, % shrub, % bare, 2001-2016 change in land cover classification), phenology and climate (peak NDVI month, SD NDVI, total annual precipitation, SD annual precipitation), and topography (slope, elevation, soil color, and proximity to major rivers). These variables were transformed (using square root and log transformations where applicable) to follow normal distributions. A multiple linear model was run with all potential variables with r 2 as the dependent variable. Significant (p-value < 0.05) variables were identified and assessed for multicollinearity. The variables which had low collinearity (r < 0.5) and jointly explained the most variance (r 2 > 0.9) were retained as the variables that best describe the spatial variation in NDVI sensitivity.
The region was stratified into unique zones based on these explanatory variables. For each zone, spatial means for each climate variable (P, T max , and VPD max ) were calculated, and linear models were fit to each zone.

Calculating Counterfactual T max
To estimate the increase in temperatures attributable to anthropogenic climate change, counterfactual T max (T max,cf ) was calculated by perturbing PRISM T max with CMIP6-derived attributable increases in T max (T max,att ) (E. L. . Note that while CMIP6 models have much coarser resolution than PRISM (∼100 vs. ∼4 km), no downscaling was used. Instead, as attributable temperature increases do not vary as significantly across the landscape compared to other variables, T max,att was calculated as a spatial mean across the full study region, and subtracted from spatially-averaged PRISM T max for each zone identified in Section 2.2 using the following procedure.
First, for each CMIP6 T max historical + SSP2-4.5 spatiotemporal cube (152 simulations across 25 models), spatial means across the full study region were calculated. Next, for each of the 25 CMIP6 models, model-specific means were taken across the simulations, and a multi-model ensemble mean was derived from the 25 individual model means. The multi-model ensemble mean was compared to PRISM to determine whether a bias correction was needed: CMIP6 T max biases remained constant throughout the time period (1895-2020), and as the counterfactual experiment relies on temperature differences, no bias correction was required.
Two methods were compared to derive T max,att from the CMIP6-based time series. In the first method, for each of the 25 models, 1850-1880 mean monthly T max were calculated across the historical + SSP2-4.5 cubes, and used as estimates of the PI T max . These PI T max values were subtracted from 10 years rolling means in the time series of each model ensemble to derive T max,att . The second method used the PI CMIP6 experiments: for each model, the monthly mean across the last 200 years of the PI experiments was subtracted from the ensemble members in 10-year moving windows. Both methods were compared and yielded similar results for T max,att . The first method was used to derive T max,att . Finally, to derive T max,cf , the range of the T max,att values was subtracted from PRISM T max .

Calculating Counterfactual VPD max
For each zone, counterfactual VPD max (VPD max,cf ) was calculated. VPD max is the difference between saturation vapor pressure (SVP max ) and actual vapor pressure (AVP max ), and SVP max is a function of T max (Equation 1, Daly et al., 2015). Actual SVP max was calculated from actual T max and was used with VPD max to derive actual AVP max . SVP max,cf was calculated from the 25 estimates of T max,cf . To isolate the effect of anthropogenically increased temperatures (and hold precipitation constant), we assume that AVP max does not change in the counterfactual as there is limited available water to evaporate over this region-indeed, observed trends in AVP for this region are small compared to SVP. Furthermore, the trend that does exist corresponds to anomalously low precipitation since 2000 (Ficklin & Novick, 2017). AVP max was subtracted from the 25 SVP max,cf estimates to derive estimates of VPD max,cf .
Furthermore, VPD max,cf was calculated for 1950-2020 in order to identify Time of Emergence (TOE), the date when anthropogenic forcing first significantly increased VPD max , after accounting for uncertainty. TOE was calculated as when actual (PRISM) VPD max is consistently above the upper bound of the CI of VPD max,cf . The CI was calculated using a bootstrapping procedure in which the annual mean and standard deviation of each year's attributable increase in VPD max (from the distribution of the 25 CMIP6 models) was used to generate a 10,000 element array for each year of potential VPD change values, and the 5th and 95th percentiles were used as the CI.

Calculating Counterfactual NDVI-Accounting for Potential Non-Linearity
Before estimating counterfactual NDVI (NDVI cf ), the potential for non-linearity between NDVI and VPD max was considered. It has been shown that, across large gradients of VPD, the relationship between VPD and plant productivity is non-linear (Grossiord et al., 2020;McDowell & Allen, 2015;Novick et al., 2016;Q. Zhang et al., 2019). Yet, across a small VPD max gradient, the relationship between measures of plant productivity and VPD max may be fairly represented with a linear regression. To test for non-linearities-specifically, whether the difference between actual and counterfactual VPD max in this study falls within the linear portion of the relationship between VPD max and NDVI-we explored several approaches. First, we fit a non-linear model, using a self-starting function which guesses its own coefficients. The model fit was comparable between the non-linear model and the best fit linear model (correlation between NDVI and NDVI est = 0.91 and 0.92, respectively, for the linear and non-linear models). We then tested for potential nonlinearity by substituting space for time (see Text S1 in Supporting Information S1). Each zone was split into cooler and warmer zones, ensuring that the mean difference in VPD max between cool and warm zones was at least as large as the maximum (2020) difference between actual VPD max and VPD max,cf . We then fit models for each zone and took the difference between the VPD max model coefficients. We then examined the CI of the coefficients, finding they overlap for each zone (Table S3 in Supporting Information S1). These results indicate that across the range of attributable differences in VPD max experienced by each zone, the influences of non-linearity are relatively small (see Text S1 in Supporting Information S1 for more information).

Calculating Counterfactual NDVI and ANPP
To calculate the decrease in NDVI attributable to human-induced warming, NDVI was first reconstructed using the model fit for each zone (NDVI est ). As this study focuses on the effect of anthropogenically-increased temperatures on NDVI, we derived counterfactual NDVI (NDVI cf ) using counterfactual VPD max while holding precipitation constant. As NDVI = f(VPD max ) is linear across the range of VPD max experienced by the study region (Section 2.5), and as the interaction term (VPD max * precip) is negligible (below, Section 3.1), we calculate NDVI cf = f(VPD max , precip) using a linear model. Therefore, NDVI cf was calculated using the same model intercept, precipitation value, precipitation coefficient, and VPD max coefficient, with the range of VPD max,cf z-scores, as: A bootstrapping procedure was used to estimate the confidence intervals for attributable changes in NDVI. Since precipitation was held constant, the annual changes in NDVI are a function of the uncertainty in the annual changes in VPD max (σ ΔVPD ) and the standard error term for the VPD slope term from Equation 4 (σ β1 ). As with T max,cf , to represent the first term, the annual mean and standard deviation of each year's attributable change in VPD max was used to generate an 18 years × 10,000 element array of potential VPD change values. Next, the NDVI-VPD slope and standard error values were used to generate an 18 × 10,000 array of potential slope values. These matrices were multiplied, and the resulting distributions used to generate confidence intervals.
Counterfactual annual net primary productivity (ANPP cf ) was then estimated for each zone by fitting models which relate NDVI to ANPP. Two ANPP products were considered-one from the Rangeland Management Production Service (RPMS) and one from GrassCast (Hartman et al., 2020;Reeves et al., 2020). RPMS includes all primary production while GrassCast is a function of solely herbaceous production (Hartman et al., 2020;Jones et al., 2021;Reeves et al., 2020). As certain livestock only eat herbaceous vegetation (e.g., cattle), while other livestock will graze on shrubs (i.e., goats), both products provide valuable information. Both products agree for low NDVI (Aug NDVI < 0.5), while GrassCast estimates were higher than RPMS estimates for high NDVI values. Therefore, as the study region includes many non-herbaceous pixels, and given the relatively strong agreement between both products at lower NDVI values (which includes all areas except high-precipitation zones in our study region), the RPMS equations were used to estimate changing ANPP. Note, RPMS provides ANPP geospatial data for download: these data were used to validate the suitability of the equations used by calculating ANPP using MODIS NDVI and comparing it to the RPMS ANPP data ( Figure S1 in Supporting Information S1). RPMS estimates ANPP (kg/ha) using two models: a quadratic and a log-based regression for high and low NDVI, respectively (Equation 3, Reeves et al., 2020).
For each zone, 2003-2020 NDVI est was first used to estimate ANPP est , then NDVI cf was substituted into the models to estimate ANPP cf . To estimate uncertainty in the ANPP cf estimates, the confidence intervals for NDVI cf estimates were carried over into the ANPP cf calculation.

Modeling NDVI
Of all the regressions, after screening for collinearity (removing predictors correlated beyond ±0.5), the model using standardized summer maximum VPD (VPD max ) and winter-summer precipitation as predictors was identified as the top model for predicting August NDVI (r 2 = 0.87) (Equation 4). The correlation between the two predictor variables was 0.48, and the VIF was 1.3, indicating low collinearity. Moreover, when run with an interaction term (VPD max * P), the interaction term was insignificant and the VPD max and precip coefficients remained largely consistent. The model indicates that NDVI responds strongly to winter-through-summer moisture supply (Jan-Aug precipitation), modulated by summer evaporative demand (Jul-Aug VPD max ).
August NDVI = −0.5 * Jul.Aug VPDmax + 0.6 * Jan.Aug Precip These slope terms in Equation 4 were highly significant (VPD max p-value = 1.6e −4 ; precipitation p-value = 2.5e −5 ; standard error for both coefficients = 0.1). While very simple, Equation 4 has important implications: since the magnitude of the slope coefficients are very similar, a one-sigma standardized anomaly shift in VPD max will have an impact similar to a one-sigma standardized anomaly shift in precipitation.

Examining Spatial Variation in Model Fit and Climatic, Topographic, Land Cover Variables Across the Region
While Equation 4 best describes August NDVI for the spatial aggregate of the region, the sensitivity of NDVI to climate variables will vary across the study region due to the influences of topography, climate, land cover type, and other variables (Figure 2). The study region encompasses much of the Colorado Plateau, and is climatically and topographically heterogeneous, ranging from grasslands to forests (Figure 2b), across elevations (Figure 2c), and variable temperature (and by extension VPD) and precipitation regimes (Figures 2e and 2f). As such, there may be areas that are more sensitive to increasing VPD than others. Figure 3 shows the spatial patterns in the r 2 , VPD max p-values, and VPD max and precipitation coefficients from the pixel-wise OLS using Equation 4. The model fit (r 2 ) is strongest in the middle and southeast of the study region ( Figure 3a), that is, in New Mexico. These areas are mostly shrublands and grasslands which receive relatively low annual precipitation and experience mid-to-high VPD max (Figures 2b, 2e, and 2f).
Visually examining areas with the strongest r 2 values, two sub-regions emerge. Significant VPD max p-values (p ≤ 0.05, Figure 3b) accompanied by strong r 2 values (>0.5, Figure 3a) are mostly located in the southeastern region, concentrated in New Mexico. These areas are also characterized by highly negative VPD max coefficients ( Figure 3d). This indicates that the highest-risk areas for increases in summer VPD max are concentrated in the southeast of the region. Interestingly, these locations tend to be relatively cooler, mid-elevation areas (Figures 2c and 2e) covered with shrublands ( Figure 2b). Conversely, the model fit is weakest in the northeastern, eastern, and far southwestern parts of the region. These poor model fit areas often correspond with high precipitation, forested areas. Note that while we find numerous significant VPD max p-values in the forested areas of Colorado (upper-right) and Arizona (bottom-left), these areas have lower r 2 values.
Moreover, VPD max is not a significant predictor for the northwestern and central parts of the region, concentrated in the Colorado River Basin and along the San Juan River. These are the hottest and driest areas of the study region (Figures 2e and 2f), in major river basins, and with low climatological NDVI values (Figure 2a). Moreover, these regions have the largest precipitation coefficients, demonstrating the dominant importance of precipitation on NDVI. This indicates that for the center of the study region, relatively low precipitation may mask any effect of VPD on NDVI. Moreover, at very low values of NDVI (Figure 2a), satellite retrievals of NDVI are also likely to be influenced by bare soil reflectance, especially in areas with red soils, which is not uncommon in the Four Corners area.

Examining Controls of Spatial Variations and Creating Unique Zones
Spatial variation in model skill (Figure 3) appears to be primarily linked to four variables: the phenology of vegetation (the seasonality of "greenup" or peak NDVI month), total annual precipitation, land cover type, and whether the pixel has experienced a change in land cover type. While NDVI peaks in August for most of the region, in the northwest, NDVI peaks in winter or spring months (Figure 2d). In this region, we find frequent non-significant VPD max coefficient p-values (Figure 3b). These areas interestingly are close to major rivers, suggesting that the differences in phenology may be explained by other moisture inputs to the region, such as snowmelt and rainfall runoff from the tributaries to the Colorado River. In further analysis, since we do not account for runoff and snowfall, all non-summer-peak NDVI pixels were excluded from this study. Total annual precipitation is also related to model fit. High r 2 values are concentrated in areas with precipitation <500 mm/ year (Figure 2f), while summary statistics indicate a secondary break near 250 mm/year. The highest r 2 values tend to arise when annual total precipitation falls between 250 and 500 mm/year.
The model fit also varies between land cover types, and based on whether a pixel experienced a change in land cover type. Forests have significantly lower r 2 values than shrublands and grasslands. We therefore also exclude non-shrubland and non-grassland areas, given both the relatively poor model fit for forested areas, and because we are primarily interested in impacts on rangelands and on those communities who use rangeland resources. When land cover types are restricted to grasslands and shrublands, land cover still remains significant.
Finally, a few pixels (<1%) experience a change in major land cover classification (e.g., forest to grassland). As our study does not account for land cover change, these pixels were also excluded from the study.
Overall, these results indicate that, as expected, sensitivity to rising VPD varies between different vegetation types and with increasing moisture inputs (Grossiord et al., 2020;Novick et al., 2016;Rao et al., 2022). After selecting all summer peak-NDVI rangeland pixels, the study region was stratified into six unique zones based on land cover type (grasslands vs. shrublands) and total annual precipitation (<250, 250-500, and >500 mm total annual precipitation) (Figure 4). Importantly, these zones exhibit similar phenology ( Figure S2 in Supporting Information S1). Moreover, the series of potential models (Section 2.2) were fit again to these unique six zones, and Equation 4 was again in the top three models (defined by r 2 values).

Counterfactual Temperature and Vapor Pressure Deficits
For each zone depicted in Figure 4, counterfactual T max (T max,cf ) were estimated ( Figure 5). The counterfactual values represent estimates of what T max would have been without human-induced warming. The calculated ensemble average attributable increase for Jul-Aug T max for 2010-2020 was ∼1.5°C compared to the preindustrial baseline, comparable to the results of A. P.  for this region. This attributable increase in warming is substantially larger than the interannual standard deviation of Jul-Aug T max in this region.
Next, counterfactual VPD max (VPD max,cf ) for 1950-2020 was derived using T max,cf ( Figure 6). The TOE of attributable increases in VPD max was estimated as 2005 for all zones. TOE is when actual VPD max is greater than the upper bound of the CI of the range of estimates of VPD max,cf , or, as shown in red in Figure 6, when the CI of ΔVPD max is above zero. These results imply that due to human-induced warming, after accounting for uncertainty from the range of estimates from the 25 CMIP6 models, since 2005, VPD max has been higher than it would have been without that warming. Moreover, since 2000, the attributable increases in VPD max have been particularly pronounced-specifically, on average, they have been greater than 1SD of pre-1980 actual VPD max (blue dotted line in Figure 6). Therefore, for most of the study period (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020), actual VPD max was consistently and substantially higher than estimated VPD max,cf , including uncertainty from model fit (CI shown in ribbons). At the end of the study period, the largest differences between actual VPD max and VPD max,cf are about 5 hPa in the low precipitation zones, as these zones are warmer and the relationship between temperature and SVP is non-linear/exponential. Furthermore there will be less latent heat flux and more sensible heat flux in low precipitation zones as compared to other regions. Conversely, attributable differences for high precipitation zones are about 3 hPa. These results imply that modest increases in air temperatures can have very large impacts on VPD, with the influence of climate change typically being about twice the inter-annual standard deviation by the end of the study period.

Estimating Attributable Reductions in NDVI and ANPP
To calculate the decrease in NDVI attributable to anthropogenic warming, NDVI was first reconstructed for each zone (NDVI est ) using Equation 4. All six models performed well, with r 2 values ranging from 0.71 to 0.91, and significant VPD max and precipitation slope coefficients. For low and mid-precipitation zones, the model fit was strong (r 2 0.84-0.90), yet was lower for high precipitation zones (r 2 0.71-0.77). Note the sample sizes for high precipitation zones are relatively small (number of pixels: n = 199 grassland; n = 390 shrubland). Moreover, the middle precipitation zones (n = 635 grassland; n = 3,182 shrubland) have the largest VPD max coefficients (−0.43/−0.47), indicating an important temperature control over NDVI. Conversely, the low and high precipitation zones have the largest precipitation coefficients. This indicates different climatic controls over NDVI for each precipitation regime.
Next, counterfactual NDVI (NDVI cf ) was derived for each zone (Figure 7). The results imply that for the full study period, including uncertainty in the linear model fit and from the suite of CMIP6-based VPD cf estimates, human-induced temperature increases have led to reductions in NDVI. Notably, the attributable reductions in NDVI increased during the study period, particularly during the first decade and after 2018, when there was both observed substantially elevated VPD and low precipitation. The mean estimated attributable reduction in NDVI in 2020 ranges from 0.35 to 0.85 standard deviations of the interannual variability in NDVI. There are patterns in NDVI cf across precipitation regimes-the strongest reductions in NDVI corresponding with increases in VPD are in mid precipitation shrublands and grasslands, with mean estimates of 0.78-0.85 SD reductions in shrublands and grasslands, respectively. These areas with the strongest responses correspond to pixels with the greatest VPD max coefficients (Figure 3d).
While NDVI is a useful proxy for vegetation productivity, ANPP is a more direct and more commonly used measure of forage abundance and overall plant productivity. As with the attributable reductions in NDVI, largely by design of the equation relating NDVI to ANPP, the largest standardized reductions in ANPP are found for mid-precipitation shrublands and grasslands (Table 1). In these regions, the magnitude of attributable reductions in ANPP is notable: ∼1 SD at the end of the study period, or around 50% of the observed anomaly (Table 1). Across the full study period, for these mid precipitation regimes, these reductions range from 37 to 75 lbs/acre of ANPP in shrublands and 49-95 lbs/acre of ANPP in grasslands.
These reductions were further estimated at the census tract-level. Figure 8 depicts the reductions in ANPP attributable to human-induced warming in terms of standard deviations of ANPP per census tract. The largest attributable reductions in ANPP are mostly found in New Mexico, with some census tracts reaching greater than 1SD of reductions in ANPP.

Discussion
In many monsoonal areas, including the greater Four Corners region, vegetation greenness typically peaks during the hottest months corresponding with monsoonal rains. Our results indicate that peak greenness and overall productivity is related to winter-summer precipitation and summer VPD. Previous studies have identified winter precipitation as important for plants in this region, including grasses and shrubs, as it influences available soil moisture for later months (from Garfin et al., 2013, pp. 154-156). Furthermore, for the Four Corners region, hot and dry conditions in the early summer (June-July), are characterized by extremely high VPD. It is during this time, shortly before peak greenness, as well as in the summer (July-August) if the monsoonal rains are delayed or fail to materialize, that there is particularly high potential sensitivity to high VPD (Weiss et al., 2009;. This sensitivity is magnified by the non-linear relationship between warm air temperatures and SVP (Equation 1). From year-to-year, summer T max and precipitation also covary: therefore, low precipitation years can further increase temperatures due to land-atmospheric coupling and surface energy balances. As such, the combined influences from failed monsoonal rains and high temperatures can significantly reduce vegetation productivity on rangelands.
The decades of our study period have experienced low precipitationaverage Jan-Aug precipitation from 2003 to 2020 was −0.4 standard  Note. 2020 ANPP for each zone are depicted (column 1), including anomalies (from 2003 to 2020 mean) (column 2). ANPP decreases attributable to human-induced warming ("Attributable Decreases") are shown in terms of raw reductions, as a percent of the observed anomaly, and in standard deviations of ANPP. Confidence intervals for the attributable decreases are shown in parentheses, representing the uncertainty from the NDVI cf estimates shown in Figure 7. deviations below the long-term mean ( Figure S4 in Supporting Information S1). Even in this lower precipitation era, the results presented here indicate that increased temperatures still have had substantial influences on VPD max (Figure 6), NDVI (Figure 7), and ANPP (Figure 8). At regional scales, the influence of VPD is similar in magnitude to that of precipitation (Equation 4). Hence, increased VPD levels are having a persistent and damaging influence. Moreover, the detrimental effect of increased VPD on NDVI is detectable despite elevated CO 2 : while elevated atmospheric CO 2 concentrations increase plant water use efficiency and may offset negative effects of increasing VPD in certain sites (e.g., Gonsamo et al., 2021;Yang et al., 2018), studies have indicated that in many water-limited sites, decreases in available moisture can counteract gains from increased CO 2 , resulting in decreased vegetation productivity (e.g., Dannenberg et al., 2022;Ficklin & Novick, 2017;Jiao et al., 2021;Wang et al., 2022;W. Yuan et al., 2019). Increasing water use efficiency may be offset by reductions in soil moisture due to increased evapotranspiration, and evaporation (not transpiration) will be important in arid regions. Overall, these findings correspond with studies that have identified recent increases in temperature and VPD as significant drivers of reduced vegetation productivity in the Southwest (e.g., Dannenberg et al., 2022;Jiao et al., 2021;Rao et al., 2022).
Significant attributable reductions in NDVI are identified for all regions and years in the study. While there is substantial variation in the magnitude of the impacts across zones, the overall trend from 2003 to 2020 indicates a sizable increase in attributable reductions over the past two decades. Comparing zones, the driest areas in the study region-largely found in Arizona-have the smallest decrease in NDVI attributable to human-induced warming. These are also the areas with the hottest temperatures and lowest climatological NDVI: the average peak NDVI for these zones is ∼0.22, while peak NDVI for bare ground can be ∼0.1 (Huete et al., 2002). These zones are highly water limited. It is likely that for these zones, precipitation, and therefore soil moisture, is low enough that, barring increases in precipitation, increased temperatures only have a marginal effect. Indeed, studies indicate that in the most arid areas, soil moisture/precipitation is the dominant control over vegetation growth, while VPD is the strongest control in higher soil moisture areas (Famiglietti et al., 2021;Novick et al., 2016). The areas with the largest attributable decrease in NDVI are concentrated in New Mexico, at low-to-mid elevations, which exhibit slightly cooler temperatures and higher precipitation than the driest zones. These areas may receive enough precipitation to have a more pronounced sensitivity to increased VPD, and are thus the highest-risk areas to increasing VPD. These findings correspond with Novick et al., 2016, in which VPD was found to have a stronger control over stomatal conductance in wetter sites than in drier sites. These areas have also been identified as having high plant-water sensitivity due to both plant and soil traits, in which increasing VPD corresponds with decreased live fuel moisture content and higher wildfire risk (Rao et al., 2022).
Our results also indicate that attributable increases in VPD in the mid-precipitation areas correspond with significant reductions in ANPP. In mid-precipitation areas, we find that reductions in NDVI attributable to human-induced warming roughly correspond to −25 to 60 lbs/acre reduction of ANPP in shrublands and −40 to 100 lbs/acre of ANPP in grasslands, accounting for ∼50% of ANPP deficits by 2020. To put this estimate in context, in 2018, exceptional drought covered much of the greater Four Corners region and soil moisture conditions were "very poor" (USDM, 2018). The average ANPP in 2018 for the mid-precipitation zones was ∼−72 lbs/ acre below the 2003-2020 mean, while our counterfactual ANPP analysis suggests that without human-induced warming, ANPP in 2018 would have been only around −12 lbs/acre below the 2003-2020 mean. In other words, human-induced warming contributed to approximately three-fourths of the observed ANPP reductions in 2018 (i.e., 60 lbs/acre). This corresponds to the findings of Dannenberg et al., 2022, who showed that while half of the 2020 GPP anomaly was due to observed reduced soil moisture, nearly half was due to increased VPD. The results presented here, however, add important context to Dannenberg et al. 2022-namely that human-induced warming has been largely responsible for such impacts.
The reductions in ANPP on both grasslands and shrublands are significant in terms of fodder-while cattle primarily consume grasses, sheep will eat forbes, some shrubs, and broad-leaved plants (Milchunas, 2006). In particular, where these large reductions in NDVI/ANPP coincide with high density livestock, these areas may be experiencing impacts on livestock ( Figure S3 in Supporting Information S1). Our counterfactual ANPP analysis can be further understood in the following context: if a rangeland is at full stocking capacity (i.e., the livestock consume the maximum amount of forage they can without damaging the productivity of forage), and assuming that half of ANPP is edible by cattle, the reductions in ANPP attributable to anthropogenic climate change will correspond to either additional fodder purchases or decreases in herd size. To sustain the same herd size, taking the price of alfalfa at $240/ton, the 2020 reductions in ANPP in mid-precipitation grasslands result in an additional ∼$6/acre. Alternatively, if no additional fodder is purchased, to maintain the same grazing demand, the herd size would need to be reduced by ∼8%. Our counterfactual ANPP analysis provides an indication of the effect of anthropogenic warming on vegetation productivity in units that are of greater relevance to land users, as compared to indices like NDVI. While ANPP does not perfectly correspond to actual forage availability, it is a closer estimate than NDVI (Jones et al., 2021;Zimmer et al., 2021).
While a precise economic impact assessment is beyond the scope of this study, our results indicate that human-induced warming accounted for a large portion of recent reductions in four corners vegetation production. In 2020, these arid conditions contributed to a $5.1 billion dollar disaster in the western United States (NOAA, 2022). According to NOAA "There were considerable crop and livestock impacts across the West and Central states from both the persistent heat and increasingly dry conditions. The combined drought and heat also assisted in drying out vegetation across the West that contributed to the Western wildfire potential and severity." In 2018, according to NOAA, extreme drought conditions across the Four Corners area led to $3.6 billion dollars in impacts, with "the agriculture sector impacted across the affected states …. Ranchers have also be [en] forced to sell-off livestock early in some regions due to high feeding costs." The research presented here suggests that a substantial portion of these billions of dollars of damages were almost certainly linked to large human-induced increases in VPD and air temperature.
Moreover, since we do not attribute further increases in temperature due to land-atmospheric coupling, our study may yield conservative estimates of attributable increases in VPD max . While the experiments in many land surface models in CMIP experiments account for land-atmospheric coupling, some studies have indicated that the coupling in the models is too weak, or the latent heat flux is too high (Mueller & Seneviratne, 2014;K. Yuan et al., 2022). This would provide conservative estimates of heating, contributing to a conservative estimate of attributable increases in VPD max . The strong VPD sensitivity indicated by our results-namely the model fit between actual NDVI, VPD, and precipitation-supports the possibility that there could be stronger-than-modeled positive feedbacks between drought, heat, and latent heat flux/evapotranspiration from the land surface.
Finally, the results of this study may be of interest to regions outside of the Southwestern U.S. There are many other semi-arid regions of the world with similar climates. Some areas, such as East Africa and southern Madagascar, have experienced similar repetitive droughts, extreme temperatures, and anomalously low NDVI. Since this study only requires precipitation, VPD, and NDVI, this methodology may be applicable to other semiarid regions dominated by grasslands and shrublands with more limited data availability, and could also be leveraged in forecasting applications, including to help identify areas at greatest risk for further desiccation due to increasing temperatures.
In conclusion, our results emphasize how a modest 1-2°C warming can dramatically increase VPD and SVP, especially in hot areas as the SVP equation is non-linear (Equation 1). Furthermore, areas or seasons with low average AVP will see little relief in Clausius-Clapeyron-type responses-i.e., if relative humidity remains the same, warming air temperatures will do little to increase AVP, because a ∼7% per °C influence will have little influence when and where average AVP is low. When and where VPD and precipitation have similar levels of influence (Equation 4), we will likely see persistent and substantial reductions in vegetation productivity.

Data Availability Statement
CMIP6 model output for pre-industrial, historical, and SSP2-4.5 experiments were accessed from the Lawrence Livermore National Laboratory (LLNL) ESGF node (https://esgf-node.llnl.gov/search/cmip6/). MODIS NDVI data are available from Didan and Huete (2015). PRISM temperature, VPD max , and precipitation data is available from the PRISM Climate Group at Oregon State University (PRISM, n.d.). Land cover data is available from the Multi-Resolution Land Characteristics (MRLC) Consortium (Dewitz and USGS, 2021;Rigge et al., 2020). Soil color and characteristics were retrieved from the Soil Survey from the USDA Natural Resources Conservation Service (USDA NRCS, https://nrcs.app.box.com/v/soils). Finally, data and scripts used to create figures in this manuscript may be accessed via GitHub (E. L. Williams, 2023).