Regional differences in the response of California’s rangeland production to climate and future projection

Rangelands support many important ecosystem services and are highly sensitive to climate change. Understanding temporal dynamics in rangeland gross primary production (GPP) and how it may change under projected future climate, including more frequent and severe droughts, is critical for ranching communities to cope with future changes. Herein, we examined how climate regulates the interannual variability of GPP in California’s diverse annual rangeland, based on the contemporary records of satellite derived GPP at 500 m resolution since 2001. We built Gradient Boosted Regression Tree models for 23 ecoregion subsections, relating annual GPP with 30 climatic variables, to disentangle the partial dependence of GPP on each climate variable. The machine learning results showed that GPP was most sensitive to growing season (GS) precipitation, with a reduction in GPP up to 200 g cm−2 yr−1 when GS precipitation decreased from 400 to 100 mm yr−1 in one of the driest subsections. We also found that years with more evenly distributed GS precipitation had higher GPP. Warmer winter minimum air temperature enhanced GPP in approximately two-thirds of the subsections. In contrast, average GS air temperatures showed a negative relationship with annual GPP. When the pre-trained models were forced by downscaled future climate projections, changes in the predicted rangeland productivity by mid- and end of century were more remarkable at the ecoregion subsection scale than at the state level. Our machine learning-based analysis highlights key regional differences in GPP vulnerability to climate and provides insights on the intertwining and potentially counteracting effects of seasonal temperature and precipitation regimes. This work demonstrates the potential of using remote sensing to enhance field-based rangeland monitoring and, combined with machine learning, to inform adaptive management and conservation within the context of weather extremes and climate change.


