A near-global, high resolution land surface parameter dataset for the variable infiltration capacity model

Hydrologic models predict the spatial and temporal distribution of water and energy at the land surface. Currently, parameter availability limits global-scale hydrologic modelling to very coarse resolution, hindering researchers from resolving fine-scale variability. With the aim of addressing this problem, we present a set of globally consistent soil and vegetation parameters for the Variable Infiltration Capacity (VIC) model at 1/16° resolution (approximately 6 km at the equator), with spatial coverage from 60°S to 85°N. Soil parameters derived from interpolated soil profiles and vegetation parameters estimated from space-based MODIS measurements have been compiled into input files for both the Classic and Image drivers of the VIC model, version 5. Geographical subsetting codes are provided, as well. Our dataset provides all necessary land surface parameters to run the VIC model at regional to global scale. We evaluate VICGlobal’s ability to simulate the water balance in the Upper Colorado River basin and 12 smaller basins in the CONUS, and their ability to simulate the radiation budget at six SURFRAD stations in the CONUS.

codes to subset the VICGlobal parameters to a particular domain. In addition to parameters, meteorological forcing data are required to run VIC. We do not include meteorological forcing data as part of VICGlobal. Instead, we direct readers to existing forcing datasets with near-global coverage, such as the reanalysis datasets MERRA-2 15 and GLDAS 16 , or real-time measurement-based datasets -see e.g. Xiao et al. 17 , Livneh et al. 12,18 , Bohn et al. 19 Finally, a note on the file format: The upgrade from VIC version 4 (VIC-4) to VIC version 5 (VIC-5) introduced two "drivers" for running the model. The Image driver takes NetCDF files as inputs, while the Classic driver takes ASCII text files. The VICGlobal parameter files are available in two formats: one for VIC-5 Classic, and one for VIC-5 Image.

