Changes of Streamflow Caused by Early Start of Growing Season in Nevada, United States

The fluctuation of streamflow in snowmelt-dominated watersheds may be an indicator of climate change. However, the relationship between the start of growing season (SOS) and the streamflow in snowmelt-dominated watersheds is not clear. In this study, we update the Coupled Hydro-Ecological Simulation System (CHESS) model by incorporating the Growing Season Index (GSI) module to estimate the start of the growing season. The updated CHESS model is then used to calculate the streamflow in the Cleve Creek, Incline Creek and Twin River watersheds located in Nevada in the United States from 1981 to 2017. This updated CHESS can be applied in any regions that are suitable for deciduous vegetation. The streamflow in the static and dynamic scheme in the three watersheds have been simulated between 1981 and 2017 with the NS of 0.52 and 0.80 in the Cleve Creek, 0.46 and 0.75 in the Incline Creek, and 0.42 and 0.70 in the Twin River watersheds, respectively. The results illustrate that the SOS have come around 3–5 weeks earlier during the last 37 years. The results illustrate a high correlation between the temperature and the timing of the SOS. Early SOS leads to a substantial increase in total annual transpiration. An increase in annual transpiration can reduce aquifer recharge and increase cumulative growing season soil moisture deficit. Comparing to the streamflow without vegetation, the streamflow with vegetation is smaller due to transpiration. As the SOS comes earlier, the peaks of the streamflow with vegetation also come earlier. If the shifts in SOS continue, the effects on annual rates of transpiration can be significant, which may reduce the risk of flooding during snowmelt. On the other hand, earlier SOS may cause soil moisture to decline during summer, which would increase the drought stress in trees and the risk of wildfires and insect infestation.


