Estimating climate change effects on net primary production of rangelands in the United States

The potential effects of climate change on net primary productivity (NPP) of U.S. rangelands were evaluated using estimated climate regimes from the A1B, A2 and B2 global change scenarios imposed on the biogeochemical cycling model, Biome-BGC from 2001 to 2100. Temperature, precipitation, vapor pressure deficit, day length, solar radiation, CO2 enrichment and nitrogen deposition were evaluated as drivers of NPP. Across all three scenarios, rangeland NPP increased by 0.26 % year−1 (7 kg C ha−1 year−1) but increases were not apparent until after 2030 and significant regional variation in NPP was revealed. The Desert Southwest and Southwest assessment regions exhibited declines in NPP of about 7 % by 2100, while the Northern and Southern Great Plains, Interior West and Eastern Prairies all experienced increases over 25 %. Grasslands dominated by warm season (C4 photosynthetic pathway) species showed the greatest response to temperature while cool season (C3 photosynthetic pathway) dominated regions responded most strongly to CO2 enrichment. Modeled NPP responses compared favorably with experimental results from CO2 manipulation experiments and to NPP estimates from the Moderate Resolution Imaging Spectroradiometer (MODIS). Collectively, these results indicate significant and asymmetric changes in NPP for U.S. rangelands may be expected.