Methods
This section describes how we used freely-available data to compile Classic driver input files for the VIC model. First, we created parameter files for VIC-5 Classic, then we converted them to NetCDF format for VIC-5 Image. VIC-5 Classic requires three parameter files: a soil parameter file, a vegetation parameter file, and a vegetation library file. An optional elevation band file can be provided to resolve sub-grid variability in elevation, which is important in regions with complex topography. The parameters are arranged as a relational database: each grid cell has a unique identifier, called a grid cell number, in the soil parameter file, that VIC uses to find the corresponding rows of data in the vegetation parameter and elevation band files. The Image driver uses a different setup, with all parameters stored in a single NetCDF file.
Soil parameters. The soil parameter file for VIC-5 Classic is an ASCII text file that includes soil parameters such as hydraulic conductivity and porosity, but also other kinds of static parameters, such as average precipitation and time zone offset from GMT. Each row of the soil parameter file represents one grid cell, and each column represents a different variable. We compiled the soil parameter file using MERIT 20 elevation data, soil texture data from the FAO HWSD, pedotransfer tables relating soil texture to other soil properties, and interpolated weather station data (WorldClim 21 ). Any remaining parameters were set to suggested values from the VIC model's documentation 2 . The following sections describe the estimation of each variable in the soil parameter file, summarized in Table 1. elevation and land mask. The VICGlobal soil parameter file uses the Multi-Error-Removed Improved-Terrain (MERIT 20 ) digital elevation model (DEM) to define the elevations, latitudes, and longitudes of each land grid cell. The MERIT DEM is an error-corrected and extended version of the SRTM DEM, with 3 arc-second resolution and coverage from 60°S to 85°N and 180°W to 180°E. Specifically, MERIT is a combination of the SRTM, AW3D, and Viewfinder Panoramas' DEMs, corrected for striping, speckle, absolute bias, and tree height bias. We used bilinear interpolation to aggregate MERIT to 1/16° resolution and derive a 1/16° MERIT-based land mask and DEM ( Figure S1). Soil texture data. Soil texture (percent sand, silt, and clay) and bulk density were obtained from the FAO HWSD, a gridded soil parameter dataset derived from in-situ measurements of the soil column. We used a 0.05° resolution NetCDF dataset converted from the original HWSD Microsoft Access database by Wieder et al. 22 . We resampled the HWSD soil data from 0.05° to 1/16° resolution using bilinear interpolation with the MATLAB ® function griddedInterpolant. While HWSD has near global coverage, there are missing data in some places around the world, notably Greenland and northern Africa. We filled in these missing data using inpainting, a gap-filling method from the field of image processing. We used the MATLAB ® function inpaintnans 23 , which uses a partial differential equation method to fill in missing data, to fill gaps in the HWSD data over the MERIT land mask. Figure 1 shows the HWSD bulk density data before and after inpainting.
The HWSD data are divided into "topsoil" and "subsoil" parameters. The first 30 cm of the soil column are considered topsoil and the lower 70 cm subsoil. VIC is typically run with three soil layers, so we created a three-layer soil parameter file by breaking up the 30 cm HWSD topsoil layer into two soil layers: one of 10 cm and one of 20 cm, so the final soil parameter file has three layers, with thicknesses of 10 cm, 20 cm and 70 cm, from top to bottom of the soil column. Ten centimeters has been a common choice for the uppermost layer soil depth in VIC modeling applications since its use by Liang et al. 24 . Soil layer depths are typically used as calibration parameters. VICGlobal values should be considered a starting estimate.
www.nature.com/scientificdata www.nature.com/scientificdata/ Calculating soil parameter values based on soil textures. Pedotransfer functions (e.g. Cosby et al. 25 ) relate readily available soil properties, such as soil texture, to less easily-observed properties, such as hydraulic conductivity. After resampling the HWSD data from 1/4° to 1/16° resolution, we estimated soil parameters by classifying each grid cell's USDA soil texture class and assigning physical soil properties based on a lookup table included with the VIC documentation 2,26 . The lookup table (Table 2) relates the 12 USDA soil texture classes to bulk density, field capacity, wilting point, porosity, saturated hydraulic conductivity, and slope of the soil water retention curve in Campbell's equation. We classified soil textures using the USDA soil texture triangle, as implemented by the MATLAB ® function soil_classification 27 . Figure 2 shows the derived USDA soil texture map.
We used these along with the lookup table to estimate saturated hydraulic conductivity (K sat ), the exponent in Campbell's equation for hydraulic conductivity (expt), fractional soil moisture at the critical point (wcr fract ), where the critical point is about 70% of field capacity, fractional soil moisture at the wilting point (wpwp fract ), quartz content, and porosity for each soil layer. The lookup table 26 did not include quartz content, so we supplemented it with the soil texture-quartz content lookup table from Peters-Lidard et al. 28 .
We set the variable infiltration capacity parameter b 0 2 infilt = . , the maximum baseflow fraction threshold d 0 001 s = . , and maximum soil moisture threshold w 0 9 s = . , their suggested values from the VIC documentation. These parameters, along with maximum baseflow velocity (dsmax) and soil depth, are typically calibrated. We set the baseflow curve exponent c = 2, the soil thermal damping depth dp = 4 m, soil density = 2685 kg/m 3 , surface roughness = 0.001 m, and snow roughness = 0.0005 m, also based on guidance from the VIC documentation. The soil moisture diffusion parameter phi s is not used in the current version of VIC, so we set it to the no-data value (−999). The final few soil parameters -dsmax, initial soil moisture (initm), and bubbling pressure (bubble)were calculated using the following equations, based on guidance from the VIC documentation.  www.nature.com/scientificdata www.nature.com/scientificdata/ content at the critical point. Equation (3) calculates bubbling pressure as a function of expt, based on linear regression of bubbling pressure vs. expt 30 . Figures S2-S9 in the Supplementary Information show maps of each soil parameter. We assumed residual soil moisture, the amount of soil moisture that cannot be removed from the soil by drainage or evapotranspiration, was zero. Elevation bands. VIC uses an elevation band file (also called a snow band file) to account for subgrid heterogeneity in grid cell elevations. The assumption of uniform elevation over an entire grid cell can lead to modeling errors in mountainous regions, where higher topography is associated with cooler temperatures and higher precipitation rates. The elevation band file accounts for subgrid variability in topography by dividing each grid cell into a number of elevation bands, each of which is simulated separately. VIC adjusts temperature, pressure, and precipitation depending on the elevation in each band. We prepared an elevation band file with five elevation bands by comparing the 1/16° DEM used for the soil parameter file with a 30 arc-second DEM. Both DEMs were derived by aggregating MERIT data. For simplicity, we assumed precipitation was evenly distributed among elevation bands within a grid cell. The elevation band file is provided with the caveat that using elevation bands requires more computing power; users may wish to turn elevation bands on or off (via the VIC global parameter file) depending on their needs. Vegetation parameters. VIC-5 Classic uses a vegetation parameter file to define the fractional cover of different vegetation types within each grid cell and some of their physical properties. Other vegetation parameters are stored in the "vegetation library" file. (VIC-5 Image simply stores all parameters in a single "parameter" file.) The VIC-5 Classic vegetation parameter file consists of information about fractional cover of each land cover type in each grid cell, and their corresponding root zone depths and root fractions within each root zone. The vegetation parameter file can optionally include time-varying LAI, fractional canopy cover, and albedo data, but it is simpler to specify these in the vegetation library (at the cost of not representing some spatial heterogeneity).
We used MODIS land cover data from the 0.05° MODIS MCD12C1 Collection 6 data product 31 to assign fractional land cover values to each grid cell by calculating the average land cover for MCD12C1 observations over the 2017 calendar year. We chose 2017 because it was the most recent year with data in all the MODIS-based datasets used for this study, and there is very low interannual variability of land cover 32 in MCD12C1 Collection 6. Figure 3 shows majority land cover types from the 2017 MCD12C1 observations. Like all global land cover data products, MCD12C1 makes classification errors. Sulla-Menashe et al. 32 reported 67% overall IGBP classification accuracy for 2001 land cover. Classification errors are more common in the "mixed" land covers, such as cropland/natural vegetation mosaic, shrublands, grasslands, and savannas. Fortunately for our purposes, the vegetation parameters for commonly-confused land covers tend to be fairly similar themselves, which reduces the impact of misclassification on land surface modelling results. For example, the LAI of open shrubland is not too different from the LAI of closed shrubland. www.nature.com/scientificdata www.nature.com/scientificdata/ We calculated root fraction as a function of land cover class following the method of Zeng 33 , who defined the following formula (Eq. 4) for use in parameterizing land surface models: where Y = cumulative root fraction, d = depth, and a and b are empirical parameters defined by Zeng 33 for each International Geosphere-Biosphere Programme (IGBP) land cover type, based on a rooting depth database compiled from more than 200 field surveys. We used this formula with depths of 0.1 m, 0.7 m, and dr, corresponding to three root zones. The value of dr, the maximum rooting depth for each IGBP land cover type, was taken from Zeng 33 . This method assumes that the depth and distribution of roots depends only on the land cover type; we assume that land cover type is the primary control on root characteristics. Table 3 shows root fractions and root zone depths for each IGBP land cover type. Like previous large-scale VIC vegetation cover datasets, our vegetation parameter file neglects land cover change over time. However, it does have a few other advantages over past vegetation parameter datasets. The land cover classification used in the N2001 and L2013 VIC parameter sets is referred to as "UMD-NLDAS" because it is a modified version of the AVHRR-based University of Maryland (UMD) land cover product 34 . The UMD-NLDAS classification was modified for the North American Land Data Assimilation project (NLDAS 35 ) to exclude open water, urban, and snow and ice land cover classes (see BV2019). VICGlobal uses 17 IGBP land cover classes, including urban, barren, perennial snow and ice, and inland water bodies, permitting better description of land cover variability than the 11 UMD-NLDAS classification.
Vegetation library file. The vegetation library maps each land cover type to a set of vegetation parameters (Table 4). We adapted the LDAS vegetation library 36 for use with the 17 IGBP land cover classes, taking monthly average LAI, fractional canopy cover (fcanopy), and albedo values obtained from recent MODIS data products. We set architectural resistance (r 0 ) and minimum stomatal resistance (r min ) to values from literature (described below). The rest of the parameters, which are described in the N2001 paper, were left to their original LDAS vegetation library values. This section describes how we estimated LAI, fcanopy, albedo, r 0 , and r min , and how we transferred the remaining parameters from the 11 UMD-NLDAS land cover classes to the 17 IGBP land cover classes.
We used MODIS observations from the year 2017 to calculate monthly average LAI, fcanopy, and albedo for each IGBP land cover type. We calculated LAI and albedo from the MODIS-based Global LAnd Surface Satellite dataset (GLASS [37][38][39] ) and fcanopy from NDVI observations (MCD13C1 40 ) The expression used for fcanopy follows BV2019: where NDVI min and NDVI max are the minimum and maximum values of NDVI observed for that month. Monthly LAI, fcanopy, and albedo values were calculated by averaging over all grid cells of the same land cover type, counting only cells that were at least 90% homogenous, to avoid noise from grid cells with multiple land covers. Excepting perennial snow and ice land cover, the vegetation parameters in the VIC vegetation library should describe snow-free vegetation. Therefore, before calculating LAI, fcanopy, and albedo for each land cover class, we used fractional snow cover data from MOD10CM 41 , a global 0.05 degree monthly snow cover dataset, to exclude grid cells with more than 90% snow cover. Additionally, we set albedo to 0.05 for open water, and we set LAI and fcanopy to 0 for open water and perennial snow and ice. The resistances r min and r 0 play a role in determining how much plant transpiration occurs. Higher resistance means less transpiration. Stomatal resistance is resistance to the release of water through the plant stomata, and For the other parameters in the vegetation library file (displacement height, roughness length, etc.), we assigned values using the existing LDAS vegetation library. Since there are 17 IGBP land cover classes, and only 11 UMD-NLDAS land cover classes in the LDAS vegetation library, we re-assigned some IGBP land cover classes to take the parameters of UMD-NLDAS land cover classes. We remapped barren land, permanent wetlands, snow and ice, urban land, and water bodies to take the parameters of "grasslands" from the LDAS vegetation parameter file. While the characteristics of the barren, snow and ice, urban, and water land cover types clearly differ from those of grasslands, their low LAI and fcanopy values, corresponding to sparse vegetation, essentially "turns off " the other vegetation parameters in the VIC model, as pointed out by BV2019. The other remappings were more straightforward. Croplands and croplands/natural vegetation mosaics inherited values from "croplands, " savannas became "wooded grasslands, " and woody savannas became "woodlands. " We were thus able to assign vegetation parameter values to the each of the 17 IGBP land cover classes.
To calculate global average time series of seasonally-varying vegetation parameters would be of limited interest as the seasonal cycle would average out across the equator. Therefore, we calculated average monthly fcanopy, LAI, and albedo for each vegetation type in each hemisphere, and we developed two separate vegetation library   www.nature.com/scientificdata www.nature.com/scientificdata/ files: one for the northern hemisphere and one for the southern hemisphere. Maps of January and July LAI, fcanopy, and albedo are shown in Fig. 4. For illustrative purposes, the parameter values in this figure have been averaged over the 17 IGBP land cover classes using area-based weighting. Figures S14-S19 show maps of the remaining vegetation parameters. Figures S20-S22 show the cycle of LAI, fractional canopy cover, and albedo for each vegetation type, averaged separately over each hemisphere.

