Impacts of Climate Change on Hydroclimatic Conditions of U.S. National Forests and Grasslands

: The conterminous United States includes national forests and grasslands that provide ecological, social, economic, recreational, and aesthetic services. Future climate change can alter long-term hydroclimatic conditions of national forests and grasslands and lead to negative consequences. This study characterizes shifts in hydroclimatology and basin characteristics of US National Forests (NFs) and National Grasslands (NGs) in response to climate change over the 21st century under the DRY, MIDDLE, and WET climate models with the representative concentration pathway (RCP) 8.5 emission scenario. Climatic projections for three climate models ranging from the driest to wettest conditions were obtained from the Multivariate Adaptive Constructed Analogs (MACA) dataset. Then, the variable inﬁltration capacity (VIC) model was used to model hydrological responses of the selected future climates. Changes in regional hydroclimatic conditions of NFs and NGs were assessed by the magnitude and direction of movements in the Budyko space. The Fu’s equation was applied to estimate changes in basin characteristics. The results indicate that NFs and NGs are likely to experience larger changes in basin characteristics compared to the average of the United States. In general, across the conterminous US, the NFs in mountainous regions are likely to have larger changes in hydroclimatic variables than NFs with lower elevation and NGs. Comparing Forest Service regions, Paciﬁc Northwest, Intermountain, and Northern regions may have a less arid climate with lower freshwater availability. The Southwestern, Northern, Intermountain, and Rocky Mountain regions are likely to experience higher shifts in their basin characteristics. This study can help environmental scientists, and land and water managers improve future land management plans. Author Contributions: Conceptualization, H.H., M.A.; methodology, M.A. T.C.B.; software, validation, analysis,


