On the value of satellite remote sensing to reduce uncertainties of regional simulations of the Colorado River

As the major water resource in the southwestern United States, the Colorado River is experiencing decreases in naturalized streamflow and is predicted to face severe challenges under future climate scenarios. To better quantify these hydroclimatic changes, it is crucial that the scientific community establishes a reasonably accurate understanding of the spatial patterns associated with the basin hydrologic response. In this study, we employed remotely sensed land surface temperature (LST) and snow cover fraction (SCF) data from the Moderate Resolution Imaging Spectroradiometer (MODIS) to assess a regional hydrological model applied over the Colorado River Basin between 2003 and 2018. Based on the comparison between simulated and observed LST and SCF spatiotemporal patterns, a stepwise strategy was implemented to enhance the model performance. Specifically, we corrected the forcing temperature data, updated the time-varying vegetation parameters, and upgraded the snow-related process physics. Simulated nighttime LST errors were mainly controlled by the forcing temperature, while updated vegetation parameters reduced errors in daytime LST. Snow-related changes produced a good spatial representation of SCF that was consistent with MODIS but degraded the overall streamflow performance. This effort highlights the value of Earth observing satellites and provides a roadmap for building confidence in the spatiotemporal simulations from regional models for assessing the sensitivity of the Colorado River to climate change.

1 Introduction

Physically based numerical models of the coupled water–energy cycle have emerged as powerful tools to address critical societal needs (Fatichi et al., 2016), including flood forecasting (Maidment, 2017), irrigation operation (Gibson et al., 2017), weather and climate prediction (Baker et al., 2017; Senatore et al., 2015), and evaluations of water scarcity (Zhou et al., 2016). Over the last 3 decades, several types of hydrologic models have been developed with different levels of conceptualization that often change with the domain size due to computational constraints. One class of models, denoted as regional or macroscale models, was originally designed to serve as land surface scheme of atmospheric models and is routinely used to simulate hydrologic processes in continental basins (>10⁵ km²) at spatial resolutions of 10 to 25 km (e.g., Lawrence et al., 2011; Liang et al., 1994; Niu et al., 2011). These processes include infiltration, evapotranspiration, runoff production, and snow accumulation and ablation that are typically simulated in a regular grid without considering lateral transfers across cells (Clark et al., 2015). In recent years, the National Water Model has combined a regional hydrologic model applied at the unprecedented resolution of 1 km with routing schemes to generate operational hydrologic predictions over the continental United States (Lahmers et al., 2019, 2021).

In many cases, hydrologic models are applied under prescribed meteorological forcings using an optimal set of parameters that are calibrated by minimizing differences between simulated streamflow and observations at one or more locations (e.g., Gou et al., 2021; Li et al., 2019; Nijssen et al., 1997; Xiao et al., 2018; Yun et al., 2020; Zhang et al., 2017). While widely used, this approach has two important limitations. First, input and structural uncertainties are often not taken into account (Gupta and Govindaraju, 2019), causing an inflation of parametric uncertainty that can exacerbate the problem of equifinality (Beven and Binley, 1992). Second, this calibration method relies only on an aggregated measure of the hydrologic response and does not consider the model's ability to capture the spatially variable internal processes (Becker et al., 2019; Ajami et al., 2004). As a result of these two limitations, this calibration approach could cause the undesirable outcome that the model provides the right answer for the wrong physical reasons (Rajib et al., 2018; Tobin and Bennett, 2017), which can in turn induce wrong conclusions when the model is applied under nonstationary conditions due to changes in land cover and/or climate.

Satellite remote sensors provide spatially distributed estimates of hydrologic states and fluxes, including soil moisture (Entekhabi et al., 2010; Njoku et al., 2003; Kerr et al., 2001), land surface temperature (LST; Shi and Bates, 2011; Wan and Dozier, 1996), snow cover fraction (SCF, Painter et al., 2009), evapotranspiration (Boschetti et al., 2019; Fisher et al., 2020), and changes in water storage (Tapley et al., 2004). These products can reduce parametric, structural, and input uncertainties in hydrologic models by including additional constraints in the calibration process (Wood et al., 2011; Fatichi et al., 2016; Ko et al., 2019). Despite this potential, the use of remote sensing products to reduce hydrologic simulation uncertainty has been explored in only a few studies. For instance, in studies by Corbari and Mancini (2014), Crow et al. (2003), and Zink et al. (2018), satellite LST was used with river discharge to calibrate model parameters, and it was found that including LST in the process improved the simulation of evapotranspiration as estimated by eddy covariance towers or other satellite products. This outcome was also found by Gutmann and Small (2010), who applied a regional model at 14 flux towers and showed that incorporating remotely sensed LST estimates in the calibration allowed the achievement of two-thirds of the improvements gained by ingesting more accurate ground LST data. In other efforts, satellite LST products have been used to verify the performance of hydrologic models, as done by Koch et al. (2016), with the North America Land Data Assimilation System (NLDAS), Xiang et al. (2014), with the TIN-based Real-time Integrated Basin Simulator (tRIBS), Xiang et al. (2017), with the Weather Research and Forecasting (WRF)-Hydro model, and Wang et al. (2021), with the Variable Infiltration Capacity (VIC) model. Finally, a few studies have enhanced streamflow simulations (Bennett et al., 2019; Bergeron et al., 2014; Tekeli et al., 2005) by improving the timing of snowmelt using remotely sensed snow cover fields.

The Colorado River basin (CRB) is a regional watershed where hydrologic simulations are needed to support short- and long-term water management decisions. Its water resources are used by almost 40 million people in seven states of southwestern U.S. (Arizona, California, Colorado, Nevada, New Mexico, Utah, and Wyoming) to irrigate ∼22 000 km² of land and to generate over 4200 MW of hydroelectric power (USBR, 2012). The mean annual discharge of the CRB is 20.2 km³, with high interannual variability resulting from large variations in climatic forcings (Christensen et al., 2004; Gautam and Mascaro, 2018). Until 2021, the CRB was able to meet the demand of all users by storing runoff in a large system of dams, mainly operated by the U.S. Bureau of Reclamation (USBR), and transporting water through canals and aqueducts, including the Central Arizona Project. However, declines in the mean flow observed over the last 2 decades (Hoerling et al., 2019; Udall and Overpeck, 2017), combined with increasing demands, led to the first-ever declaration of water shortages in the CRB in January 2022. The water cuts affecting users in Arizona and Nevada (CAP, 2021) are expected to become more severe in the near future and impact the agricultural sector (Mitchell et al., 2022; Norton et al., 2021).

In previous studies on the hydrologic responses of the CRB using the VIC model, confidence in the model results was built mainly through comparisons against estimates of naturalized flow (e.g., Christensen et al., 2004; Vano et al., 2012, 2014; Xiao et al., 2018). The CRB is characterized by a marked difference between the colder and wetter Upper Basin, where more than 90 % of streamflow is generated (Li et al., 2017), and the warmer and drier Lower Basin, with reduced runoff production due to low precipitation, high evaporative demand, and channel transmission losses (Rajagopalan et al., 2009). As a result of this large contrast, limiting the calibration of VIC to the use of naturalized flow in the Upper Basin may lead to uncertainty in its ability to simulate the spatiotemporal hydrologic response.