Introduction
The dominant form of annual precipitation in high-elevation mountains in the western United States is winter and spring snow [1][2][3]. Most of the snow accumulation in seasonal snowpack represents a natural storage reservoir [4,5]. Because most of the precipitation is concentrated in spring and winter, the western United States relies heavily on natural and man-made storage for sufficient water for agriculture, industry, and drinking during the dry summer and fall seasons [6]. Snow is a critical component of the hydrologic cycle for vegetation in the western United States and always makes the largest contribution for conifers which are photosynthetically active in the spring and early summer when meltwater is available [7]. Forest transpiration and growth are partially determined by the start of growing season (SOS) since it is limited by water availability. Hydrological models incorporating biophysical controls on forest transpiration need an accurate prediction of SOS. Process models represent snowpack as a state variable containing a mass of water that accumulates during the winter months and releases during the spring months. Since the meltwater is a vital source for generating spring runoff and re-saturating the soil to provide water for transpiration and growth, it is important to investigate the mechanisms of flow generation in snow-dominated mountain watersheds for predicting the impacts of climate change on the availability of water resources.
Previous research demonstrates the average annual temperature has increased 1-2 • C since the 1940s over western America, especially during the winter and spring seasons [8][9][10][11]. These changes have a complicated influence of climate variability on terrestrial ecosystems [12]. The worldwide warming directly influences the hydrological cycle since evapotranspiration (ET) is a critical water loss in the water balance of ecosystems [13], especially in semi-arid and arid regions [14]. Ecosystems determine streamflow quantity [15] and timing [16] by altering the ET and groundwater recharge processes [17]. Generally speaking, deforestation decreases ET and thus increases streamflow [18] and peak flow rates [19], while afforestation increases ET [20], raises shallow groundwater tables and reduces streamflow quantities [21].
It is generally assumed that warming will make SOS occur earlier in boreal and temperate ecosystems, where wintertime temperature is metabolically limiting. A simple assumption of increased primary productivity is made from the fact that more days are available for carbon assimilation and biomass growth. Several previous models have demonstrated a similar trend of the relationship between the SOS and net productivity in deciduous forests. For instance, Fu et al. [22] propose that warmer spring temperatures lead to an earlier start of net carbon uptake activity and higher spring and annual net ecosystem productivity (NEP) values in a deciduous broadleaf forest (DBF) and evergreen forest (EF). Hmimina et al. [23] illustrate that the early SOS and increasing photosynthetic activity results in a predominant greening trend over 42.0% of the vegetated area and translate to a 20.9% gain in gross primary productivity during last three decades.
However, other studies have suggested somewhat contradictory results. For example, Dunn et al. [24] find that ecosystem carbon exchange responds strongly to air temperature, moisture status, potential evapotranspiration, and summertime solar radiation. The seasonal cycle of ecosystem respiration significantly lags that of photosynthesis, limits by the rate of soil thaw and the slow drainage of the soil column. Net uptake of carbon is enhanced and respiration is inhibited by multiple years of rainfall over evaporative demand. Winchell [25] find that earlier snowmelt associated with climate warming, counterintuitively, lead to colder atmospheric temperatures during the snow ablation period and concomitantly reduce rates of net carbon uptake. Winchell concludes that when the snowmelt period occurs one month earlier than average, the forest experiences an ablation period mean air temperature of 2.7 • C, approximately 2 • C colder than an ablation period that occurs one month later than average. As most climate models project peak snow water equivalent (SWE) to occur earlier under various warming scenarios, they expect to see a decrease in total growing season carbon uptake if the post-snowmelt period is unable to compensate for the decrease in ablation period carbon uptake. Wohlfahrt et al. [26] show that reductions in the quantity and duration of daylight associated with earlier snowmelts are responsible for diminishing returns, in terms of carbon gain from early SOS.
Another uncertainty is how SOS heterogeneity affects streamflow. Streamflow may or may not increase as the regional climate warms. A common approach simulating hydrologic budgets is used for quantifying these uncertainties by dividing a heterogeneous area into smaller, more homogeneous subsites, and executing the individual model on each sub-site [27,28]. The Regional Hydro-Ecological Simulation System (RHESSys) model is built to compute a spatially varying snowpack for regional applications and is able to explicitly simulate the effects of topographic and subsurface heterogeneities on downslope redistribution of water [29]. Additionally, this model is widely used to identify saturated areas that produce runoff [30], evaluate irrigation systems [31], and examine flood potential associated with disturbances such as deforestation [32]. Based on the RHESSys model, Tang et al. [33,34] develop Coupled Hydro-Ecological Simulation System (CHESS) model by building a new model structure simulating the integrated water, carbon, nutrient dynamics and vegetation growth at regional and watershed scales. For land-type cells, CHESS explicitly simulates the routing of the surface and subsurface flows through space. For channel-type cells, the surface and subsurface flows are first channelized and then routed. However, both the RHESSys and CHESS models do not consider the influences of SOS on the terrestrial ecosystem.
The object of this study is to investigate the impact of climate warming-induced shifts in phenology on the hydrology and ecology of the vegetation, particularly focusing on hydrological changes. We first update the previous CHESS model by incorporating the daily Growing Season Index (GSI) module, which can be used as a tool to estimate SOS in a year. The calculated SOS is then validated by the Normalized Difference Vegetation Index (NDVI) obtained from the satellite. Two cases with the application of the updated CHESS model are performed on the three watersheds (Cleve Creek, Incline Creek, and Twin River) located in Nevada. Note that both SOS and temperature increase have impact on streamflow and temperature increase has a strong relationship with SOS. The impact on streamflow caused by temperature increase should be isolated. Therefore, the first case is performed in the watersheds without vegetation under two different temperatures and the second case is performed with vegetation under the same conditions.

Study Regions
Three mountain watersheds (Cleve Creek (Figure 1a), Incline Creek ( Figure 1b) and Twin River (Figure 1c)) located in Nevada, U.S., are investigated for the following reasons: (1) these mountain watersheds have warm summers and cold winters with annual mean temperatures ranging from 1.9 to 7.1 • C [35]; (2) the annual mean temperature of these watersheds have an overall increasing trend in recent years [10,36] as shown in Figure 2. The annual mean temperature approximately increases 1.5 • C in Cleve Creek watershed and Incline Creek watershed, individually. At Twin River watershed, the annual mean daily temperature increases around 4 • C; (3) long-term historical streamflow observations recorded by the US Geological Survey gauge stations can be employed to calibrate and evaluate model simulations; (4) these watersheds covered by multiple types of vegetation such as evergreen, deciduous, shrub and grass are suitable for investigating the connections between ecological and hydrological processes.