Introduction
Net primary productivity (NPP), or the rate of assimilation of CO 2 through photosynthesis, is the fundamental link between the atmosphere and the biosphere. Human activities emit large quantities of CO 2 into the atmosphere and alter climate patterns, which directly impacts NPP (Greer et al. 1995;Nakićenović et al. 2000). Furthermore, humans already consume almost a quarter of potential NPP (Haberl et al. 2007).
Rangelands provide a multitude of biological and social benefits that are a key component of global sustainability (Reynolds et al. 2007;Munasinghe 2009). The future of goods and services derived from rangelands, such as fuelwood and protein, however, is uncertain as are impacts expected from climate change. In general, across the extent of U.S. rangelands, climate change is being expressed as warmer temperatures and variable precipitation which are expected to continue into the foreseeable future (IPCC 2013). Low and variable precipitation, extreme temperatures and high evaporative demands render most rangelands inherently vulnerable to degradation under climate change (Reynolds et al. 2007;Korner 2000;Le Houérou 1996;Sala et al. 2000). Thus, projections of NPP in rangelands are fundamental to understanding climate change impacts on global ecosystems and sustainability of goods and services.
In rangelands, productivity is mainly determined by the distribution of precipitation and resultant effects on soil water availability (Campbell et al. 1997;Knapp et al. 2001;Izaurralde et al. 2011). Soil water availability is subsequently dependent on other environmental factors such as temperature, vapor pressure deficit, soil properties and CO 2 concentration via effects on stomatal conductance (Izaurralde et al. 2011;Morgan et al. 2004bMorgan et al. , 2011. Elevated CO 2 may increase NPP in semi-arid regions by increasing plant water use efficiency (Fay et al. 2003;Izaurralde et al. 2011). Furthermore, plant species vary in their response to these factors and alteration of NPP can be expected in the future as species respond to climate change through range shifts or local population dynamics (Morgan et al. 2007;Polley et al. 2012).
Much of our present knowledge, empirical evidence, and expectations of NPP response to climate change are deduced from localized experiments, which have primarily focused on forests or crops. The limited applicability of most experimental studies, however, prohibits confident extrapolation to a comprehensive regional picture of future trends of NPP (Norby and Luo 2004). Therefore, regional analyses should consider interactions of biophysical factors, such as nitrogen deposition and elevated CO 2 , (Hungate et al. 2003;Milchunas et al. 2005;Maestre et al. 2012;Wouter et al. 2012) to improve understanding of spatial variability in response of NPP to climate change and to identify areas where adaptation measures may be required (Fussel and Klein 2006;Cleland et al. 2007;Joyce et al. 2013). In response to this need, a regional analysis of NPP, projected to 2100, for rangelands of the coterminous U.S. was conducted using the ecosystem process model, Biome-BGC (Running and Hunt 1993). Dominant metabolic pathway was used to broadly indicate influence of current species composition (Epstein et al. 2002;Zavaleta et al. 2003;Polley et al. 2012). Current species composition, coupled with Biome-BGC, and climate projections from four general circulation models (GCM) under three emissions scenarios were used to model NPP at a spatial scale of 5 arc minutes (about 9.3 by 7.1 km at 40°north latitude). In addition, climate and dominant drivers of NPP were characterized, temporally and spatially, across U.S. rangelands for each GCM and scenario. Finally, spatial and temporal patterns of changes in NPP, from a contemporary baseline without regard for possible future land cover changes or potential disturbance regimes, were quantified.

Climate data
A spinup simulation was performed to create steady state conditions for the ecosystems from which future projections began. Biome-BGC requires daily estimates of maximum and minimum temperature, precipitation, daylength, solar radiation and vapor pressure deficit. Temperature and precipitation data for the spinup period were obtained from the PRISM (Parameter-elevation Regressions on Independent Slopes Model) project (Daly et al. 2000). Vapor pressure deficit, solar radiation and day length are not produced by PRISM. These variables were created from PRISM data using the MT-CLIM algorithms (Kimball et al. 1997). The resulting climatological data for the spinup were developed at a daily temporal resolution and matched the spatial resolution of the GCM data. The delta method (Mote and Salathé, 2009) was used to temporally downscale the monthly PRISM Data, from 1940-2000, to a daily temporal resolution which were used for the spinup period (Supplementary Material (SM) 1).
The projection covered three emissions scenarios from 2001-2100 derived from the Third Intergovernmental Panel on Climate Change (IPCC) Assessment Special Report on Emission Scenarios (Nakićenović et al. 2000). These emissions scenarios encompass different assumptions regarding socioeconomic drivers such as gross domestic product, population growth, technological innovation and greenhouse gas emissions (SM 2). For the present study, the IPCC scenarios corresponded to the A1B, A2, and B2 scenarios. Climate data processed by Coulson et al. (2010a, b) Coulson et al. (2010a, b) spatially downscaled data from a selected suite of GCM's used in the Fourth IPCC Assessment resulting in estimates of monthly averaged maximum and minimum temperature and precipitation at the 5 arc minute spatial resolution. The GCM data were spatially downscaled by Coulson et al. (2010a, b) using the ANUSPLIN software to fit a twodimensional spline function at each month's for each climate variable. These procedures are described at length in Joyce et al. (2011Joyce et al. ( , 2014. As with the PRISM data monthly estimates of maximum and minimum temperature and precipitation from the GCM's were temporally downscaled to a daily time-step to accommodate the needs of Biome-BGC using the delta method (SM 1). From these daily data, estimates of daylength, solar radiation and vapor pressure deficit were derived using the MT-CLIM algorithms.

Spatial domain of analysis
Defining the extent and composition of the area for estimating potential future rangeland NPP required two steps (SM 1). For the first step, the spatial extent of rangelands in the coterminous U.S. was identified from Reeves and Mitchell (2011) (SM 3). In the second step, the relative abundance of C3, C4, and shrub plant functional types were quantified at each pixel in the rangeland domain from Step 1. Cross-hatched regions depicted in SM 3 were used for validation and comparison purposes but not for describing NPP response to climate change, because, although shrub or herb dominated, these regions are outside the rangelands study area as defined from Reeves and Mitchell (2011). All modeling of NPP projections was conducted at the 5 arc minute pixel level. Results, however, were often aggregated to larger spatial extents including the six assessment regions depicted in SM 3.

Nitrogen deposition, carbon dioxide concentrations, and soils
The emission scenarios included decadal, aspatial estimates of NO x emissions (SM4). Decadal values were linearly interpolated to estimate annual values. Further, nitrogen deposition was assumed to increase at the same rate as nitrogen emissions, therefore, permitting the rate of increase of nitrogen emissions as a proxy for the rate of increase of nitrogen deposition. Since these data were tabular, data from Holland et al. (2005), describing current patterns of nitrogen deposition, were used to apportion projections of deposition across the landscape for each scenario. This ensured that the spatial pattern and relative magnitude of nitrogen deposition was held constant throughout the projection, although total deposition was based on the temporal trends found in estimated nitrogen emissions for each scenario (SM 4).
As with nitrogen emissions, decadal estimates of CO 2 concentrations accompanied each IPCC emissions scenario and are represented in SM 4, which shows the trajectory of CO 2 concentrations across the projection. Estimates of decadal CO 2 concentration were aspatial and applied symmetrically to every pixel as CO 2 is assumed to be distributed homogenously across the study area.
To provide linkages with the forthcoming Fifth IPCC Assessment Report (AR5) (Moss et al. 2010;Rogelj et al. 2012) (http://www.ipcc.ch/activities/activities.shtml#.Ug6Z0UG1Hh4), suggested linkages between the CO 2 concentrations found in the Representative Concentration Pathways and the socioeconomic scenarios used in the present study is provided in SM 4. Our linkages are based on similarities of CO 2 concentration by 2100, and differ, for the B2 scenario, from the linkages conducted by Rogelj et al. (2012).
The final parameters needed to simulate future NPP with Biome-BGC included soil data, which were derived from the State Soil Geographic database (STATSGO) (U.S. Department of Agriculture and Natural Resources Conservation Service NRCS 1994). The percentage of sand, silt and clay and depth to bedrock were quantified for the study area. These data were gridded to the 5 arc minute resolution to match all other spatially explicit simulation parameters. These soils data are used by Biome-BGC to control belowground processes such as decomposition and water balance.

Simulating vegetation productivity
Using the present vegetation distribution, climate, CO 2 , NO x , and soils data as described above, estimates of annual NPP on U.S. rangelands were produced by simulating ecosystem dynamics through application of Biome-BGC across the entire projection for every GCM and scenario combination (SM 5). The contemporary distribution of plant functional types and photosynthetic pathways were held constant throughout the projection. This was done because it isn't feasible to make spatially-explicit, realistic assumptions about vegetation change over a century time span over the study area. Therefore, results only reflect changes in CO 2 concentration, nitrogen deposition, and climate variables, which was the focus of the study. Change in NPP at each pixel was characterized by the slope of the linear trend with respect to time (2001-2100).

Analysis of NPP and drivers of change
Though the Biome-BGC simulations were conducted across the entire projection for every GCM and scenario combination, NPP data were averaged across each GCM to produce scenario-level results. Biome-BGC was run for each emission scenario using climate projections from each of several GCMs. NPP derived from different GCMs was then averaged for each emission scenario. The spatial differences between the various GCM's are intensively explored in Joyce et al. (2011) and are therefore omitted here. Spearman correlation coefficients were computed at each pixel to identify which of the aforementioned drivers was most closely related to patterns of NPP through time for each scenario. The bioclimatic factor with the highest correlation to estimated NPP was retained and considered the dominant driver of change in productivity. In a similar fashion, temporal trends in NPP were described using a linear trend to determine the direction and magnitude of change. Precipitation and temperature were also evaluated using a linear trend enabling the direction and magnitude of change relative to present day conditions to be quantified. For NPP and climate data, the response of the three emissions scenarios were averaged together and differences among them were represented as standard deviations about the mean.
2.6 Quantifying model agreement Direct validation of projected estimates of NPP was not possible because no spatially extensive field-referenced estimates of NPP currently exist for coterminous U.S. rangelands. There are data, however, from remote sensing instruments that can be used to characterize agreement among NPP estimates. Net primary productivity derived from the Moderate Resolution Imaging Spectroradiometer (MODIS) satellite has been used for evaluating NPP and biomass dynamics of rangeland biomes (Heinsch et al. 2006;Reeves et al. 2006;Zhang et al. 2010). The MODIS-derived NPP estimates are available globally at 1 km 2 spatial resolution and were compared with NPP derived using Biome-BGC within the regions identified in SM 3. This analysis does not constitute an accuracy assessment but provides insight into locations where model results may be questionable or, conversely, where they are more reliable.

Projected climatic trends
Average annual temperature increased (Fig. 1a) in each of the three scenarios but varied considerably by region (Fig. 1b). Near 2018, the average annual temperature was projected to surpass the upper limit of historic variation measured since 1940 and continued to rise steadily throughout the projection. The Interior West and Eastern Prairie assessment regions exhibited the greatest amount of warming (0.032°C±0.002 and 0.003 years −1 , respectively). Some ecological subsections (Bailey and Hogg 1986) in the Interior West, at the upper altitude limits of the coterminous U.S., experienced significant increases in temperature, often exceeding 4°C by the end of the projection. The Northern Great Plains, Southern Great Plains, and Southwest all warmed at about the same rate (~0.03°C±0.002, 0.0001, 0.001 year −1 , respectively), while the Desert Southwest warmed at the slowest rate (0.02°C±0.002 years −1 ). Despite these generalities, significant sub-regional variation in temperature trends was exhibited (Fig. 1b).
In contrast to temperature, projected changes in precipitation were not as dramatic, but more variable, in terms of scenario differences. From a national perspective, precipitation remained close to or within the historical range of variability of precipitation for all three scenarios across the entire projection (Fig. 1c). Overall, the mean precipitation difference from 2001 to 2100 across all three scenarios was −1.5 % with a standard deviation of 5.2 %. This standard deviation represents the variability inherent among precipitation projections associated with the three scenarios. Among the assessment regions, the Northern Great Plains had the greatest rate of precipitation increase (Fig. 1d) but was also highly variable (0.06 %±0.06 years −1 ) between scenarios. Annual precipitation in the Eastern Prairies and Southern Great Plains declined at an average rate of −0.07 and −0.06 % year −1 and exhibited less variability than in the northern Great Plains. The Southwest assessment region exhibited the greatest decrease in annual precipitation, albeit with high variability between scenarios. The mean precipitation response was positive in the Desert Southwest, but the region exhibited extreme variability (standard deviation of 4 times the mean).

Vegetation productivity
The overall NPP trends suggest increases of 16, 19 and 15 %, for the A1B, A2, and B2 scenarios, respectively, for all coterminous U.S. rangelands by 2100 (Fig. 2). Although the trend for the entire projection was strongly positive with an average annual increase of 0.26 %, NPP did not begin moving in a positive direction until approximately 2030 (Fig. 2). The earlier stages of the projection suggest decreasing NPP from present day to the 2030s while tracking the trend in precipitation (Fig. 3). After this initial lag, NPP begins to increase steadily and diverge from the trend in precipitation (Fig. 3). The divergence of the NPP response from precipitation is related to CO 2 fertilization and ever increasing temperatures (Fig. 3). Lengthened growing season increases opportunity for CO 2 fixation if moisture is sufficiently abundant (Christensen et al. 2004). Thus, losses due to moisture limitations in the southern extent of the study area begin to be offset and surpassed by increases in NPP due to increased growing season length in northern latitudes beginning in the 2030s. The amount of disagreement among scenarios reached its apex in the last decade of the simulation period (Fig. 2). These general trends, however, bely patterns across rangelands in the western United States (Fig. 4, SM 6).
The NPP of the Southern Great Plains, Northern Great Plains, and Eastern Prairies increased at roughly the same rate of 3.6 % increase per decade while the Interior West increased the most at 4.5 % per year. NPP of the Desert Southwest and Southwest regions decreased slightly during the projection at −0.67 and −1.05 % per decade, respectively. While average regional responses act as a general guide, within each assessment region there was significant variation in terms of the NPP response. For example, declines in productivity for southern and southwestern Texas and central Arizona often exceeded 11 % during the projection (Fig. 4).
Interpretation of information portrayed in Fig. 4 requires evaluation of the deviation or disagreement among the scenarios in terms of NPP response. In many instances, parts of the Desert Southwest and Southwest experienced significant declines in NPP, yet exhibited high disagreement among scenarios suggesting that the ultimate fate of NPP of these regions is more uncertain. In contrast, the variation in NPP response across the northern and southern Great Plains among the three scenarios was quite low.  positively biased), the Biome-BGC-derived NPP was considerably higher. Though Biome-BGC-derived NPP compared reasonably well with MODIS and experimental results, multiple sources of uncertainty and several caveats should be considered.

Drivers of change
Relations between bioclimatic drivers and projected NPP trends are shown in Fig. 5. The Northern and Southern Great Plains and Eastern Prairies were characterized as grassdominated expanses comprised of mixtures of C4 (warm season) species (SM 3) and C3 (cool season) species in greatly varying proportions. Data portrayed in SM 3 indicate the proportional abundance of differing photosynthetic pathways as many areas are mixtures of species. Net primary productivity of these regions tended to be most influenced by temperature (Fig. 5), but the southern end of the Southern Great Plains and eastern edge of Northern Great Plains experienced change most related to precipitation. In contrast, the Interior West, dominated by C3 grasses or shrubs, tended to be most highly correlated with CO 2 in terms of NPP through time. As expected, NPP trends of the southerly rangeland extent were most highly correlated with precipitation and water relations. Nitrogen deposition (not shown; SM8) was rarely found to be a greater determinant of productivity than climatic or CO 2 concentration factors, which may reflect relatively low levels of projected change compared with other factors.

Discussion
Overall, the balance of changing NPP across U.S. rangelands was positive with gains outpacing losses by the 2030s. Trends of increasing NPP were observed for northern and central regions while decreasing NPP was apparent across the southwestern U.S. (Fig. 4). Additionally, NPP was projected to increase proportionally more at high elevations, Fig. 5 Bioclimatic driver with the highest correlation to estimated NPP trends particularly in the Interior West (Fig. 4). The proportionally larger increase in NPP at higher elevations was similarly noted by Riedo et al. (2001) in alpine rangeland landscapes at double ambient CO 2 conditions. Drivers of changing productivity followed a geographic pattern resembling that of NPP with the primary driver being precipitation in southern regions, temperature in northern regions, and CO 2 in the Interior West. Nitrogen deposition had relatively little influence at these scales.
Similar to our results, Polley et al. (2013) suggested that NPP will decrease in southern rangelands in response to warmer temperatures and declining precipitation, and will increase on northern rangelands as a result of warmer temperatures and greater precipitation. Model results could not be directly validated with NPP measures across U.S. rangelands, but several studies provide comparable localized results. Results from the Prairie Heating and CO 2 Enrichment (PHACE) study ) for a mixed-grass prairie in Wyoming compare favorably with our results. Morgan et al. (2011) reported an increase of aboveground productivity by an average of 33 % over 3 years, which substantiates our estimate of a 28 % increase in productivity for the Northern Great Plains by 2100 in response to a similar total increase in CO 2 . Similarly in a Colorado shortgrass steppe, enrichment of CO 2 to double ambient levels positively affected the growth of C3 grasses and significantly increased aboveground biomass (Morgan et al. 2007). For tallgrass prairie, Owensby et al. (1999) found that elevated CO 2 could increase productivity of above-and below-ground biomass, but response was dependent on water stress. Aboveground NPP increased with CO 2 enrichment only during very wet years in the Mojave desert, but was unchanged over a 10 year period suggesting CO 2 has little influence in a perennial plant community characterized by low plant cover and rare recruitment events (Newingham et al. 2013;Newingham et al. 2014). These findings are consistent with results from our model and suggest that desiccation effects of increased temperature can be offset to some extent by CO 2 enrichment via reduced transpirational demand (Leakey 2009;Morgan et al. 2004bMorgan et al. , 2011Woodward and Kelly 2008), but responses in arid regions are less certain. Within the range evaluated here, CO 2 concentration ultimately reduces moisture limitations via greater water use efficiency (Bachelet et al. 2001;Christensen et al. 2004;Morgan et al. 2008Morgan et al. , 2011Polley et al. 2003). NPP may decrease, however, if CO 2 becomes sufficiently abundant to reduce N mineralization rates, through greater production of structural components, and N availability to plants (King et al. 2004;Morgan et al. 2001Morgan et al. , 2004aMorgan et al. , 2007Milchunas et al. 2005).
It is worth noting that all models are a simplification of reality and interpretation of model results needs to consider uncertainty, inputs, and model assumptions. Differences among the GCM's used here have been extensively examined (Randall et al. 2007, Tabor andWilliams 2010). Of importance for interpreting model output for the present study is the increasing disparity among models as time progresses and the more pronounced variation in arid regions where precipitation is the primary driver of NPP change. Projections of productivity in desert areas, omitted here, are additionally problematic, because Biome-BGC does not evaluate response of Crassulacean Acid Metabolism (CAM) plants (Running and Hunt 1993). Although more uncertain than projections for northern regions, we note that any decrease in productivity in arid regions is likely to have large ecological impacts because NPP is already low (Tietjen et al. 2010).
Although Biome-BGC models the complex interaction of water, carbon, and nitrogen between the soil and atmosphere, it does not model all factors determining rangeland productivity and vegetation type is held constant. Changes in species assemblages, disturbance regimes, and land use will undoubtedly affect future NPP in complex and unexpected ways (Morgan et al. 2007;Suttle et al. 2007;Walther 2010). Thus, projections presented here represent a potential for NPP based on climate and soils in the absence of other change. NPP derived by Biome-BGC from 2000-2012 compared reasonably well with MODIS-derived NPP from the same time period, although the Biome-BGC-derived NPP was considerably higher (SM 7). This is not necessarily a reflection of poor model performance because Biome-BGCderived NPP represents a kind of potential rather than actual NPP. Satellite data, such as MODIS, reflects land management effects on vegetation productivity and also relies on models to estimate NPP. In contrast, Biome-BGC models NPP based on climate and other biophysical drivers without an implicit link to land management practices such as livestock grazing.
In this vein, the present work is novel given the spatially explicit application of Biome-BGC accounting for multiple and interacting factors simultaneously applied across a large region representing diverse climate futures. Results demonstrate the importance of regional or national analyses aimed at understanding the response of rangeland NPP to climate change, which has important implications for sustaining livelihoods and food security (Izaurralde et al. 2011;Polley et al. 2013). For example, on native pastures in Queensland, Australia, Hall et al. (1998) estimated that the combined effects of warmer temperature and CO 2 enrichment could raise carrying capacity of livestock by 3-45 % depending on the rainfall scenario. Others have predicted a reduction in livestock production even with increasing NPP, due to concurrent reductions in plant nitrogen content, higher production variability, and lower animal performance (Hanson et al. 1993). Although the overall increase in potential NPP for rangelands is potentially important, localized benefits should not be used as evidence to suggest that mitigation aimed at reducing greenhouse gas concentrations does not need to be pursued. On the contrary, the complexity of estimating future rangeland NPP indicates a strong need for improved forecasting tools and coordinated field monitoring (Luo et al. 2011).