The objective of this study is to improve the physical reliability of VIC simulations in the CRB by incorporating remotely sensed fields of LST and SCF obtained from the Moderate Resolution Imaging Spectroradiometer (MODIS). LST is an important variable that impacts the coupled water–energy balance, while SCF provides information on snow conditions which are crucial for quantifying runoff generation. We start from a parameterization of VIC that led to good estimates of monthly discharge in the period 2003–2018. We then apply a stepwise procedure to reduce uncertainties in model forcings, parameters, and structure based on comparisons of simulated and remotely sensed LST and SCF fields. While based on VIC, the methods proposed here can provide guidance to refine the calibration and reduce uncertainties in other physically based hydrologic models and to identify areas for structural improvement.

Figure 1 (a) Digital elevation model of the CRB. (b) Channel network and eight subbasins analyzed in this study. The red circle marks Imperial Dam. (c) Dominant vegetation type in each pixel indicated in the legend. (d) Time-averaged vegetation fraction, f_v. (e) Total soil depth. All maps are at 0.0625^∘ (∼6 km) spatial resolution. Values of f_v and soil depth are from the baseline simulation.

Table 1 Spatially averaged mean annual precipitation (P), snow water equivalent (SWE), runoff (Q), and runoff ratio ($Q/P$), along with area, mean elevation, mean soil depth, and percentage of trees in the CRB and its subbasins.

2 Study area and datasets

2.1 Study basin

The CRB has a total area of approximately 630 000 km², covering seven states in United States and a small portion of Mexico. Here, we considered the drainage area above Imperial Dam, plus the Gila River (Fig. 1). The Colorado River Compact of 1922 divides the CRB into the Upper and Lower basins. As revealed by the land cover map shown in Fig. 1c, most of the basin is covered by shrub or scrub ecosystems (∼60 %), followed by various forest types (∼24 %). Table 1 summarizes the mean hydroclimatic and land surface features of the subbasins. The Upper Basin consists of the Green, Upper Colorado, Glen Canyon, and San Juan River subbasins. These higher-elevation subbasins (except Glen Canyon) receive more snowfall than the rest of the CRB, resulting in the presence of a significant snowpack (mean annual snow water equivalent, or SWE, ranges from 13.7 to 58.8 mm) that eventually leads to the generation of ∼90 % of the CRB runoff. While the Lower Basin receives about 60 % of the mean annual precipitation of the subbasins in the Upper Basin per unit area, its runoff ratio (i.e., the fraction of annual precipitation becoming runoff) is 3 times smaller than that of the Upper Basin.

2.2 Remote sensing and ground-based datasets

We integrated different remotely sensed and ground-based data. Meteorological forcings were obtained from the gridded (0.0625^∘ or ∼6 km) daily datasets of Livneh et al. (2013) and Su et al. (2021) for precipitation, maximum temperature, minimum temperature, and wind speed. We also used the Parameter-elevation Regressions on Independent Slopes Model (PRISM) 30-year normal (Di Luzio et al., 2008) for temperature corrections. For assessing streamflow performance, we used monthly naturalized flow records from USBR at four interior locations of the Upper Basin. Note that this is the largest available resolution for the reconstructed naturalized flow, since the river is highly regulated. To improve the simulation of spatial patterns, we used two products from the Aqua MODIS sensor, i.e., daily LST (MYD11A1) and monthly SCF (MYD10CM). The LST product is available at 1 km resolution twice a day at about 13:00 LT (local time; daytime) and 01:00 LT (nighttime; Wan et al., 2021). The percent of missing data, largely due to cloud cover, varies from 42 % to 95 %, with larger values in the winter season and July (Fig. S1 in the Supplement). Monthly SCF is provided at 0.05^∘ (∼5 km) resolution as the average of SCF for days with a prescribed level of sky clearness (Hall and Riggs, 2021). Both MODIS products were aggregated to the 0.0625^∘ scale used in the model. We also validated simulated and remotely sensed LST using measurements at 14 eddy covariance towers (Baldocchi et al., 2001) selected based on available data (>300 d during 2003–2018). The station locations are shown in Fig. S2, with 12 located in the Lower Basin at elevations from 987 to 2618 m. Five stations were forested, and the remaining were covered by a short canopy. We extracted records of observed longwave radiation at the stations and used them to compute LST, following Wang et al. (2021). We also used the National Land Cover Database (NLCD) Multi-Resolution Land Characteristics (MRLC) rangeland and tree canopy cover products, which contain canopy cover fraction at 30 m resolution for forests and shrublands (Homer et al., 2020).

Figure 2 Schematic explaining how LST is computed in VIC (LST_V) as compared to MODIS (LST_M) in a pixel covered by short vegetation (tile A) and tall trees (tile B). f_v is the vegetation fraction, T_air is the air temperature, T_s, T_f, and T_c are simulated temperatures for the surface, canopy, and canopy air, LW_d,v is the downward longwave radiation from the canopy, and LW_d is the downward longwave radiation from the atmosphere. The subscripts A and B refer to the variables in each tile.

3 Methods

3.1 The variable infiltration capacity model

We used the VIC model version 5.0 (Hamman et al., 2018) to simulate the hydrologic response of the CRB from 2003 to 2018 at an hourly time step and 0.0625^∘ resolution. VIC is a macroscale, physically based model that solves the water and energy balance on a regular grid. Land surface heterogeneity in each cell is modeled through land cover tiles, each with a single vegetation class on top of a three-layer soil column. The model requires meteorological forcings as inputs and returns outputs over the grid. Fluxes and state variables simulated at grid cells are calculated as the areal-weighted average of separate computations of the water and energy balances for each land cover tile. Here, we adopted the VIC version with the clumped vegetation scheme proposed by Bohn and Vivoni (2016), where the vegetation fraction (f_v) accounts for spacing between plants in each tile. This modification allows simulating the energy balance with a higher fidelity, as shown by Bohn and Vivoni (2016), through the comparison with ground estimates of evapotranspiration in the southwestern U.S. and northwestern Mexico.

Table 2 List of spatially variable forcings, vegetation and soil parameters, and state variables involved in the computation of the energy balance (symbols defined in the main text and Appendix A). Forcings and state variables vary each hour. Parameters are either constant in time or vary each month (denoted with an asterisk^*).