Weather Forcing Data
By incorporating the module of estimating the daily GSI to the updated CHESS model, the time series of daily maximum and minimum temperature ( • C), as well as total precipitation (mm) are required as inputs. The daily weather records from the nearest SNOwpack TELemetry (SNOTEL) stations are employed to derive the weather forcing data since there is no weather station located in these watersheds. Because elevation differences exist between the SNOTEL station and the corresponding watershed, we apply the local lapse rate for adjusting the effect of elevation on daily maximum temperature, daily minimum temperature and daily total precipitation, individually. Table 1 shows the parameters we used in this study.     The separation of precipitation into snow and rainfall is necessary for a water balance calculation in a cold region. The separation step is essential to determine whether water is available for runoff and soil infiltration, or stored as snow. The snow proportion is a function of daily minimum and maximum threshold temperatures. If the average temperature is greater than the maximum threshold temperature, then all precipitation is considered as rainfall. If the average temperature is smaller than the minimum threshold temperature, then the precipitation is considered to fall as snow. When the average temperature is between the minimum and maximum threshold temperature, the rainfall proportion of the mixed precipitation is calculated by dividing the difference between the maximum threshold and average temperature by the difference between the minimum and maximum threshold temperatures. The other portion of the precipitation is considered as snow. The equation is presented as, where P r represents the portion of rainfall. T Ave is the daily average temperature. T TMin and T TMax is the minimum and maximum threshold temperature in degree Celsius, respectively. For all tests, T TMin = −0.1 • C and T TMax = 0.1 • C [36].

