Observations of surface energy fluxes and meteorology in the seasonally snow-covered high-elevation East River watershed during SPLASH, 2021–2023

From autumn 2021 through summer 2023, scientists from the National Oceanic and Atmospheric Administration (NOAA) and partners conducted the Study of Precipitation, the Lower Atmosphere, and Surface for Hydrometeorology (SPLASH) campaign in the East River watershed of Colorado. One objective of SPLASH was to observe the transfer of energy between the atmosphere and the surface, which was done at several locations. Two remote sites were chosen that did not have access to power utilities. These were along the valley floor near the East River in the vicinity of the unincorporated town of Gothic, Colorado. Energy balance measurements were made at these locations using autonomous, single-level flux towers referred to as atmospheric surface flux stations (ASFSs). The ASFSs were deployed on 28 September 2021 at the Kettle Ponds Annex site and on 12 October 2021 at the Avery Picnic site and operated until 19 July and 21 June 2023, respectively. Measurements included basic meteorology; upward and downward longwave and shortwave radiative fluxes and subsurface conductive flux, each at 1 min resolution; 3-D winds from a sonic anemometer and ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}/{\mathrm{CO}}_{\mathrm{2}}$ from an open-path gas analyzer, both at 20 Hz from which sensible, latent heat, and CO₂ fluxes were derived; and profiles of soil properties in the upper 0.5 m (both sites) and temperature profiles through the snow (at Avery Picnic), each reported between 10 min and 6 h. The system uptime was 97 % (Kettle Ponds) and 90 % (Avery Picnic), and collectively 1184 d of data was obtained between the stations. The purpose of this article is to document the ASFS deployment at SPLASH, to document the data acquisition and post-processing of measurements, and to serve as a guide for interested users of the data sets, which are archived at Zenodo (https://doi.org/10.5281/zenodo.10313363, Cox et al., 2023b; https://doi.org/10.5281/zenodo.10327409, Cox et al., 2023c; https://doi.org/10.5281/zenodo.10313894, Cox et al., 2023d; https://doi.org/10.5281/zenodo.10307825, Cox et al., 2023e; https://doi.org/10.5281/zenodo.10310520, Cox et al., 2023f) with the Creative Commons Attribution 4.0 International license.

1 Introduction

The Upper Colorado River Basin (UCRB), which spans parts of Colorado, Utah, Wyoming, Arizona, and New Mexico, is a snow-dominated watershed. The water stored there in the mountain snowpack supplies reservoirs used to support communities and industry throughout the southwest United States. In turn, management of these reservoirs of water depends on the fidelity of measurements and predictions of hydrological inputs and timing of discharge within the riverine system. To support advancement of the tools used for these predictive services, the National Oceanic and Atmospheric Administration (NOAA) supported a campaign to make detailed observations pertaining to snowpack, soil, and precipitation in the UCRB from 2021 to 2023, termed the Study of Precipitation, the Lower Atmosphere, and Surface for Hydrometeorology (SPLASH) (de Boer et al., 2023). SPLASH was carried out in the East River watershed, an approximately 300 km² high-altitude area along the western slope of Colorado's Rocky Mountains near Crested Butte. The valley is additionally the focus of a U.S. Department of Energy (DoE) Watershed Function Project Science Focus Area (Hubbard et al., 2018) managed by the Lawrence Berkeley National Laboratory (LBNL) and is home to the Rocky Mountain Biological Laboratory (RMBL).

The goals of the 2-year SPLASH campaign were to collect data suitable for analyzing subgrid-scale variability and physical process representation in models such as the National Water Model (NWM) and the High-Resolution Rapid Refresh (HRRR) forecast model. Of particular interest was the seasonal water cycle in complex terrain from the lower atmosphere through the upper soil. SPLASH included observations made by a network of radars, aircraft, and in situ stations and was coordinated jointly with other complementary programs, including the Surface Atmosphere Integrated field Laboratory (SAIL), which was operated by the DoE Atmospheric Radiation Measurement (ARM) program (Feldman et al., 2023) and the Sublimation of Snow (SOS) experiment (Lundquist et al., 2023), supported by the National Science Foundation. The general scope and coordinating activities amongst the programs and with SPLASH's host institution, RMBL, are described by the referenced studies.

The purpose of this article is to provide details on the deployment, data acquisition, post-processing, data organization, and availability of one component of SPLASH: the atmospheric surface flux stations (ASFSs). The ASFS contribution was supported by the NOAA Physical Sciences Laboratory (PSL) and the Cooperative Institute for Research in Environmental Sciences (CIRES). The ASFSs are autonomous, single-level flux towers designed to observe energy exchanges between the atmosphere and surface at sufficiently high resolution for supporting process-based research on the coupling between these layers. The ASFS observations have been post-processed and quality-controlled and have been archived at Zenodo with the Creative Commons Attribution 4.0 International license. We begin by contextualizing the observations within the surface energy budget equation in Sect. 2. In Sect. 3, we describe the ASFS system and deployment, followed in Sect. 4 by a description of the deployment sites. Section 5 contains details on data post-processing and uncertainties. Section 6 provides information about data availability and formatting and conclusions are in Sect. 7. Appendix A contains a list of abbreviations and symbols defined in the text.

2 Surface energy balance

The SPLASH observations described here focus on the terms necessary to solve the net surface heat flux equation, which represents the balance of vertical radiative, conductive, and turbulent energy transfer across the infinitely thin plane of the surface such that

$$\begin{array}{}\text{(1)}& F=Q-H_{\mathrm{s}}-H_{\mathrm{l}}+C,\end{array}$$

where H_s is the turbulent sensible heat flux, H_l is the turbulent latent heat flux, C is the soil conductive flux at the surface, and Q is the net surface radiation, defined as follows:

$$\begin{array}{}\text{(2)}& Q=({\text{SW}}_{\mathrm{D}}-{\text{SW}}_{\mathrm{U}})+({\text{LW}}_{\mathrm{D}}-{\text{LW}}_{\mathrm{U}}),\end{array}$$

where D and U refer to downwelling and upwelling components of the shortwave (SW) and longwave (LW) radiative fluxes. Q is defined as positive warming of the surface, while H_s, H_l, and C are positive cooling of the surface, as is the typical convention (these are also the conventions used in the data set).

In practice, observations of Eq. (1) frequently do not balance (e.g., Foken et al., 2006), and thus it is more practical to formulate Eq. (1) including a residual, R, such that

$$\begin{array}{}\text{(3)}& F=Q-H_{\mathrm{s}}-H_{\mathrm{l}}+C+R.\end{array}$$

Grachev et al. (2020) provide a thorough analysis and review of R and its implications for terrestrial applications in balancing Eq. (3). We encourage readers to consult that resource for more details, while our purpose here is to contextualize the observations we made during SPLASH. In keeping with Grachev et al. (2020),

$$\begin{array}{}\text{(4)}& R=T+S+X,\end{array}$$

where T represents unaccounted horizontal and vertical transports (e.g., from spatial heterogeneity); S is the storage in the vertical column defined by the instrument locations relative to the surface, which in practice does not meet the infinitely thin requirement of Eq. (1); and X represents a variety of sources of uncertainty in the terms comprising Eq. (3). S can be further subdivided into sources. Again, following Grachev et al. (2020),

$$\begin{array}{}\text{(5)}& S=S_{\mathrm{a}}+S_{\mathrm{g}}+S_{\mathrm{p}}+S_{\mathrm{c}}+S_x,\end{array}$$

where S_a represents flux divergence within the air column between the surface and the actual measurement heights; S_g is storage of energy in the ground between the surface and the subsurface conductive flux measurement; S_p is the photosynthetic heat storage; S_c is a vegetation biomass heat storage term; and S_x is all other terms, such as the heat flux from precipitation (e.g., He et al., 2023) and the latent energy associated with snowmelt. For the purposes of the present study, we further subdivide S_g into components of the soil (S_soil) and snow (S_snow) layers thusly:

$$\begin{array}{}\text{(6)}& S_{\mathrm{g}}=S_{\text{soil}}+S_{\text{snow}}.\end{array}$$

While the purpose of the present study is not to provide a discourse on these many details, it is practical to provide the equations because it permits us to be more precise in defining the measurements that were made and their associated processing, as well as to highlight those terms that have not yet been constrained or for which relevant measurements were not made.

3 ASFS

The ASFSs were initially designed and used for observing the energy exchange over sea ice in the central Arctic Ocean at unstaffed locations accessed only periodically during the Multidisciplinary drifting Observatory for the Study of Arctic Climate (MOSAiC) (Cox et al., 2023a; Shupe et al., 2022). The device (Figs. 1e, f and 2d, e) is an aluminum pipe structure approximately 3 m long and 2 m high, set on skis similar to a rigid Nansen sled, on which an electronics box is mounted. From the main structure protrudes a horizontal boom 3 m outward from the front at a height of 2 m from the surface and a vertical mast 3 m high, both for mounting instrumentation. The station is powered using an EFOY Pro Duo 2400 direct methanol fuel cell supplied by 120 L of methanol, which generates the 60–70 W required for typical use for up to 2 months between refueling. The instrumentation suite (Table 1) includes basic meteorology, sensors for detecting surface temperature and height (for snow depth), a sonic anemometer, an open-path gas analyzer, a four-component broadband radiation suite, a GPS, soil heat flux plates, and a soil temperature and moisture probe. The data acquisition and communications are managed by a Campbell Scientific CR1000X data logger.

Figure 1 The Kettle Ponds Annex site (KPS-A). (a) © Google Earth (Landsat/Copernicus) looking northward (upvalley). Panels (b) and (c) are slope and aspect maps of the KPS-A vicinity (USGS, 2023). Panel (d) is a photo of the valley taken facing east from the road depicted in (a). Black dots in (a)–(d) denote the locations of ASFS-30. Panels (d), (e), and (f) are photos of ASFS-30 in February 2023, June 2023, and September 2021, respectively.

Figure 2 The Avery Picnic site (AYP). Panels (a)–(c) are as in Fig. 1. Black dots denote the locations of ASFS-50. Panels (d) and (e) are photos of ASFS-50 in January 2022 and October 2021, respectively. Image in (a) is from © Google Earth (Landsat/Copernicus), and (b) and (c) are from USGS (2023). Photo credit for panel (d) is Michael R. Gallagher.

Table 1 Instrumentation, sampling, and uncertainties.

The details of ASFS use and post-processing of data during MOSAiC are described in Cox et al. (2023a). Two ASFSs built for MOSAiC were repurposed for SPLASH. The stations are optimized for cold climates and perform well, maintaining ice-free instrumentation (Cox et al., 2021, 2023a). There are, however, several aspects of the UCRB alpine environment that differ from Arctic sea ice in ways that were significant for the SPLASH deployment. While the systems were largely unchanged for SPLASH, some adjustments were made:

1. The stations were originally designed for a snow depth of <0.5 m, typical of the Arctic. For SPLASH, the stations needed to be raised to accommodate up to 2 m of snow. Instead of a costly redesign of the structure, we opted to position the sleds on top of wooden picnic tables, which increased the height of the boom to 2.75 m. While simplistic, this arrangement was found to be robust: we did not need to reseat any of the equipment, and the level measured at the mast was maintained within 1°. Due to winter accumulation of snow, the effective measurement height of the instrumentation changed over time.

2. The stations are capable of direct peer-to-peer radio communications (2.4 GHz) as a primary form of data transfer to a base station. At SPLASH, without line of sight in the mountainous terrain, the data were stored locally on microSD cards and retrieved manually during refueling visits, approximately monthly. Only infrequent (3 h) samples and diagnostics were transferred daily to the base station for monitoring purposes (displayed on the NOAA/PSL website, https://psl.noaa.gov/splash/, last access: 4 April 2025) using the station's backup (Iridium dial-up data) satellite link communications.

3. At MOSAiC, complementary subsurface measurements were collected using buoys deployed by collaborators (Nicolaus et al., 2022; Rabe et al., 2022). To fill this gap for SPLASH, a 0.5 m soil probe was added to the ASFS instrument suite to collect soil temperature and water content information. In year 2, an experimental thermistor string for measuring temperature profiles through snow was also added to one station. Details of these sensors are provided later in the article.

4. The EFOY is most efficient at sea level and specified to operate up to 1500 m in altitude. Additionally, while it cannot be started from a frozen state, it performs best in the cold and can overheat easily in its insulated box. The SPLASH domain has warm, snow-free summers and is near 3000 m elevation. Thus, we anticipated a decrease in efficiency in both harvesting energy and ventilation. While we do not have objective measures, our experience was that the system's fuel burn rate was similar to that at sea level. We improved the ventilation but still experienced several instances of overheating in summer. Overall, the lifetime of the fuel cells was shorter than expected, and some instruments had to be turned off to conserve power and extend the life of the system (described later).

4 Sites

The East River is located north of Gunnison, Colorado, along the western slopes of the Rocky Mountains. It flows southeast from Emerald Lake (3186 m) ∼62 km to a confluence (2446 m) with the Taylor River near Almont, forming the Gunnison River, a tributary of the Colorado River. Gothic (2895 m), an unincorporated town that is home to RMBL, is located near the headwaters of the East River at the base of Gothic Mountain (3850 m). The slopes of the subalpine valley near Gothic include mixed stands of aspen, subalpine fir, and Engelmann spruce, while few trees are found in the valley floor where fescue grasses and sagebrush dominate (Johnston et al., 2001). The area around Gothic is transiently ranchland during the shoulder seasons as cattle are moved between the high-elevation and down-valley rangelands. The ASFSs were installed in this valley floor environment. Once installed, the stations operated largely autonomously. Regular maintenance visits were made every 4–6 weeks to retrieve data and refuel, and irregular visits were made as needed to conduct repairs and collect samples. Maps and a photograph of the valley annotated with the locations of all SPLASH sites can be found in Fig. 1 from de Boer et al. (2023). More detailed information about the sites specific to the present article can be found in Figs. 1 and 2.

4.1 Kettle Ponds Annex

A high density of measurements from SPLASH, SAIL, and SOS were located at the Kettle Ponds site (“KPS” in de Boer et al., 2023), ∼2.25 km southeast of Gothic (“GTH” in de Boer et al., 2023). Kettle Ponds Annex site (KPS-A) was situated on a parcel of land owned by RMBL approximately 400 m south (and down-valley) of KPS. ASFS station 30 (“ASFS-30”), installed at 38°56.3686^′ N, 106°58.1781^′ W on 28 September 2021, was the only equipment stationed at KPS-A. ASFS-30 operated until 19 July 2023. KPS-A lies within a wide, gently sloping portion of the valley, locally 4–7° (USGS, 2023) (Fig. 1a–c). The elevation at the height of the station boom was 2856 ± 0.6 m. In addition to the slope southward along the main trajectory of the valley, there is significant local component of the slope westward (Fig. 1c) toward the East River, which is ∼60 m lower in elevation than KPS-A. The nearest of the kettle holes that form the namesake of the site is located ∼180 m upslope. Photographs of the station can be found in Fig. 1d–f.

ASFS-30 was installed on two picnic tables with the boom facing approximately south to reduce shading of the radiation sensors at the end of the boom. Note that the direction the boom points is defined as “station north” with respect to the relative wind direction. Conductive flux plates were installed approximately 3.2 m offset from the centerline of the boom on either side of the station at 5 cm depth in the soil. A soil probe was installed approximately 4 m off centerline of the boom to the east. The boom height over the snow-free surface as measured by snow height sensor was 288 cm. Fencing was installed during the summer at a radius of approximately 8 m from the end of the boom to protect the instrumentation from cattle.

Surface soil samples were collected, and soil properties were measured from the immediate area around the station (Table 2). A texture analysis was performed by the Colorado State University Soil, Water and Plant Testing Laboratory. These data were also obtained for KPS. The dry soil density of the samples at KPS-A and Avery Picnic site (AYP) was also calculated after the samples were dried at 205 °C for 48 h. We obtained measurements of thermal conductivity using a handheld METER Tempos probe in the field in the vicinity of the sample locations. These measurements were converted to dry thermal conductivity using coincident, co-located measurements of soil volume water content (VWC) following the method of Peng et al. (2017). The values found in Table 2 represent means of values obtained within the upper 10 cm of the soil on 12 June and 29 August 2023.

Table 2 Soil properties obtained from Avery Picnic (AYP), Kettle Ponds Annex (KPS-A), and Kettle Ponds (KPS).

4.2 Avery Picnic

The Avery Picnic (AYP) site is approximately 1.8 km upvalley from Gothic along Colorado State Highway 317. The site is adjacent to the road on its west side, just a few meters above the East River, which is ∼60 m further west. At AYP, ASFS station 50 (“ASFS-50”) was installed at 38°58.3455^′ N, 106°59.8113^′ W on 12 October 2021 and operated until 21 June 2023. The elevation at boom height was 2933 ± 0.7 m. The terrain surrounding AYP is generally flat, and while not sloping like KPS-A, the surface is undulating in the vicinity of the meandering East River (Fig. 2a and c). The valley floor is narrower than at KPS-A, and the ground rises steeply to the east and west of the station at distances of ∼150–250 m (Fig. 2b). ASFS-50 and a separate snow-profiling thermistor string (described later) were the only equipment installed at AYP. Photographs of the station can be found in Fig. 2d and e.

The deployment approach closely mimicked that of ASFS-30 at KPS-A with the station positioned on two picnic tables, the boom pointing south, and flux plates positioned 3.5 m off either side of boom centerline. The 0.5 m soil probe was deployed approximately 4 m east of the rear of the station instead of near the boom because of difficultly encountered with a layer impenetrable to the auger near 0.45 m depth. The location that was chosen also had this layer slightly deeper, and the probe was offset from being flush with the mean surface by 3 cm, but because of local variability in the surface the uppermost sensor was covered by soil and not exposed to air. The boom height over the snow-free surface as measured by the snow height sensor was 302 cm. Fencing was also installed around ASFS-50 at AYP in summer. The same procedure for soil sampling conducted at KPS-A was also carried out at AYP (Table 2). Additional spatial sampling of near-surface VWC using a handheld Fieldscout time domain reflectometry (TDR) soil moisture meter was done periodically in the vicinity of both KPS-A and AYP, and the data may be obtained from Zenodo (Intrieri et al., 2023).

5 Data processing

In this section, we provide details relevant for user interpretation of the data set pertaining to the observing strategy, post-processing, and uncertainties of the measurements. Data gaps occurred both from downtime for entire stations and individual sensors. Instrument uptimes during the experiment are displayed in Fig. 3.

Figure 3 Time series of daily uptime (% operational) for subsets of the station sensors at (a) Kettle Ponds Annex (KPS-A) and (b) Avery Picnic (AYP). White=0 % uptime, and dark green=100 %. Summary values on the right are relative to the deployment periods of individual sensors.

5.1 Meteorology

Air temperature, pressure, and humidity were measured at both sites using a Vaisala PTU307 probe. The sensor was affixed to the boom, and the precise instrument heights (relative to snow-free ground) are found in Table 1. A detailed description of this sensor and intercomparisons can be found in Cox et al. (2023a). Post-processing of the data collected during SPLASH is identical to that described in the earlier study. Figure 4a shows time series of air temperature as an example of the available meteorological data.

Figure 4 Sample measurement time series from the ASFS at SPLASH: (a) 10 min means of air temperature; (b) soil temperature (colors; color map, left axis) and snow depth (black; right axis) at Kettle Ponds Annex (KPS-A); (c) as in (b) for Avery Picnic (AYP); (d) net radiation, Q, at both sites; and (e) temperatures measured at AYP by the SIMBA thermistor string. Values above the black contour are in air, values below are in snow, and values at negative depth (below the horizontal dashed line) are in soil. The gap in the contour from 22 March to 4 April 2023 is the period when the top of the thermistor string was covered by snow. Dates are MM/DD beginning in September 2021 and continuing through July 2023.

Downtime was only experienced for these systems when an entire station was down (Fig. 3), which was due to overheating (summer) or (likely) ice-obstructed exhaust lines (winter). This occurred from 19 May through 2 June 2022 at KPS-A and 20 January through 21 February, 7–14 March, and 26 May through 12 June 2023 at AYP. The total uptime for the systems relative to the deployment beginning and end dates was ∼97 % and ∼90 % for KPS-A and AYP, respectively. Some individual measurements experienced more downtime than the parent system (see Fig. 3).

5.2 Snow depth

Snow depth was calculated using Campbell Scientific SR-50A acoustic rangers. The sensors were mounted on the boom and averaged 12 samples of distances to the target (surface) each minute. These distances were corrected for temperature dependencies on the speed of sound and then converted to snow depth by subtracting the measured value from the snow-free reference value reported in Sect. 3. The data were quality-controlled using an internal diagnostic reported by the sensor for signal strength and also manually by using a despiking routine. Snow depths are shown in Fig. 4b and c (black lines). There is a gap in the AYP record prior to 1 February 2022 because of a failure in the originally installed sensor (Figs. 3b and 4c). Data collected during the snow-free period (not displayed in Fig. 4) are preserved in the data set with a caution flag applied because of noise associated with a poorly defined reflecting surface caused by vegetation.

5.3 Radiation

Broadband shortwave (285–3000 nm; SW_D, SW_U) and longwave (4.5–40 µm; LW_D, LW_U) radiation were measured using Hukseflux SR30-D1 pyranometers and IR20-T pyrgeometers, respectively. The sensors were fixed to a common mounting plate that could be manually adjusted to achieve level, which was measured continuously by inclinometers embedded in the SR30s. The mean level of the upward-facing pyranometers for the duration of the experiment was 0.5 ± 0.28° (KPS-A) and 0.11 ± 0.28° (AYP). The sensors were heated and ventilated. They have been shown previously to perform well in self-maintenance for ice mitigation (Cox et al., 2021, 2023a), and icing was not observed in the SPLASH data except for occasional cases of snow accumulation on the upward-facing dome during active precipitation. The SR30s directional response is <1 % (Cox et al., 2023a), and the sensor's infrared loss characteristics are <3 W m⁻² (Wang et al., 2018). Calibrations were factory (Table 1). The SR30s are digital RS485 protocol, and the IR20s are analog; a 100 Ω 0.01 % reference resistor was used to derive case temperature from the IR20s for calibration. The quality control procedure included both manual review and automated detection of spurious data using the “QCRAD” methodology (Long and Shi, 2008). The time series of net radiation, Q, is shown in Fig. 4d. The uptime for the radiative fluxes is similar to that of the meteorology (Fig. 3).

To supplement the broadband longwave data, Apogee SI-400 infrared thermometers were mounted to the boom to help characterize the surface temperature. The sensor at AYP performed as expected, but unfortunately the one at KPS-A shows values more characteristic of the sensor temperature than the target temperature, a feature most obvious when observing warm temperatures over the melting snow surface in spring when the target is ∼0 °C. These data from KPS-A are preserved in the data set with a cautionary flag, and we consider them suspect.

The surface skin temperature was derived from the pyrgeometers, accounting for emissivity and reflections, using the following equation:

$$\begin{array}{}\text{(7)}& T_{\mathrm{s}}={\left\{{\displaystyle \frac{{\text{LW}}_{\mathrm{U}}-(\mathrm{1}-\mathit{\epsilon }){\text{LW}}_{\mathrm{D}}}{\mathit{\epsilon }\mathit{\sigma }}}\right\}}^{\mathrm{0.25}},\end{array}$$

where ε is emissivity, and σ is the Stefan–Boltzmann constant. The emissivity was set to 0.985 when the surface was covered by snow (Warren, 1982) and to 0.975 for grassland (Hu et al., 2019) during the snow-free period. No further attempt to constrain the emissivity was conducted for the present data set, and doing so is considered ongoing work by the SPLASH team and collaborators.

Occasionally in spring very large gradients in near-surface temperature (>10 °C between T_s and 2 m) were observed in association with above-freezing air being advected over the snow surface. Under these conditions, the air in the path between the LW_U sensor (mounted to the boom) and the surface (target) was sufficiently warmer than the target to bias the measurement high due to emission at wavelengths >14 µm where the IR20-T is sensitive and the water vapor and CO₂ in the atmosphere are highly absorptive. These factors are contributions to S_a from Eq. (5), which we accounted for in LW_U so that the T_s could be more accurately derived. To estimate this influence, we ran a series of broadband longwave flux calculations using the Rapid Radiative Transfer Model (RRTM) (Mlawer et al., 1997) with input from representative radiosonde profiles launched by SAIL (Keeler et al., 2021) and specifying a range of near-surface temperature gradients between 0 and 14 °C. From the results, we calculated the flux divergence between the surface and a nominal instrument height of 2 m over snow as a function of the gradient, which is approximately linear in the ranges tested and of similar scale to the calibration uncertainty. A simple linear fit is used to apply a first-order correction to the LW_U based on the initially observed gradient:

$$\begin{array}{}\text{(8)}& S_{\text{a,LW}}=\mathrm{0.29193}\cdot \mathrm{\Delta }T-\mathrm{0.021642},\end{array}$$

where ΔT is the air temperature minus an initial calculation of the surface skin temperature. Analogous contributions to S_a in SW_D and SW_U are considered negligible, and we did not adjust LW_D in the data set.

Figure 5 presents a cross-validation between the radiative fluxes (1 h means) observed by ASFS-30 at KPS-A and those made ∼400 m upvalley at KPS by the NOAA Global Monitoring Laboratory (GML) (Soldo et al., 2023), also part of SPLASH. The comparison is independent: the equipment, calibration procedures, maintenance, mounting, and post-processing all differ between the GML and ASFS stations. We expect SW_D to be comparable when averaged over 1 h across the distance between the stations with the potential exception of the timing of mountain shadows near sunrise and sunset. The results in Fig. 5a demonstrate excellent reproducibility between the systems, with a combined bias of <1.9 %, which is better than the (combined) target uncertainty of 2.8 % (McArthur, 2005). We also expect LW_D to be comparable. Indeed, 69 % of the observations fall within the expected (combined) 2.8 W m⁻² 1σ difference (McArthur, 2005; Cox et al., 2021), which compares well to a 68 % nominal expectation of a normal distribution. However, approximately 9 % of the observations comprise a long tail of warm values at KPS-A relative to KPS. We do not have a definitive explanation for these differences. The two stations were at different heights above the surface, but the differences are uncorrelated with the near-surface temperature gradient. The differences tend to occur when the atmospheric boundary layer is well-mixed beneath thick cloud cover. We can rule out icing of the KPS sensor (which would result in a cold bias) because the occurrences are equally likely when the temperature is above 0 °C, as below. One speculative possibility is that during precipitation events, the case heating of the IR20-T at KPS-A increases the temperature of droplets temporarily occupying the dome surface, thereby enhancing the measured signal. Differences in LW_U (Fig. 5c) and SW_U (Fig. 5d) are not easily separable from actual heterogeneity in the surface but are shown for completeness. The comparison of SW_U shows a systematic difference of $\sim +\mathrm{8}$ % at KPS-A relative to KPS, which appears only in winter (thus are not related to differences in vegetation). We speculate that the difference is explainable by reflection of sunlight from snow on the nearby mountain slopes. The GML station featured a housing that shades the ∼5° of the hemisphere closest to the horizon. The intention is to prevent the direct beam from appearing in SW_U near sunrise/sunset due to minor errors in level, but in the East River valley it would have reduced the sunlight reflected from neighboring slopes as well. On the ASFS, the full hemisphere was observed. Note that neither approach produces a spurious result but rather a different perspective. The GML system likely produces a better estimate of the surface albedo (${\text{SW}}_{\mathrm{U}}/{\text{SW}}_{\mathrm{D}}$), a property of the surface, whereas the ASFS station likely provides a more complete estimate of the local net SW (SW_D−SW_U), a property of the local radiation field. Regardless, the differences are overall small and could also be explained in part by spatial heterogeneity in snow morphology and surface slope.

Figure 5 Comparison of 1 h mean radiative fluxes between measurements made by ASFS-30 at Kettle Ponds Annex (KPS-A) and those made at Kettle Ponds (KPS) by the NOAA Global Monitoring Laboratory (GML), approximately 400 m upvalley from KPS.

5.4 Soil properties and soil conductive flux

The soil probe was a 0.5 m length Campbell Scientific SoilVUE10 TDR sensor that measures VWC, temperature, electrical conductivity, and permittivity at 5, 10, 20, 30, 40, and 50 cm depth. In addition to periods when the stations were down, the probe at KPS-A suffered a severed cable in June 2022 that could not be repaired until September (Fig. 3a). As an example of data from the probes, soil temperatures are plotted in Fig. 4b and c (colors). In addition to the probe, conductive flux plates (Hukseflux HFP01) were installed in the soil 5 cm below the surface. Two plates were installed at each site (Sect. 4), but unfortunately one of the flux plates failed at KPS-A. There were two problems with this sensor: first, it was inadvertently initially a single-ended measurement (floating ground), and, second, in June 2022 the sensor wire was severed, and repairs were not possible. We have retained the observation prior to June 2022, which is noisy and exhibited an offset (the offset was found using periods of isothermal soil temperatures and subtracted off), but the data are given a cautionary flag and considered suspect.

Soil conductive flux was measured in two ways. The first method uses the flux plates, which made direct measurements in situ and thus represent the flux at 5 cm depth. Consequently, relative to the flux at the surface interface, the flux plates have not been adjusted for storage and are best defined as C+S_g (Eqs. 3 and 6) where the relative values of S_soil and S_snow may be dramatically different between winter and summer and may vary substantially throughout the day as well. We were not able to calculate S_soil due to a lack of co-located temperature data, which we believe may have been spatially variable due to the fact that a comparison of the flux plates at AYP suggests local spatial heterogeneity of at least up to a factor of 2. Calculations of S_snow are considered ongoing work (see Sect. 5.5).

We corrected for deflection error in each of the flux plates (Sauer et al., 2007; Morgensen, 1970), which arises from differences between the conductivity of the plate (0.76 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{K}}^{-\mathrm{1}}$) and the surrounding matrix. The first step is to calculate the thermal conductivity of the soil. We began with the measured dry thermal conductivity of the soil (Sect. 4.1). Then, a time-variant conductivity was calculated based on the observations of VWC following Peng et al. (2017). For times when the soil probe VWC was unavailable, a mean value was used, and the dependent variables in the data set are given a caution flag. The soil properties necessary to make the calculation are found in Table 2.

The second estimate of conductive flux comes from the soil probe. We calculated a layer-averaged C (i.e., C−S_soil) with reference to the upper 10 cm (KPS-A) and 7 cm (AYP) and also S_soil for the same layer. The equations (e.g., Morris, 2018) for these calculations are

$$\begin{array}{}\text{(9)}& C-S_{\text{soil}}=k_{\text{eff}}{\displaystyle \frac{T_{\mathrm{s}}-T_{\text{soil}}}{z}},\end{array}$$

where T_soil refers to the temperature at depth z, T_s is the surface temperature from Eq. (7), and k_eff is effective soil thermal conductivity, and

$$\begin{array}{}\text{(10)}& S_{\text{soil}}=C_{v}z{\displaystyle \frac{T_{\mathrm{s}}^{n+\mathrm{1}}-T_{\mathrm{s}}^{n-\mathrm{1}}+T_{\text{soil}}^{n+\mathrm{1}}-T_{\text{soil}}^{n-\mathrm{1}}}{\mathrm{2}(t_{n+\mathrm{1}}-t_{n-\mathrm{1}})}},\end{array}$$

where t is time at index n (references to n in the equation are indices, not computations), and C_v is the soil heat capacity, which was found using Eq. (10) from Abu-Hamdeh (2003) using inputs of the soil probe's VWC and values found in Table 2.

5.5 Snow temperature profiles

During the second year of SPLASH (Fig. 3b), a stand-alone thermistor string was added to the suite of instruments at AYP. The sensor is manufactured by SAMS Enterprise as the Snow Ice Mass Balance Apparatus (SIMBA) (Jackson et al., 2013). It was originally intended for use as a buoy in polar sea ice (e.g., Lei et al., 2022), but it was reconfigured by the manufacturer for terrestrial applications for the Scottish Avalanche Information Service. The terrestrial version, used for SPLASH, is 1.9 m in length with thermistors at 2 cm spacing affixed to a ridged plastic support structure and suspended from a tripod. In addition to ambient temperature, it also features a low-power heating cycle applied to the thermistors for which cooling rates are reported. These rates are correlated with the thermal environment, adding information for distinguishing between air and snow (Jackson et al., 2013). We installed the SIMBA prior to snowpack onset on 1 November 2022 approximately 3 m east of ASFS-50, and the snow was permitted to accumulate and ablate around the sensor throughout the winter and spring (Fig. 4e).

There were two unanticipated issues with this system. Initially, we planned to obtain profiles every 10 min, but a loss of the configuration file on 23 December 2022 was compounded by a firmware issue on 1 January 2023 that prevented us from accessing the data until after the snow melted in spring, resulting in profiles being collected instead every 6 h after the configuration changed. This same firmware problem also resulted in a loss of clock synchronization (with GPS), and consequently from January through June the clock drifted. We estimate that the drift was a maximum of 80 s but was likely much less. The second problem is that the snow depth in spring 2023 was much deeper than anticipated. Consequently, from 22 March to 4 April, the snow surface was over the top of the sensor. It is also notable that there was a gradient in snow depth with deeper snow at the ASFS and shallower snow east of the ASFS, where the SIMBA was located. The difference near peak snow depth was ∼25 cm.

The archived data set includes the temperatures and associated thermistor cooling rates. The SIMBA data may also be useful in deriving profiles of conduction, density, and k_eff through the snow (Sledd et al., 2024), allowing for C and S_snow to be obtainable at AYP during snow-covered period in 2022–2023. This work is experimental and ongoing, and, if successful, higher-order metrics will be made available as part of a future update to the data set.

5.6 Turbulent fluxes

Turbulent sensible (H_s) and latent (H_l) heat fluxes were derived using eddy covariance methodology. The fundamental measurements were made using a METEK uSonic-3 Cage Multipath sonic anemometer and a co-located LI-COR LI-7500DS gas analyzer. The raw data were acquired at 20 Hz, after which they were aggregated to 10 Hz by averaging, rotated into earth coordinates, and then quality-controlled (despiking, physically possible limits, sensor diagnostics). Horizontal wind speed and direction (1, 10 min means) were calculated from these processed components. For the fluxes, integration windows used in the provided data set are reported at 10 min intervals, though in practice the calculation is performed for the nearest power-of-2 number of samples, so the actual integration period is ∼13.65 min, centered; missing data within the window were filled with the median value of the window. To perform the flux calculation, u, v, and w were rotated using the planar-fit method (Wilczak et al., 2001), which puts u and v into the streamline of the horizontal wind direction (U,V) while rotating vertical (w→W) velocity to the long-term mean, where the mean in this case was determined using the observations from October 2021 through June 2022. The approach is thought to provide a more robust estimate of fluxes in sloping terrain, in particular for momentum fluxes. The coefficients (Fig. 6) for the method at both sites were an acceptable fit, but (as expected) it is only over the sloping meadow at KPS-A where the streamlines are substantially different from the level anemometer. Nevertheless, for consistency, the planar fit was used for both sites.

Figure 6 Planar fits (Wilczak et al., 2001) for Kettle Ponds Annex (KPS-A) (a) and Avery Picnic (AYP) (b). Panels (c) and (d) show for KPS-A and AYP the earth coordinate w^′ (blue) and the rotated value, W^′, respectively, after planar fit is applied (red). The targets for the distribution are for the mean and 1σ variability of the red distribution to be < the sensor uncertainty of 0.15 m s⁻¹.

Power (Welch) spectral densities and cross-spectral densities were calculated using the Python “SciPy” module; the fast Fourier transforms were Hamming-filtered and linearly detrended, and the resulting spectra were smoothed before integration (Grachev et al., 2018; Cox et al., 2023a). We also applied a correction for undesirable high-frequency (>1 Hz) correlated noise between W^′ and T^′ that is unique to the sonics used for these ASFSs and that we devised during an earlier study (Cox et al., 2023a). Shadow corrections for the sensor head geometry are performed internally by the sensor following Wyngaard and Zhang (1985). Buoyancy fluxes were then calculated from the cospectra and corrected for moisture fluctuations to yield H_s. While it is possible to make the correction directly using the humidity measurements from the co-located LI-COR, the LI-COR suffered more downtime than the sonic anemometer, so instead we use the method of Schotanus et al. (1983), which is based on the Bowen ratio and is deemed necessary for environments where the Bowen ratio <1 (∼36 % of the time at SPLASH). The main downtime in the LI-COR during SPLASH occurred at KPS-A (Fig. 3a), where the sensor was turned off after May 2022 to conserve power due to an aging EFOY generator (see Sect. 3). Secondarily, LI-COR data were lost when the sensor windows were covered by water, ice, or dust.

H_l and CO₂ mass fluxes were calculated using the gas measurements from the LI-COR. The mounting of the sonic anemometer and LI-COR on the mast follows the recommendations of Kristensen et al. (1997), both in its horizontal (being in line with the dominant wind direction) and vertical offsets, to minimize losses in covariance signal at high frequencies associated with the sensor separation. We thus do not perform a sensor separation correction in post-processing and assume the mounting geometry minimizes the error to <2 %, as found by Kristensen et al. (1997) (relative to a strictly horizontal separation and a random wind direction). LI-COR data were checked for quality using the recommended diagnostics, and data were also only retained when the optics were determined to be clean and unobstructed, as evidenced by the CO₂ reference signal strength exceeding 94 %. The Webb density correction (Webb et al., 1980) was applied to both gas fluxes. H_l was calculated using the latent heat of sublimation over dry-snow conditions and the latent heat of vaporization during snow-free conditions and when the surface snow temperature $>-\mathrm{0.2}$ °C, indicating the likelihood of wet or melting snow (when the liquid is assumed to evaporate preferentially to the sublimation of ice in the slurry). In addition to reporting the CO₂ mass flux ($\mathrm{mg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$), we also report an inferred photosynthetic storage value (S_p) in W m⁻² following the formulation used by Grachev et al. (2020) over grassland in the Columbia River basin. Because the quality of these CO₂ fluxes is not easily verified (e.g., Kohsiek, 2000), nor were CO₂ fluxes the primary focus of our measurement strategy, we have conservatively marked all CO₂ data with a caution flag in the files. In addition to the heat fluxes, other conventional variables related to stress (e.g., friction velocity, u^∗, drag coefficient, C_d) and the Monin–Obukhov stability parameter ($z/L$) are provided. We also provide diagnostics in the files (e.g., angle of attack) as well as the spectra and co-spectra from which the fluxes are calculated.

To supplement eddy covariance data during gaps, we also calculated turbulent fluxes using the bulk aerodynamic method of Andreas et al. (2010). The Andreas et al. (2010) method was designed for snow-covered sea ice and has a similar form to the “COARE” algorithm widely used over global oceans (Fairall et al., 2003). The algorithm operates iteratively to constrain a solution satisfying co-dependence between Obukhov length, L, and u^∗. We found that the method adapted well to snow-covered periods at SPLASH. The primary assumptions are that air in pore spaces of surface snow is saturated (for calculating H_l) and that surface roughness length (z₀) can be reasonably approximated. We calculated a median wintertime z₀ for each station using sonic anemometer data (Andreas et al., 2010; Grachev et al., 2007; Paulson, 1970) and found it to be close to that reported by Andreas et al. (2010) for wintertime snow in the Beaufort Sea ($\mathrm{2.3}\times {\mathrm{10}}^{-\mathrm{4}}$ m), so we used this published value for the snow-covered period.

Figure 7 Comparisons of 1 h means of H_s in (a) and (b) and H_l in (c) and (d) between eddy covariance (“ec”) and the bulk calculation (“bulk”). Outliers >4σ were removed from the comparison.

Table 3 List of file names and DOIs for level-1 and level-2 data sets.

We further adapted this method for snow-free periods. To do so, we first specified z₀ calculated as previously described for summer months ($\mathrm{2.4}\times {\mathrm{10}}^{-\mathrm{3}}$ at AYP and $\mathrm{3.6}\times {\mathrm{10}}^{-\mathrm{3}}$ at KPS-A). This was sufficient to produce useful estimates of H_s (Fig. 7a and b). Through much of summer, the soil surface layer is dry in the East River valley, and the assumption of saturation at the surface valid for snow fails for calculation of H_l. To calculate bulk H_l for snow-free conditions, we began by assuming that the surface layer relative humidity (RH) gradient tends toward 0 % and thus that H_l is dominated by diabatic processes. Kim et al. (2021) demonstrate that this is an appropriate approach based on their reformulation of the Penman–Monteith equation in terms of RH. This diabatic assumption applies best when the surface and near-surface air are in thermodynamic equilibrium, and the adiabatic term is negligible. This “surface flux equilibrium” has been shown to be a sufficient approximation to estimate evapotranspiration for daily-to-monthly means in inland continental locations using only single-level atmospheric measurements (McColl and Rigden, 2020). Its utility at subdaily scales when nonstationarity can introduce larger adiabatic influences on fluxes remains an open question (Kim et al., 2021). The results, shown in Fig. 7c and d, exhibit bias and correlative relationships to the eddy covariance measurements comparable to those reported by Grachev et al. (2022) over more homogeneous terrain in the Columbia River basin of Oregon. These statistics include large scatter amongst individual samples, which arises in part from random errors in the eddy covariance measurements, commonly 10 %–25 %, and sometimes larger (Kessomkiat et al., 2013).

Our implementation of concepts described by Kim et al. (2021) and McColl and Rigden (2020) makes use of the same Andreas et al. (2010) bulk methodology adapted to include the conventional “alpha” (α) approach for estimating surface specific humidity (Lee and Pielke, 1992):

$$\begin{array}{}\text{(11)}& q_{\text{sfc}}=\mathit{\alpha }{q}_{\text{sat}}\left(T_{\mathrm{s}}\right).\end{array}$$

We calculate α following Kondo et al. (1990):

$$\begin{array}{}\text{(12)}& \mathit{\alpha }={\displaystyle \frac{q}{q_{\text{sat}}\left(T_{\mathrm{s}}\right)}}+\mathit{\beta }\left(\mathrm{1}-{\displaystyle \frac{q}{q_{\text{sat}}\left(T_{\mathrm{s}}\right)}}\right),\end{array}$$

where $q/q_{\text{sat}}\left(T_{\mathrm{s}}\right)\sim {\text{RH}}_{\text{2m}}$ in accordance with Kim et al. (2021) and McColl and Rigden (2020). β is a wetness parameter that is calculated using soil VWC (Eq. 7 in Lee and Pielke, 1992) and helps to account for periods when the equilibration assumption is not met due to soil moistening. Except for the period immediately after snowmelt and peak saturation during rain events, the soil is dry and $\mathit{\beta }\ll \frac{q}{q_{\text{sat}}\left(T_{\mathrm{s}}\right)}$. Therefore, for periods when soil moisture data are unavailable, we set β=0 (consistent with McColl and Rigden, 2020); inclusion of β increases correlation (r) between eddy covariance H_l and analogous bulk estimates over snow-free surfaces by <0.05.

6 Data availability

The data are made available in daily netCDF4 formatted files and may be downloaded from Zenodo (Creative Commons Attribution 4.0 International license). Several types of files are provided. Level-1 files (Cox et al., 2023e: https://doi.org/10.5281/zenodo.10307825; Cox et al., 2023f: https://doi.org/10.5281/zenodo.10310520) are raw data with only technical corrections (e.g., sign conventions) and metadata applied. These level-1 files include “slow” and “fast” versions, which refer respectively to data with a native resolution of 1 min (most data) and 20 Hz (sonic anemometer and gas analyzer). These data are largely for archival purposes as it is recommended that most users use level-2 data. A separate level-1 data set was also created and made available for the SIMBA data (https://doi.org/10.5281/zenodo.10327409, Cox et al., 2023c).

A processed level-2 data set (Cox et al., 2023d: https://doi.org/10.5281/zenodo.10313894; Cox et al., 2023b: https://doi.org/10.5281/zenodo.10313363), which includes quality control, corrections, and derived variables as outlined in this article, is the data set recommended for most scientific purposes. There are three types of level-2 files: “sledmet*”, “sledseb*”, and “sledwind10Hz*”. Sledmet contains the 1 min means of all relevant variables, and sledseb contains 10 min means of the same variables, as well as all the bulk and eddy covariance calculations. Sledwind10Hz contains the quality-controlled, aggregated, and rotated earth coordinate sonic anemometer and gas analyzer data to facilitate users who may wish to recalculate turbulent fluxes using other approaches. Flags containing quality control information are found in the files, denoted *_qc with the coding for the flags encapsulated in the file's global attributes. Table 3 summarizes the file names and DOIs. These data are found within the SPLASH Zenodo “community” amongst complementary data from other contributing SPLASH projects. Complementary data from SAIL can be found at the ARM archive, https://arm.gov (US Department of Energy, 2025), and complementary data from SOS can be obtained through the University of Washington, https://depts.washington.edu/mtnhydr/Pages/NewDataPages/SOSdata.html (University of Washington, 2025).

7 Conclusions

The data set described here was collected by two single-level, autonomous atmospheric surface flux stations (ASFSs) during the Study of Precipitation, the Lower Atmosphere, and Surface for Hydrometeorology (SPLASH) campaign (de Boer et al., 2023) in the East River watershed of Colorado. Collectively between the stations, 1184 d of data was obtained between 2021 and 2023 that include meteorology, surface energy budget variables, and (sub)surface properties.

Appendix A: List of abbreviations and symbols defined in the text

Organizations and campaigns
ARM	Atmospheric Radiation Measurement program
CIRES	Cooperative Institute for Research in Environmental Sciences
DoE	Department of Energy
GML	NOAA Global Monitoring Laboratory
LBNL	Lawrence Berkeley National Laboratory
MOSAiC	Multidisciplinary drifting Observatory for the Study of Arctic Climate
NOAA	National Oceanic and Atmospheric Administration
PSL	NOAA Physical Sciences Laboratory
RMBL	Rocky Mountain Biological Laboratory
SAIL	Surface Atmosphere Integrated field Laboratory
SOS	Sublimation of Snow experiment
SPLASH	Study of Precipitation, the Lower Atmosphere, and Surface for Hydrometeorology
Locations
AYP	Avery Picnic site
GTH	Gothic
KPS	Kettle Ponds site
KPS-A	Kettle Ponds Annex site
UCRB	Upper Colorado River Basin
Models and sensors
ASFS-30	Atmospheric surface flux station 30
ASFS-50	atmospheric surface flux station 50
HRRR	High-Resolution Rapid Refresh model
NWM	National Water Model
RRTM	Rapid Radiative Transfer Model
SIMBA	Snow Ice Mass Balance Apparatus
TDR	Time domain reflectometry
Variables and parameters
C	Conductive heat flux
F	Net surface heat flux
Hl	Latent heat flux, eddy covariance measurement
Hlb	Latent heat flux, bulk model calculation
Hs	Sensible heat flux, eddy covariance measurement
Hsb	Sensible heat flux, bulk model calculation
keff	Effective thermal conductivity
LWD	Downwelling longwave radiative flux
LWU	Upwelling longwave radiative flux

Variables and parameters
Q	Net surface radiative flux
qsat	Saturation specific humidity at the surface
qsfc	Soil surface specific humidity
R	Surface energy balance residual
RH	Relative humidity
S	Energy storage in the vertical column between the sensor and the surface
Sa	Flux divergence in air
Sc	Vegetation biomass heat storage
Sg	Storage in the subsurface above the measurement of C
Sp	Photosynthetic heat storage
Ssnow	Portion of Sg in snow
Ssoil	Portion of Sg in soil
SWD	Downwelling shortwave radiative flux
SWU	Upwelling shortwave radiative flux
Sx	Unspecified contributions to S
T	Sources of residual in the surface energy budget due to heterogeneity
Ts	Surface temperature
Tsoil	Soil temperature
VWC	Volume water content
X	Sources of uncertainty contributing to R
α	Lee and Pielke (1992) wetness coefficient
β	Lee and Pielke (1992) wetness coefficient
ε	Surface emissivity
σ	Stefan–Boltzmann constant (Eq. 7), standard deviation (text)