Background
Rangelands cover more than 30% of the global land area and are biologically diverse landscapes that include grasslands, shrublands, woodlands, wetlands, and deserts (Lund 2007). Changes in rangeland production affect many important ecosystem services, including carbon storage, nutrient cycling, water resource protection, biodiversity, and wildlife habitat (Li 2012, Roche et al 2015. Gross primary production (GPP) is a widely used metric for quantifying photosynthesis rate at the ecosystem scale. Recent studies have shown increasing interannual variability in GPP globally, with 57% of the variability contributed by rangeland-dominated ecosystems (Zhang et al 2016). In semi-arid rangelands, precipitation is, in most cases, the primary controller of GPP fluctuations (Jin and Goulden 2014, Zhang et al 2016. With projected increasing precipitation variability, the interannual variation in GPP is also predicted to increase (Zscheischler et al 2014). Therefore, monitoring GPP and understanding its variability may be a useful benchmark indicator of ecosystem response to climate change.
In California, annual rangelands, which include open grasslands and woodlands dominated by an understory of herbaceous annual plants, are the primary forage source for the State's $3 billion livestock industry (California Department of Food & Agriculture 2015). These annual rangelands comprise approximately 4 million hectares, encompassing portions of the Central Valley, the Coast Ranges and Sierra Foothill Region (FRAP 2017). California rangelands are highly variable, in part, because the plant community phenology, consisting of annual grasses and forbs, is dynamically linked to seasonal weather patterns (Stromberg et al 2007). For example, the growing season (GS) is initialized by fall rains (September-December) exceeding 1.25-2.5 cm during a seven-day period, which promotes germination. Plant growth is slow in winter (December-February), due to low air temperatures and low solar radiation. In late winter/early spring (February-April), rapid growth begins as air temperature increases (>7.2 • C) and solar radiation intensifies. Finally, in mid-late spring (April-June), plant production reaches peak biomass once root zone moisture is depleted. Thus, California rangeland GPP varies significantly seasonally and from year-to-year (Jin and Goulden 2014), in response to short-term changes in climate and solar radiation. The large gradient of topography and long term climate across California's diverse landscape also shapes the spatial patterns of rangeland GPP (Liu et al 2021) and potentially a differentiating sensitivity to climate change.
Over the past decades, researchers have relied on regression and/or correlation analyses to examine drivers regulating variability in annual GPP and plant biomass production on rangelands at different scales. Taking advantage of the high correlation between absorbed photosynthetically active radiation and GPP, Jin and Goulden (2014), for example, reported GPP was positively related to annual precipitation, with higher sensitivity in grasslands than woody shrublands in California. A few other studies also reported positive response of rangeland GPP to annual precipitation amount in Northern China (Guo et al 2015), North America (Wu and Chen 2012), and across the entire globe (Huang et al 2016).
Annual GPP can be further regulated by other climate variables. Seasonal distribution of rainfall events has a significant impact on rangeland GPP by regulating the active GS (Xu and Baldocchi 2004). Warming air temperature was found to reduce GPP in northern Eurasia's southwestern grasslands, although it was the primary reason for recent GPP increases in northern Eurasia's shrublands (Dass et al 2016). In California, rangeland production is predicted to increase as air temperatures become warmer in the San Francisco Bay Area as early as mid-century (Chaplin-Kramer and George 2013); however, another study predicted a 14%-58% decrease in rangeland production across California using a precipitation-based model (Shaw et al 2011). Previous small scale studies showed that variations in GPP were related to the combined effects of various climatic variables, and the relationship between GPP and climatic variables may vary across different locations. Thus, to fully understand variations in GPP and predict future changes, it is necessary to build individual models with respect to different regions and consider a suite of relevant climate variables, including precipitation amount, seasonal precipitation distribution, air temperature, and solar radiation.
A long-term record of annual GPP at moderate resolution is critical for a more robust analysis of interannual variability in GPP and for building data-driven models. Advances in both remote sensing technologies and data analysis are providing novel, cost-effective opportunities to obtain data at management-relevant scales for terrestrial-ecosystem research, which has otherwise relied on traditional statistical methods and field-based monitoring that require substantial time and labor resources (Briske et al 2017). The Moderate Resolution Imaging Spectroradiometer (MODIS) has been collecting terrestrial observations for more than 20 years on a daily basis since the launch of Terra satellite in 1999 (Savtchenko et al 2004). GPP products derived from MODIS satellite images at 500 m spatial resolution are available at an eight-day interval for the entire globe (Running et al 2019). Machine learning models have also become more interpretable while keeping the powerful non-linear, multi-parameter model fitting functionalities (Molnar 2019). These advances are key because traditional statistical models lack the ability to disentangle the interactions among variables (Lee et al 2020, Liu et al 2021.
In this study, we aimed to improve our understanding of temporal variability in rangeland annual GPP and assess how its sensitivity to changes in weather and climate may vary across different ecoregions in California. The machine learning models were built to learn the complex climate impacts on GPP, by leveraging 20 years of contemporary records on weather and remote sensing based GPP. Using the developed models, we further predicted and evaluated how rangeland annual GPP would change under different climate projections by mid and end of century over California's diverse landscape. We hypothesized that the primary climatic variable regulating California annual rangeland GPP varies by regions. This spatial variability results from the influence of other climatic and environmental variables that become more pronounced as precipitation increases, especially when water is no longer the primary limiting resource.

Study area
Our study focused on California annual rangelands, spanning 33 • -41 • N and 118 • -124 • W, with a wide range of microclimates resulting from complex geographic and topographic diversity. The study area encompasses four major ecoregions as delineated by the Ecological Classification and Mapping Task Team (Cleland et al 2007): Northern California Coast Ranges (NCCR), Northern California Interior Coast Ranges (NCICR), Central Coast Ranges (CCR), and Sierra Nevada Foothills (SNF) ( figure 1(a)). Annual rangelands feature a Mediterranean climate with hot, dry summers and mild, moderately wet winters. The active GS of annual rangelands coincides with the rainy season: germination typically ensues in late fall (November), and senescence occurs after the rainy season ceases in late spring or early summer (April-May).
Precipitation varies considerably across ecoregions, with long-term-means (LTMs) ranging from 160 mm yr −1 in the CCR to 2290 mm yr −1 in the NCCR ( figure 1(a)). Coastal regions are relatively cooler in the summer and receive higher amounts of precipitation than inland regions, especially those in rain shadows. Across the study area, air temperature increased from the coast to inland, north to south, and from valleys to mountains. Due to the significant heterogeneity of the landscape, the four major ecoregions are further divided into 23 ecoregion subsections based on similar ecological potentials, determined by many biotic and environmental factors such as climate, physiography, soils, hydrology, and potential natural communities (figure 1(a)) (Cleland et al 2007). We hereafter refer to ecoregion subsections by their parent ecoregion abbreviation followed by a unique number as specified in figure 1(a).

Contemporary satellite and climate data
We obtained the version 6 gap-filled GPP product (MYD17A2HGF) at 500 m resolution every eight days from MODIS Terra and Aqua during 2000-2019 (https://earthdata.nasa.gov/). MODIS provides frequent multispectral moderate-resolution (250-500 m) observations, revisiting the globe once to twice per day since 2000. We calculated water year (October-September) (U.S.Geological Survey 2016) cumulative GPP from the 8 day product (Running et al 2019). A LTM annual GPP layer (820 ± 283 g cm −2 yr −1 ) was derived from 2001-2019 as a measure of spatial variation in forage production (figure 1(c)). Annual rangeland GPP varied significantly from year to year as shown in figure 1(d), with a statewide mean ranging from 1017 g cm −2 yr −1 in 2005-750 g cm −2 yr −1 in 2008 (see detailed discussion in supplemental materials section 1).
Meteorological data during 1986-2019 came from Daymet hosted at the Oak Ridge National Laboratory Distributed Active Archive Center (Thornton et al 2016). The gridded daily surface weather parameters at a 1 km resolution were interpolated from more than 8000 meteorological stations, based on a digital elevation model (Thornton et al 2016). We aggregated daily precipitation, min/max/mean air temperatures, and solar radiation to monthly and seasonal averages, resulting in a total of 30 climatic variables to be explored (as detailed in table S1). To quantify seasonal precipitation variation between the driest and wettest months during the GS, a Precipitation Concentration Index (PCI) was calculated as the ratio of the sum of squares of the wet season monthly precipitation (November-May) over squared water year precipitation multiplied by a scaling factor of 58 (Oliver 1980, Sloat et al 2018.

Future climate projection data
For future climate, we obtained daily climate projections downscaled for California from Cal-Adapt Data Server (http://albers.cnr.berkeley.edu/data/scripps/ loca/). This dataset is bias-corrected and downscaled to a resolution of 1/16 • (∼6 km) from 32 global climate models (GCM) using the localized constructed analogues statistical method (Pierce et al 2018).
The data include both the historical simulations during 1950-2005 and future climate projections during 2006-2099 under two emission scenarios that have been widely used (e.g. by the California's climate change assessments) for climate related policy decisions (IPCC 2013). Representative Concentration Pathway (RCP) 4.5 (Thomson et al 2011) is an intermediate scenario where greenhouse gas emissions peak around 2040 and then decline (IPCC 2013). A worst-case scenario, RCP 8.5, represents rising emissions throughout the 21st century (Thorne et al 2017). In this study, we used projections from four GCMs that have demonstrated good performance in reproducing California's historical climate (Pierce et al 2018): HadGEM2-ES, MIROC5, CNRM-CM5, and CanESM2.
For each model and scenario projection, we acquired daily mean precipitation, air temperature (minimum, maximum, and mean), and solar radiation layers. To match the contemporary record, a bias adjustment was applied to the GCM data, based on 20 years of overlapping time periods between Daymet and GCM historical simulations. We first calculated the differences between Daymet and GCM data derived for each climate variable (see table S1) for each overlapping year during 1986-2005. LTM offsets were then calculated, and applied to the original GCM projection to derive bias-adjusted future climate data.

Machine learning GPP model development
Using the contemporary records of satellite derived GPP and Daymet weather during 2000-2019, we built Gradient Boosted Regression Trees (GBRT) models for each ecoregion subsection to investigate the complex associations between the interannual variability in annual rangeland GPP and climatic drivers. The annual GPP layers were resampled to 1-km resolution with a linear method to match the Daymet grid size. The LTM GPP at 1 km was included as an additional independent variable (table S2: LTMGPP, top one or two variables), in order to account for spatial variability in GPP caused by topography and soils.
GBRT is one of the most popular ensemble machine learning methods that combine outputs from many individual trees to overcome the potential poor performance and reduce the risk of overfitting faced by single learner methods (Friedman 2001). Instead of building trees independently, it fits regression trees sequentially, with each tree aimed at minimizing prediction residuals of its predecessors (Friedman 2001). GBRT is a powerful tool for discovering and representing nonlinearities and complex interactions among predictors (Elith et al 2008), because the hierarchical tree structure inherently models the dependence of the response variable to one predictor given the values of inputs higher in the tree. For each ecoregion subsection, a GBRT model was trained with 70% randomly selected pixel-years, using key variables selected from table S1.
Feature selection was conducted to remove redundant climatic variables and reduce potential overfitting issues. We first iteratively eliminated variables with respect to their correlation to other variables (i.e. when two variables are highly correlated (r > 0.8), removing the variable having the lower distance correlation (Székely et al 2007) with annual rangeland GPP. We further removed variables that had low predictive power using the Recursive Feature Elimination with cross-validation technique (Guyon et al 2002). Finally, a permutation importance-based feature selection was performed to further reduce the number of variables. Permutation importance is computed as a decrease in the model performance (for example, measured by R 2 or Root Mean Squared Error (RMSE)) when a variable is permuted. We sequentially permuted variables with the least permutation importance score and fitted models with the most important three to nine variables. For the purpose of balancing model accuracy and complexity, we kept the top seven variables to build ecoregion subsection-specific machine learning models (figure S1).
We investigated how each of the key climatic variables affects interannual rangeland GPP variation using partial dependence plots from the GBRT models (Friedman 2001, Friedman andMeulman 2003). Partial dependence function shows the response to one target predictor after accounting for the average effects from all other predictors, by marginalizing the model over the distribution of other independent variables (Molnar 2019). It is therefore used here to visualize and disentangle the interactive response of GPP to a suite of climate variables. The y-axis of partial dependence plots represents the difference between the marginalized prediction and mean GPP for each value of the predictor on the x-axis. For example, a wide range on the y-axis indicates strong sensitivity of GPP to the target predictor when other confounding factors are excluded.

Future GPP projection and attribution
To predict future rangeland GPP, we applied the trained GBRT models, built specifically for each ecoregion-subsection, to bias-adjusted climate projections from 2040-2099 as described in 2.3. For a consistent comparison, GPP estimations by the same models driven by the Daymet climate data during 2001-2019 were used as a contemporary baseline. The results were summarized by taking the means over three-time windows: contemporary (2001-2019), mid-century (2041-2059), and end of century (2081-2099), under two emission scenarios.
We further quantified changes in GPP by midcentury and end of the century caused by each type of climatic variable (e.g. precipitation amount, precipitation distribution, minimum air temperature, mean air temperature, maximum air temperature, and solar radiation). For each type of variable, the attribution was performed by comparing the difference between GPPs predicted for the contemporary period and GPPs predicted with the future climate simulation for the target type of variables while keeping all other variables as the contemporary climate. Differences between this newly predicted GPP and the original contemporary GPP indicate the extent to which changes in GPP are attributed to changes in precipitation, air temperature, and solar radiation.

Sensitivity of interannual rangeland GPP to climate
For each of the 23 subsections, we chose the top seven climate variables (see table S2) to model annual rangeland productivity at 6-km, based on GBRT modelbased feature selection as detailed in section 2.3. Predicted GPP was shown in good agreement with the MODIS GPP product when tested against the validation dataset, with an R 2 of 0.97, RMSE of 30.7 g cm −2 yr −1 , and a mean absolute percentage error of 4.1% ( figure S2).
Overall, precipitation-related variables were found important for 23 subsections for capturing interannual variation in rangeland productivity (table S2 and figure 2). Specifically, at least one variable representing precipitation amount was selected for all subsections, either for individual winter and/or spring months or seasonal cumulative amount, depending on the region (figure 2(a) and table S2). Over eight subsections, precipitation amount during the GS enhanced annual GPP before reaching a plateau at 300-700 mm yr −1 , depending on the subsection (figure 3(c)); production in the drier subsections of the central coast were most sensitive, with a decrease by more than 200 g cm −2 yr −1 when GS precipitation varied from 400 to 100 mm yr −1 .
PCI was critical in 21 subsections, except for the two driest subsections on the central coast. Partial dependence plots from subsection-specific GBRT models showed that annual GPP decreased rapidly before the curves flattened at a PCI equal to 13 ( figure 3(a)). This suggests that years with more evenly distributed precipitation (i.e. lower PCI values) had higher GPP than years with less evenly distributed precipitation (i.e. higher PCI values), when everything else was held constant.
Winter monthly minimum air temperature (November-February), was the third important variable category, followed by mean air temperature, selected by 19 and 9 subsections, respectively (figure 2). T max and solar radiation during GS were important in less than five subsections, mostly in north coast interior and central coast (table S2).
Air temperatures (T min , T mean , and T max ) showed different controls on GPP. Warmer January minimum air temperature enhanced annual GPP as it warmed from 0 • C to 6 • C ( figure 3(b)). In contrast, Figure 2. Prevalence of selected (a) climate variables and (b) climate variable groups, summarized among the 23-subsection specific models. For example, an appearance in 21 subsections for the precipitation concentration index (PCI) was selected as one of the seven top variables in 21 subsection models. Tmin, tave/tmean, and tmax stand for minimum, average, and maximum air temperature. Cumuprcp stands for cumulative precipitation. Srad stands for shortwave solar radiation. See table S1 for a full list of variable abbreviation definitions. mean November air temperature and GS maximum air temperature showed a negative impact on annual GPP (figures 3(d)-(f)). For example, as November mean air temperature increased from 7 • C to 14 • C, GPP decreased from 10∼40 g cm −2 yr −1 above the LTM (during 2001-2019) to 5∼40 g cm −2 yr −1 below Table 1. Projected GPP (g cm −2 yr −1 ) for the entire study area (lower and upper 25% quantiles, mean and standard deviation) during contemporary (2001-19), mid-century (2041-59), and end of century (2081-99) periods under two different emission scenarios, as compared to the contemporary time period. Results summarized from the predicted GPP ensemble at 6 km by machine learning models driven by future climate projections from four GCMs. Areas are expressed as percentage of the entire study area. the LTM in subsections SNF-3, SNF-4, NCCR-3, and NCCR-4 (figure 3(d)). Maximum air temperature during the GS showed a more substantial negative impact than mean air temperature. Annual GPP decreased more than 100 g cm −2 yr −1 as GS maximum air temperature increased from 17 • C in a colder year to 21 • C in a warmer year in CCR-2, CCR-3, NCICR-1, and NCICR-3 ( figure 3(f)). Despite the large latitudinal difference, all subsections had higher annual GPP in years with a lower amount of solar radiation (figure 3(e)).

Changes in GPP under future climate projections
When driven by future climate projections from four GCMs, the GBRT model predictions showed that, cross the entire study area, mean annual GPP was projected to decrease from 816 ± 271 g cm −2 yr −1 to 809 ± 260 g cm −2 yr −1 by mid-century and to 812 ± 264 g cm −2 yr −1 by the end of the century (2081-2099) under RCP4.5 (table 1). Although mean reductions in GPP were projected to be small, these changes were not spatially uniform within the study area (figures 4 and 5). Both the lower and upper quantiles of the predicted mid-century GPP showed more significant reductions (e.g. by 20 and 26 g cm −2 yr −1 , respectively). Predicted GPP decreased in over 59% of the study area by 27 ± 21 g c m −2 yr −1 by the end of the century, which was partly counteracted by an increase over 41% of the study area (23 ± 18 g cm −2 yr −1 ) compared to the contemporary period (2001-2019) (table 1). Under the RCP8.5 emission scenario, the model predicted GPP changes similar to the RCP4.5 model over our study area; however, compared to the RCP4.5 model, a smaller extent of the study area would exhibit GPP reductions (i.e. only 56% and 47% by mid-and end of century, respectively) (table 1).
The magnitude of changes in GPP was larger though than under RCP4.5, especially by the end of the century, i.e. by −37 ± 32 g cm −2 yr −1 (table 1).

Spatial patterns of future GPP change and climate change contributions
Future GPP reductions were predicted to occur mostly in moderate to wetter areas, such as NCCR, central and southern Sierra Nevada, relatively wetter parts of the Central Coast, and to a lesser degree in NCICR (figures 5 and 6(a)). Less evenly distributed precipitation during the rangeland GS, as indicated by the projected increase in PCI (Figures. S4c and S5c), contributed the most to the GPP decrease by both mid-and end of the century, especially over regions with intermediate rainfall ( figure 6(b)). For example, GPP declined by as much as 15-30 g cm −2 yr −1 in subsections CCR-11, SNF-3, and SNF-4 (more than 29% of the entire study area) where the impact of PCI was dominant. In wetter regions, reduction in GPP was further caused by warmer T mean during GS such as in CCR and SNF (figures S4(e) and S5(e)). Warmer Feburary T min (figures S4(d) and S5(d)) was also predicted to induce a GPP reduction in costal subsections NCCR-1 and CCR-11, except in the wettest costal subsection NCCR-4.
In contrast, increases in predicted future GPP were found mostly in drier and less productive subsections, as shown in the 6 km map (figure 5) and summarized at the ecoregion subsection level (lefthand side of figures S4-S6), with a few exceptions. GPP in drier subsections that received less than 513∼522 mm yr −1 precipitation, mostly in the CCR, benefited from a predicted increase in precipitation ( figure S4(b)). Increasing minimum air temperature (figures S4(d) and S5(d)) further enhances rangeland productivity in these areas ( figure 6). The projected increase in maximum air temperature under RCP 4.5, however, caused a significant decline in GPP in the southern drier end of the central coast and northern interior coast, by 45 and 12 g cm −2 yr −1 in CCR-2 and CCR-3, and 13 g cm −2 yr −1 in NCICR-3, respectively.

Impacts of climate variables on rangeland GPP
This study leveraged 20 years of remote sensing records and machine learning methods to quantify and disentangle the interactive response of annual GPP to a suite of climate variables across California's annual rangelands. We found a high degree of spatial variability in responses across the 23 ecologically distinct subsections, as a result of varying sensitivity to climate change. While the relative importance of climatic variables varied in controlling interannual GPP variability across ecoregion subsections, variables related to precipitation, in general, were major drivers of rangeland productivity responses. Moving from drier to wetter areas, the driving climatic variable for GPP shifted from total GS precipitation to the relative distribution of GS rainfall (PCI). This is consistent with previous findings from field observations-that is, water is likely no longer the limiting factor to plant growth in wet regions in most years. Instead, timing and monthly distribution of precipitation become more important because, as intervals between rainfall events increase, the duration and severity of soil water stress is expected to increase (Knapp et al 2008). Furthermore, we found greater total precipitation in wetter areas reduced GPP, likely in part because oversaturated soils can create anaerobic conditions that impair growth and development of upland (i.e. non-wetland species) plants (Wilson and Livingston 1932). Wetter years may also have colder temperatures and lower solar radiation resulting from greater cloud cover.
Air temperatures (minimum, mean, and maximum) showed different controls on GPP. We found a positive relationship between GPP and winter minimum air temperature in all but four subsections (SNF-1, SNF-2, SNF-3, and NCCR-3), likely due to warmer temperatures in the early season promoting germination and early growth . In contrast, Mean and maximum air temperatures had a negative control on GPP in several subsections across the study area, likely by increasing water demand and thus causing water stress (Devine et al 2019). Higher solar radiation during GS was also shown to cause GPP reduction in three subsections, probably via enhanced drying processes. Further work is needed to better understand coupling effects of climate variables, such as the combined impacts of cold temperatures and high precipitation on GPP, as well as reasons for spatial differences in climate controls on GPP, such as why GPP is not as responsive to air temperature in some areas as it is in other areas.

Implications of future changes in rangeland GPP
While our models predicted small changes in future rangeland productivity for the study area as a whole (table 1), impacts at the ecoregion subsection scale were more marked and suggest rangeland productivity responses to climate change will be highly variable at the local level (figure 4). We found greater magnitude of GPP changes (mostly increases) in many drier subsections (figure 6(a)) in response to slightly higher predicted precipitation and higher GS minimum temperatures over time. Considering the entire study area, predicted increases in GPP from favorable future precipitation and temperature conditions in mostly drier subsections (41% of the study area) were counterbalanced by predicted decreases in GPP from unfavorable changes in temperature, rainfall distribution, and solar radiation in wetter subsections (59% of the study area)-resulting in the overall small decrease in GPP.
However, it is worth noting that several uncertainties in GPP projections exist due to inconsistencies among GCMs, especially in precipitation simulations. For example, although four GCM projections agreed on larger variation in monthly precipitation (less evenly distributed) in the future, the amount of precipitation change, however, varied among GCMs. A drier climate was projected by MIROC5 versus a wetter climate by CNRM-CM5 (figures 4(c) and (d)), resulting in large discrepancies between GPP projections: by mid-century under RCP4.5, GPP projections with MIRCO5 (dry model) versus CNRM-CM5 (wet model) differed by as much as 36 g cm −2 yr −1 (figures 4(c) and (d); see detailed discussion in supplemental materials section 2). The sign of projected precipitation changes also varied by regions, in contrast to the warming trend across the state by all model simulations.
Climate change threatens California's annual rangelands, which are largely rainfall-dependent and therefore highly vulnerable to weather extremes (Macon et al 2016, Roche 2016). More extreme weather events and more frequent year-to-year weather swings in the future will lead to greater interannual rangeland productivity, potentially threatening the sustainability of forage supply and, ultimately, the resilience of ranches and rangelands to climate change. Other modeling approaches have similarly projected declines in productivity in various rangeland systems around the globe, with associated decreases in ecosystem goods and services, including livestock production (Boone et al 2018). The high degree of spatial variability in productivity impacts projected by our models means local, climate-informed decision-making will be increasingly critical to identifying adaptation strategies for mitigating climate risks.
This work can be leveraged to create nearreal time rangeland productivity tracking and forecasting tools to support adaptive decision-making across California's extensive and highly variable rangelands. Future research and development will need to combine new data and emerging technologies with land management expertise to create locally-calibrated decision support systems. Working with land managers to co-develop decision-support tools will ensure on-the-ground relevance to sustainability and resilience at local scales and, therefore, greatly increase manager adoption. Broader benefits potentially include producing locally-trusted tools and data to help inform policy and management guidance (e.g. drought and disaster early warning systems) at regional and national scales.

Conclusion
Climate change is expected to impact forage production across California's annual rangelands, threatening the economic viability and environmental sustainability of these highly valued systems. This study projected the degree to which GPP may respond to changes in different climatic variables by mid-century and end of century under RCP4.5 (intermediate emissions) and RCP8.5 (worst-case emissions) scenarios using remote sensing and machine learning techniques. Our machine learning-based analysis highlighted regional differences in GPP vulnerability to climate and provided insights on the intertwining and potentially counteracting effects of seasonal temperature and precipitation regimes. While total annual precipitation dominates the variations in annual GPP in drier ecoregions, changes in seasonal precipitation distribution affects the interannual GPP variability, especially in the wetter areas. We found that warmer winter season minimum air temperature promotes GPP for the majority of ecoregions; however, the maximum temperature during the GS reduced production in the CCR and NCICR. These results can inform rangeland management and conservation in the context of climate adaptation during an era of climate change.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.