Estimating Start of Growing Season (SOS)
There are numerous studies to investigate phenology shift such as species-level observations, satellite remote sensing of ecosystem production and atmospheric monitoring of carbon dioxide concentrations [37][38][39]. Previous studies have shown that minimum temperature, vapor pressure deficit (VPD) and photoperiod are the common set of meteorological variables leading to the variation in the seasonal phenology [40,41]. The product of these three indicators develops a combined model of the GSI that is calculated daily and integrated as a 21-day running average to avoid reaction to short-term changes in environmental conditions [42].
Minimum temperature: Many of the biochemical processes of plants are sensitive to low temperatures [43]. Although ambient air temperatures certainly have an impact on plant growth, constraints on phenology are more closely related to restrictions on water uptake by roots when soil temperatures are suboptimal [44]. The minimum temperature is a strong indicator for predicting the onset of plant growth. To incorporate a range of species, we choose a range encompassed by a lower minimum temperature threshold of −2 • C (T MMin ) and an upper threshold of 5 • C (T MMax ) [43]. To eliminate the effect of the unit, the minimum temperature index (iT min ) is presented as follows [42]: where iT Min is the daily indicator for minimum temperature and is bounded between 0 and 1 and T Min is the observed daily minimum temperature in degrees Celsius. For all tests, T MMin = −2 • C and T MMax = 5 • C [42]. Vapor pressure deficit (VPD): water stress causes partial to complete stomatal closure, reduces leaf development rate, induces the shedding of leaves and slows or halts cell division. Since it is complicated to describe all the processes, the VPD of the atmosphere is selected as an index of the evaporative demand as a surrogate in this study. At low values, latent heat losses are not going to exceed available water. At high values, photosynthesis and growth are significantly restricted. Previous studies recommend that VPD less than 900 Pa should exert little effect on stomata whereas values greater than 4100 Pa generally are sufficient to force complete stomatal closure, even when the soils are moist [45,46]. Although these restrictions are various by locations and species, we select a widely used set of parameters (VPD Min = 900 Pa and VPD Max = 4100 Pa) in this research [41]. The VPD index (iVPD) is provided as follows: where iVPD is the daily indicator for VPD and is bounded between 0 and 1 and VPD is the observed daily VPD in Pascals. For all tests, VPD Min = 900 Pa and VPD Max = 4100 Pa. Photoperiod: photoperiod is a credible annual climatic cue for plants because it can be kept within a small range in different years at a given location. Studies have shown that photoperiod is critical for both leaf flush and leaf senescence throughout the world [47]. The photoperiod index (iPhoto) is given as follows: where iPhoto is the daily photoperiod indicator and Photo is the daily photoperiod in seconds. For all tests, Photo Min = 36,000 s (10 h) and Photo Max = 39,600 s (11 h) [41]. GSI is a daily indicator of the relative constraints to foliar canopy development or maintenance due to climatic limits. The product of the three indices quantifies the GSI that is calculated daily and integrated as a 21-day running average to avoid reaction to short-term changes in environmental conditions [42]. The daily metric iGSI is presented as follows: where iGSI is the daily GSI, iT min , iVPD and iPhoto are the minimum temperature indicator, the VPD indicator and the photoperiod indicator, respectively. The daily GSI is then calculated as the 21-day moving average of the daily indicator, iGSI. The moving average serves to avoid single extreme events from prematurely triggering canopy changes.
The iGSI values are computed daily at each watershed using the parameters in Equations (2)-(5). The satellite-derived NDVI is used for validating the mean iGSI for the corresponding 8-day satellite data composite period. Phenology observations are averaged over all species for each observation date. The average date of leaf flushing is defined when the average canopy leaf area exceeds 10% of the seasonal maximum. Conversely, the average date of leaf coloration is defined as the date when leaf coloration exceeded 90%. The SOS is defined as the time when iGSI exceeds 0.5 in the spring and the EOS (end of growing season) is defined as the time when iGSI drops to 0.5 in the fall [41].
To find the impact of SOS on streamflow, the influences caused by temperature increase should be isolated. Two cases with the application of the updated CHESS model in these three watersheds are proposed. In each case, two scenarios at different temperatures are performed. One scenario is under the normal temperature condition. The other scenario is simulated under the condition that the temperature is 2 • C higher than the normal temperature throughout the whole year. The simulation is performed in the entire period from 1981 to 2017. Since the precipitation is abundant and the difference between the normal temperature and high temperature scenario is obvious for the three watersheds in 1995, we select the year 1995 to illustrate the specific variations of soil water storage (SWS), evaporation, transpiration, and streamflow when the temperature increases. The first case simulates the evaporation, SWS, and streamflow without vegetation in 1995 and is used as a baseline. Because the plant is removed in this case, the transpiration is 0 in this simulation. The second case is used for investigating the impact of varying SOS on evaporation, SWS, and streamflow with vegetation in 1995. The other conditions are the same as the previous case. Figure 3 shows the percentage of snow in precipitation in the recent 37 years. To reduce the high-frequency variability, we applied a 10-year moving average to remove statistical uncertainties and smooth the annual value of each phenological event. Although snow occupies the main fraction in several years such as 2001-2003 in these three watersheds, a shift from snow towards rainfall becomes the dominating trend related to the potential global warming. The percentage of snow in Cleve Creek, Incline Creek and Twin River watersheds decreases by more than 5%, 10%, and 12% as the annual temperature rises around 1.5 • C, 1.5 • C and 4.0 • C, respectively.