Data Records
Soil and vegetation parameters for the VIC model are available for download at Zenodo 48 in NetCDF format for version 5 of the VIC model. The files are stored as zip archives. parameters_classic.zip contains ASCII text files with soil parameters, vegetation parameters, elevation bands, and two vegetation library files -one of the northern hemisphere and one for the southern hemisphere -for VIC-5 Classic. parameters_global.zip contains a NetCDF "parameter" file with all the soil and vegetation parameters described above and a NetCDF "domain" file describing the VICGlobal domain (all land mass between 60°S and 85°N) for VIC-5 Image. MATLAB ® codes for subsetting either set of parameters from the entire VICGlobal extent to a subregion of interest. Additionally, parameter and domain files pre-subsetted to North America, South America, Africa, Eurasia, and Oceania are available for download.

technical Validation
Streamflow and snow-water equivalent in the Upper Colorado Basin. Having created input files for the VIC model, we tested the parameters in a large, well-studied river basin. We used the VICGlobal parameters to run VIC in water balance mode over the Upper Colorado River Basin (UCRB), a 293,600 km 2 basin in the western United States. We ran VIC once using the VICGlobal parameters and once using the L2013 parameters. Both simulations used the meteorological forcing data from L2013, at a six-hourly timestep, in water balance mode, for the 6-year period from Oct. 1, 2005 to Sept. 30, 2011.
We compared estimated streamflow from the VICGlobal and L2013 simulations with naturalized streamflow estimates from the U.S. Bureau of Reclamation 49 (USBR) at Lees Ferry, Arizona (Fig. 5). Naturalized flow is measured streamflow adjusted for the effects of reservoir storage and management and consumptive uses such as irrigation. We compared our VIC model outputs with naturalized streamflow because our VIC implementation does not simulate consumptive water use or reservoir storage. Due to differences in soil and vegetation parameters between the two sets of input files, there are notable differences in the hydrographs from each simulation. Relative to the L2013 results, the uncalibrated VICGlobal peak flows' timing is too early and their magnitude is too high. This is expected given that the L2013 parameters have been calibrated to get a good match to gauge data.
To understand the cause of this mismatch, we examine seven commonly calibrated soil parameters, which are difficult or impossible to estimate from measurements: ds, ws, dsmax, b infilt , and the thicknesses of each soil layer (t 1 , t 2 , t 3 ). Taking a closer look at these soil parameters in the UCRB (Fig. 6), we see that the infiltration capacity parameter b infilt is considerably higher in the VICGlobal parameters than it is in the L2013 parameters, which www.nature.com/scientificdata www.nature.com/scientificdata/ would tend to cause higher runoff rates. ds is lower for VICGlobal than for L2013, so nonlinear baseflow occurs at a lower fraction of dsmax, tending to make baseflow peaks occur earlier. dsmax is considerably higher for VICGlobal than for L2013 in much of the UCRB, so the maximum baseflow rate is higher for VICGlobal. Finally, the thicker soil layers in L2013 mean that more water can infiltrate into the soil before baseflow occurs. Table 5 describes each of the seven calibration parameters and their influence on VIC model outputs. VIC users seeking more guidance on calibrating soil parameters should consult the VIC model documentation and relevant literature [50][51][52] . We calibrated the VICGlobal parameters b infilt , dsmax, and t 3 , the same parameters calibrated by L2013, to get a good match between predicted and observed (USBR naturalized) streamflow. However, decreasing b infilt on its own was not enough to reduce the high runoff estimates produced by the VICGlobal parameters (VIC's sensitivity to b infilt depends on the water-holding capacity of the upper two soil layers). By introducing the thickness of the second soil layer t 2 as a fourth calibration parameter, we were able to reduce runoff and increase transpiration. The final set of calibrated parameters was b 0 038 infilt = . , dsmax = 0.60 mm⁄day, t 2 = 1.5m, and t 3 = 1.6m. For this analysis, we used manual calibration because the model run time made automated methods, which require hundreds to thousands of model runs, impractical. We used a custom MATLAB ® application -a graphical user interface for running the VIC model, tuning its parameters, and displaying its outputs -to assist with manual calibration. We assumed the calibrated parameters were uniform over the basin. In addition to calibration by trial and error, VICGlobal users may also wish to explore automated calibration methods such as the Shuffled Complex Evolutionary algorithm 53 (SCE-UA) or Dynamically-Dimensioned Search 54 (DDS) when practical.
The calibrated VICGlobal simulation outperformed the L2013 simulation, with a Kling-Gupta efficiency 55 (KGE) of 0.24, compared to −0.26 for L2013 and −1.7 for the uncalibrated VICGlobal simulation. We also compared simulated snow-water equivalent (SWE) between the VICGlobal and L2013 simulations. Spatial patterns of simulated SWE were consistent between the two simulations ( Fig. 7a-c), as were patterns of snow accumulation and melt (Fig. 7d). VICGlobal SWE was 3 mm higher than L2013 SWE, on average. Snow sublimation, including canopy sublimation, was higher for L2013 than for VICGlobal, which helps explain the slight overestimate of SWE by VICGlobal relative to L2013. Overall, VICGlobal is able reproduce the timing, magnitude, and spatial pattern of L2013-simulated SWE in the UCRB, with no need for parameter calibration.  www.nature.com/scientificdata www.nature.com/scientificdata/ Water balance in 12 unregulated CONUS basins. Beyond the Upper Colorado Basin, we evaluated the VICGlobal parameters' potential for modelling the water and energy balance in 12 basins, ranging from 1500-25000 km 2 , chosen for good spatial coverage of the CONUS. Modelled discharge was compared with monthly observations at USGS reference stream gages 56 at each basin outlet; we used the DDS method to calibrate b infilt , dsmax, t 2 , and t 3 . (These basins are small enough for automatic calibration to be practical.) The calibration was performed with L2013 meteorological forcing data, for calendar year 1993, with the VIC model run from 1990 -1992 as spin-up. After 500 model evaluations (50 for the Clearwater River), the average discharge calibration KGE was 0.47, with a maximum of 0.83 for the Clearwater River in Idaho and a minimum of 0.01 for the White River in Arkansas. Using the calibrated parameters, we performed a validation run from 1994-2011. Table 6 shows goodness of fit between modelled and measured discharge for the 18-year validation run. Figure 8 shows discharge plots for each basin over the validation period. Good matches can be seen for the Clearwater, New, Homochitto, Mattawamkeag, Gasconade, Trinity, and Little Fork rivers, while the White, Sheyenne, Brazos, and San Simon rivers did not respond well to calibration, suggesting that parameters other than the four calibrated here are to blame. See e.g. Demaria et al. 50 for more insight on VIC calibration.  www.nature.com/scientificdata www.nature.com/scientificdata/

Surface radiation budget validation with SURFRAD. To evaluate how well the (uncalibrated)
VICGlobal parameters simulate the surface radiation balance, we ran VIC over six SURFRAD 57 sites, using soil parameters from the 1/16° grid cell containing the SURFRAD sites. We ran hourly simulations in energy balance mode from 1995-2011, with 1994 as a spin-up year. Meteorological inputs taken from meteorological stations at the sites provided input data for the model, except for precipitation, which we took from L2013.
The first row of Fig. 9 shows modelled and observed net radiation, upwelling longwave, downwelling longwave, upwelling shortwave, and downwelling shortwave radiation averaged over each day from 1995-2011 for six SURFRAD sites in the CONUS. Downwelling shortwave and longwave radiation predictions similar to the observations because the SURFRAD data were used as inputs for the VIC model (but not identical because precipitation inputs were taken from L2013 due to lack of ground measurements at the sites). There is a positive bias for net radiation resulting from a slight low bias for upwelling shortwave radiation. Overall, the bias is small.
The second and third rows of Fig. 9 show scatterplots of predicted vs. observed upwelling shortwave and longwave radiation, with one-to-one lines shown in black. The correlation between predicted and observed upwelling shortwave and longwave radiation is close to 1 (ranges from 0.96-0.99 for all sites). Running VIC with VICGlobal parameters allows simulation of upwelling longwave radiation with an RMSE of 25 W/m 2 and RMSE of 15 W/ m 2 for upwelling shortwave radiation. In Fig. 9, we have excluded data at times when snow covers the ground to address the scale-issue -the spatial scale of the VIC simulations (a 1/16° grid cell) is much larger than that of the SURFRAD measurement -because of snow's large role in determining upwelling solar radiation, we excluded   www.nature.com/scientificdata www.nature.com/scientificdata/ times when either the VIC model or SURFRAD measurements had snow on the ground using an albedo threshold of 0.4; none of the VICGlobal albedos for non-snowy land surfaces are this large.

Usage Notes
We have described VICGlobal, a globally-consistent 1/16° VIC parameter dataset with soil and vegetation parameters derived from the latest satellite-based remote sensing datasets (MODIS and MERIT, which is based on SRTM data) and in-situ soil data from the FAO HWSD. In addition to its higher resolution, VICGlobal has an advantage over previous global VIC setups due to its inclusion of seasonally-varying fractional canopy cover, LAI, and albedo, and because it explicitly accounts for barren, wetland, open water, and perennial snow and ice land covers. VICGlobal is provided in geographic coordinates, referenced to the WGS84 ellipsoid and datum.
VICGlobal has a few limitations. Its parameters are uncalibrated, so users must calibrate sensitive yet hard-to-measure parameters such as soil depth and the variable infiltration capacity parameter to get a good match between simulated and observed discharge. Several of the vegetation parameters, such as roughness length and displacement height, are assumed constant in time, even though realistically these parameters change as vegetation blooms and senesces throughout the year. And while we believe our monthly, hemisphere-average fractional canopy cover, LAI, and albedo are a major improvement over past global datasets, the most realistic parameter set would have them vary from grid cell to grid cell, even for the same vegetation type. Despite its limitations, we hope that VICGlobal, with its relatively high spatial resolution, wide coverage, and easy availability will be a valuable resource for VIC users.

Code availability
The scripts used to create the VICGlobal data set can be found on the corresponding author's Github page (https://github.com/jschap1/vicglobal-prep and https://github.com/jschap1/vegpar). The VICGlobal parameters were subset to the Upper Colorado River Basin using the subsetting codes included with the VICGlobal dataset, archived on Zenodo 48 .

Fig. 9
First row: Hourly modelled and observed net radiation, upwelling longwave, downwelling longwave, upwelling shortwave, and downwelling shortwave radiation from 1995-2011 for six SURFRAD sites in the CONUS. Dashed lines show predictions, while solid lines show observations. Second and third rows: scatterplots of predicted vs. observed upwelling shortwave and longwave radiation, with one-to-one lines shown in black. Key: Black = net radiation, red = downwelling shortwave, cyan = upwelling shortwave, green = downwelling longwave, blue = upwelling longwave.