The Reynolds Creek Long-Term Dataset: A long-term meteorological dataset derived from station observations in the Reynolds Creek Experimental Watershed

AbstractA forty-year gridded meteorological forcing dataset spanning the water years 1984 to 2023 (October 1st to September 30th) has been compiled for the Reynolds Creek Experimental Watershed (RCEW) in southwest Idaho, USA. This Reynold Creek Long-Term (RCLT) dataset consists of hourly, 10-meter resolution grids of air temperature, vapor pressure, precipitation mass and phase, incoming shortwave and longwave radiation, and wind speed and direction. These variables are foundational for many ecological and hydrological Land Surface Models (LSMs) used in research and operational applications and have been interpolated and calculated from hourly measurements from the dense meteorological station network within the RCEW. The elevation relief of the 240 km2 RCEW spans the historical winter rain-to-snow transition, so a research application example is presented in which we show how the snow-dominated area of the basin has evolved over the forty-year data record. This dataset, stored in cloud-optimized Zarr format, enables future model development, benchmarking, and uncertainty analyses of existing models, independent validation of gridded atmospheric reanalysis datasets, and novel investigations of hydroclimatic variability across snow-dominated semi-arid environments.
1 Site DescriptionIn 1960, the United States Congress allocated funding for an experimental research watershed to advance hydrologic research in western U.S. rangelands. Since then, the USDA Agricultural Research Service (ARS) Northwest Watershed Research Center (NWRC) has managed the scientific infrastructure in the Reynolds Creek Experimental Watershed (RCEW) in Southwest Idaho (43.205°, -116.75°). At approximately 240 km2, the RCEW is characteristic of semiarid snow-dominated environments found throughout the Great Basin and spans an elevation gradient of 1100 to 2244 meters above sea level. Characteristic of mountain climatology, a significant elevational and directional precipitation gradient exists due to the prevailing northeast-trending storms during the winter and spring, when most of the annual precipitation occurs (Hanson, 2001). Since 1984, average annual precipitation ranged from 228 mm at site RC.057 site in the low northeast elevations to 1,086 mm at site RC.163 in the highest southwest elevations of the watershed. The NWRC has a long history of publishing station-based hydrometeorological datasets, including Slaughter et al. (2001), Reba et al. (2011) and Godsey et al. (2018). This dataset updates the gridded temperature, humidity, and precipitation dataset reported by Kormos et al. (2018), includes additional wind and radiation data, and appends nine more recent years of conditions (2015-2023) that have experienced a wider range of weather variability (Monteiro & Morin, 2023).
2 Instrumentation and Variable DistributionThe Spatial Modeling for Resources Framework (SMRF; Havens et al. (2017)) was employed to distribute each of the ten land-surface meteorological variables to a 10-meter regular grid. Each forcing variable has either been empirically derived or directly interpolated from hourly station measurements across the catchment domain. Owing to the considerable length of time encompassed by the dataset, many different sensors have been deployed in the watershed over the forty-year data record with differing levels of accuracy, which are not reported here. Invalid data were preliminarily removed for all measured variables besides precipitation (which necessitated a unique approach detailed in subsection 2.6), and temporal interpolation was performed for data gaps of two hours or less. Gaps lasting longer than two hours were left empty with the foreknowledge that spatial interpolation from nearby sites in the high-density network would act as data surrogates. It is worth noting that the NWRC has always employed a full-time staff of technicians tasked with the calibration and servicing of each sensor deployed in the RCEW. The gridded interpolation methods for all ten modeled variables from the six measured variables are described in Hedrick et al. (2018) and further elaborated upon in the following subsections.
3.1 Air TemperatureHourly measurements of air temperature (Ta) were made at 35 individual sites over the 40-year period (Figure 2), with the number of sites significantly increasing after water year 2000. Measurements are currently made using various incarnations of the widely used Vaisala HMP series of temperature and humidity sensors with ventilated radiation shields.
A modified inverse distance weighting (IDW) approach was used to distribute Ta across the 10-meter grid. In this process, the elevational trend is calculated at each time step, constrained to be negative due to the general relationship between elevation and Ta, then subtracted from the station measurements to produce a temperature residual. These residuals are distributed using standard IDW and added to each grid cell’s position on the elevation gradient slope line. This approach for distributing Ta improves station representativeness for areas with complex local topography.
2.2 Vapor PressureGridded actual vapor pressure (ea) values were interpolated from measurements of relative humidity (RH) at the same 35 sites as Ta over the 40-year data record, using the same Vaisala HMP instruments referred to in subsection 2.1. The empirical Tetens equation was used for deriving ea from RH and Ta:
ea = RH × 0.6108 × e ^ (17.3 × Ta / 237.3 + Ta). (1)
Station-derived ea was then distributed to the 10-meter grid using the same modified IDW detrending approach described in subsection 2.1. Dew point temperature was also calculated for the wet bulb temperature calculation (subsection 2.6.2 below) but was not stored in the RCLT dataset because it can be calculated from the vapor pressure and air temperature using the existing empirical relationships (e.g., the Clausius-Clapeyron equation).
2.3 WindWind speed (us) and direction (udir) were measured at a total of 29 sites over the data record, though only three sites were available prior to 1994 and four sites prior to 2002. The sparseness of wind measurements through the early years of this dataset is likely a source of uncertainty in the distributed wind grids for that period, though we should note that the pre-2002 measurements captured the full elevation gradient in RCEW at low (RC.076), mid (RC.127), and high elevation sites (RC.176).
Station measurements of wind were distributed to the 10-meter grid using the maximum upwind slope (maxus) terrain parameter described in Winstral, et al. (2002) and Winstral et al. (2009). In short, the underlying digital elevation model (DEM) is used to calculate a maxus value (in degrees) over a user-defined upwind distance (here 300 meters) for all possible upwind directions (0˚ to 360˚) in 5˚ increments. The resulting 72 layers of maxus grids are stored in a lookup library. Then, for each station the measured wind speed is adjusted to simulate what the wind speed would have been on a flat surface (‘flatwind’) using the maxus value for the measured wind direction at the site. Once the adjusted ‘flatwind’ speeds and wind direction components have been distributed across the entire grid using standard IDW, the distributed wind directions are used to find the maxus value for each grid cell and the distributed ‘flatwind’ speeds are converted back to actual wind speeds.
For the gridded dataset, wind speed and direction were converted into U- and V-components to match the conventions of NWP models such as the WRF and HRRR models. The U-component represents the East-West wind speed, with positive values indicating wind out of the west, while the V-component represents the North-South wind speed, positive value indicating wind out of the south.
2.4 Shortwave RadiationGridded incoming shortwave radiation (Sin) has been measured at 23 sites across the RCEW but cannot be directly spatially interpolated from measurements due to the complex terrain and the variable vegetation canopy present across the catchment. Instead, a three-step process produced the hourly gridded Sin product.
Hourly clear sky atmospheric Sin was modeled, then corrected for surrounding terrain in each 10-meter grid cell following Dozier (1980) and Dubayah (1994), respectively.Station measurements of Sin were divided by clear sky radiation to derive a cloud factor (Cfac) at each station pixel (Cfac =1 represents cloud-free conditions), which was then distributed using standard IDW across the domain.Canopy-corrected Sin values were estimated using empirical relationships presented by Link and Marks (1999), where direct beam shortwave radiation under canopy (Rb) can be represented by:Rb = Sb,in × e ^ (-μh / cos(⁡θ)).
In this equation, Sb,in is the above canopy cloud corrected direct beam radiation, μ is a canopy extinction coefficient, h is the height of the canopy, and θ is the solar zenith angle. Diffuse shortwave radiation under canopy (Rd) is computed by adjusting the above canopy cloud corrected diffuse radiation (Sd,in) by the canopy optical transmissivity (τ):
Rd= τ × Sd,in.
The terrain-, cloud-, and canopy-corrected Sin presented in this dataset is then the sum of Rb and Rd.
When a comparison was performed between modeled and measured Sin across the 40-year record, we discovered that the values for τ and μ presented in Link and Marks (1999), which were derived in the Canadian Boreal forests, led to an overestimation of Sin at open sites by up to 40%, and an underestimation of Sin at forested sites by 30% or more. Therefore, a manual adjustment of the τ, μ, and height parameters was performed to produce Sin values that more closely matched the station observations. Importantly, the scale difference between point measurements of Sin and averaged Sin across a 10-meter by 10-meter grid cell precludes direct agreement since shortwave radiation at the ground surface varies over very short length scales. However, the general trends in modeled Sin magnitude as a function of cloud cover were well-represented in the spatial dataset.
Calculating net shortwave radiation from the RCLT incoming solar product requires an estimate of land surface reflectance, or albedo, but there are many ways to derive an albedo product from both models and remote sensing products. Users of the RCLT dataset can use their own methods, but for simplicity we also include here visible and infrared bands of snow albedo for when snow is present. These albedo estimates use a time decay approach to capture albedo change as a function of springtime snow metamorphism, terrain factors, and solar zenith angle (Marshall and Warren, 1987). Importantly, these albedo estimates do not apply for snow-free conditions and should be masked for applications involving land surface models.
3.5 Longwave RadiationDue to a lack of continuous upward looking pyrgeometers, or sensors measuring incoming longwave radiation (Lin), across the dataset time domain, incoming longwave radiation was not directly interpolated from measurements. Instead, clear sky Lin was empirically estimated from distributed air temperature and humidity using the methods of Brutsaert (1975), then corrected for surrounding terrain using equations presented in Marks and Dozier (1979). Next, cloud-corrected Lin was computed using the empirical relationship described by Garen and Marks (2005). The final step was to adjust the cloud-corrected Lin using vegetation maps derived from the LANDFIRE 2016 dataset (LANDFIRE, 2016) and empirically derived transmissivity (τ) reported by Link and Marks (1999) then adjusted as described in subsection 2.5.
2.6 PrecipitationTo satisfy the standard input requirements of an energy and mass balance snow model, The RCLT dataset contains four distinct variables related to precipitation in the basin. These variables of precipitation mass, temperature of the falling hydrometeor, initial density of newly fallen snow, and the snow proportion of precipitation are described in the following subsections.
2.6.1 Precipitation MassDespite being the foundation upon which hydrologic models rely, precipitation measurements are often the largest source of predictive hydrologic uncertainty (Bárdossy et al., 2022). Across the wide spectrum of snow-dominated watersheds in the Western U.S., the majority of in situ measurements are made with weighing buckets fitted with alter shields that reduce wind speeds above the bucket orifice and thus increase the gauge catch efficiency (CE), or the ratio of measured precipitation to a “true” value (Thériault et al., 2021). However, many sites lack co-located wind speed measurements for applying the necessary World Meteorological Organization (WMO) transfer functions for undercatch correction (Kochendorfer et al., 2018), which can lead to low biases in regional and basin estimates of precipitation. Additionally, precipitation exhibits high spatial heterogeneity in complex terrain, which cannot be captured by a single measurement site in a large mountain basin.
To overcome the issue of spatial representativeness, the RCEW measurement network was initially planned to contain one gauge for every square mile of the watershed (n=110). By the beginning of this dataset in water year 1984, the number of sites had been reduced to the 25 stations (~one measurement per 10 km2) used here to produce the hourly gridded precipitation fields.
To address the undercatch issue, a unique dual-gauge approach (Hamon, 1970) was used for all but one of the sites in the RCEW. This method requires two co-located Belfort-type weighing buckets with one existing unshielded and the other fitted with a single-alter shield. The combination of the shielded and unshielded measurements allows an empirical extrapolation of more accurate ‘actual’ precipitation data compared with single shielded gauges employing a WMO transfer function (Hanson et al., 2004). Site RC.124B is the only single shielded gauge in the basin and was here corrected using the WMO transfer function. Hourly measurements of precipitation mass were distributed across the 10-meter grid using a Detrended Kriging interpolation method (Garen, 1995) identical to the approach by Kormos et al. (2018).
3.6.2 Precipitation TemperatureThe precipitation temperature variable is represented by the hourly computed ice or wet bulb temperature calculated with a widely used Newton-Raphson iterative solution to the psychrometric equation (Campbell and Norman, 1998). This approach requires air temperature (subsection 2.1), dew point temperature (from calculated vapor pressure in subsection 2.2) and estimated atmospheric pressure from elevation. Prior work has demonstrated that wet bulb temperature is the most suitable method for partitioning snow from rain in a semiarid watershed such as the RCEW (Marks et al., 2013).
3.6.3 Density of New SnowHourly estimated of new accumulated snow density is included in the RCLT dataset for energy and mass balance snow models that may require it. For this long-term application, we computed new snow density from a lookup table based on previous work (Susong et al., 1999) using the calculated precipitation temperature (see subsection 2.6.2) and precipitation mass (see subsection 2.6.1) in each grid cell.
3.6.4 Snow Fraction of PrecipitationIn addition to new snow density, the lookup table from Susong et al. (1999) was also used to determine the amount of precipitation that fell as snow across the catchment. When time step precipitation temperatures fell between -0.5°C and +0.5°C, the precipitation was defined as mixed phase, while colder temperatures resulted in 100% snowfall, and warmer temperatures were 100% rainfall.