Model Validation
The calculated SWE from 1981 to 2017 are validated with the observed data in Figure 4. We notice that the updated CHESS model can simulate the variation of SWE in the three watersheds well. However, there are still several discrepancies between the simulated results and the observed data at the peaks (e.g., Cleve Creek watersheds in 2005, Incline Creek watersheds in 1991 and Twin River watersheds in 1993). The discrepancies may be caused by the spatial and temporal variation in the physical characteristics of the snow pack (e.g., the snow density, structures, and snow grain size) [48]. The algorithm may overestimate or underestimate the parameter in the updated CHESS model based on the hypothesis that the snow grain size is uniform both temporally and spatially.
Time-series plots of iGSI and NDVI values for these three watersheds are also calculated from 1981 to 2017. Since the results of every year in this period are similar, we only choose one year (2001) to display as an example in Figure 5. The NDVI obtained from MODerate-resolution Imaging Spectroradiometer (MODIS) are used for validating the calculated iGSI. The image of NDVI contains 16 days of data, with the image's data occurring on the 9th day of the 16-day period. This product combines the data choosing the best representative pixel from the 16-day period [49]. In Figure 5, the iGSI-predicted dynamics of growing season for the three watersheds agree with the observed NDVI well especially the start date and the end date. In these three cases, the updated CHESS model can accurately predict a rise at around the 130th day and a drop at around the 300th day. In other words, the updated CHESS model can precisely find SOS and EOS in a year. However, obvious difference can be observed within the growing season. This is because the updated CHESS uses the average parameter to calculate the iGSI for the diversity of plants.             The interplay between precipitation and streamflow in the dynamic and static scheme is plotted in Figure 8. The relationship between the annual precipitation and streamflow is similar for the three watersheds. The increment speed of annual streamflow in the static scheme is faster than that in the dynamic scheme as the annual precipitation increases. The difference between the static and dynamic scheme is small in the low precipitation. The percentage of streamflow distributes within 70-87% and 31-38% in the static and dynamic scheme, respectively, while the percentage of the streamflow distributes within 62-76% and 25-33% in the high precipitation, individually. The results illustrate that the streamflow is overestimated in the static scheme because the static scheme underestimates the water consumption caused by the extended growing season. When the precipitation increases and provides sufficient water to the vegetation, the difference between these two schemes become large, which means the transpiration stands a larger proportion of precipitation.   Creek, and (c) Twin River watersheds between the static and dynamic scheme. Figure 9 presents the influence of the temperature on the evaporation. When the temperature becomes higher, the evaporation increases more. Both scenarios have the maximum evaporations around July in the three watersheds. Most of the time, the evaporation in the high-temperature scenario is greater than the evaporation in the low-temperature scenario due to more snowmelt and rainfall occurring in the high-temperature scenario.  Figure 9 presents the influence of the temperature on the evaporation. When the temperature becomes higher, the evaporation increases more. Both scenarios have the maximum evaporations around July in the three watersheds. Most of the time, the evaporation in the high-temperature scenario is greater than the evaporation in the low-temperature scenario due to more snowmelt and rainfall occurring in the high-temperature scenario.  Most of the SWS in the high-temperature scenario before May is greater than that in the low-temperature scenario. This is due to higher infiltration caused by higher temperature. When precipitation in summer and autumn decreases, the soil storage, however, remains almost the same. The SWS has a significant increase in the high-temperature scenario after October 1 st in the Incline Creek and Twin River watersheds because new snow falls.  Most of the SWS in the high-temperature scenario before May is greater than that in the low-temperature scenario. This is due to higher infiltration caused by higher temperature. When precipitation in summer and autumn decreases, the soil storage, however, remains almost the same. The SWS has a significant increase in the high-temperature scenario after October 1st in the Incline Creek and Twin River watersheds because new snow falls. The comparison of streamflow in these two scenarios is plotted in Figure 10d-f. In the Cleve Creek watershed, the variations of the stream flows are quite similar under these two scenarios where the maximum streamflow occurs in June. Because the evaporation and the soil storage are greater in the high-temperature scenario, the maximum streamflow in the high-temperature scenario is smaller than those in the low-temperature scenario. In the Incline Creek and Twin River watersheds, higher temperature causes the maximum streamflow occur much earlier. The streamflow in these two scenarios gradually converges after the snow melts. Figure 11a-c plot the results of evaporation with vegetation under two different temperature scenarios. The evaporation in the high-temperature scenario is much greater than that in the low-temperature scenario before May but the difference between these two scenarios gradually decreases after May. On the other hand, there is a short period in the three watersheds when the evaporation in the low-temperature scenario is higher than that in the high-temperature. The main reason is that all the snow melts more quickly in the high-temperature scenario while a portion of the snow is still melting in the low-temperature scenario. We note that this period does not occur in the simulation without vegetation because the consumption of water from transpiration is much higher than the The comparison of streamflow in these two scenarios is plotted in Figure 10d-f. In the Cleve Creek watershed, the variations of the stream flows are quite similar under these two scenarios where the maximum streamflow occurs in June. Because the evaporation and the soil storage are greater in the high-temperature scenario, the maximum streamflow in the high-temperature scenario is smaller than those in the low-temperature scenario. In the Incline Creek and Twin River watersheds, higher temperature causes the maximum streamflow occur much earlier. The streamflow in these two scenarios gradually converges after the snow melts. Figure 11a-c plot the results of evaporation with vegetation under two different temperature scenarios. The evaporation in the high-temperature scenario is much greater than that in the low-temperature scenario before May but the difference between these two scenarios gradually decreases after May. On the other hand, there is a short period in the three watersheds when the evaporation in the low-temperature scenario is higher than that in the high-temperature. The main reason is that all the snow melts more quickly in the high-temperature scenario while a portion of the snow is still melting in the lowtemperature scenario. We note that this period does not occur in the simulation without vegetation because the consumption of water from transpiration is much higher than the evaporation and the limited water is consumed at a faster rate in the vegetated simulation. If the SOS comes earlier, this period would last longer. Figure 11d-f compares the transpiration in these two scenarios. The maximum transpiration occurs around July in the high-temperature scenario that is one month earlier than the normal temperature scenario. The transpiration is also much higher in the hightemperature scenario than that in the low-temperature scenario which shows the transpiration starts early and increases as the temperature increases. However, the difference of the transpiration between these two scenarios gradually decreases after July which means the EOS in various temperature does not significantly affect the transpiration.