Since our adjustment strategy is based on the comparison of simulated and remotely sensed LST and SCF, we describe how these variables are simulated using the schematic in Fig. 2. The governing equations are reported in Appendix A, while the most influential parameters are in Table 2. In our simulations, 16 vegetation classes are used, which include four types of tall trees, namely deciduous forest, evergreen forest, mixed forest, and woody wetlands. For other canopy types (e.g., tile A in Fig. 2), the energy balance is solved over a control volume that combines the fractions of vegetation (f_v,A) and bare soil ($\mathrm{1}-f_{\mathrm{v},A}$) using a weighted aerodynamic resistance. A single surface temperature (T_s,A) is computed and assumed to be uniform over the tile and equal to the foliage temperature ($T_{\mathrm{f},A}=T_{\mathrm{s},A}$). For tall trees (e.g., tile B in Fig. 2), a vegetated overstory and an understory without vegetation are introduced. If snow is absent, the overstory foliage temperature is assumed equal to air temperature ($T_{\mathrm{f},B}=T_{\mathrm{air}}$), and a single T_s,B in the understory is calculated with the scheme described above. When snow is present, T_s,B is calculated by solving the energy balance in the overstory, understory, and the atmosphere surrounding the canopy. Since the satellite sensor observes the top of the surface, the simulated LST by VIC (LST_V) that is compared against MODIS (LST_M) is the weighted average of foliage temperature in tiles with tall trees and the ground temperature in other tiles. In the case of Fig. 2, this leads to the following:

$$\begin{array}{}\text{(1)}& {\mathrm{LST}}_{\mathrm{V}}={\displaystyle \frac{A_A\cdot T_{\mathrm{s},A}+A_B\left[f_{\mathrm{v},B}\cdot T_{\mathrm{f},B}+\left(\mathrm{1}-f_{\mathrm{v},B}\right)\cdot T_{\mathrm{s},B}\right]}{A_A+A_B}},\end{array}$$

where A_A and A_B are the areas of tiles A and B, respectively.

To compute SCF in the grid cells, VIC allows the subdivision of each tile into elevation bands to capture changes in forcing temperature due to terrain heterogeneity. Elevation bands are the same for all tiles in a grid cell and limited typically to three bands in total. Given the mean elevation of each elevation band, the air temperature forcing is adjusted using a lapse rate of −6.5 ^∘C km⁻¹ and is then used to solve the energy balance within each tile. Depending on temperature and precipitation, snow may be simulated within a tile, and SWE is calculated. When SWE > 0, SCF is assumed to be 100 %, such that a tile within that elevation band is fully covered with snow; otherwise, SCF is 0, and the elevation band within the tile is snow-free (i.e., a binary outcome). SCF in the grid cell is the area-weighted average of the SCFs from all tiles and elevation bands.

3.2 Baseline simulation

We created a first model parameterization, labeled as baseline, based on applications by Xiao et al. (2018) and Bohn and Vivoni (2019). Hourly gridded meteorological forcings were generated from the daily grids of Livneh et al. (2013) and Su et al. (2021), using MetSim (Bennett et al., 2020; Bohn et al., 2013; Bohn and Vivoni, 2019). Model parameters were obtained from Livneh et al. (2015), with a few updates as follows. Land surface parameters were based on MODIS and NLCD products from Bohn and Vivoni (2019), which include a land cover classification and climatological monthly means of leaf area index (LAI), f_v, and albedo. We replaced the elevation data used in prior VIC studies with the 30 m United States Geological Survey (USGS) National Elevation Dataset (USGS, 2017). The model was tested against monthly naturalized streamflow records, by manually adjusting seven soil parameters that affect runoff production and the parameters controlling the relation between snow albedo with snow age. As shown in Fig. S3, under the baseline simulation, VIC captured the monthly streamflow in key subbasins of the Upper Basin, where most runoff is produced, and at the basin outlet well with a Nash–Sutcliffe efficiency (NSE) > 0.9.

Figure 3 Flowchart of the stepwise calibration procedure.

3.3 Model improvements with remote sensing products: overview of the stepwise calibration strategy

The baseline simulation was aimed at reproducing the streamflow response and did not consider the model's ability to capture spatial patterns of hydrologic variables. We designed a stepwise strategy aimed at reducing the three main sources of uncertainty in the simulation of LST and SCF. A schematic of the procedure is reported in Fig. 3; here, we provide an overview of the steps and describe the details of each step in the corresponding sections in Sect. 4. In the first step (“Forcing-adj” or forcing adjustment), we targeted input uncertainty and modified air temperature to reduce errors in nighttime LST. In the second step (“Veg-adj” or vegetation adjustment), we focused on modifying spatially variable vegetation parameters affecting daytime LST identified among those reported in Table 2. The first two steps were guided by metrics quantifying the agreement between simulated and remotely sensed LST, including the correlation coefficient (CC), root mean squared error (RMSE), and bias (mean LST_V–mean LST_M) between the (1) time series of daily LST_V and LST_M at each grid cell and (2) daily spatial maps. These metrics were obtained for both daytime and nighttime through comparisons at the MODIS overpass time. To further quantify the improvements in our calibration approach, for each step we computed the structural similarity index measure (SSIM; Wang and Bovik, 2002) and the SPAtial EFficiency (SPAEF) metric (Koch et al., 2018) between spatial maps of observed and simulated long-term climatological mean LST; these two metrics were chosen since they have been specifically designed to compare spatial patterns.

After improving LST, we reduced structural uncertainty by modifying the computation of the snow energy balance in a step labeled as “Snow-adj” (or snow adjustment). As described above, when snow exists in tiles covered by tall trees, the downward longwave radiation into the understory (or ground) snowpack is assumed to originate from the overstory (indicated as LW_d,v in Fig. 2; tile B). For areas without tall trees, the downward longwave radiation reaching the understory comes from the atmosphere (indicated as LW_d). To account for this in the clumped canopy scheme, we modified the downward longwave radiation as the weighted average of [$f_{\mathrm{v}}\cdot {\mathrm{LW}}_{\mathrm{d},\mathrm{v}}+(\mathrm{1}-f_{\mathrm{v}})\cdot {\mathrm{LW}}_{\mathrm{d}}$]. In addition, we adjusted the empirical relation controlling the change in albedo during snowmelt to reduce the bias between VIC and MODIS SCF. All modifications of the model parameters were performed via manual tuning.

4 Results

4.1 Comparison of VIC and MODIS LST with ground observations

First, we provide an overview of the comparison among the time series of LST that were (1) observed at the 14 eddy covariance stations, (2) simulated by VIC, and (3) retrieved from MODIS at the co-located 6 km pixel. The error metrics for the 14 stations are summarized through box plots in Fig. 4a–c, while the time series of LST at a representative site for daytime and nighttime are shown in Fig. 4d and e. Station values and VIC simulations at the overpass times were extracted for comparison with MODIS. Dates with missing data in MODIS and station records were not considered. We find MODIS LST to be very strongly correlated with ground measurements (CC > 0.91) and characterized by RMSE from ∼1.5 to 5.3 ^∘C. Bias is slightly positive (negative) during the daytime (nighttime), with a median of 0.3 ^∘C (−1.6 ^∘C). The error metrics for VIC reveal that performance degrades moderately with larger variability across the stations, where CC ranges from 0.70 to 0.95, the median RMSE is 6.3 ^∘C (5.8 ^∘C) for daytime (nighttime), and the median bias is 1.1 ^∘C (−3.3 ^∘C) for daytime (nighttime). The error metrics against ground data provide a reference for evaluating the model improvements, as discussed next.

Figure 4 (a–c) Box plots of CC, RMSE, and bias comparing VIC and MODIS LST to observations at 14 sites. Time series of daytime (d) and nighttime LST (e) at one site (station ID “FuF”, whose location is shown in Fig. S2).