Introduction
The United States National Forest System includes National Forests (NFs) and National Grasslands (NGs), which are divided into eight regions in the conterminous United States (CONUS) managed by the U.S. Forest Service (USFS) [1]. The USFS provides a wide range of services for present and future generations that have broadened to include hydrological, ecological, social, economic, recreational, and aesthetic services [2][3][4]. A challenge for the agency is that climate change may lead to shifts in hydroclimatic conditions of NFs and NGs, such as a decrease in freshwater availability and may cause changes in the structure and composition of forests and grasslands at various spatial and temporal scales [2,[5][6][7][8][9][10]. Assessments of potential shifts in hydroclimatic conditions of NFs and NGs can help decision-makers to mitigate the negative consequences of deforestation [11].
This study characterizes changes in regional hydroclimatology and basin characteristics of NFs and NGs across the CONUS. Climate change has already affected forests and grasslands throughout the CONUS via changes in wildfire frequency, forest pathogens, insect populations, and timber resources [12][13][14][15]. Although significant progress has been

Materials and Methods
The downscaled Multivariate Adaptive Constructed Analogs (MACA) datasets [30] were used to characterize changes in future climate variables under three climate models representing future wet, dry, and intermediate conditions. Then, the variable infiltration capacity (VIC) model [31] was applied to simulate future hydrological variables driven with the future climate scenarios. Finally, the Budyko framework was used to project shifts in hydroclimatic conditions of NFs and NGs, and Fu's equation was used to help assess changes in integrative basin characteristics. In this study, the 30-year average of hydroclimatic variables from 1986 to 2015 represents current conditions and the 30-year average of hydroclimatic variables from 2070 to 2099 describes future conditions.

U.S. National Forests and Grasslands
The CONUS is divided into eight Forest Service Regions (Figure 1). Each region encompasses NFs and NGs with diverse landscapes, ecosystems, fauna, and flora [2,3]. Although the Southern region is the largest region, spanning from Texas to Virginia, most NFs and NGs are located in the western United States. The total NF + NG area is about 114 million ha that NFs account for nearly 98% of this total NF + NG area in the CONUS.

Hydroclimatic Projections
Current climate variables  including minimum and maximum temperature and precipitation at the grid size of~4 km (1/24 degree) were obtained from a combination of the daily surface weather and climatological summaries (Daymet) [32] and the Parameterelevation Regressions on Independent Slopes Model (PRISM) datasets [33]. The daily precipitation and temperature were derived from Daymet and biased corrected with PRISM at the monthly scale. The North American Regional Reanalysis (NARR) dataset [34] was used to obtain the historical wind speed data. Readers are referred to Oubeidillah et al. (2014) [35], Naz et al. (2016) [36], and Heidari et al. (2020) [20] for the detailed description of the historical climate dataset.
Future climate variables (2070-2099) were obtained from the downscaled MACA datasets [30]. The MACA climate dataset provides CONUS-wide downscaled (to the grid size of~4 km, or 1/24 degree) climate projections for the RCP 4.5 and RCP 8.5 emission scenarios for twenty climate models. Three downscaled climate models were chosen to represent a range of possible climate scenarios, on average, ranging from wet to dry, plus one model that represents the middle of the range (Table 1) [20,37]. The WET, DRY, and MIDDLE climate models under RCP 8.5 emission scenario were selected based on a range of changes in precipitation from current to future conditions. The use of these climate models allows to characterize the estimated range in projected annual precipitation across all projections for entire the CONUS. Note, opposite patterns in individual ensemble members can originate from internal variability in the large-scale circulation [38][39][40]. A good representation of current climate is indeed a necessary condition required to realistically simulate future climate changes. Figure S1 compares the 30-year average of annual precipitation for each HUC08 watershed over the historical period . The 30-year average values for each HUC08 watershed are approximately the same meaning that the variation between climate models is not significant at the 30-year average scale. Additionally, 30-year average of climate models over the historical period is close to the selected baseline historical model that is the combination of Daymet and PRISM.
Therefore, for the consistency over the historical period, we decided to use the combination of Daymet and PRISM models. The three climate models were able to reproduce current climate conditions with the same statistics. Although the 30-year average values for the selected baseline is close to the historical MACA climate models, the performance and skills of each selected global climate model (GCM)are not specifically evaluated in this study. Readers are referred to Naz et al. (2016) [36], and Joyce and Coulson, (2020) [37] for the detailed information about the climate models and historical projections.
The projected climatic variables were then used as inputs to the VIC (version 4.1) model [31] at the grid size of~4 km (1/24 degree) to project precipitation, potential evapotranspiration, and water yield across the CONUS over the 21st century at a HUC8 level scale. The VIC model is a macroscale semi-distributed hydrological model to simulate land-atmosphere fluxes and the water and energy balances at the land surface [41]. The VIC model has been widely used to simulate streamflow over a number of large river basins in North America [42].
The VIC model uses the variable infiltration capacity curve to solve full water and energy balances and estimate infiltration and surface runoff. The model takes into account snow processes. The land surface can be modeled as a grid of uniform cells at a daily or sub-daily time step. Each grid cell can represent the spatial variability of precipitation, topography, and vegetation [41]. Several key assumptions have been made in the VIC model, including that the atmosphere is the only source of incoming water for a grid cell, and that grid cells are independent of each other meaning that there are no horizontal water and energy exchanges between grid cells [43].
Topography, soil characteristics, vegetation, land surface classification, and meteorological forcing are key hydrological inputs to the VIC model. Organized and calibrated VIC input data for the US Geological Survey (USGS) eight-digit hydrologic subbasins across the entire CONUS were obtained from Oubeidillah et al. (2014) [35]. The VIC model uses the Penman-Monteith equation to estimate potential evapotranspiration.
Aggregated monthly runoff data from USGS National Water Information System gauge observations (WaterWatch dataset) [44] was used to calibrate the VIC model for each HUC8 basin. Calibrating the VIC model using the WaterWatch monthly runoff data allows the homogenous application of the VIC model to all relevant grid cells. To calibrate the VIC model, the simulated monthly total streamflow (surface runoff plus baseflow) of each HUC8 basin was matched with the monthly runoff from the USGS WaterWatch runoff dataset. Readers are referred to Oubeidillah et al. (2014) [35] and Naz et al. (2016) [36] for the detailed description of the VIC model set up, calibration, and evaluation. Figure S2 compares the observed versus simulated annual water yield for each NFs and NGs from 1986 to 2015 period. The VIC model shows a strong linear correlation (0.9799) between observed and simulated mean annual water yield.

Movements in the Budyko Space
The combined changes in hydroclimatic variables can be characterized by movements in the Budyko space [45,46]. The x-axis in the Budyko space shows the aridity index defined as: Aridity Index = PET P (1) and the y-axis in the Budyko space characterizes the evaporative index defined as: where PET, P, and Q are respectively potential evapotranspiration, precipitation, and water yield ( Figure 2). P-Q in Equation (2) can be simplified to actual evapotranspiration [20].

Movements in the Budyko Space
The combined changes in hydroclimatic variables can be characterized by movements in the Budyko space [45,46]. The x-axis in the Budyko space shows the aridity index defined as:

=
(1) and the y-axis in the Budyko space characterizes the evaporative index defined as: where , , and are respectively potential evapotranspiration, precipitation, and water yield ( Figure 2). P-Q in Equation (2) can be simplified to actual evapotranspiration [20].  Movements in the Budyko space over time from current conditions to future conditions can be a combination of shifts in aridity index (∆x) and evaporative index (∆y). Improved understanding of the relative contribution of such drivers is a challenge to comprehensively characterize changes in future hydroclimatic conditions and basin characteristics in response to climate change. The movement can be described by a direction (D) and magnitude (M) defined as [20]: where the direction of movements identifies regional differentiation and magnitude of movements describes sensitivity in response to climate change [17]. A river basin may move toward the right (D = 0) in the Budyko space over time by increase in the aridity index, meaning that the climatic condition gets more arid. Conversely, a river basin may move toward the left (D = 180) by decrease in aridity index, indicating that the climatic regime becomes less arid. Additionally, a river basin may move upward (D = 90) over time by increasing evaporative index indicating lower river discharge, or downward (D = 270) indicating higher freshwater availability and lower evaporation. Direction of movement in the Budyko space can be representative of regional differentiation in response to climate change.
Magnitude of change in the Budyko space can be representative of sensitivity of a river basin to climate change. Basins with high magnitude of change are more prone to experience prolonged drought or long-term wetting period which can considerably affect their water and natural resources. Readers are also referred to Heidari et al. (2020) [20] for further explanation of movements in the Budyko space. Changes in long-term hydroclimatic conditions can alter the frequency of extreme events such as drought [47], and flood [48][49][50]. In this study, we obtained the current and future potential evapotranspiration (PET), P, and Q from the VIC model and applied Equations (1)-(4) to estimate shifts in the Budyko space of NFs and NGs across the CONUS from current to future periods. Note, the PET obtained from the VIC model is potential evapotranspiration from open water. It was assumed that there is sufficient open water supply so that the PET values are the maximum (potential) evapotranspiration capacity.

Changes in Basin Characteristics
Previous studies have proposed several analytical equations describing relationships between the aridity and evaporative indices, including the Fu's one-parameter equation [29]. Fu's equation accounts for influencing factors such as basin size, seasonal variability, and soil and vegetation characteristics that can affect the relationship. Fu's equation is defined as: where ω is a free parameter that has no physical meaning ( Figure 2). The ω free parameter can represent an integrative property of the catchment that globally ranges between 1.3 to 4.6 with a median value at 1.8 [51]. Previous studies reported that differences in ω are highly correlated with differences in land cover, vegetation cover, basin slope, and area [52][53][54][55].
A region with a sufficiently high ω is energy-limited, while a region with a sufficiently small ω is water-limited. In a water-limited region, forests have deeper and larger root systems that help them access more soil water [56]. The basin and vegetation characteristics of a region with small ω (i.e., a water-limited region) can be more affected by hydroclimatic changes [57]. Characterizing change in ω provides an understanding of how long-term climate and hydrological changes interactively affect forests and grasslands [56]. In this study, shifts in ω of NFs and NGs from current to future conditions was used to represent changes in integrative properties of NFs and NGs such as land cover, vegetation cover, and climate conditions. Under stationary climatic conditions, the parameter omega in Fu's equation is primarily governed by surface vegetation. However, long-term changes in climate normal (e.g., 30-year average annual temperature, precipitation, and relative humidity) beat changes in hydroclimatic indices in the Budyko space. Thus, the current study characterized these climatic influences while surface vegetation is kept the same as the control period.

Results
Hydroclimatic variables of NGs and NFs may differently respond to climate change. This section summarizes changes in regional hydroclimatic conditions from the past period  to future period (2070-2099) for groups of NFs and NGs in eight National Forest service regions. Changes in hydroclimatic conditions of NFs and NGs were assessed by characterizing the direction and magnitude of shifts in the Budyko space. Finally, Fu's equation was applied to determine hotspot regions with the highest change in their basin characteristics. The results for NFs and NGs reported below are precisely for HUC8 watersheds containing all or part of NFs and NGs. Table 2 summarizes average precipitation (PCP), water yield (YIELD), and PET for NFs and NGs and the entire CONUS for the current (1986-2015) period from the combination of DAYMET and PRISM dataset. Values for HUC8s that include any proportion of NFs and NGs lands are area-weighted by fractional HUC area to make a comparison. The 30-year average of precipitation and water yield of NFs are above the CONUS average. However, the amounts of precipitation and water yield in NGs are considerably lower than the CONUS average. Although most NGs are in areas with low precipitation and water yield compared to the NFs, the 30-year averages of potential evapotranspiration of NFs and NGs are relatively close (Table 2). A higher proportion of precipitation was converted to water yield in the NFs than in the CONUS as a whole, confirming that NFs tend to be an important source of water yield.   Projected changes in water yield among the three scenarios largely reflect the changes precipitation, with increases projected for the MIDDLE and WET scenarios in the CONUS as well as in NFs and NGs, and decreases projected for the DRY scenario for CONUS and NGs. The exception is NFs under the DRY scenario, where average water yield decreases despite a small average increase in precipitation.

Changes in Hydroclimatic Variables
The percentage change in average potential evapotranspiration is small compared to changes in the other hydroclimatic variables. While potential evapotranspiration increases under the MIDDLE and DRY scenarios for all NFs and NGs, it consistently decreases under the WET scenario. Although NGs have higher changes in the average potential evapotranspiration and water yield under MIDDLE and DRY scenarios, NFs see greater changes in potential evapotranspiration and water yield under the WET scenario. Figures 3 and 4 illustrate spatial changes from the current to future time period in precipitation, potential evapotranspiration, and water yield of USFS land (NFs and NGs) of the eight USF regions of the CONUS. As shown in Figure 3, the direction and magnitude of changes in all three variables of most regions highly depends on the selected future climate model. Regarding the direction in average precipitation, for instance, in the Southwestern, Rocky Mountain, and Southern regions precipitation increases for the WET scenario and decreases under the DRY scenario, but consistently increases in the Northern, Intermountain, Pacific Southwest, Pacific Northwest, and Eastern regions.       Despite the variability across future scenario, at least two findings stand out. First, precipitation is projected to increase in most regions under all three scenarios but decrease substantially in two regions (Southwestern and Southern) under the DRY scenario. Second, water yield is projected to increase under the WET and MIDDLE scenarios in all but one case (Southwestern under the MIDDLE scenario) and to decrease in five regions under the DRY scenario, with large percentage decreases in the Southwestern and Southern regions. Note, the DRY, MIDDLE, and WET climate models reflect conditions in general, not necessarily in every region. The NFs in high mountainous regions will experience higher changes in hydroclimatic variables than NGs and NFs at lower elevations ( Figure 4). However, the direction of change varies based on the selected future climate model. Precipitation, PET, and water yield for the high mountainous NFs are likely to experience the largest changes in response to climate change. Table S1 provides percentage changes in precipitation, potential evapotranspiration, and water yield of U.S. national forests and grasslands under WET, MIDDLE, and DRY climate scenarios in details.
The Wilcoxon signed-ranked test was used to assess the significance of the future changes in PCP, PET, and YIELD under DRY, MIDDLE, and WET scenarios at the 5% significance level (Table S2). Under the WET scenario, change in PCP is significant in most NFs and NGs within the Pacific Southwest, Intermountain, Southwestern, Eastern, and Southern regions. Change in PET is statistically significant for most NFs and NGs located in the Rocky Mountain and Northern regions. Changes in YIELD is significant in most regions within the Pacific Southwest and Pacific Northwest. Under the MIDDLE scenario, change in PCP is statistically significant in most NFs and NGs within the Northern and Intermountain regions. Change in PET is significant for most NFs and NGs located in the Southwestern, Southern, and Eastern regions. Changes in YIELD is significant in most regions within the Southwestern and Pacific Northwest regions. Under the DRY scenario, change in PCP is significant in most NFs and NGs within the Southern, Rocky Mountain, and Southwestern regions. Change in PET is statistically significant for most NFs and NGs located in the Eastern, Southern, Rocky Mountain, Southern, and Intermountain regions. Changes in YIELD is significant in most regions within the Southern, Southwestern, and Pacific Northwest regions.  The NFs in high mountainous regions will experience higher changes in hydroclimatic variables than NGs and NFs at lower elevations ( Figure 4). However, the direction of change varies based on the selected future climate model. Precipitation, PET, and water yield for the high mountainous NFs are likely to experience the largest changes in response to climate change. Table S1 provides percentage changes in precipitation, potential evapotranspiration, and water yield of U.S. national forests and grasslands under WET, MIDDLE, and DRY climate scenarios in details.

Movements in the Budyko Space
The Wilcoxon signed-ranked test was used to assess the significance of the future changes in PCP, PET, and YIELD under DRY, MIDDLE, and WET scenarios at the 5% significance level (Table S2). Under the WET scenario, change in PCP is significant in most NFs and NGs within the Pacific Southwest, Intermountain, Southwestern, Eastern, and Southern regions. Change in PET is statistically significant for most NFs and NGs located in the Rocky Mountain and Northern regions. Changes in YIELD is significant in most regions within the Pacific Southwest and Pacific Northwest. Under the MIDDLE scenario, change in PCP is statistically significant in most NFs and NGs within the Northern and Intermountain regions. Change in PET is significant for most NFs and NGs located in the Southwestern, Southern, and Eastern regions. Changes in YIELD is significant in most regions within the Southwestern and Pacific Northwest regions. Under the DRY scenario, change in PCP is significant in most NFs and NGs within the Southern, Rocky Mountain, and Southwestern regions. Change in PET is statistically significant for most NFs and NGs located in the Eastern, Southern, Rocky Mountain, Southern, and Intermountain regions. Changes in YIELD is significant in most regions within the Southern, Southwestern, and Pacific Northwest regions. Figure 5 shows the current aridity and evaporative indices of the NGs and NFs, and Figure 6 presents the averages of the current aridity and evaporative indices of the eight USFS regions. The NFs and NGs within the Southwestern, Rocky Mountain, and Intermountain regions currently have the most arid conditions along with the highest evaporative index. In contrast, the Southern, Eastern, and Pacific Northwest regions have the lowest aridity index and, along with the Pacific Southwest region, the lowest evaporative index, meaning that USFS and in these regions currently have generally a less arid climate with higher freshwater availability.     Although the USFS lands of the Northern, Pacific Northwest, Southern, and Eastern regions, which are least arid under current conditions, are projected to incur only minor changes in aridity index under all three scenarios, the Southwestern, Intermountain, and Pacific Southwest regions are more likely to experience higher change in aridity index, especially under selected scenarios. The Intermountain, Pacific Southwest, and Pacific Northwest regions consistently are projected to experience a decrease in aridity index, meaning that they are likely to be less arid in the future. However, the Southwestern region under the DRY climate scenario shows the highest increase in aridity index, of approximately 5% meaning that this USFS region is vulnerable to a more highly arid condition in the future.

Movements in the Budyko Space
An increase in the evaporative index indicates that more of the available precipitation is evaporated than under current conditions. The generally increasing evaporative index     Although the USFS lands of the Northern, Pacific Northwest, Southern, and Eastern regions, which are least arid under current conditions, are projected to incur only minor changes in aridity index under all three scenarios, the Southwestern, Intermountain, and Pacific Southwest regions are more likely to experience higher change in aridity index, especially under selected scenarios. The Intermountain, Pacific Southwest, and Pacific Northwest regions consistently are projected to experience a decrease in aridity index, meaning that they are likely to be less arid in the future. However, the Southwestern region under the DRY climate scenario shows the highest increase in aridity index, of approximately 5% meaning that this USFS region is vulnerable to a more highly arid condition in the future.
An increase in the evaporative index indicates that more of the available precipitation is evaporated than under current conditions. The generally increasing evaporative index  among the regions and scenarios reflects the overall increase in projected temperature and thus energy for evaporation (plus in some cases the decrease in precipitation). The significance of the future changes in the aridity and evaporative indices was assessed using the Wilcoxon signed-ranked test at the 5% significance level under DRY, MIDDLE, and WET scenarios (Table S3)  Although the USFS lands of the Northern, Pacific Northwest, Southern, and Eastern regions, which are least arid under current conditions, are projected to incur only minor changes in aridity index under all three scenarios, the Southwestern, Intermountain, and Pacific Southwest regions are more likely to experience higher change in aridity index, especially under selected scenarios. The Intermountain, Pacific Southwest, and Pacific Northwest regions consistently are projected to experience a decrease in aridity index, meaning that they are likely to be less arid in the future. However, the Southwestern region under the DRY climate scenario shows the highest increase in aridity index, of approximately 5% meaning that this USFS region is vulnerable to a more highly arid condition in the future.
An increase in the evaporative index indicates that more of the available precipitation is evaporated than under current conditions. The generally increasing evaporative index among the regions and scenarios reflects the overall increase in projected temperature and thus energy for evaporation (plus in some cases the decrease in precipitation).
The significance of the future changes in the aridity and evaporative indices was assessed using the Wilcoxon signed-ranked test at the 5% significance level under DRY, MIDDLE, and WET scenarios (Table S3). Under the WET scenario, change in aridity index is significant in most NFs and NGs within the Rocky Mountain, Intermountain, Pacific Southwest, and Northern regions. Change in evaporative index is statistically significant for most NFs and NGs located in the Pacific Southwest, Pacific Northwest, and Northern regions. Under the MIDDLE scenario, change in aridity index is statistically significant in most NFs and NGs within the Northern, Intermountain, and Pacific Southwest regions. Change in evaporative index is significant for most NFs and NGs located in Rocky Mountain, Southwestern, and Northern regions. Under the DRY scenario, change in aridity index is statistically significant in most NFs and NGs within the Southwestern, Southern, Rocky Mountain regions. Change in evaporative index is statistically significant for most NFs and NGs located in the Southern, Eastern, Rocky Mountain, Northern, Southwestern, and Pacific Northwest regions. Figure 8 (left panel) shows in the Budyko space the changes depicted in Figure 8 (right panel) using the Forest Service region numbers. The hydroclimatology of the Northern, Intermountain, Pacific Southwest, and Pacific Northwest regions show the most consistency across the three climate scenarios. While the Northern, Intermountain, and Pacific Northwest regions are projected to experience less arid condition but with lower freshwater availability, the Pacific Southwest region is projected (under the WET and DRY scenarios) to move to the lower-left, meaning less arid conditions and higher freshwater availability. Long-term changes in the hydroclimatology of the Rocky Mountain, Southwestern, Southern, and Eastern regions vary in aridity index across the three scenarios, with aridity consistently decreasing under the WET scenario and increasing under the DRY scenario, with variable results under the MIDDLE scenario.
Under the WET scenario, all forest service regions move to the left meaning less arid conditions in the future. However, the USFS regions located in the Central and Western United States have a higher magnitude of changes indicating that these regions are more sensitive to climate change.
Under the MIDDLE scenario, similar to the WET scenario, most USFS regions move to the left. However, the magnitude of movements is smaller than in the WET scenario. Besides, the NFs and NGs within the Southwestern region behave differently under the MIDDLE scenario. Movements within the Southwestern region are more likely to the right indicating that this region may experience drier hydroclimatic conditions in the future.
Under the DRY scenario, the Eastern, Southern, Rocky Mountain, and Southwestern regions consistently show drier conditions with less water yield with the highest magnitude in the Southwestern and Southern regions, respectively. NFs and NGs within the Northern, Pacific Northwest and Intermountain regions behave differently under the DRY scenario. However, movements within the Pacific Southwest region are to the left indicating that this region is more likely to experience less arid hydroclimatic conditions in the future.
Our findings are sensitive to the choice of climate model. The results presented in this study under DRY, MIDDLE, and WET climate models are indicative of the effect the uncertainty associated with future climate change impacts on the hydroclimatology of U.S. megaregions. Under all three scenarios, the Southwestern region has a high magnitude of change meaning that its hydroclimatic conditions have high sensitivity to future climate change. While the Southwestern region is getting warmer and drier with a high magnitude of change under the DRY and MIDDLE scenarios, it is more likely to experience wetter and less arid hydroclimate conditions with high magnitude under the WET scenario.
In addition to the average magnitude and direction of each USFS region, we assessed the distribution of changes in frequency, magnitude, and direction of each region using the wind rose diagram (Figures S3-S5). The wind rose diagram visualizes the summary of movements in the Budyko space for the HUC8 river basins within NFs and NGs. This type of diagram has been used in global hydroclimatic change assessments [58]. Figures S3-S5 provide wind rose diagrams of movements in the Budyko space under WET, MIDDLE, and DRY scenarios, respectively. Figure Table S4 provides changes in the direction and magnitude of each U.S. NFs and NGs in the Budyko space under the WET, MIDDLE, and DRY scenarios.

Changes in Basin Characteristics
Changes in integrative properties of basins in response to hydroclimatic shifts can be represented by changes in ω of Fu's equation (Equation (5)). Change over time in ω of a HUC8 could demonstrate changes in integrative properties of NFs and NGs such as differences in land cover, vegetation cover, and climate. Table 4 represents current and future average ω of NFs and NGs across the CONUS. The change in ω from the current to the future period is small under all three climate scenarios. Under all three climate scenarios, NFs have larger changes in ω than NGs, meaning that NFs are more likely to experience higher changes in their basin characteristics compared to NGs assuming that a specific amount of change in omega has the same effect on likelihood of shift in characteristics in all USFS regions and both forests and grasslands. Additionally, under the WET and MIDDLE climate scenarios, an increase in ω for NFs and NGs is above the CONUS average, indicating that they may experience higher shifts in their properties compared to other regions. Under the DRY scenario, although NFs have higher changes in ω, change in ω of NGs is below the CONUS average. NFs under all three climate scenarios and NGs under the WET, and MIDDLE climate scenarios experience increasing ω while NGs and CONUS average have minor decrease in the ω under the DRY climate scenario.  Table 5 provides current and future ω of Fu's equation for the eight USFS regions under the WET, MIDDLE, and DRY scenarios. The ω for eight regions is nearly always increasing under three climate scenarios, the main exception being the decrease in ω of the Southwestern region under the DRY climate scenario. The ω in the Eastern, Southern, Pacific Northwest, and Pacific Southwest regions changes only slightly, indicating that these regions are more likely to experience only minor shifts in the basin characteristics.
However, the Southwestern, Northern, Intermountain, and Rocky Mountain regions are more likely to experience higher shifts in their physiographic and ecological characteristics. In this study, we assessed changes in ω from current conditions to future conditions to find which NFs and NGs are more prone to be affected by long-term hydroclimatic changes. The NFs and NGs with larger changes in ω are likely to experience higher moving in the direction of a shift in integrative properties of basin. Thus, we figured out which US Forest Service regions have the highest sensitivity to experience considerable shifts in the structure and composition of forests and grasslands.

Discussion
Future hydroclimatic changes may have various consequences for forests. For example, while a temperature rise is likely to decrease growth in lower elevations, it might increase growth in higher elevations [16]. Thus, future studies are needed to determine the relationships between hydroclimataology of US river basins and the composition, distribution, and growth of various types of forests and grasslands.
The results are subject to several sources of uncertainty, including those introduced by the choice of the climate models and the hydrological models we used. Climate models are uncertain because they highly depend on future anthropogenic and natural forcing scenarios. Besides, the future climate projections can be affected by the existence of internal climate variability and incomplete understanding and imprecise climate models [59,60]. Using ensemble simulation rather than specific climate scenarios can be a prospect for this study to assess decrease in climate model variability.
In addition to the uncertainty in the choice of climate models, some uncertainties remain, related to the hydrological model. The VIC model may not capture all physical basin characteristics, water management regulations, landcover changes [36]. Besides, there are some uncertainties associated with parameters of the VIC model and structural deficiencies of the model simulation [61][62][63]. Using multiple hydrologic models can be a prospect for this study to evaluate the modeling uncertainty associated with the role of hydrological models in hydroclimatic changes assessment [35]. Furthermore, estimated changes in ω can originate from the assumptions associated to Fu's equation.
Despite these limitations, comparisons of differences in long-term hydroclimatic variables from the current period  to the future period (2070-2099) are more likely to be accurately estimated than are the absolute amounts from which the differences are computed. The structural deficiencies and assumptions in the VIC model can be addressed by aggregating outputs over longer time periods such as 30-year average [61].
Additionally, we focused on a range of possible future climate conditions, from driest to wettest, to account for current uncertainty about long-term future climatic conditions. Although variability in the projection of future climate conditions can be affected by some sources of uncertainty and variability in climate change scenarios, projected shifts in hydroclimatic conditions of some NFs and NGs showed some consistency across climate change models in terms of the direction and magnitude of changes. The results suggest a need for flexible and regionalized forest management strategies under these scenariosensitive, spatially-heterogeneous projections of future hydroclimatic conditions.

Summary and Conclusions
Understanding change in hydroclimatic conditions of NFs and NGs, and how those changes will impact natural resources, have risen in importance as our understanding of climate change has improved [15]. This study focuses on shifts in long-term regional climate and hydrologic trends of two important sources of natural resources in the United States, National Forests and National Grasslands.
Changes in current climatic conditions (1986-2015) were first derived from a combination of Daymet and PRISM with wind speed data from the NARR dataset at the 4 km by 4 km spatial scale. Changes in future climatic variables (2070-2099) were obtained from the MACA climate dataset at 4 km by 4 km spatial scale for three climate models, yielding projections for three quite different climate scenarios. The forcing parameters from that data were then used as inputs to the VIC hydrological model to obtain hydrological estimates at a daily timestep. The output hydroclimatic variables were aggregated at the HUC8 river basin scale and annual timestep. The aggregated hydroclimatic variables were then used to estimate current and future hydroclimatic conditions of NFs and NGs across the CONUS. Long-term changes in those conditions were depicted as movement in the Budyko space. Finally, the Fu's equation was implemented to assess shifts in integrative basin characteristics of the NFs and NGs.
The response of basins to climate change varies from NFs to NGs and from one climate scenario to another. However, there was some significant consistency in regional long-term hydroclimatic changes on NFs and NGs. Both NFs and NGs are projected to experience less arid hydroclimatic conditions under both the WET and MIDDLE climate scenarios, in contrast to generally more arid conditions under the DRY scenario. Under the DRY climate scenario, the Northern, Eastern, Southern, Rocky Mountain, and Southwestern regions all show drier conditions with less water yield. Under all three climate scenarios, NFs are likely to experience larger changes in the basin characteristics than are NGs. The hydroclimatology of the Southwestern Forest Service region is found to have the highest sensitivity to climate changes across the CONUS, with especially high aridity under the DRY climate scenario. The dramatic differences in results between the WET and DRY scenarios demonstrates current uncertainty about the future impacts of climate change.
Long-term hydroclimatic shifts may fundamentally change conditions of US forests and grasslands. This finding highlights the need to incorporate climate change impacts into forest and grassland resource management and planning. This study can be used as a roadmap for land and water managers to develop options for adapting to future climate change, such as increasing landscape diversity and facilitating biological diversity [16]. This study can help environmental scientists and water managers to mitigate the negative economic, social, and environmental consequences of climate change on US NFs and NGs.
Supplementary Materials: The following are available online at https://www.mdpi.com/1999-4 907/12/2/139/s1, Figure S1: 30-year average of precipitation (1986-2015) for HUC08 watersheds; Figure S2: Comparison of simulated and observed annual water yield for 1986-2015 period for each NFs and NGs; Figure S3: Wind rose diagrams of movements in the Budyko under the WET climate model; Figure S4: Wind rose diagrams of movements in the Budyko under the MIDDLE climate model; Figure S5: Wind rose diagrams of movements in the Budyko under the DRY climate model; Figure S6: Wind rose diagrams of movements in the Budyko for NFs and NGs under DRY, MIDDLE and WET climate models; Table S1: Changes in precipitation, PET and water yield of U.S. national forests and grasslands; Table S2: P-values of Wilcoxon signed-ranked test at the 5% level (PCP, PET, Yield); Table S3: P-values of Wilcoxon signed-ranked test at the 5% level (Aridity, Evaporative); Table  S4: Changes in the direction and magnitude of the Budyko space for the U.S. national forests and grasslands.