The Impact of Varying SOS
Water 2021, 13, x FOR PEER REVIEW 17 of 23 evaporation and the limited water is consumed at a faster rate in the vegetated simulation. If the SOS comes earlier, this period would last longer. Figure 11d-f compares the transpiration in these two scenarios. The maximum transpiration occurs around July in the high-temperature scenario that is one month earlier than the normal temperature scenario. The transpiration is also much higher in the hightemperature scenario than that in the low-temperature scenario which shows the transpiration starts early and increases as the temperature increases. However, the difference of the transpiration between these two scenarios gradually decreases after July which means the EOS in various temperature does not significantly affect the transpiration.  12a-c shows the results of SWS under these two temperature conditions. The SWS has a large variation before June and in general is greater in the high-temperature scenario than in the low-temperature scenario. We note that between May and October the SWS in the low-temperature scenario is greater than that in the high-temperature scenario in the Incline Creek and Twin River watershed because snow melting can last longer in the low-temperature scenario. When the precipitation occurs in November, the soil storage increases in the high-temperature scenario. Figure 12d-f describes the streamflow under these two different temperature conditions. The streamflow in the low-temperature scenario is generally greater than that in the  12a-c shows the results of SWS under these two temperature conditions. The SWS has a large variation before June and in general is greater in the high-temperature scenario than in the low-temperature scenario. We note that between May and October the SWS in the low-temperature scenario is greater than that in the high-temperature scenario in the Incline Creek and Twin River watershed because snow melting can last longer in the low-temperature scenario. When the precipitation occurs in November, the soil storage increases in the high-temperature scenario.
Figure 12d-f describes the streamflow under these two different temperature conditions. The streamflow in the low-temperature scenario is generally greater than that in the high temperature scenario in the Cleve Creek watershed. In the Incline Creek and Twin River watersheds, the streamflow in the high-temperature scenario is greater at the beginning of the year while the streamflow in the low-temperature scenario is greater after May. Compare to the streamflow without vegetation in Figure 10d-f, the streamflow under both various temperature scenarios with vegetation is smaller due to transpiration. As the temperature increases, the streamflow has an obvious increase because more snow melts and rainfall dominate the precipitation. However, the transpiration consumes a portion of the water which leads to the smaller streamflow. As the SOS becomes earlier, the peak of the stream with vegetation would come earlier.
Water 2021, 13, x FOR PEER REVIEW 18 of 23 high temperature scenario in the Cleve Creek watershed. In the Incline Creek and Twin River watersheds, the streamflow in the high-temperature scenario is greater at the beginning of the year while the streamflow in the low-temperature scenario is greater after May.
Compare to the streamflow without vegetation in Figure 10d-f, the streamflow under both various temperature scenarios with vegetation is smaller due to transpiration. As the temperature increases, the streamflow has an obvious increase because more snow melts and rainfall dominate the precipitation. However, the transpiration consumes a portion of the water which leads to the smaller streamflow. As the SOS becomes earlier, the peak of the stream with vegetation would come earlier.