4.2 Errors in the simulation of LST in the baseline simulation and their controls

Figure 5 shows maps of CC, RMSE, and bias of the time series of LST_V and LST_M at each pixel for daytime and nighttime periods over the entire simulation from 2003 to 2018. To help the interpretation, box plots of the metrics in the grid cells within the CRB and three subbasins are presented in Fig. 6. Results for other subbasins are reported in Figs. S4–S6 and Table S1 in the Supplement.

Figure 5 Spatial maps of (a, d) CC, (b, e) RMSE, and (c, f) bias between time series of LST_V and LST_M over 2003–2018 at each pixel. The top (bottom) row presents daytime (nighttime) comparisons.

Figure 6 Box plots of (a–d) CC, (e–h) RMSE, and (i–l) bias between time series of LST_V and LST_M in pixels of the entire CRB and of three representative subbasins. Box plots show the median, with 50 % and 90 % confidence intervals. Different simulations are plotted in different colors.

Overall, CC is high (>0.8) throughout the CRB, with values like those found against station data. CC is relatively higher for daytime than nighttime. On the other hand, RMSE maps show that simulated LST matches better with MODIS during the nighttime, with values largely consistent with those found for stations. For both times of the day, RMSE is slightly larger in the Upper Basin. Results for RMSE suggest that model performance for LST is relatively better at nighttime without solar radiation forcing and tends to be better in drier and hotter regions in the Lower Basin. Bias maps reveal that simulations of LST during daytime (nighttime) are warmer (cooler) than MODIS observations in most of the CRB, with a median bias of 1.2 ^∘C (−0.7 ^∘C). These findings are largely consistent across the subbasins and with the station observations.

Spatial patterns of the metrics are complex, suggesting that LST simulation errors are impacted by several model parameters and forcings. To gain insights into these controls, we computed the correlation coefficient between the maps of error metrics between the time series and key parameters or forcings involved in the energy balance. Model parameter maps were created by calculating the area-weighted averages within each grid cell. For monthly LAI, albedo, and f_v, we computed the annual mean map. For T_air, we calculated the mean across the entire study period. Figure 7 summarizes the results in each subbasin for RMSE and bias using heatmaps (also see Fig. S7 for CC). For daytime LST, the key factors change across the subbasins, while results are more spatially uniform for nighttime LST. During daytime, we found that the Green and Upper Colorado subbasins dominated by snow and evergreen forests exhibit different controls compared to the other subbasins. Here, RMSE is highly correlated to f_v and LAI, while bias is mainly controlled by T_air. In the other subbasins, albedo and, to a lesser extent, T_air are the dominant factors related to daytime RMSE. Different parameters affect the patterns of bias, including albedo in all subbasins, most vegetation parameters, root depth in the San Juan and Little Colorado, and T_air in the Little Colorado. Considering nighttime LST, T_air and, to a lower degree, soil depth are the main factors related to RMSE at all sites. Interestingly, nearly all parameters and T_air are linked to nighttime bias. This is explained by considering that T_air is correlated with elevation, and elevation is correlated with all other parameters (Fig. S8).

Figure 7 Heatmaps showing the Pearson correlation coefficient between (1) the spatial map of T_air or key soil and vegetation parameters involved in the energy balance and (2) the spatial map of the error metrics – (a, c) RMSE and (b, d) bias – between the time series of LST_M and LST_V for the baseline simulation. The correlation coefficients are computed for each subbasin. Symbols are explained in Table 2. Top (bottom) row is for daytime (nighttime) LST.

Figure 8 presents the intra-annual variability in the error metrics between daily pairs of LST_V vs. LST_M fields, shown as monthly averages. As found previously, CC is high for both times of the day and relatively higher for daytime, while RMSE is larger during the daytime. VIC simulations during the daytime are positively biased throughout the year, while bias changes sign for nighttime LST, being positive in winter and negative from April to July. In addition, both the RMSE and bias of daytime LST are higher from April to July. This indicates that simulated daytime LST degrades when incoming solar radiation is high, especially during snowmelt events after peak SWE, typically around the end of March. To corroborate this, we repeated the analyses in snow-dominated grid cells (mean annual maximum SWE > 30 mm) and for all other cells, finding higher daytime RMSE in April for snow-dominated cells than other cells, indicating that the LST during the ablation process is also more difficult to capture.

Figure 8 Time series of multiyear monthly average CC, RMSE, and bias between VIC and MODIS daily LST fields for the baseline simulation and each adjustment step.

4.3 Stepwise reduction in uncertainty in the simulation of LST and SCF

4.3.1 Forcing adjustment

We first focused on the improvement in simulated LST at nighttime. Figure 7 indicates that T_air is a key input affecting the energy balance at nighttime. Alder and Hostetler (2019) compared two air temperature datasets, finding that Livneh et al. (2013) products tend to be colder than PRISM in the mountain areas of the CRB. Based on this, we adjusted the daily minimum and maximum T_air in Livneh et al. (2013) and Su et al. (2021) to match the climatological (1981–2010) monthly means from PRISM. If $T_{\mathrm{air},\mathrm{d},\mathrm{m}}^{\mathrm{L}}$ is the maximum or minimum daily T_air on day d and month m, then the bias-corrected value, $T_{\mathrm{air},\mathrm{d},\mathrm{m}}^{\mathrm{L},\mathrm{BC}}$, was obtained as follows:

$$\begin{array}{}\text{(2)}& T_{\mathrm{air},\mathrm{d},\mathrm{m}}^{\mathrm{L},\mathrm{BC}}=T_{\mathrm{air},\mathrm{d},\mathrm{m}}^{\mathrm{L}}-\left({\stackrel{\mathrm{‾}}{T}}_{\mathrm{air},\mathrm{m}}^{\mathrm{P}}-{\stackrel{\mathrm{‾}}{T}}_{\mathrm{air},\mathrm{m}}^{\mathrm{L}}\right),\end{array}$$

where ${\stackrel{\mathrm{‾}}{T}}_{\mathrm{air},\mathrm{m}}^{\mathrm{P}}$ and ${\stackrel{\mathrm{‾}}{T}}_{\mathrm{air},\mathrm{m}}^{\mathrm{L}}$ are the climatological monthly means of maximum or minimum T_air from PRISM and Livneh et al. (2013), respectively. Once we bias-corrected T_air, we regenerated the hourly forcings using MetSim. As shown in Fig. 9, the Forcing-adj simulations improved bias, which was reduced in most subbasins. The nighttime RMSE also slightly decreased throughout the basin. These outcomes are reflected in the time series of Fig. 8 that also show that improvements (lower RMSE and bias) occur largely in the warm season. On the other hand, the Forcing-adj simulations did not improve VIC performance during the daytime, only yielding a slight increase in bias (Figs. 6 and 8) that was fixed in the next steps.

Figure 9 Spatial maps of the RMSE and bias between time series of nighttime LST_V and LST_M during 2003–2018 at each pixel for all steps. Top (bottom) row presents results of RMSE (bias).

4.3.2 Vegetation parameter adjustment

