On the Sensitivity of Future Hydrology in the Colorado River to the Selection of the Precipitation Partitioning Method

The accurate partitioning of precipitation is important for snowpack and streamflow simulations in mountainous regions. However, most hydrologic models employ a partitioning method based on near‐surface air temperature (Ta) and the sensitivity of hydrologic simulations to the selection of the precipitation partitioning method (PPM) under climate change has been underexplored. To address this, we compared simulations using two PPMs, a conventional Ta‐based scheme and a wet‐bulb temperature (Tw)‐based approach (TA and TW) with the Variable Infiltration Capacity model. Our study focused on the Colorado River Basin (CRB) which is vulnerable to future snowfall and streamflow reductions due to climate warming. Historical simulations demonstrated the improved performance of the TW scheme in simulating snowfall fraction (SF) and snow water equivalent (SWE) in relation to in situ observations and a gridded SWE product. Subsequently, we evaluated the influence of the PPM selection on snowpack and streamflow projections from eight climate models under two emissions scenarios. The ensemble median results revealed more substantial decline in annual SF with the TW scheme as compared to the TA method (−12% vs. −9%) from the historical (1976–2005) to the far future (2066–2095) period, especially at lower elevations. Hydrologic simulations showed that the TA scheme underestimated reductions in SWE and streamflow (Q) under warming as compared to the TW scheme. This finding has important implications on future projections for the CRB and in other mountainous basins where climate warming can shift conditions from snowfall to rainfall.