Drivers of the Vegetation Phenological Variations
The timing of SOS and EOS is controlled by different biological and environmental factors. In this study, we investigate several critical factors that have a great impact on vegetation phenology such as temperature, photoperiod and water availability. Temperature significantly affects vegetation phenology [50]. The temperature increase may lead to earlier spring and later autumn while temperature decrease may delay the timing of spring and advances of autumn [37]. A warming environment enhances the meristem cells

Drivers of the Vegetation Phenological Variations
The timing of SOS and EOS is controlled by different biological and environmental factors. In this study, we investigate several critical factors that have a great impact on vegetation phenology such as temperature, photoperiod and water availability. Temperature significantly affects vegetation phenology [50]. The temperature increase may lead to earlier spring and later autumn while temperature decrease may delay the timing of spring and advances of autumn [37]. A warming environment enhances the meristem cells to be stimulated and plant cell elongation to accelerate [51]. In this study, the annual temperature of the three watersheds approximately increases by 2 • C. Photoperiods also play a key role to control autumn phenological events such as leaf senescence [52]. The shortened photoperiod induces bud set and senescence when the photoperiod is smaller than the threshold of the growing-permitting period [53]. Water availability greatly affects plant Phenology in the vegetated watersheds [54]. Two impact factors, the timing of snow cover and snowmelt, notably alter plant phenology, particularly for shrubs and grasses at high altitudes and latitudes [55]. The varying snow effect also can alter the interaction between plant phenology, soil water content, and soil temperature. A portion of snow meltwater infiltrates into the frozen soil when the snow melts in early spring in a warming environment, which stimulates root activities even before the air temperature rises above zero and brings forward the timing of spring phenology [56].

The Effect of Vegetation Phenology on Evapotranspiration (ET) and Streamflow
We find that the ET values have increased greatly under the vegetation phenology scheme in the three vegetated watersheds. Transpiration plays the largest role in the total evapotranspiration during the growing season. This is consistent with the previous research that the percentage of transpiration, forest floor evaporation, and interception evaporation is 65%, 15%, and 20%, respectively [57]. The reason for the increase of ET is mainly attributed to the extended growing season lengths. Within the extended growing season, soil evaporation is in the opposite direction to transpiration and canopy evaporation due to the canopy interception which results in reducing soil water and soil evaporation.
Previous studies have demonstrated that climate warming is the main driver for streamflow decrease [58,59]. Consistent with these studies, both SOS and EOS are strongly linked with streamflow. Vegetation growth in the long-term plays a more vital role. When the SOS comes earlier and EOS terminates later, extra water consumed by vegetation may result in an aggravating decrease of soil water and groundwater. Therefore, water discharge in streams tends to decrease because precipitation would firstly supply soil deficit and recharge groundwater. The earlier SOS and delayed EOS accelerate the process of water transfer from land surface to the atmosphere dramatically because transpiration is the largest proportion of evapotranspiration [60,61].
We note that the vegetation phenology dynamics significantly affect the streamflow of these three watersheds during the last three decades and interact with long-term increase and interannual variability in growing season length. In this study, the streamflow decreases the dynamic ecosystem response to varying climate being integrated, which is consistent with the result that the increase of growing season lengths is attributed to daily minimum temperatures in spring and autumn [62]. This illustrates that the extended growing season increases vegetation water use and leaf interception in the watersheds. Vegetation first intercepts the precipitation and infiltrated water is mainly consumed via transpiration [63,64], which results in the decrease of water discharge to the stream.