Figure 7 shows that both static and time-varying vegetation parameters affect the error metrics of LST. In the Veg-adj step, we modified a set of influential parameters by incorporating new datasets. We first replaced the climatological mean monthly values of LAI, albedo, and f_v with yearly varying monthly estimates from MODIS. Second, we updated f_v using new products from MRLC. In the baseline simulation, f_v was derived from the normalized difference vegetation index (NDVI) retrieved from MODIS (Bohn and Vivoni, 2016, 2019). MRLC released 30 m grids of mean annual f_v for major vegetation types in the CRB that were used to linearly rescale values of f_v in the shrub and trees classes to match the annual climatology of MRLC as follows:

$$\begin{array}{}\text{(3)}& f_{\mathrm{v},\mathrm{m}}^{\mathrm{Resc}}=f_{\mathrm{v},\mathrm{m}}^{\mathrm{b}}{\displaystyle \frac{{\stackrel{\mathrm{‾}}{f}}_{\mathrm{v}}^{\mathrm{MRLC}}}{{\stackrel{\mathrm{‾}}{f}}_{\mathrm{v}}^{\mathrm{b}}}},\end{array}$$

where $f_{\mathrm{v},\mathrm{m}}^{\mathrm{b}}$ is f_v in month m used in the baseline simulation, $f_{\mathrm{v},\mathrm{m}}^{\mathrm{Resc}}$ is the rescaled value, and ${\stackrel{\mathrm{‾}}{f}}_{\mathrm{v}}^{\mathrm{MRLC}}$ and ${\stackrel{\mathrm{‾}}{f}}_{\mathrm{v}}^{\mathrm{b}}$ are long-term mean annual values of MRLC and the baseline parameters.

Figure 7 indicates that r_min, r_arc, d₀, and z₀ affect errors in the simulation of LST, especially in the Green and Upper Colorado subbasins. Distributed estimates for these parameters are not currently available. Thus, we adjusted their values to reduce the bias between daytime LST_V and LST_M, guided by the process equations reported in Appendix A. Reducing z₀ and d₀ leads to lower aerodynamic resistance and higher sensible heat flux and, in turn, lower LST_V. Increases in r_min and r_arc lead to lower values of latent heat flux and higher LST_V. Adjusting z₀ has a greater impact than modifying the other parameters, such that the iteratively scaling of z₀ in each pixel was performed at 25 %, 50 %, 150 %, or 250 %, depending on the daytime LST bias (Fig. 10). Changes were limited within physically plausible ranges. Next, we applied the same method to update d₀, r_min, and r_arc, but variations for these three parameters were minimal, as documented in Fig. S9.

Figure 10 Same as Fig. 9 but for daytime LST.

The Veg-adj simulation did not lead to significant changes in model performance at nighttime, confirming that the dominating factor affecting nighttime LST was T_air. On the other hand, improvements in the simulation of daytime LST were remarkable. Figure 6 shows that both RMSE and bias were reduced at all locations, both in terms of the median (∼0.9 ^∘C) and variability in each subbasin (lower width of the confidence intervals), with values slightly higher than those found between MODIS and station observations (Fig. 4). These improvements were even more apparent in the maps of Fig. 10, which also showed that the complex spatial patterns of the errors in the baseline simulation have been replaced by more uniform and smoother patterns. The Veg-adj simulation also decreased large errors in the simulation of daytime LST from April to July, with lower RMSE, higher CC, and bias close to 0 ^∘C throughout the year (Fig. 8).

4.3.3 Adjustment of snow dynamics

The Snow-adj step was aimed at improving the simulation of SCF. We first modified the computation of longwave radiation for tall trees, which improved the simulation of SCF during the snow accumulation season. Next, a parameter of the relation controlling the decay of snow albedo was modified from 0.92 to 0.80, leading to an enhanced simulation of SCF in the ablation season. Figure 11 presents bias maps between simulated and observed mean monthly SCF and seasonality of SCF in snow-dominated cells for the baseline, Veg-adj, and Snow-adj simulations. Time series of SCF in 2 pixels are also shown to visualize differences in regions with positive and negative bias. In the baseline simulation, SCF bias was positive, which occurs mainly during May through October. Forcing corrections reduced SCF, as T_air was increased in mountainous areas. Adjustments in the Snow-adj step reduced bias in most locations during the accumulation and ablation seasons. When averaged over time and in the CRB, SCF bias was relatively small. When focusing on single pixels, however, the bias magnitude was larger, with differences in seasonality depending on location. For example, bias reached +20 % (Fig. 11c) from April to December and −20 % (Fig. 11d) from November to March. As expected, Snow-adj changes mainly impacted LST simulations in mountains, while a marginal influence occurred in the rest of the CRB. Overall, the daytime LST bias map improved, while RMSE in mountain regions for both daytime and nighttime remained similar. To complete the model performance assessment, we reported, in Figs. S10 and S11, the maps of simulated and observed long-term climatology of monthly SCF in the snow season and LST, respectively, over 2003–2018. Error metrics between the maps are presented in Table S2, which shows that the overall trend of the metrics specifically designed to compare spatial patterns, SSIM and SPAEF, are in line with the changes in RMSE and bias that have been used in the rest of the paper.

Figure 11 (a) Spatial maps of bias between mean monthly SCF (VIC minus MODIS). Circles indicate the locations of two grid cells with positive and negative bias. (b) Time series of multiyear mean monthly SCF (in %) for snow-dominated cells. RMSE and bias from monthly SCF comparisons are reported. (c, d) Same as panel (b) but for sites with positive and negative bias, respectively.

Figure 12 Monthly time series of naturalized streamflow (NFL) and streamflow from baseline, Forcing-adj, Veg-adj, and Snow-adj simulations at (a) Green, (b) Upper Colorado, (c) San Juan, and (d) Upper Basin for 2003–2013. NSE values are also reported.

Figure 13 (a) Climatological monthly mean of the water balance components for the baseline simulations in the Upper Basin. P is precipitation, ET is evapotranspiration and sublimation, RO is surface and underground runoff, and ΔSM (ΔSWE) is the differences between soil moisture (snow water equivalent) at the end and beginning of the month. (b–e) Difference between each variable for the Forcing-adj, Veg-adj, and Snow-adj simulations and the baseline simulations.

4.4 Impacts on VIC streamflow performance and water balance

As shown previously (Corbari and Mancini, 2014; Crow et al., 2003), improving the simulation of hydrologic spatial patterns could affect streamflow performance since structural limitations and different degrees of conceptualization require further tuning. We investigated this in Fig. 12, using time series of monthly runoff in the Green and San Juan subbasins and the Upper Basin. Model performance is very good for baseline simulations, since its calibration was tailored to naturalized streamflow records. Forcing and vegetation parameter adjustments slightly lowered performance (changes in NSE ≤ 0.05), whereas changes for the snow adjustment led to streamflow overestimation in May in all subbasins, especially in the Green subbasin (NSE reduced to 0.57). Overall, simulated streamflow performance here is consistent with Tang and Lettenmaier (2010), who found that incorporating MODIS snow cover degrades streamflow metrics. We attribute this degradation in performance to a number of reasons. First, remotely sensed spatiotemporal data of SCF have limitations in its usefulness for tracking SWE, which is the modeled state variable more directly impacting streamflow. Second, VIC uses a binary scheme for depicting SCF in elevation bands within each time of each grid cell, limiting its accuracy in representing topographic variations. To address these limitations, enhancements are needed in both the 1simulation of snow physics and remote sensing of the spatial variation in snow depth or SWE at high spatiotemporal resolutions.