The phase of precipitation, whether liquid, solid, or a combination, has a first-order control on snowpack dynamics and subsequent hydrologic processes (Harpold, Kaplan, et al., 2017;Vionnet et al., 2022).Precipitation partitioning directly influences the simulation of snow accumulation and ablation (Cho et al., 2022;Jennings & Molotch, 2019;Pan et al., 2003), thereby impacting the relative amounts of snow storage, sublimation, and meltinduced runoff (Harder & Pomeroy, 2014;Liu et al., 2022;López-Moreno et al., 2020;Mizukami et al., 2013;Sprenger et al., 2022).Moreover, the impact of precipitation partitioning will propagate into the dynamics of streamflow production, especially in snowmelt dominated regions.Given that hydrologic models are typically calibrated using streamflow records from a limited number of gauging stations, errors in the precipitation phase partitioning may lead to biases in the model parameterization.More importantly, this may affect the model sensitivity in simulating the impact of climate change, as the precipitation regime shifts towards a lower-to-no snow future.Despite its importance, hydrologic models typically divide precipitation into snowfall and rainfall based solely on near-surface air temperature (T a ).The simplest approach uses a threshold such that snowfall occurs when T a < 0°C (Motoyama, 1990).Other schemes involve a higher T amax (maximum T a at which snow forms), ranging from +0.5 to +2.5°C (Kim et al., 2021), but generally assume that the snowfall fraction (SF) is a function of T a and that T amax is close to 0°C.In practice, the threshold temperature can be calibrated to improve hydrologic model performance.For example, the North American Land Data Assimilation System uses a higher T amax (+1.5°C) to overcome the warm bias identified in the forcing product (Xia et al., 2018).To match hydrologic simulations to snow water equivalent (SWE) data, Sun et al. ( 2019) used higher T amax values which showed spatial variations in mountainous regions of the western U.S. (e.g., +3.2, +4.2, and +2.9°C in Southern Rockies, Wasatch and Uinta, and Arizona/New Mexico).
Recent studies have identified that accounting for near-surface atmospheric humidity can improve the accuracy of precipitation partitioning methods (PPMs, Jennings et al., 2018).The underlying physical mechanism for this improvement is as follows.When air is subsaturated, snow sublimation and the evaporation of coating droplets on snow provide cooling effects which increase the likelihood that snow remains in a solid state under air temperatures that are higher than 0°C (Heymsfield et al., 2021;Mitra et al., 1990).Differing humidity levels may partially explain the spatial variations of T amax noted in prior studies (Feiccabrino et al., 2015;Rajagopal & Harpold, 2016;Sun et al., 2019), suggesting that T a alone is insufficient to determine the precipitation phase.To address this, new formulations using the wet-bulb temperature (T w ) or dewpoint temperature (T d ), which account for humidity, have improved precipitation partitioning (Ding et al., 2014;Sims & Liu, 2015;Tamang et al., 2020).For example, Wang et al. (2019, hereafter W19) developed a PPM to estimate SF from T w over a full range ( 5 to +5°C in T w , and 0 to 100 in SF), which improved snowpack simulations relative to SWE observations.Despite these advancements, the influence of the selection of the PPM on streamflow projections under climate change has not been addressed.Importantly, integrating the humidity level into a PPM could lead to different sensitivities in response to climate warming.In addition, there is large uncertainty with respect to how relative humidity (RH) levels will change in arid regions such as the CRB (e.g., Byrne & O'Gorman, 2016;Harpold, Rajagopal, et al., 2017;Held & Soden, 2006;Simpson et al., 2024;Willett et al., 2007).In this study, we aim to address this gap and answer the following questions: (a) Does a T w -based PPM more accurately represent SF in mountainous regions?, and (b) What are the implications of the selection of the PPM on future projections of snowfall, snowpack dynamics, and streamflow?We selected the Colorado River Basin (CRB) for our study due to its arid and semiarid climate and mountainous terrain.As an important regional watershed, the CRB is facing historical drought conditions that have impacted water supplies (Bass et al., 2023;Xiao et al., 2018).Processbased modeling has identified the role of a warming-induced shift from snowfall to rainfall in future streamflow declines (Whitney, Vivoni, Bohn, et al., 2023), supporting prior efforts with the Budyko framework (Berghuijs, Woods, & Hrachowitz, 2014).We employed the Variable Infiltration Capacity (VIC) model to test two PPMs based on T a and T w and conducted point simulations at selected SNOTEL sites in the CRB and basinwide simulations using observed meteorological forcings and outputs from statistically-downscaled General Circulation Models (GCMs) over historical and future periods.

Precipitation Partitioning Method
We compared two PPMs that estimate SF.The first method partitions total precipitation (P) into snowfall (P S ) and rainfall (P R ) using near-surface T a .SF is defined as P S /P.The method uses two thresholds: T amin (minimum T a at which rain can fall) and T amax .We used values of 0.5°C (T amin ) and +0.5°C (T amax ), consistent with the default settings in VIC version 5.0: . (1) Earlier versions of VIC had only one T a threshold (usually set as 0°C) to partition snowfall and rainfall (e.g., Christensen & Lettenmaier, 2007;Nijssen et al., 2001).The second PPM, proposed by W19, uses T w to compute SF using a sigmoid function as: where the parameters a (6.99 × 10 5 ), b (2.0°K 1 ), and c (3.97°K) were optimized by W19 to fit the observed relation between SF and T w from global weather stations located at high latitudes (beyond 35°N and S), as described in Behrangi et al. (2018).We calculated T w as a function of T a and RH using Stull (2011), as shown in Appendix A.
Figure 1 compares the SF derived from the two PPMs.When SF only varies with T a (SF TA ), 50% of snowfall occurs at a constant T a = 0°C independent of RH (T a50 , Figure 1a).In comparison, the second PPM estimates a higher SF (SF TW ) when T a is above freezing and T a50 increases with lower RH (i.e., T amax is higher for drier atmospheric conditions).In mountainous areas with low humidity that are common in the CRB, T a50 is often higher than the freezing point by several degrees (e.g., T a50 = +5°C when RH = 50% in Figure 1b).Differences in SF between two schemes (SF TW SF TA ) can be divided into three regions (All Snow, Mixed, and All Rain, Figure 1c).SF TW equals SF TA in regions of All Snow (SF = 100%) and All Rain (SF = 0%), such that the PPM selection plays no role in SF and the subsequent snowpack dynamics.However, SF TW can be much higher than SF TA in the Mixed region, typically between T a = 0-6°C depending on RH, leading to the potential for large differences in snowpack dynamics and hydrologic processes.

Hydrologic Model Framework
Hydrologic simulations were carried out using VIC, a regional model that solves the water and energy balance on a regular grid, where each cell is divided into land cover tiles atop a three-layer soil column (Liang et al., 1994).
The snow module in VIC uses an energy balance approach to simulate the accumulation and ablation of snow and considers the effect of sub-grid topography by dividing grid cells into different elevation bands (Andreadis et al., 2009).As shown in Figure 2, each elevation band has a unique fractional area (A f ) and elevation (H).At each timestep, VIC calculates an adjusted T a using a lapse rate of 6.5°C/km and the mean elevation for each band.VIC also allows a user to specify the fraction of the total precipitation received by each band (set to be equal for simplicity in this study).Band-specific T a and P values are used to estimate P S and P R .VIC then solves the snow energy balance in each elevation band and calculates area-weighted average fluxes and storages.We ran VIC at 1/16th degree (∼6 km) resolution and an hourly timestep.We used VIC version 5.0 (Hamman et al., 2018) that incorporates a clumped vegetation scheme and its parameters from remote sensing products (Bohn & Vivoni, 2016, 2019;Xiao et al., 2022).Soil parameters were calibrated using monthly naturalized streamflow estimates at key gauges in the CRB (Whitney, Vivoni, Bohn, et al., 2023).In addition to T a and P, we used MetSim to generate other required daily meteorological forcings (wind speed, longwave and shortwave radiation, atmospheric pressure, and vapor pressure) and to disaggregate all forcings to an hourly timestep (Bennett et al., 2020;Bohn et al., 2013).We used the hybrid scheme of Bohn et al. ( 2019) to (a) disaggregate P uniformly during the 24 hr of a day when the daily minimum T a was below 0°C to depict snowfall (nearly all days with T w < 0°C have daily minimum T a < 0°C), and (b) apply a triangle scheme to obtain hourly P in cases with rainfall to represent a storm with a climatological parameters of monthly peak timing and duration derived from meteorological data.Streamflow is calculated as the sum of runoff and baseflow from contributing areas since routing has a small role in determining monthly and annual totals (Cortés-Salazar et al., 2023).

Study Basin and Its Hydroclimatic Regime
Hydrologic simulations were conducted in the CRB which has a total area of ∼630,000 km 2 and is divided into the Upper and Lower Basin (UCRB and LCRB, Figure 3).Despite the vast spatial extent, 85% of the annual streamflow of the CRB comes from 15% of the area located in mountainous regions (Christensen & Lettenmaier, 2007).Furthermore, most of these highly productive headwaters are in the UCRB, which produces more than 90% of the total streamflow at Imperial Dam (Lukas & Payton, 2020).It is noteworthy that there are several mountainous areas in the LCRB that are also important in the region (e.g., Gila).In the CRB, snowmelt from mountain headwaters is the primary source of streamflow, where the streamflow efficiency is highest.Li et al. (2017) estimated that snowmelt contributes around 70% of streamflow, with a higher fraction (>80%) in high elevation subbasins that have more snowpack (Figure 3c).We used the delineation from Whitney, Vivoni, Bohn, et al. (2023) to obtain eight analysis subbasins for which the spatial average and standard deviations of elevation and mean annual values of T a , RH, and T w are shown in Table 1.Whitney, Vivoni, Bohn, et al. (2023) reported on the subbasin selection process, stream gauging locations, and the naturalized streamflow records (Bureau of Reclamation, 2020) used for model calibration and testing.Overall, the UCRB has lower T a and higher RH, in particular in the Green and Upper Colorado subbasins that have higher elevation, as compared to the LCRB.

Snow Data Sets and Meteorological Forcings
We obtained daily meteorological data from the Global Historical Climatology Network (GHCN) to estimate SF.Among the GHCN stations in the CRB, we selected a subset of stations (77 in total) with at least 10 years of highquality data (missing <15 days within one year) from 2001 to 2018 (Figure S1 in Supporting Information S1).We followed the method described in Knowles et al. (2006) to derive SF, which categorizes any recorded precipitation (P) as solid if accompanied by a positive snow depth measurement for the calculation of snowfall water  equivalent (SFE).SF is then calculated as SFE/P.Despite uncertainties on days with a mixture of snow and rain, this method avoids overreliance on the reported snowfall amounts, which can be unreliable due to observer bias (Knowles et al., 2006).The SF derived from the GHCN records were used as a benchmark to evaluate the two PPMs.We also obtained daily meteorological measurements and SWE from SNOTEL sites that underwent qualitycontrol procedures (Sun et al., 2019;Yan et al., 2018).These procedures eliminated erroneous values and corrected known measurement issues, including snowfall undercatch (Livneh et al., 2014) and a warm bias at cold temperatures (Oyler et al., 2015).We selected a subset of 124 SNOTEL stations in the CRB with high quality data (<15 days of missing data per year for at least 20 years from 1971 to 2018, Figure 3b).In addition, a daily gridded SWE product at 4 km resolution from 1981 to 2021 was obtained from Broxton et al. (2019, hereafter UA product).The UA product was generated by assimilating measurements of SWE and snow depth data with gridded PRISM (Parameter Regression on Independent Slopes Model) precipitation and temperature data sets (Daly et al., 1994).Details of the methodology and its evaluation have been reported in Zeng et al. (2018).The product has also been used previously for hydrologic model assessments (e.g., Hao et al., 2023;W19).
For model simulations, two sets of daily gridded (6-km) meteorological forcings were processed: (a) an interpolated product from Livneh et al. (2015) for historical evaluations against the SNOTEL observations, UA product, and naturalized streamflow data, and (b) statistically-downscaled products of Pierce et al. (2014) from 8 models under Representative Concentration Pathway 4.5 and 8.5 emissions scenarios (hereafter R45 and R85) for climate change experiments.The GCMs (Table S1 in Supporting Information S1) were selected based on their performance in reproducing historical climate (P and T a ) in the CRB (Gautam & Mascaro, 2018).Figure 4 illustrates T a and RH conditions of each GCM during different seasons (Fall: October-December; Winter: January-March; Spring: April-June; Summer: July-September).These are shown as seasonal averages for different elevation ranges to illustrate the effect of mountainous terrains.When averaged across all GCMs and elevations, the far future period (2066-2095) has a higher T a (+3.2 and +5.2°C) but similar RH (+0.04 and 0.66% under R45 and R85, respectively) as compared to the historical period .Individual GCMs show large variations in far future scenarios as depicted by the spread among the symbols (cf., GCM variability analysis in Whitney, Vivoni, Bohn, et al., 2023).GCM forcings are also compared to the three regions of the RH-T a space (All Snow, Mixed, and All Rain).GCMs whose conditions remain in regions All Snow and All Rain under warming are insensitive to the PPM selection.For instance, all GCMs at all elevations are in All Rain region during Summer (Figure 4d).GCMs with conditions in Mixed region, however, are likely to be sensitive to the PPM, for instance, at low elevations (1,000-2,000 m) in the Fall and Winter (Figures 4a and 4b) such that the T abased PPM underestimates SF.Moreover, the PPM selection is critical when warming leads to conditions that cross different regions (e.g., shift from All Snow to Mixed region in Winter at 2,000-3,000 m, Figure 4b).

Numerical Experiments
We conducted climate change experiments with VIC using the statistically-downscaled GCM forcings from 1950 to 2095 under the two emissions scenarios (R45 and R85).We then selected two 30-year periods for analysis: historical  and far future (2066-2095) periods.We carried out three numerical experiments with the two PPMs: (a) TA refers to using the T a -based PPM, (b) TW refers to using the T w -based PPM, and (c) TWC is similar to TW but involves recalibrating soil parameters such that the streamflow from TW matched that obtained from TA in the historical period.The purpose of TWC was to obtain a similar streamflow performance to the calibrated model of Whitney, Vivoni, Bohn, et al. (2023) who employed the T a -based PPM.In the following analysis, we used two notation systems to represent different experiments (TA, TW, TWC) and time periods (historical, far future).Using SF as an example, differences between PPMs for the same forcing are labeled as ΔSF and calculated as: (SF TW SF TA ).In contrast, differences between historical and far future periods for the same PPM are referred to as dSF, calculated as: (far future-historical).Following Whitney, Vivoni, Bohn, et al. (2023), we carried out spatiotemporal analyses for each GCM and then obtained ensemble-median statistics.

Snowfall Under Historical and Far Future Periods
We first compared SF from the GHCN stations to the T a and T w -based VIC simulations using observed historical forcings from Livneh et al. (2015) over the period 2001-2018.GHCN stations at elevations from 1,079 to 2,928 m reported an annual-averaged SF of 31.5%.In contrast, the annual-averaged SF at co-located VIC cells for the T abased simulation was underestimated at 22.9%.The T w -based method significantly reduced this bias, resulting in an annual-averaged SF of 33.7%.To illustrate this, Figure 5 shows the relation between SF and elevation from GHCN observations (OBS) and simulations with the two PPMs, with the inset in Figure 5a showing the cumulative distribution functions of SF.Using the T w -based scheme improved the accuracy of annual SF, which was systematically underestimated in the T a -based PPM across most elevations.When averaged for two major elevation bands with GHCN stations (Figure 5b), the T w -based PPM also improves annual-averaged SF relative to the observations.
Next, we inspected the differences in SF between the two PPMs for the historical and far future periods forced by the GCMs (Figure 6, Tables S2 and S3 in Supporting Information S1).As expected from the above, SF in the TW simulation is larger than that derived from TA.This leads to mean annual differences in SF (ΔSF) from +9.0 to +10.0% across GCMs (median of +9.6%) in the historical period (Table S2 in Supporting Information S1).Low variations in ΔSF among GCMs indicate ensemble median values are robust.The UCRB had larger differences and lower spatial variations in ΔSF (12.5% ± 2.0%) as compared to the LCRB (7.3% ± 5.1%).While this comparison suggests that higher elevations have larger differences in ΔSF, the trend does not consistently hold for subbasins with elevations ≥2,000 m.For instance, the Green subbasin, despite having a lower average H (2,215 m) than the Upper Colorado subbasin (2,542 m), exhibits a slightly larger ΔSF (13.4% vs. 12.7%).To explain this pattern, we evaluated relations between ΔSF and elevation (Figure S2a in Supporting Information S1).In the historical period, ΔSF increases with H up to ∼2,000 m, after which a slight decrease is noted.In comparison, locations at higher and lower elevations typically fall into All Snow and All Rain regions, respectively, where SF TW is close to SF TA .The trend primarily reflects an air temperature effect, as the selection of the PPM is most sensitive when T a is close to the freezing point (i.e., Mixed region in the RH-T a space).Relations between ΔSF and T a in Winter and Spring seasons (Figure S2b in Supporting Information S1) further confirmed this pattern as ΔSF peaks around 0°C.In the far future period, ΔSF decreases across all elevations due to: (a) a significant shift from Mixed region to All Rain at low to mid elevations (H < 2,000 m) such that sensitivity to the PPM is greatly reduced (e.g., region All Rain has SF TW = SF TA ), and (b) a smaller shift from region All Snow to All Rain at high elevations due to warming which leads to an asymptotic value of ΔSF of 11%.The spatial pattern of ΔSF also exhibits changes in the far future period relative to the historical conditions, including a sharper decline noted from north to south (Table S3 in Supporting Information S1).For example, under the R85 scenario, the Green and Upper Colorado subbasins show large differences in snow fraction (ΔSF = ∼12%).In contrast, the Gila and Lower Colorado subbasins exhibit a smaller ΔSF (1.7% and 0.7%, respectively) as compared to the historical period (4.8% and 3.6%).
Seasonal variations in ΔSF are shown in Figure 6 for the historical and far future periods.Overall, higher values of ΔSF are noted in Winter, followed by Fall and Spring during the historical period (Table S3 in Supporting Information S1).In Summer, most areas exhibit low ΔSF, as T a > 0°C, as expected.Spatial variations of ΔSF among the seasons become less pronounced in the far future period due to an overall decrease in ΔSF, especially in the LCRB.For instance, the far future period (R85) in the CRB has ΔSF values in Winter, Fall, and Spring that are below the historical period (Table S3 in Supporting Information S1).Clearly, the PPM selection leads to differences in SF that vary spatially, and which evolve under warming, especially for sites and times crossing regions in the RH-T a space.To pursue this further, Figure 7 shows monthly variations of ΔSF as a function of elevation in relation to similar depictions of T a and RH.In the distribution of T a and RH (Figures 7a-7d), conditions corresponding to Mixed region exist.Those conditions are clearly shown where ΔSF > 0, which includes a large fraction of the CRB elevations above 1,000 m (Figure 7f).Altitudes with the most pronounced ΔSF vary progressively during the year.In the historical period, lower elevation sites (1,000-2,000 m) have the largest ΔSF in the Winter, whereas large values of ΔSF shift toward higher elevations in the Spring.In the far future period, ΔSF has an overall lower magnitude, with a near disappearance in the Summer (ΔSF = 0).Interestingly, increases in T a in the far future period led to a shift of high ΔSF values toward elevations above 2,000 m and a reduction in ΔSF below 2,000 m during most of the year.Elevation controls on ΔSF indicate that important sensitivities on the selection of the PPM are possible under warming.
As detailed in Table S4 in Supporting Information S1, the TW simulations showed a greater decline in SF (dSF) across time (far future-historical) than TA in all subbasins in absolute terms, with dSF TW ranging from 11.9% to 7.4%, as compared to 8.5% to 5.2% for dSF TA for the CRB in R45 and R85.We also included the relative difference (Diff.column in Table S4 in Supporting Information S1) to confirm this finding.This indicates that the T w -based scheme yields more declines in snowfall under climate warming.Furthermore, there are important differences among subbasins that are linked to their position in the RH-T a space.For instance, higher elevation subbasins such as the Green and Upper Colorado, where warming impacts on SF are high (dSF < 7.0% across all cases), have a relatively small sensitivity to the selection of the PPM (dSF TW dSF TA from 0.6% to 1.2%) as compared to lower elevation subbasins.Overall, the UCRB shows a greater decline in dSF than LCRB for TA ( 11.8% vs. 5.7%, R85), yet the difference among the schemes in UCRB is lower ( 2.4% vs. 4.2%, R85).This suggests that lower elevations have a smaller decline in SF under warming, but the sensitivity to the PPM is more pronounced.Given the influence of the selection of the PPM on SF, an important sensitivity is expected for snowpack dynamics, as explored next.

Snowpack Dynamics Under Historical and Far Future Periods
We compared SWE simulations using the two PPMs to the SNOTEL observations and UA product over the period 2001 to 2018 using Livneh et al. (2015) as forcing (Figure S3, Tables S5 and S6 in Supporting Information S1).
The T a -based approach yielded an overall mean annual bias of 12.0 mm, correlation coefficient of 0.85, and root mean squared error (RMSE) of 31.7 mm with respect to all SNOTEL sites.A persistent negative bias in SWE is observed below 3,500 m, ranging from 39.9 to 4.1 mm for SNOTEL and 25.4 to 0.1 mm for the UA product.SWE underestimation was reduced when applying the T w PPM, which had a low mean annual bias of 2.6 and 0.8 mm for SNOTEL and UA, respectively.Spatial maps of SWE differences between the simulations and the UA product show an overall bias improvement by the T w PPM, but also indicate regions of overestimation.For instance, for elevations from 2,500 to 3,000 m, the T w -based PPM reduced the bias and RMSE as compared to T a , but the comparison worsened at elevations above 3,500 m.While SWE simulations could be further improved through parameter calibration (e.g., Houle et al., 2017;Sun et al., 2019), our focus here was on assessing the influence of the PPM selection.
Figure 8 presents the spatiotemporal variations of mean annual SWE differences (ΔSWE) between the TA and TW schemes for the historical and far future periods.Annual values of SWE were larger in the TW simulation (ΔSWE = +2.5 mm in historical period) as compared to TA, with the UCRB exhibiting a higher ΔSWE (+4.6 mm) as compared to the LCRB (+0.5 mm), as shown in Table S7 in Supporting Information S1.The impacts of the PPM were most evident in northern subbasins where there is more snowpack (ΔSWE of +4.7 and +7.9 mm in Green and Upper Colorado), whereas ΔSWE was lower in other subbasins, particularly in the LCRB.In the far future period, annual values of ΔSWE over the CRB narrowed as compared to the historical period (+1.9 and +1.6 mm for R45 and R85, Table S7 in Supporting Information S1) with similar responses for all subbasins and the two emissions scenarios.In contrast to annual values, the Winter season ΔSWE showed increases in the far future period for the Green and Upper Colorado subbasins, indicating heightened sensitivity to the PPM in line with ΔSF.In Spring and Fall seasons, ΔSWE in the far future period was reduced throughout the CRB, in particular for the LCRB, which led to the overall annual trends noted in ΔSWE.Overall, warming reduced the impact of the PPM selection on the snowpack dynamics in the far future period, but yielded additional sensitivity in the Winter season in those subbasins whose conditions crossed regions of the RH-T a space (All Snow to Mixed).

Water Resources Research
10.1029/2023WR035801 WANG ET AL.

Hydrologic Response Under Historical and Far Future Periods
Next, we evaluated the impact of the PPM selection on the hydrologic responses in the CRB (Figures S4 and S5, Tables S8 and S9 in Supporting Information S1) using annual averaged quantities.Overall, the TW simulation produced higher snowfall (P S ) and snowmelt (M), especially in the high elevation regions (ΔP S = +48 mm, ΔM = +42 mm in UCRB).Although ΔSF peaks ∼2,000 m, ΔP S continues to increase with elevation, driven by the consistently higher precipitation.The spatial pattern of ΔM is similar to ΔP S (Figure S4 in Supporting Information S1), except in regions below 2,000 m where ΔSWE is lower despite increases in P S (i.e., snowfall melts before enhancing the snowpack).Increases in M augment SM, although a spatial contrast was noted in how the additional M from TW affects hydrologic partitioning.In the LCRB, where elevations are predominantly below 2,500 m, the additional M is not sufficient to increase SM to soil field capacity, leading to almost no changes in total runoff and baseflow (RBF, +2 mm) and ET ( 1 mm).Components of ET exhibited different responses between TA and TW experiments.As expected, sublimation (E s ) increased due to higher SWE in the TW simulation (ΔE s = +2 mm).Plant transpiration (Tv) also increased due to higher SM.Other ET components decreased.
In contrast, the PPM scheme had a greater impact at higher elevations above 2,000 m, corresponding to the approximate freezing temperature.This elevation marks an increase in ΔSWE as well as a shift in the hydrologic partitioning between ET and RBF ( 9 mm and +10 mm for TW in UCRB, respectively).The higher snowpack increases the amount of water storage and creates a larger snowmelt pulse in the warm season that saturates the soil column, leading to a greater hydrologic partitioning to streamflow for the TW simulation.Spatial variations in SM further support the analysis, as ΔSM peaks at ∼2,500 m and stays relative stable at higher elevations.This finding is consistent with Hale et al. (2023), who identified that regions with prolonged and more extensive snow water storage tend to have more streamflow production in western U.S. Differences among the PPMs propagate in a similar fashion in the far future period (R85, Table S9 in Supporting Information S1), although reductions in P S and SWE due to warming lead to smaller differences among the schemes that affected all hydrologic variables.Clearly, there is a cascading set of impacts from the selection of the PPM influencing snowpack, SM, and streamflow dynamics, which depend on the specific location in the CRB with different sensitivities.To further explore this, Figure 9 presents the spatiotemporal variations in the differences between TW and TA simulations (ΔSWE, ΔSM, and ΔRBF) for historical and far future (R85) periods.In comparison to ΔSF, the distribution of ΔSWE values shifts toward higher elevations (H > 2,000 m) and into the Spring season when snowmelt starts.This shift in space and time highlights the importance of snowpack dynamics in modifying the effects of the PPM selection, with additional dampening occurring in the far future period.Differences in ΔSM and ΔSWE patterns indicate the role of snowmelt in increasing SM.Higher P S and M in the TW simulation led to widespread positive ΔSM above 2,000 m, except for H > 3,000 m in Spring, corresponding to the highest ΔSWE.ΔSM propagates to ΔRBF which showed positive values in most locations and times, except for high elevations (H > 3,000 m) in the Spring.This pattern can be possibly explained as Spring snowmelt occurs relatively earlier in the TA simulation, leading to greater streamflow production than in TW (Barnhart et al., 2016).Later in the Spring and Summer seasons, ΔRBF is mostly positive in regions where ΔSM values are positive.In the far future period, high values of ΔSWE are slightly reduced and shift earlier toward the Winter season.As a result, lower differences in SM occur earlier (around February) and over broader regions (H > 2,500 m).ΔRBF exhibits a 1-month advance and its magnitude becomes smaller as compared to the historical period.In essence, warming (lower SF) and changes in PPM (higher SF in TW simulation) exhibits a complex interplay in determining the hydrologic partitioning and the dynamics in hydrologic variables are heavily influenced by the snowpack condition across different regions in the CRB.

Streamflow Generation Under Historical and Far Future Periods
We next assessed the impact of the selection of the PPM on monthly streamflow (Q) at the CRB outlet for the historical and far future periods (Figure 10, Table S10 in Supporting Information S1).As noted previously, the TW simulation estimated higher SWE than TA in the historical period, which matched more closely with the reference UA product.Meanwhile, the soil parameters in VIC were calibrated to be more efficient in generating streamflow to compensate the underestimation of SF in TA.As a result, Q TW is higher than naturalized streamflow.The propagation of errors in the precipitation phase into streamflow simulations was also noted in Harpold, Kaplan, et al. (2017).To correct for this, we recalibrated soil parameters in the TWC simulation to match the naturalized streamflow data over 1981 to 2016 at six key locations (Green, Upper Colorado, San Juan, UCRB, LCRB, and CRB, Figure S6 in Supporting Information S1).We followed previous model calibration efforts by increasing the second layer soil depth from 0.7 to 1.3 m throughout the CRB.Although Q TWC does not perfectly match Q TA for all months, the annual streamflow agrees between TWC and TA.
During the historical period, water year Q TW was 17% larger (+5.3 mm) than Q TA , whereas Q TWC was similar to Q TA (+0.6 mm, Table S10 in Supporting Information S1).As expected, all simulations showed a streamflow decline in the far future periods.For instance, TA exhibited a decrease of Q (dQ = 22.5% or 7 mm) for R85, while TW had a more pronounced decline ( 26% or 10 mm), due to more significant losses in SF and SWE.By calibrating the streamflow efficiency, TWC had an even larger relative decrease in water year Q (dQ = 30% or 10 mm for R85).As a result, ΔQ values, which were positive for TW in all warming conditions, decreased for TWC in all months becoming negative in Winter and Spring.Interestingly, the sensitivity of Q varied during the year in response to the diminishing and shifting SWE in the far future periods.Most changes occur in the Spring and Summer, where a 1-month earlier shift was obtained in the peak streamflow.We also note that streamflow efficiency played a critical role in the sensitivity to the PPM selection.As a result, shifts in ΔQ for TWC result from both a lower amount of streamflow generation due to soil parameter calibration and a lower sensitivity to the PPM selection under warming.We further analyzed the impact of the selection of the PPM on SWE and Q for the UCRB (Figure S7 in Supporting Information S1), and found similar patterns of declines in SWE and Q from historical to far future periods.

Tracking the Hydrologic Impacts of the PPM and Model Calibration
Precipitation phase plays an important role in controlling the hydrological processes in watersheds with seasonal snow coverage.Alterations to the PPM cascade through the water and energy balances and affect the hydrologic response to climate change.This is especially pertinent in regions where the PPM selection can lead to changes in the watershed sensitivity (e.g., Mixed region in the RH-T a space).Figure 11 summarizes the sensitivity of the watershed processes to the selection of the PPM when averaged over the CRB under historical and far future periods.We included the TWC simulation that captured the effect of the PPM but generally increased the soil water storage and reduced the streamflow generation relative to TW.An improved representation of snowfall in the T w -based scheme led to higher amounts of P S (+27 to +34 mm/yr), a larger snowpack (ΔSWE = +2 mm), and increases in sublimation (ΔE s = +2 to +3 mm/yr) and snowmelt (ΔM = +32 to +40 mm/yr).As a result, an enhanced SM is simulated for the TWC scheme.ET is similar between two experiments in the historical period, due to reductions in soil evaporation (E soil , 5 mm/yr) and canopy evaporation (E c , 4 mm/yr) balanced with increases in plant transpiration (T v , +5 mm/yr) and sublimation.The increases in baseflow (BF, +2 mm/yr) were also balanced by the decrease in runoff (R, 1 mm/yr), to similar streamflow production (RBF) after the calibration of soil parameters.In far future periods, RBF is lower in the TWC scheme ( 2 to 3 mm/yr) due to more ET (+2 to +3 mm/yr).Differences in hydrologic variables between two PPMs remain relatively consistent in the far future periods, albeit with a reduction of differences in SF driven by the transition to warmer conditions with less snow.These trends are attributed to the stationarity of conditions at high elevations (H > 3,000 m) in the All Snow region during the Winter where differences between the PPMs are negligible.
Our analysis also reveals large spatial and seasonal variabilities in the hydrologic impacts of using the T w -based PPM in the CRB.In regions where temperatures are close to freezing point during the cold season, the PPM selection can lead to diverse snowpack responses.These variations influence the hydrologic partitioning between ET and Q, as the amount of snowpack storage affects the temporal alignment between total surface water inputs (M and P R ) and the atmospheric water demand in the warm season (Berghuijs, Sivapalan, et al., 2014;Foester et al., 2016;Hale et al., 2023).Moreover, beyond the sub-basin variability revealed in our regional study, hydrologic processes should also vary considerably at the hillslope scale.Transitions in energy or water limitations can occur in adjacent slopes within a single VIC grid cell.Such sub-grid variability affects the local and overall hydrologic response and underscores the need for hyper-resolution modeling.Future efforts at these scales could reveal the effects of climate-terrain interactions on hydrological processes in mountainous regions and improve our mechanistic understanding of PPM impacts (e.g., Ko et al., 2019;Mazzotti et al., 2021).

Sensitivity to the Selection of the T a Threshold
In this study, we used a T a threshold of 0°C as the maximum air temperature of snowfall due to its frequent application in VIC and other hydrologic models (e.g., Christensen & Lettenmaier, 2007;Hale et al., 2023;Nijssen et al., 2001) as well as its common use as a benchmark for analyses of PPM schemes (e.g., Harder & Pomeroy, 2014;Harpold, Kaplan, et al., 2017;Leroux et al., 2023;Vionnet et al., 2022).The value of 0°C, however, should be considered as a lower threshold since T amax values are likely above freezing in the CRB, especially at higher elevations (Jennings et al., 2018).Using a more realistic T amax is expected to reduce differences in the simulated hydrologic responses between the T a and T w -based PPM.For instance, Jennings et al. (2018) provided maps of threshold values of air temperature at which 50% of the precipitation is snow (T a50 ).Within CRB, T a50 from this product varies from 2.0 to 3.0°C.Similarly, Harpold, Rajagopal, et al. (2017) used the method of Dai (2008) to quantify SF in the western United States, finding a ∼20% larger value on average as compared to a PPM based on a T a threshold of 0°C.Harpold, Rajagopal, et al. (2017) also identified that the method of Dai (2008) yielded similar SF values to a humidity-based PPM over 1980-2015 across many watersheds (see Figure 2 in Harpold, Rajagopal, et al., 2017).
Given the limitations in assuming a T a threshold of 0°C for 50% SF, we conducted a VIC model sensitivity analysis that applied a different temperature threshold within the PPM scheme (Jordan, 1991).In this approach, an SF of 60% is obtained over the interval of 2.0°C < T a ≤ 2.5°C.We found that these simulations over the CRB yielded a higher SF (22.7%) as compared to the 18.7% derived from 0°C in the T a -based PPM and closer to the SF obtained for the T w -based PPM (28.3%) during the historical period .In addition, the approach of Jordan (1991) led to SF changes from the historical to the far future periods (dSF = 6.0% and 9.8% for R45 and R85, respectively) that were intermediate to the TA and TW simulations (see Table S4 in Supporting Information S1), but closer to the values of TA (dSF = 5.2% and 8.5% for R45 and R85).This comparison suggests that increasing the T a threshold alone from 0 to 2.5°C does not fully match the outcomes of incorporating the humidity effects on the PPM through the wet-bulb temperature.As a result, the differences between the TA and TW simulations cannot be solely attributed to the assumption of a T a threshold of 0°C, but the use of a different threshold, such as the Jordan (1991) scheme, would reduce the disparities between the T a and T w -based PPMs.

Impacts of PPM Selection to Physically-Based and Surrogate Hydrologic Simulations
The use of PPMs varies considerably among the meteorological and hydrologic communities at least in part due to challenges in measuring solid precipitation (Rasmussen et al., 2012) and the absence of data sets with explicit rain and snowfall components (Sun et al., 2019).While reanalysis products often contain the components, these are typically not archived independently for use in hydrologic models.As shown here, the selection of the PPM used in physically-based hydrologic models can cause significant discrepancies in simulating snowfall, snowpack, and streamflow sensitivities to warming for certain regions of the RH-T a space.Given that mountainous areas may experience a low-to-no snow future (Siirila-Woodburn et al., 2021), it is important to adopt more realistic PPMs, such as T w -based approaches, that can more properly account for the sensitivity to warming in hydrologic modeling studies, including in snow drought characterization (Hale et al., 2023), snowmelt-originated runoff estimation (Li et al., 2017), and rain-on-snow flood risk assessments (Li et al., 2019), among others.
In recent years, machine learning (ML) methods have been used to derive surrogates of physically-based models in order to reach higher computational efficiency (e.g., Ivanov et al., 2021).Current ML models often face challenges in predicting streamflow changes under warming as future climatic conditions are an extrapolation beyond the training data set (Wi & Steinschneider, 2022).As a result, improvements related to precipitation partitioning, such as explicitly incorporating rainfall and snowfall components in the training data sets, could possibly enhance surrogate models for predicting climate change impacts, particularly for snow-dominated basins.For instance, Xu et al. (2022) used snowmelt and rainfall to train an ML model, which showed good accuracy in simulating streamflow.Zhong et al. (2023) also highlighted the effectiveness of a physics-informed ML model in capturing runoff changes under warming.These examples underscore the potential for bridging ML and process-based modeling to improve our understanding of hydroclimatic change (Razavi, 2021).

Implications to Water Management in the Colorado River Basin
Water management in the CRB relies on hydrologic predictions over short time scales (1-24 months) and on scenario modeling of future conditions in the 21st century (Bureau of Reclamation, 2012;National Research Councils, 2007).Prior efforts to incorporate climate change have typically assumed a T a -based partitioning method (e.g., Cayan et al., 2010;Christensen & Lettenmaier, 2007;Whitney, Vivoni, Bohn, et al., 2023) that likely has underestimated the sensitivity of streamflow to future warming.To illustrate this, Figure 12 presents the streamflow changes (Q Diff , far future-historical) and relative differences (dQ Diff ) with respect to TA for eight subbasins in million acre-feet (MAF), a common unit used in CRB water management (1 MAF = 1.23 km 3 ).The TA simulation shows a total reduction of 1.4 and 3.5 MAF annually under R45 and R85 for the CRB, while TW exhibits larger declines ( 2.1 and 4.6 MAF) indicative of the sensitivity to the PPM selection.The soil parameter calibration in TWC had a small overall effect on Q Diff ( 2.4 and 4.7 MAF) and further reduced dQ Diff by an additional 5%-10% relative to the TA simulation.Most streamflow reductions occur in the UCRB, where the Upper Colorado is the most sensitive subbasin, whereas the smaller overall reductions in the LCRB are driven primarily by the Gila subbasin.Note that both the Gila and Green subbasins have a significant sensitivity to the PPM ( 0.2 MAF difference between TA and TWC), which are a large fraction of the total declines (Q Diff of 0.6 MAF for TA under R85).

Conclusions
In this study, we evaluated the sensitivity of the simulated watershed processes under climate change in the CRB to the selection of the PPM, including snowfall, snowpack dynamics, and streamflow generation.Two approaches based on air and wet-bulb temperatures were compared within the VIC model.Due to its arid and semiarid mountainous terrain, a wide range of air temperature and humidity levels are present in the CRB, which are critical for determining the precipitation phase.As a result, this large watershed is representative of conditions experienced in the western U.S. and in other midlatitude continental regions where snowfall can occur above the freezing point, especially when low humidity levels are present.We found that analyses of the sensitivity to the PPM were aided by inspecting regions in the two-dimensional space defined by RH and air temperature.Results from this study indicate: 1.The PPM based on the wet-bulb temperature (T w ) yielded higher estimates of SF that better matched weather station observations as compared to the more commonly used method based on near-surface air temperature (T a ).Higher SF for the T w -based PPM increased SWE and influenced other watershed processes, including a 17% increase in annual streamflow in the historical period, which required recalibration of soil parameters in the model.2. Both methods yielded significant declines in SWE and streamflow by the end of 21st century for the ensemble median conditions for the two emissions scenarios.The T a -based PPM, however, underestimated future declines in streamflow for two reasons: (a) underestimation of the shift from snowfall to rainfall, and (b) overestimation of streamflow efficiency due to the original soil parameter calibration.As a result, prior T a - based simulations underestimate future streamflow declines in the CRB (e.g., Whitney, Vivoni, Bohn, et al., 2023,Whitney, Vivoni, Wang, et al., 2023).3. The selection of the PPM has hydrologic impacts that vary spatially throughout the CRB and evolve under future warming.In higher elevation regions, the impact of the PPM selection is more substantial as measured by the amount of snowfall and snowpack, but the lower Winter T a in these areas provides a buffer against warming.In contrast, lower elevations in the LCRB, with more ephemeral snowpacks, exhibited heightened reductions in SWE and streamflow in a low-to-no snow future.
In summary, our findings highlight the critical role of the selection and proper use of the PPM when simulating watershed processes in the CRB and evaluating its susceptibility to warming conditions.As warming shifts precipitation to less snowfall and more rainfall, the latter will be poised to become a more significant determinant of overall streamflow generation (Ban et al., 2023) and can interact with other changes, such as forest disturbances (Whitney, Vivoni, Wang, et al., 2023), to impact the overall streamflow production of the basin.Importantly, watersheds that currently have ephemeral snowpacks will be more sensitive to the warming effect and could lead to conditions that transform the watershed processes.As an example, prior efforts in the CRB using VIC and the air temperature threshold of 0°C likely underestimated the effects of warming on the streamflow response.Future studies in the CRB should either employ a PPM based on the wet-bulb temperature or apply the air temperature scheme with a spatially variable threshold derived from station data (e.g., Jennings et al., 2018).It is also important that studies include a description of the selected PPM scheme as this will allow comparisons across different approaches.In addition, the use of the two-dimensional framework of RH and air temperature presented in this study can be applied to screen for those sites and time periods where the selection of the PPM is critical.This can be extended, for example, to other PPM schemes and to outputs from more recent downscaled climate change projections (e.g., Pierce et al., 2023).Furthermore, this framework provides an effective tool for communicating these important findings to water managers and to support a more comprehensive understanding of long-term hydrologic changes that are possible in the CRB.

Figure 1 .
Figure 1.Snowfall fraction (SF) as a function of T a and relative humidity.(a) SF based on T a .(b) SF based on T w with inset showing relation between SF and T w .(c) Differences in SF (SF TW SF TA ) with three regions (All Snow, Mixed, and All Rain).White dashed lines represent 50% SF (or T a50 ).

Figure 2 .
Figure 2. Schematic of different elevation bands (B 1 , B 2 , B 3 , …) with areas A f and heights H in each grid cell and its effects on T a , the partitioning of P, and snow water equivalent.

Figure 3 .
Figure 3. (a) Location of the Colorado River Basin with a simplified stream network.(b) Elevation (H) distribution from USGS (2016) with SNOTEL sites.(c) Mean annual snow water equivalent (mm) from a 4-km gridded product (Broxton et al., 2019) from 2004 to 2018 overlayed with subbasins.

Figure 4 .
Figure 4. T a and relative humidity from eight General Circulation Models (shown with different markers, full names in TableS1in Supporting Information S1) averaged over different elevations (1,000-2,000, 2,000-3,000, and 3,000-4,000 m, labeled in (c)) for (a) Fall, (b) Winter, (c) Spring, and (d) Summer seasons in relation to the precipitation partitioning method regions (All Snow, Mixed, and All Rain, Figure1).Shaded area denotes where SF TW > SF TA .

Figure 5 .
Figure 5.Comparison of co-located average annual snowfall fraction over period 2001-2018 from simulations that employ the T a and T w -based precipitation partitioning methods and the Global Historical Climatology Network data (OBS) for: (a) all sites as a function of elevation, and (b) averaged for sites within the two major elevation bands.

Figure 7 .
Figure 7. Spatiotemporal variations of T a (top), relative humidity (middle), and ΔSF (bottom) under historical (left) and far future (R85, right) periods, showing the ensemble median.Inset in (f) depicts the elevation distribution in the Colorado River Basin.Δ signs represent TW-TA.

Figure 10 .
Figure 10.Mean monthly snow water equivalent (SWE) and Q in (a, d) historical and far future periods under (b, e) R45 and (c, f) R85 scenarios for the entire Colorado River Basin organized by water year (October 1-September 30).Results for the historical and far future periods are based on ensemble median values from the statistically-downscaled General Circulation Models.Dashed lines in (a, d) represent reference SWE and Q from UA product and USBR naturalized streamflow, respectively.Note that the UA product is only available for a portion of the historical period (1982-2005).

Figure 11 .
Figure 11.Conceptual diagram of annual-averaged water balance components over the Colorado River Basin (units in mm/ yr).All values represent the ensemble median for TWC during historical (a) and far future periods for (b) R45 and (c) R85 scenarios.Parentheses are differences (TWC-TA).

Figure 12 .
Figure 12.Subbasin streamflow changes including absolute values (Q Diff , million acre-feet, bars) and relative differences with respect to TA (dQ Diff , %, symbols) in the TA, TW, and TWC simulations for: (a) R45 and (b) R85.Results are ensemble median values.

Table 1
Spatially-Averaged Elevation (H) and Mean Annual T a , Relative Humidity (RH), and T w in Subbasins ObtainedOver 1981-2018 From Livneh  et al. (2015)Note.Parentheses show spatial standard deviations.
WANG ET AL.