Long-Term Influence of Early SOS
The extension of the growing season has been associated with recent global warming. Changes in the SOS not only have far-reaching consequences for plant and animal ecosystems but also may lead to long-term decreases in streamflow and changes in vegetation cover which affect the climate system. If the shifts in SOS observed in this study continue, they have critical implications for reservoir operations. Snowmelt occurs earlier, and the growing season may increase in length caused by earlier SOS, which may reduce the risk of flooding during snowmelt. On the other hand, flood risk might increase if the temperature persistently increases. Earlier SOS may cause soil moisture to decline during summer, increasing drought stress in trees, leading to high risk of wildfires and insect infestation.
Land surface roughness refers to the configuration or microrelief. Flow resistance determines the speed with which surface runoff responds. By affecting the time by which surface runoff remains on the land surface, the roughness also influences the volume of water that can infiltrate into the soil. In many parts of the northern hemisphere, studies indicated that the onset of growing season comes at least a week to 10 days earlier over the 20th century [65]. If this pattern of early SOS continues, as expected in response to continued climate warming, the effects on annual rates of transpiration could be significant. Increases in early spring transpiration would likely result in earlier declines in soil moisture content, declines in summer seasonal streamflow, lower summer base flow, and longer periods of summer low flow conditions. Watersheds are strongly correlated with total runoff because vegetated soil has greater roughness and higher capacity for infiltration [66]. The magnitude of the peak for the streamflow will become smaller and earlier, leading to flattened hydrographs.

Future Work
This study presents the drivers, mechanisms and the hydro-ecological response to the warming climate. However, several challenges need to be improved in future studies. The current model employs empirical parameters via a statistical relationship to describe plant growth. Lacking vegetational mechanisms may cause uncertainties when predicting phenology feedback in future climate. Inverse modeling methods can help to reduce the uncertainties by incorporating eco-physiological and morphological responses [56]. More climate manipulation experiments are needed to investigate the processes controlling plant phenological dynamics and to improve the CHESS model. Root production contributes to 33-67% of the terrestrial net primary production [67], and the response of root phenology to climate change is quite different from the aboveground phenology [68]. Maintaining species coexistence in multispecies plant communities is essential in plant phenology. Ongoing climate change alters the timing and length of phenological events, which may desynchronize seasonal interactions among species, resulting in altered biodiversity and ecosystem primary productivity [69]. Our model has the limitation that the species of the vegetation are constant in varying climate conditions. Filling such knowledge gaps can help us to better understand and predict plant phenological variations in response to the ongoing climate change.

Conclusions
The simulated results show that the updated CHESS model can accurately predict the variation of SWE in the three watersheds. But discrepancies between the simulated results and the observed data occur at the peaks due to the spatial and temporal variation in the physical characteristics of the snowpack. The static algorithm may overestimate or underestimate the parameters in the updated CHESS model based on the hypothesis that the snow grain size is uniform both temporally and spatially. By incorporating the GSI module, the updated CHESS model can be used to examine the timing of SOS and EOS. The simulated iGSI has been validated by NDVI from MODIS containing 16 days of data. The results indicate there is a high relationship between the temperature and the timing of the SOS and EOS. The SOS starts earlier and EOS ends later in a year when the watershed has a greater temperature increase. The streamflow in the static and dynamic scheme in the three watersheds are simulated between 1981 and 2017 with the NS of 0.52 and 0.80 in Cleve Creek, 0.46 and 0.75 in Incline Creek, and 0.42 and 0.70 in the Twin River watersheds, respectively.
Compared to the streamflow without vegetation, the streamflow under the two temperature scenarios with vegetation is smaller. As the temperature increases, the streamflow has an obvious increase with more snow melts and more rainfall-dominated precipitation. However, the transpiration leads to flattened streamflow hydrographs. As the SOS becomes earlier, the peak of the stream with vegetation occurs earlier. If the shifts in SOS observed in this study continue, the effects on annual rates of transpiration could be large and the growing season may increase in length, which may reduce the risk of flooding during snowmelt. On the other hand, flood risk may increase if temperature persistently increases.

Data Availability Statement:
The data displayed in this study is available in the current manuscript.