In addition to streamflow, we explored the impacts of each calibration step on the water balance. For this aim, we computed the climatological monthly mean of the water balance components for the Upper Basin, where most runoff is generated. Results are presented in Fig. 13, which shows fluxes (P, ET, and RO; see the Fig. 13a caption for their definition) and changes in state variables (ΔSM and ΔSWE) for the baseline simulations (Fig. 13a) and the difference between a given variable simulated in each calibration step and the variable from the baseline simulation (Fig. 13b–d). The Forcing-adj and Veg-adj steps lead to small changes in ET and RO, with a decrease in both fluxes in the summer months and an increase in the other months. The modification of these fluxes is due to a change in the storage components, with (1) lower SWE (i.e., negative ΔSWE) and higher SM from November to February and (2) higher SWE and lower SM from March to July. The Snow-adj step modifies the seasonality of SWE compared to the baseline by increasing this storage component in February and March and reducing it in April and May. This, in turn, leads to an opposite behavior for SM, which is ultimately translated into a positive (negative) change in RO in May and June (July and August). In all cases, the changes in runoff occurred in a similar way for both the surface and underground components.

5 Summary and conclusions

In this study, we made improvement to a regional hydrologic model in the Colorado River basin using MODIS observations of land surface temperature and snow cover. Based on the remotely sensed data, we corrected the meteorological forcings, updated the vegetation parameters, and revised snow-related processes to enhance the model performance. The adjustments increased the consistency between VIC and MODIS LST and SCF fields, thus enhancing credibility of the spatial simulations. Our conclusions are summarized as follows:

1. MODIS products provided spatiotemporal information that can be used to identify uncertainties in a hydrologic model calibrated with streamflow records at a few locations. Although baseline simulation performance for LST was high (mostly CC > 0.8), spatial errors within the CRB were non-negligible. The baseline simulation had lower RMSE of LST for nighttime and cold season conditions. Baseline model discrepancies were primarily associated with energy exchanges at land surface during periods of higher solar radiation.

2. Simulated nighttime LST values were dominated by the initial air temperature, such that improvements were obtained from forcing corrections. This led to a reduction in nighttime LST bias from −7 to 6 ^∘C in the baseline case to −5 to 5 ^∘C in the Forcing-adj simulation. Vegetation adjustments led to large improvements in daytime LST, with RMSE reductions from 7.5 to 2.5 ^∘C, but were less effective at night. In addition, the range of daytime RMSE of LST was reduced from 4 to 10 ^∘C in the baseline case to 2.5 to 3.5 ^∘C in the Veg-adj simulation.

3. Updated snow physics reduced the negative bias in SCF during the accumulation season. We further adjusted melting snow albedo to improve performance in the ablation period. Unlike other modifications, runoff was substantially impacted by the lower snow albedo. Thus, the consistency between VIC and MODIS snow cover did not ensure an improved streamflow simulation, demonstrating the limitations of the regional application in accurately capturing the variation in SWE in mountainous areas. A possible solution to improve the spatial credibility of the hydrologic model without degrading streamflow performance is by incorporating satellite products and ground observations into a multi-objective calibration.

Our work complements and expands efforts on validating physically based hydrologic simulations through remote sensing products. The adjustment steps led to the improvements in simulated LST that are in line with studies using hydrologic models with various levels of sophistication. For instance, simulations of Xiang et al. (2017) in a semiarid basin in northern Mexico found LST RMSE of 4.3 ^∘C for the daytime and 1.9 ^∘C during the nighttime compared to MODIS. The hyperresolution (∼80 m) simulations of Ko et al. (2019) in the same basin resulted in bias of −1.4 ^∘C and CC of 0.87, and the high-resolution simulations with VIC in central Arizona by Wang et al. (2021) yielded LST biases between −1.5 and 3.6 ^∘C. To our knowledge, this study is the first to improve the simulated spatial patterns of hydrologic variables in the CRB using remote sensing products. By increasing the credibility of the spatial model outputs, this effort builds confidence in using regional hydrologic models for water resources predictions and decision-making under the ongoing megadrought in the Colorado River. Finally, we identified several future research avenues to further improve the fidelity of hydrologic models through the incorporation of remote sensing products. First, once the key parameters involved in the physical equations simulating a variable observed by satellite sensors have been identified as done here, a robust multiparameter sensitivity analysis could be conducted to investigate possible interactions among the parameters; this effort will help further refine the calibration. Second, automatic calibration strategies could be designed and applied to simultaneously target the simulation of multiple variables (here, LST and SCF).

Appendix A

We describe the solution of the energy balance in VIC, which leads to the computation of ground surface temperature (T_s) and canopy foliage temperature (T_f) used to compute the land surface temperature variable, LST_V, that is compared against the MODIS estimate, LST_M. We emphasize the main parameters and variables involved in the computation of these state variables. More detailed descriptions can be found in previous publications (Andreadis et al., 2009; Bohn and Vivoni, 2016; Cherkauer et al., 2003; Cherkauer and Lettenmaier, 1999; Liang and Lettenmaier, 1994). We first illustrate the original algorithm introduced in the first version of VIC (Liang and Lettenmaier, 1994), then the snow overstory scheme introduced by Cherkauer and Lettenmaier (2003), and finally the clumped-canopy scheme implemented by Bohn and Vivoni (2016).

A1 Original scheme from Liang and Lettenmaier (1994)

In Liang and Lettenmaier (1994), the minimal unit of simulation is the tile with a homogeneous land cover, i.e., the big-leaf approach. The energy balance equation for the tile can be expressed as follows:

$$\begin{array}{}\text{(A1)}& R_{\mathrm{n}}=\mathrm{LH}+\mathrm{SH}+\mathrm{GH},\end{array}$$

where R_n is net radiation, SH is sensible heat flux, LH is latent heat flux, and GH is ground heat flux. The parameters and variables involved in the computation of each term are summarized in Table 2. Net radiation is determined by the following:

$$\begin{array}{}\text{(A2)}& R_{\mathrm{n}}=(\mathrm{1}-\mathit{\alpha })\cdot R_{\mathrm{s}}+\mathit{\epsilon }\cdot \left(R_{\mathrm{L}}-\mathit{\sigma }\cdot T_{\mathrm{s}}^{\mathrm{4}}\right),\end{array}$$

where R_S and R_L are downward shortwave and longwave radiation, α is albedo, ε is surface emissivity (0.98 for water; 0.97 for other conditions), and σ is the Stefan–Boltzmann constant.

The latent heat is computed as follows:

$$\begin{array}{}\text{(A3)}& \mathrm{LH}={\mathit{\rho }}_{\mathrm{w}}\cdot {\mathit{\lambda }}_{\mathrm{v}}\cdot \left(E_{\mathrm{c}}+E_{\mathrm{t}}+E_{\mathrm{b}}\right),\end{array}$$

where ρ_w is the density of liquid water, λ_v is the latent heat of vaporization, E_c is evaporation from wet canopy, E_t is plant transpiration, and E_b is evaporation from surface soil moisture. For any given time, the maximum value of E_c, denoted as E_c,max, is calculated as follows:

$$\begin{array}{}\text{(A4)}& E_{\mathrm{c},\max }={\left({\displaystyle \frac{W}{W_{\max }}}\right)}^{\mathrm{2}/\mathrm{3}}\cdot E_{\mathrm{p}}\cdot \left({\displaystyle \frac{r_{\mathrm{a}}}{r_{\mathrm{a}}+r_{\mathrm{arc}}}}\right),\end{array}$$

where W is the amount of canopy interception at a given time, W_max is the maximum amount of water that the canopy can intercept (computed as 0.2 ⋅ LAI), r_arc is the canopy architectural resistance, r_a is the aerodynamic resistance, and E_p is the potential evaporation derived from the Penman–Monteith equation with a canopy resistance set to zero. This is calculated as follows:

$$\begin{array}{}\text{(A5)}& E_{\mathrm{p}}={\displaystyle \frac{\mathrm{\Delta }{R}_{\mathrm{n}}+{\mathit{\rho }}_{\mathrm{a}}\cdot c_{\mathrm{p}}\cdot \mathit{\delta }e\cdot \frac{\mathrm{1}}{r_{\mathrm{a}}}}{\left[\mathrm{\Delta }+\mathit{\gamma }\cdot \left(\mathrm{1}+\frac{r_{\mathrm{s}}}{r_{\mathrm{a}}}\right)\right]\cdot {\mathit{\lambda }}_{\mathrm{v}}}},\end{array}$$

where Δ is the slope of the saturation vapor pressure temperature relationship, ρ_a is the air density, c_p is the specific heat of air, δe is the vapor pressure deficit, γ is the psychrometric constant, and r_s is the surface resistance. The aerodynamic resistance is calculated as follows:

$$\begin{array}{}\text{(A6)}& r_{\mathrm{a}}={\displaystyle \frac{\mathrm{1}}{C_{\mathrm{w}}+u\left(z\right)}},\end{array}$$

where u(z) is the wind speed at the measurement height z, and C_w is the transfer coefficient for water. This is defined as follows:

$$\begin{array}{}\text{(A7)}& C_{\mathrm{w}}=\mathrm{1.351}\cdot {\displaystyle \frac{k^{\mathrm{2}}}{{\left[\ln \left(\frac{z}{z_{\mathrm{0}}}-\frac{d_{\mathrm{0}}}{z_{\mathrm{0}}}\right)\right]}^{\mathrm{2}}}}\cdot F\left(Ri\right),\end{array}$$

where k is the von Kármán constant, z₀ is the roughness length, d₀ is the displacement height, and F(Ri) is a function of the Richardson number, Ri, that accounts for atmospheric stability. z₀ and d₀ have different values for each vegetation type and for bare soil and snow. Ri is defined as follows:

$$\begin{array}{}\text{(A8)}& Ri={\displaystyle \frac{g\cdot \left(T_{\mathrm{air}}-T_{\mathrm{s}}\right)\cdot z}{\left(\frac{T_{\mathrm{air}}+T_{\mathrm{s}}}{\mathrm{2}}\right)\cdot u(z{)}^{\mathrm{2}}}},\end{array}$$

where g is the gravitational acceleration, and T_air is the air temperature. When $W\ge E_{\mathrm{c},\max }$, $E_{\mathrm{c}}=E_{\mathrm{c},\max }$; otherwise, E_c is a fraction of E_c,max determined as a function of precipitation and W.

The transpiration, E_t, is calculated as follows:

$$\begin{array}{}\text{(A9)}& E_{\mathrm{t}}=\left[\mathrm{1}-{\left({\displaystyle \frac{W}{W_{\max }}}\right)}^{\frac{\mathrm{2}}{\mathrm{3}}}\right]\cdot E_{\mathrm{p}}\cdot \left({\displaystyle \frac{r_{\mathrm{a}}}{r_{\mathrm{a}}+r_{\mathrm{arc}}+r_{\mathrm{c}}}}\right),\end{array}$$

where the canopy resistance, r_c, is related to the minimal stomatal resistance, r_min, via

$$\begin{array}{}\text{(A10)}& r_{\mathrm{c}}=r_{\min }\cdot {\displaystyle \frac{G_{\mathrm{sm}}}{\mathrm{LAI}}}.\end{array}$$

G_sm is the soil moisture stress factor depending on root zone water availability (depth dependent on vegetation type). Bare soil evaporation, E_b, is equal to E_p when the shallowest soil layer is saturated; otherwise, it is computed as follows:

$$\begin{array}{}\text{(A11)}& E_{\mathrm{b}}=E_{\mathrm{p}}\cdot \left[\underset{\mathrm{0}}{\overset{A_{\mathrm{s}}}{\int }}dA+\underset{A_{\mathrm{s}}}{\overset{\mathrm{1}}{\int }}{\displaystyle \frac{i_{\mathrm{0}}}{i_{\mathrm{m}}\left[\mathrm{1}-(\mathrm{1}-A{)}^{\mathrm{1}/b_i}\right]}}dA\right],\end{array}$$

where A_s is the fraction of saturated soil. This is computed as follows (Zhao et al., 1980):

$$\begin{array}{}\text{(A12)}& A_{\mathrm{s}}=\mathrm{1}-{\left(\mathrm{1}-{\displaystyle \frac{i_{\mathrm{0}}}{i_{\mathrm{m}}}}\right)}^{b_i},\end{array}$$

where b_i is the infiltration shape parameter, i₀ is the current infiltration capacity determined by water availability, and i_m is the maximum infiltration capacity computed as the product between maximum soil moisture (equal to soil depth times porosity) and (1+b_i).

The sensible heat flux, SH, is given by the following:

$$\begin{array}{}\text{(A13)}& \mathrm{SH}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{a}}\cdot c\cdot \left(T_{\mathrm{s}}-T_{\mathrm{air}}\right)}{r_{\mathrm{a}}}},\end{array}$$

where ρ_a and c are the mass density and specific heat of air at constant pressure, respectively.

The ground heat flux, GH, is calculated by the following:

$$\begin{array}{}\text{(A14)}& \mathrm{GH}={\displaystyle \frac{\mathit{\kappa }}{D_{\mathrm{1}}}}\left(T_{\mathrm{s}}-T_{\mathrm{1}}\right),\end{array}$$

where T₁ is the soil temperature at depth D₁ (0.1 m here), and κ is the soil thermal conductivity.

The equations described above are used to estimate T_s through an iterative procedure. T_s is initially set to T_air, leading to Ri=0 and F(Ri)=1; evapotranspiration is then estimated and the energy balance is solved to update T_s (Liang and Lettenmaier, 1994). Iterative solutions for T_s are repeated until the difference between initial and final values are within a tolerance. This scheme is applied to the case of tile A in Fig. 2 when $f_{\mathrm{v},A}=\mathrm{1}$.

A2 Snow overstory scheme from Cherkauer et al. (2003)

The energy balance in VIC was improved with the snow overstory scheme of Cherkauer et al. (2003). Andreadis et al. (2009) upgraded this scheme with fully balanced energy terms and the representation of snow interception. The scheme assumes a vegetated overstory (with foliage temperature T_f) and an understory without vegetation (with surface temperature T_s), as in tile B in Fig. 2 with $f_{\mathrm{v},B}=\mathrm{1}$. If snow is not present, T_f is assumed equal to T_air, and T_s is calculated with the scheme described above. When snow is present, the energy balance is solved separately in control volumes (CVs) of the overstory, understory, and the atmosphere surrounding the canopy (with temperature T_c), respectively. The algorithm involves the following steps:

1. T_c is initially assigned equal to T_air. The snow on the canopy is determined according to snowfall and maximum interception capacity, $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}\cdot L_{\mathrm{r}}\cdot \mathrm{LAI}$, where L_r is a step function of T_f from the last time step. If there is no snow on the trees, then $T_{\mathrm{f}}=T_{\mathrm{c}}=T_{\mathrm{air}}$. If there is snow on the trees, and snow is melting, then T_f=0 ^∘C. If the snow is not melting, then the energy balance of the overstory CV with snow is solved for T_f, as follows:

   $$\begin{array}{}\text{(A15)}& R_{\mathrm{n}}^{\mathrm{snow}\text{-}\mathrm{canopy}}+E_{\mathrm{A}}={\mathrm{SH}}^{\mathrm{snow}\text{-}\mathrm{canopy}}+{\mathrm{LH}}^{\mathrm{snow}\text{-}\mathrm{canopy}},\end{array}$$

   where E_A is energy advected by precipitation, SH^snow-canopy is calculated as in Eq. (A13) but with T_s, and T_air replaced by T_f and T_c. The net radiation for snow on the canopy is as follows:

   $$\begin{array}{}\text{(A16)}& \begin{array}{rl}{R}_{\mathrm{n}}^{\mathrm{snow}\text{-}\mathrm{canopy}}& =\left(\mathrm{1}-{\mathit{\alpha }}_{\mathrm{snow}}\right)\cdot R_{\mathrm{s}}+\mathit{\epsilon }\\ & \cdot \left(R_{\mathrm{L}}+\mathit{\sigma }\cdot T_{\mathrm{s}}^{\mathrm{4}}-\mathrm{2}\cdot \mathit{\sigma }\cdot T_{\mathrm{f}}^{\mathrm{4}}\right),\end{array}\end{array}$$

   with α_snow as the snow albedo. If T_s is not available, an initial value of 0 ^∘C is used in Eq. (A16). The latent heat from snow sublimation is as follows:

   $$\begin{array}{}\text{(A17)}& {\mathrm{LH}}^{\mathrm{snow}\text{-}\mathrm{canopy}}={\displaystyle \frac{\mathrm{0.622}\cdot {\mathit{\lambda }}_{\mathrm{s}}\cdot {\mathit{\rho }}_{\mathrm{a}}\cdot \mathit{\delta }e}{P_{\mathrm{a}}\cdot r_{\mathrm{a},\mathrm{snow}}}},\end{array}$$

   where λ_s is the latent heat of sublimation, P_a is atmospheric pressure, and r_a,snow is the aerodynamic resistance near the snow surface.

2. The energy balance is then applied to the understory CV. Due to the presence of a tall tree, the shortwave radiation reaching the ground surface is reduced (by 50 %) due to the shading effect. The incoming longwave radiation is computed only as a function of T_f, while the contribution from the atmosphere is assumed negligible. T_s is then calculated by solving the energy balance. In this case, sensible heat is calculated using Eq. (A13) by replacing T_air with T_c and computing the aerodynamic resistance as follows:

   $$\begin{array}{}\text{(A18)}& r_{\mathrm{a},\mathrm{snow}}={\displaystyle \frac{\ln {\left(\frac{z-d_{\mathrm{s}}}{z_{\mathrm{s}}}\right)}^{\mathrm{2}}}{k^{\mathrm{2}}\cdot u\left(z\right)}},\end{array}$$

   where z_s is snow surface roughness and d_s is the snow depth. If there is no liquid water in the ground snowpack, then the latent heat is calculated with Eq. (A17). If there is liquid water, then Eq. (A17) is used with the latent heat of vaporization, i.e., λ_s is replaced by λ_v.

3. Once T_s is derived, T_c is updated by solving the energy balance at the CV that includes the atmosphere surrounding the canopy, as follows:

   $$\begin{array}{}\text{(A19)}& {\mathrm{SH}}_{T_{\mathrm{air}},T_{\mathrm{c}}}={\mathrm{SH}}_{T_{\mathrm{c}},T_{\mathrm{s}}}+{\mathrm{SH}}_{T_{\mathrm{c}},T_{\mathrm{f}}},\end{array}$$

   where SH${}_{T_{\mathrm{c}},T_{\mathrm{s}}}$ is the sensible heat into snow calculated in step 2, and SH${}_{T_{\mathrm{c}},T_{\mathrm{f}}}$ is the SH^snow-canopy calculated in step 1. T_c is compared with its estimate from the previous step (T_air in the first iteration). If the values are not included within a tolerance, then steps 1–3 are repeated.

A3 Clumped-canopy scheme from Bohn and Vivoni (2016)

The schemes described above are based on the big-leaf approach, where vegetation was assumed to cover the entire surface of the tile. Bohn and Vivoni (2016) introduced the clumped-canopy scheme to improve the simulation of bare soil evaporation from inter-canopy spaces. This scheme relies on the vegetation fraction (f_v). The aerodynamic resistance of each tile is updated to be the inverse of aerodynamic conductance, $\mathrm{1}/g_{\mathrm{a}}$, with the following:

$$\begin{array}{}\text{(A20)}& g_{\mathrm{a}}=\left(\mathrm{1}-f_{\mathrm{v}}\right)\cdot \mathrm{1}/r_{\mathrm{a},\mathrm{s}}+f_{\mathrm{v}}\cdot \mathrm{1}/r_{\mathrm{a},\mathrm{v}},\end{array}$$

where r_a,s and r_a,v are aerodynamic resistances for bare soil and vegetated area, respectively, computed using Eq. (A6). For the soil, a constant roughness height of 0.0001 m is used.

Because of the introduction of f_v, we improved the snow physics in the Snow-adj step. The version of VIC employed in our baseline simulation assumed that longwave radiation into the snowpack was received only from the canopy in the tiles covered by trees, even for the unvegetated fraction. In the clumped scheme, where a fraction (1−f_v) is unvegetated, this assumption is not reliable. Therefore, we updated the computation of the longwave radiation as the weighted average of the canopy longwave and longwave from the atmosphere (LW${}_{\mathrm{d},\mathrm{v},B}\cdot (\mathrm{1}-f_{\mathrm{v},B})$ was replaced by LW${}_{\mathrm{d},B}\cdot (\mathrm{1}-f_{\mathrm{v},B})$, as highlighted in Fig. 2b).