Estimating snow water equivalent using cosmic-ray neutron sensors from the COSMOS-UK network

The intensity of cosmic ray neutrons is inversely correlated with the amount of water present in the surrounding environment. This effect is already employed by around 50 neutron sensors in the COSMOS-UK network to provide daily estimates of soil moisture across the UK. Here, these same sensors are used to automatically provide estimates of snow water equivalent (SWE). Lying snow is typically ephemeral and of shallow depth for most parts of the UK. Moreover, soil moisture is usually high and vari-able, which acts to increase uncertainties in the SWE estimate. Nevertheless, even under such challenging conditions, both above ground and buried cosmic ray neutron sensors are still able to produce potentially useful SWE estimates. Triple collocation analysis suggests typical uncertainties of less than around 4 mm under UK snow conditions.

sors to produce such representative measurements.
Cosmic ray neutron sensors, positioned a few metres above ground level, are often used to provide estimates of soil moisture as it is inversely correlated with the neutron count rate (Andreasen, Jensen, Desilets, Franz, et al., 2017;Desilets et al., 2010;Evans et al., 2016;Zreda et al., 2008;Zreda et al., 2012). These sensors have also been employed to measure snow water equivalent (SWE), exploiting the sensitivity to the presence of water held in the snowpack as well as the soil (Bogena et al., 2020;Desilets, 2017;Desilets et al., 2010;Rasmussen et al., 2012;Schattan et al., 2017;Schattan et al., 2019;Schattan et al., 2020;Sigouin & Si, 2016). Importantly, work on deriving soil moisture estimates has shown these sensors to have footprints of a few hundred metres (Desilets & Zreda, 2013;Köhli et al., 2015;Zreda et al., 2008), which, in the context of SWE estimations, should allow representative measurements of heterogeneous snowpacks over a similar footprint. An alternative is to employ a neutron sensor buried a few centimetres below ground level (Avdyushin et al., 1982;Gugerli et al., 2019;Howat et al., 2018;Kodama, 1980;Kodama et al., 1979;Paquet et al., 2008), which can be used to measure the overlying accumulation of SWE, albeit with a much smaller footprint. Both methods of deploying neutron sensors are appealing as they may be operated automatically in remote locations, have limited disturbing effect on the snow, and are sensitive to the SWE directly as opposed, for example, to snow depth.
In this study, data are sourced from the COSMOS-UK network: a network of approximately 50 measurement sites, each equipped with a cosmic ray neutron sensor positioned approximately 2 m above ground level and referred to herein as a CRNS. A method is developed to produce SWE estimates using CRNSs in the COSMOS-UK network, which, to date, have been employed solely for the purpose of estimating soil moisture (Evans et al., 2016). In total, 874 daily CRNS SWE estimates from 254 separate snow events have been produced at 46 COSMOS-UK sites for the five winters between 2014/2015 and 2018/2019 and made freely available (Wallbank et al., 2020).
In order to separate the neutron signal due to SWE from that due to soil moisture, previous studies have found some benefit in using a pair of neutron counters (moderated and bare) with sensitivities to different neutron energy regimes (Desilets et al., 2010;Rivera Villarreyes et al., 2011). Others have either, (i) relied on the SWE signal dominating the soil moisture signal for sufficiently deep snow (Schattan et al., 2017), (ii) have used yearly calibrations of the SWE and ascribed differences between years to differences in the soil moisture (Sigouin & Si, 2016), or, (iii) use a method that takes into account an estimate of the neutron signal for snow-free ground (Desilets, 2017).
A variety of approaches were trialled by Bogena et al. (2020). Similar techniques have also been applied to quantify the amount of aboveground biomass, including the use of paired (moderated and bare) neutron counters, and methods accounting for soil moisture: a recent review of literature is provided in Vather et al. (2020).
Here the method proposed in Desilets (2017) (method (iii)) is adopted because the neutron signal will be strongly influenced by the underlying soil moisture for the shallow (typically <30 mm SWE) and ephemeral (typically shorter than one week) snowpacks found in the majority of the UK, and that the soil moisture can vary significantly between periods of lying snow. The neutron signal for snow-free ground will be estimated using the neutron count rate from immediately before each snow event and will optionally be adjusted using data from point soil moisture probes. This method may be applied in near real-time across all COSMOS-UK sites without the need for sitespecific calibration, or a requirement for paired neutron counters.
A subset of seven suitable COSMOS-UK sites additionally include a buried neutron sensor, referred to herein as a 'buried-CRNS', which are used to produce daily SWE estimates for 489 snow days from 126 separate snow events for the five winters. These have also been made freely available (Wallbank et al., 2020). SWE estimates from the CRNS and buried-CRNS are validated against alternative SWE estimates, either based on a snow depth measured automatically at certain COSMOS-UK sites, or, using the output of a snowmelt model (Moore et al., 1999). The availability of three independent SWE estimates allows the use of triple collocation (Stoffelen, 1998) to estimate the uncertainty in the different SWE estimation methods. This analysis suggests a typical uncertainty in the SWE estimate of less than around 4 mm for both CRNS and buried-CRNS. Even though this level of uncertainty is fairly modest, it is nevertheless greater than the estimated SWE for approximately 21% (for the CRNS), or 15% (for the buried-CRNS) of the recorded snow days.
For the CRNS and buried-CRNS, the uncertainty can be attributed to both the statistical uncertainty in the neutron count rate (which is reduced by using a daily average), and uncertainty in estimating the neutron count rate for snow-free conditions. It is shown that the greatest uncertainty in SWE estimates arises from estimating the count rate for snow-free conditions. Additionally, for the CRNS the absolute uncertainty is shown to increase considerably with the soil moisture, and to a lesser extent with the SWE, while much weaker dependencies are found for the buried-CRNS.

| SITES AND INSTRUMENTATION
All data used in this study were collected by instrumentation at COSMOS-UK sites (https://cosmos.ceh.ac.uk; Evans et al., 2016). The locations of the 46 sites used in this study are shown in Figure 1. They vary in altitude from 3 to 565 m with a mean of 156 m, and cover a variety of soil type and land use. All sites are equipped with the following instrumentation positioned within an approximately 36 m 2 compound.
• CRNS: Hydroinnova CRS-2000 or CRS-1000/B models positioned approximately 2 m above ground level. The measured neutron counts are used to produce the CRNS SWE estimate (Section 3).
• Point soil moisture sensor: two Acclima Digital TDT sensors at a depth of approximately 10 cm. Optionally used to correct estimated snow-free neutron counts (Section 3.1).
• Radiometer: Hukseflux four-component. Measured incoming and outgoing short-wave radiation is used to calculate albedo for snow cover detection (Section 3.3).
• Automatic weather station: Rotronic HC2A-S3 within the Gill MetPak Pro Base Station. Measured pressure and humidity are used to correct CRNS and buried-CRNS neutron counts (Sections 4 and 3, respectively). Measured temperature is used in the snowmelt model (Section 5.2).
• Phenocam: Mobotix S14 IP. Used to assess the performance of albedo-based snow detection (Section 3.3) and to manually F I G U R E 1 Location and elevation of the 46 COSMOS-UK study sites. The seven snow sites are named and the four excluded sites shown by pink circles quality-control depth-based and snowmelt model SWE estimates (Sections 5.1 and 5.2).
The seven named sites in Figure 1 are designated 'snow sites' and are typically located where higher than average amounts of snowfall are expected. These sites have the following two additional sensors.
• Buried-CRNS sensor: a Hydroinnova SnowFox buried in the ground at depths of up to 10 cm. The measured neutron counts are used to produce the buried-CRNS SWE estimate (Section 4).
• Snow depth sensor: Campbell Scientific SR50A sonic ranging sensor. The measured depth is used to produce the depth-based SWE estimate (Section 5.1).
Four sites marked with pink circles in Figure 1 (Alice Holt, Gisburn Forest, Harwood Forest, and Wytham Woods) are excluded from this study because the significant forest land cover within the CRNS footprint is expected to complicate the analysis Tian et al., 2016). The first COSMOS-UK site

| ESTIMATING SWE FROM THE CRNS COUNT RATE
The method of estimating SWE from the CRNS count rate is based on the formula of Desilets (2017): This describes the exponential decay of the count rate at time t, N(t), from an estimated value for snow-free ground, N θ (t), to the count rate that would be detected over an infinite depth of water, N wat , with a decay length Λ = 48 mm obtained by Desilets (2017) from neutron particle transport modelling. To reduce the statistical noise present in the hourly neutron count, the quantity N(t) is taken to be the 24 h mean, from 12 h before time t until 12 h after. Because of this, SWE (t) and N θ (t) should also be regarded as 24 h mean values over the same period.
To obtain N(t), hourly accumulations of the neutron count are first corrected for the influence of variation in the pressure and water vapour of the atmosphere as described in Evans et al. (2016), and a new empirical data-driven method is used to correct for variation in the background neutron intensity based on count rates measured at the Jungfraujoch monitoring station. Simple automatic quality-control checks are performed on the hourly corrected counts, and a 24 h running mean is taken. Three separate methods of estimating N θ (t) are investigated, each obtained by adjusting the count rate from immediately before the snow period (Section 3.1). In contrast, a single calibration is used to estimate N wat for all sites (Section 3.2). Snow days are identified at each site via the locally measured albedo which is used to define the start and end of snow periods (Section 3.3).
For soil moisture estimation, the CRNS footprint was found to be about 300 m for dry air (Desilets & Zreda, 2013) with significant dependence on humidity, or about 120-230 m (Köhli et al., 2015) depending on both humidity and soil moisture. A comparable footprint should be expected for SWE estimates.    3.1 | Estimation of the neutron count rate over snow-free ground Three different methods have been investigated to estimate the 24 h mean snow-free count rate, N θ (t). All of these are based on adjusting the count rate from immediately before the albedo-based detection of snow cover (Section 3.3). The first method only adjusts N θ (t) if N(t) increases above N θ (t) (which would imply a negative SWE in Equation (1)). The other two methods additionally attempt to correct for variations in soil moisture using hourly data from two local point soil moisture probes located at a depth of approximately 10 cm.
For all three methods, the estimate at the initial time-step, t 0 , is given by the measured count rate, N θ (t 0 ) = N(t 0 ), with t 0 chosen to be the first midnight of the day preceding the albedo-based detection of snow cover. This choice excludes the possibility of snow occurring within the 24 h averaging period implicit in N(t 0 ) due to the fact that the albedo-based detection of snow cover is only effective around midday. For subsequent times, t > t 0 , the estimation of N θ (t) is based on its value at the previous time-step.

Method 1 simply uses
Method 1 simply uses where δt is the time-step and taken to be 1 h.

Methods 2 and 3 use
where Δθ(t) is an estimate of the change in 24 h mean volumetric water content (VWC) of the soil since the previous time-step, and f is the one-to-one function used by COSMOS-UK to obtain the VWC from the neutron count rate (Equation (5)  The correlation between the 24 h mean VWC detected by the two point-probes and that derived from the CRNS, f(N(t)), is typically found to be reasonably good at all sites. An example comparison is shown for the Easter Bush site in Figure 3. Some of the largest discrepancies are caused by the fact that the CRNS counts are sensitive to water held in a snowpack and will therefore tend to overestimate the VWC during snow days (marked with red asterisks in Figure 3).
Removing snow days from the analysis (identified by the albedo, see Section 3.3), the site-specific correlation coefficient, r, between 24 h mean VWC derived from the CRNS and the average of the two pointprobes was found to be 0.72 on average across the COSMOS-UK snow sites.
Nevertheless, the site-specific gradient, b, of the ordinary least squares regression line (blue line in Figure 3b) is often found to be quite different from 1 (typically 0.5 < b < 1), perhaps owing to the point-probe measurements at approximately 10 cm depth not sufficiently representing the footprint volume of the CRNS.
Motivated by the considerations above, Method 2 estimates Δθ(t) using the average VWC change from the two point-probes, F I G U R E 2 Reduction of CRNS counts due to snow cover. The hourly count rate (fine grey line) and its 24 h mean, N(t), (black line) are shown for example snow events from 2018 at (a) Easter Bush and (b) Moor House. Snow-free days are shaded grey and estimates of the snow-free count rate are shown using a blue dotted (Method 1), dashed (Method 2) or dot-dashed (Method 3) line. Hourly point-probe soil moisture (volumetric water content, VWC) is shown using green lines (use right y-axis) F I G U R E 3 (a) Time series of daily VWC at Easter Bush, derived from either the CRNS neutron count rate (black circles for snow-free days, red asterisks for snow days) or the two point-probes (green triangles). Note a period of missing point-probe data. (b) The same CRNS VWC data plotted against the average of the VWC measured by the two point-probes. The blue line shows the ordinary least squares linear regression (snow days excluded) with the correlation coefficient (r) and gradient (b) stated in the legend where θ 1 (t) and θ 2 (t) are the 24 h averaged VWC measurements from point-probes 1 and 2 centred on time t, to align with the calculation of N(t) and N θ (t). Method 3 estimates Δθ(t) using the point-probe that has the smallest absolute change in VWC, Method 3 is a conservative approach designed to reduce the impact of artefacts in the point-probe soil moisture data which may effect only one of the probes at a given time (e.g., see green triangles in Figure 3a). In practice, automatic quality-control is also applied to the point-probe VWC data for both Methods 2 and 3 to remove artificial spikes such as those visible on 4 and 9 March in Figure 2a.
The estimation of the snow-free count rate is continued until the measured albedo shows snow to be absent from the 24 h averaging period. The fact that all three N θ estimates (blue lines in Figure 2) display discontinuities at this time as they revert to the measured count rate, N(t), suggests some uncertainty in all three methods. This is explored further in Section 6.1.

| Estimation of the neutron count rate over an infinite water depth
The value of the count rate over an infinite depth of water, N wat , is expected to vary slightly between COSMOS-UK sites (e.g., because of differences in the neutron counter itself, or the local topography) but not between snow events at the same site. Desilets (2017) recommends using the value N wat = 0.24N 0 where N 0 is a site-specific calibration parameter used to normalize the count rate in the calculation of soil moisture. Using this value in Equation (1) and taking N 0 from the COSMOS-UK calibration (Evans et al., 2016; was found to systematically underestimate the SWE. To rectify this, an estimation of N wat /N 0 based on a rearrangement of Equation (1) is employed: To calculate χ, independent estimates of SWE(t) are used based on either, the snow depth (available at snow sites only and excluding periods of melt, see Section 5.1), or, the output of a snowmelt model (Section 5.2). Daily mean values of N wat /N 0 estimated using this method are shown in Figure 4.
Uncertainty in the estimate of N wat arises predominantly from uncertainties in N θ and χ (with standard deviations σ Nθ and σ χ , respectively), and to a lesser extent from statistical noise in the count rate (where the neutron counts are assumed to follow the Poisson distribution and the factor of 24 comes from taking a daily mean). Fortunately, the uncertainty this produces in the estimate of N wat decays rapidly with increasing SWE depth, as described by its standard deviation: Later (Section 6.1), a triple collocation analysis will suggest σ Nθ ≈12 and, for the depth-based SWE estimate, σ χ ≈ 4 mm/Λ. Using these values, the uncertainty estimate N wat AE σ Nwat is plotted as a function of χ for average values of N θ and N (dashed lines, Figure 4), and is consistent with the reduction in scatter of the points as the estimated SWE depth increases.
Taking N wat /N 0 to be its average estimated value from Equation (6) for days with SWE > 10 mm, and using Methods 1, 2 and 3 to estimate N θ , gives N wat /N 0 = 0.386, 0.380, and 0.383, respectively.
Here, the snow depth was used as the independent estimate of SWE and hence χ. This used 193 snow days from 73 snow events across the 7 snow sites, with a largest SWE estimate of 34.6 mm. Alternatively, using the snowmelt model to estimate χ gives N wat /N 0 ≈ 0.395 for all three methods. For simplicity its value is taken as Despite the reasonably good agreement of estimations employing either the snow depth or the snowmelt model, the possibility of a systematic relative bias in the estimation of χ cannot be excluded. This could be produced by either a systematic relative bias in the SWE estimates, or, an inaccuracy in the value for Λ. Relative bias in χ by an unknown factor, a, creates a biased N wat estimate, N a wat , according to  (1) is used to eliminate N(t) (which has implicit χ dependence) from Equations (7) and (9), and N wat = 0.38N 0 and N θ = 0.46N 0 are used. The latter is the average value of N θ for the snow days shown on the plot N a Estimates of N a = 0:5 wat and N a = 2 wat shown in Figure 4 (lower and upper dotted lines, respectively) indicate that the bias in N a≠1 wat decays only slowly with increasing estimated SWE depth. Moreover, a bias in the N wat estimate could, in turn, significantly bias the CRNS-based SWE estimate. For example at SWE = 10 mm, N a = 0:5 wat = 0:32N 0 and N a = 2 wat = 0:42N 0 would typically produce biases in the CRNS SWE estimate by factors of almost 0.5 and 2, respectively. Because of this, the value for N wat /N 0 given in Equation (8) should, to some extent, be regarded as calibrated to the particular value of Λ used and any relative bias in the depth-based or snowmelt SWE estimates. Ultimately, the dependence of the N wat /N 0 estimate on these details will vanish if measurements at larger SWE depths are available because N a wat ! N wat ! N as χ ! ∞.
The SWE estimation method proposed in Desilets (2017) (Equation (1)) was also tested in Bogena et al. (2020). In that study a calibration against observed SWE values was used to obtain N wat = 0.436N 0 . Although this is slightly larger than the value found in this study (N wat = 0.38N 0 ) it is consistent in suggesting a significantly larger value being required than that given in Desilets (2017) (N wat = 0.24N 0 ).

| Method for identifying snow days based on albedo
An independent method for identifying days with snow cover is provided by the albedo, calculated as the ratio of the outgoing to incoming fluxes of short-wave radiation measured by the radiometer at each COSMOS-UK site. The method aims to identify only those days with sufficiently deep snow cover (e.g., greater than 2-3 cm depth) to have a potentially significant effect on the neutron count rate. Identifying shallower snow events would unnecessarily flag the CRNS count rate as being inappropriate for the COSMOS-UK network's primary role of estimating soil moisture, while at the same time the resulting SWE estimate would be expected to be less accurate than the value of zero assumed for days without detected snow cover.
Typically, the onset of a snow event is marked by a sharp increase in the albedo as a uniform carpet of snow is deposited, followed by a more gradual decline as the snow slowly melts and becomes increasingly patchy (see Figure 5a). A clear influence of snow on albedo is also present for complex snow events featuring multiple periods of accumulation and partial melt (see Figure 5b). Because the COSMOS-UK sites are only illuminated by the sun, clearly the albedo can only be calculated during daylight hours. Moreover, the calculated albedo becomes an unreliable measure of the ground's actual albedo near sunrise and sunset, and often displays sharp increases or occasionally decreases around these times, particularly on cloud-free winter days or for sites located on sloping ground. In addition, snow at the end of an event is often patchy yet may still retain a significant water content, which contrasts with patchy snow at the beginning of an event which typically holds negligible water content. This suggests a stricter criteria on the albedo should be used to identify the beginning of a snow event compared to its end. Following these considerations, days with snow cover are identified by first calculating the mean of the albedo between the times of 10:00 and 14:00 GMT (using 30 min average data reported by the radiometer). Then, a period of snow cover is deemed to begin if this mean albedo exceeds a threshold of 50% and end when it falls below 35%. The threshold values used in this procedure were chosen by comparing albedo values with phenocam images at multiple sites.
This method was tested for the period January to March 2018 which contained over 300 days with snow cover identified by the method across all COSMOS-UK sites. Over 100 days in this 3-month F I G U R E 5 Albedo (grey line) calculated for daylight hours for the two example snow events in 2018 at (a) Easter Bush and (b) Moor House. The albedo for times between 10:00 and 14:00 (GMT) is highlighted in black, the mean albedo between these times is shown as a green cross, and dashed lines mark albedo thresholds of 50% (blue) and 35% (orange) used to delimit snow-free days (shaded in grey). (a) Snowfall occurs on 28 and 29 February and 6 March. A dusting of snow on 27 February does not exceed the 50% threshold. (b) Snowfall occurs on 3, 6 and 9 February. A heavy dusting of snow on 1 February almost causes the mean albedo to reach the 50% threshold period, including days that were not identified as having snow cover by the method, were spot-checked using phenocam images and the presence of significant snow cover (greater than 2-3 cm) was noted.
This exposed only four days inappropriately deemed to be snow-free by the method: three due to rough ground producing a patchy snow cover even for significant (though still small) average snow depths, and one caused by snow melting from under the radiometer yet persisting elsewhere. All checked days that were identified as having snow cover by the method had lying snow in the phenocam images, although some could be shallow depths. Based on this, the albedobased method was deemed to be reliable at identifying snow days and snow-free days. It was used to identify a total of 913 snow days from 263 separate events, of which 508 occurred during 128 separate events at the seven snow sites.

| ESTIMATING SWE FROM THE BURIED-CRNS COUNT RATE
The COSMOS-UK snow sites (named on Figure 1) additionally have a buried-CRNS. Its footprint is much smaller than that of the CRNS (e.g., <1 m). Because of this, the buried-CRNS-based SWE estimate will be much more strongly affected by inhomogeneities in the snowpack, such as those caused by snow scouring and drifting by wind that are often produced within the COSMOS-UK compounds during windy conditions.
To obtain a SWE estimate from the buried-CRNS, first a normalized count rate, N * (t), relative to an estimated snow-free value, N θ (t), is calculated as Then an attenuation curve, which describes the reduction of the normalized count rate N * with increasing SWE depth, is used to estimate SWE. Here, the buried-CRNS count rate, N(t), is corrected for variations in pressure, humidity, and the background neutron intensity using the method employed for the CRNS, and a 24 h mean centred on time t is used to reduce statistical noise. The reduction of the count rate due to snow cover is shown in Figure 6 for the two example snow events.
The choice of attenuation curve depends on the ground underlying the buried-CRNS due to a boundary effect. Attenuation curves given in Kodama et al. (1979) and Delunel et al. (2014) are applicable for a sensor above dry soil. In contrast, another set of attenuation curves, given in Howat et al. (2018) and Kodama et al. (1979), are applicable for a sensor above water, either in the form of a deep pre-existing snowpack (Howat et al., 2018), or liquid water (Kodama et al., 1979). In principle, the attenuation curve for a sensor above wet soil should lie somewhere between the two sets of curves, depending on the soil moisture (Kodama et al., 1979). However, for COSMOS-UK sites both sets of attenuation curves tended to overestimate the SWE relative to the other SWE estimates. The overestimation was found to be less severe for the curves applicable to a sensor above water, so the attenuation curve of Howat et al. (2018) was used. This choice is partly justified by the fact that the VWC of soil under snow in the UK is often greater than 50%. Also note that the SWE depths measured in this study are often smaller than those for which the attenuation curves were developed. Figure 7 shows that the buried-CRNS count rate has a strong negative correlation with the average of the daily 24 h mean VWC recorded by the two point-probes. Excluding snow days, the average correlation coefficient for the COSMOS-UK snow sites is −0.87: higher in absolute terms than the value of 0.72 calculated using the CRNS for these same sites (Section 3.1). This is despite the fact that the CRNS was developed primarily for measuring soil moisture, and that a site-specific non-linear calibrated function (f in Equation (3)) is used to convert the count rate into VWC before the comparison with the point-probe data. Instead, the stronger correlation found for the buried-CRNS may arise from the greater similarity of its footprint with that of the point-probe, and suggests the strong sensitivity of buried-CRNS to soil moisture which was exploited in Kodama et al. (1985).  (3)). Here, the site-specific gradient, β, of the ordinary least-squares regression line (blue line in Figure 7) is used instead. Hence, for the buried-CRNS is used in place of Equation (3), while in Equations (4) and (5) b is set to 1.
Using these methods daily buried-CRNS SWE estimates are produced for 489 snow days from 126 separate events for the seven suitable sites with buried-CRNS.

| VALIDATION OF SWE ESTIMATES
Neutron count based SWE estimates can be compared to independent estimations based on either a snow depth sensor, or, the output from a snowmelt model. Because these three SWE estimates are expected to have independent errors, the statistical method of triple collocation (Stoffelen, 1998) can be used to produce uncertainty estimates for each method.

| SWE estimate from the sonic ranging sensor
The snow depth is reported by a sonic ranging sensor every 30 min for each of the COSMOS-UK snow sites (named on Figure 1). These data contain a significant amount of noise as well as an unknown background height to subtract that can vary with time (see Figure 8).
Because of this, only the daily median depth, for the 24-h period centred on 12:00 GMT, is used. The background height is estimated to be the minimum of the daily median depths from the 2 days preceding the snow event, its last day, and the following 2 days. Depth data from the end of the event was used to estimate the background height because at certain sites (e.g., Cochno, Figure 8) the height is often found to be lower at the end of the snow event, possibly as a result of snow-flattened vegetation. A rough manual quality-control was performed on the snow depth data, consulting phenocam images when necessary. This resulted in the removal of snow depth data from five snow events, leaving daily depth-based SWE estimates for 193 snow days from 73 separate snow events in total.
To convert the daily median snow depth into a SWE estimate, a density of 0.1 g/cm 3 is assumed for fresh snow (i.e., 1 cm of snow depth equates to 1 mm of SWE). As snow ages its density increases in a complicated manner. Because of this, the depth-based SWE estimate is excluded from numerical comparisons for days in which the snow depth has decreased by more than 20% from its maximum in that snow event.
Allowing a 20% decrease in the depth was judged to have limited impact on the accuracy of the depth-based SWE estimate, but allows additional comparisons to be made during the middle of the snow events.
Uncertainty in this SWE estimate arises from uncertainty in the density (a range of 0.07-0.15 g/cm 3 is typical for fresh snow (Garstka, 1964)), and uncertainty in the depth measurement, the estimated background measurement, and how representative it is of the immediate vicinity. As the depth sensor only has a small footprint (about 1 m), it is affected by potential inhomogeneity in the snowpack in a similar manner to the buried-CRNS.
F I G U R E 7 Comparison of the daily average buried-CRNS count rate at Easter Bush with the average of the daily VWC measured by the two point-probes. Snow days are marked with red asterisks. The blue line shows the ordinary least-squares linear regression (snow days excluded) with correlation coefficient (r) and gradient (β) stated in the legend F I G U R E 8 Snow depth measurements (including background height) at the Cochno site in 2018. Small red circles indicate 30 min mean snow depths, and large filled/hollow circles show the daily median values (depth not/has decreased by more than 20% respectively). Days without snow cover (determined from the albedo) are shaded grey. Precipitation and air temperature (lower plot), measured at the site, aid the interpretation of the snow-depth data

| Snowmelt model estimation of SWE
The snowmelt model used here is PACK (Moore et al., 1999) and is applied to all COSMOS-UK sites. This takes inputs of air temperature and precipitation, the latter being designated as either snow or rain based on a temperature threshold. The snowpack is partitioned into two components: a dry pack consisting of lying snow prior to melting, and a wet pack containing water melted from the dry pack. Water is added to the dry pack in the form of snow and transferred to the wet pack according to a temperature-excess formulation of the melting process. Water may also be added to the wet pack due to rain. Drainage out of the wet pack is initially slow, increasing once a critical liquid water content is reached.
Here, parameters for PACK are taken from the Grid-to-Grid hydrological model (Bell et al., 2009;Cole & Moore, 2008;Moore et al., 2006) as used operationally over Scotland (see Table 1). The model is initiated without snow at 13:00 GMT on the day before the albedo-based detection of snow cover, and is terminated at 11:00 GMT on the first snow-free day. PACK was applied at an hourly timestep and a 24 h mean of its SWE (centred on time t) was calculated to make it comparable to the CRNS and buried-CRNS based estimates.
Both the air temperature and precipitation inputs are measured locally using the COSMOS-UK instrumentation (Section 2). Precipitation is measured using a weighing raingauge without either windshield or heated rim. During snow conditions it may suffer from blockages, wind-induced undercatch and the effects of turbulence (Rasmussen et al., 2012). SWE estimates based on suspiciously low input precipitation (less than 2 mm accumulated precipitation for the snow event) were excluded. Snow events with large modelled SWE values (SWE >50 mm) were examined and input precipitation values checked. From the 263 snow events identified by the albedo method (Section 3.3), snowmelt model estimates could not be produced for six events due to missing precipitation and/or temperature data. Estimates for 27 event were removed due to suspiciously low input precipitation, and estimates for a further five events were either removed or partially removed due to suspiciously large input precipitation. This left 799 daily SWE estimates from 228 snow events. There is also uncertainty in the snowmelt model due to incorrect partitioning of precipitation into either rain or snow, and the simplified conceptualisation of the melting process. However, since the precipitation and temperature measurements can generally be regarded as being representative of the immediate vicinity, the modelled SWE estimate should be as well.

| Uncertainty estimation using triple collocation
The availability of three independent estimates of the SWE at COSMOS-UK snow sites allows error estimates to be produced using triple collocation (Stoffelen, 1998). In contrast to dual comparisons which rely on the availability of an accurate reference measurement, triple collocation instead relies on the three measurement methods having uncorrelated errors.
Let x, y, and z represent estimates of the mean daily SWE, based either on the snow depth (x), CRNS or buried-CRNS (y), or the snowmelt model (z). Suppose the depth-based estimate x, is an unbiased estimate of the true mean daily SWE, S, and any bias in the remaining two estimations (y and z) is relative, then where ε x , ε y and ε z are measurement errors, and a y and a z are the relative biases. Under the assumption that the error terms have zero expectation, and are uncorrelated either with each other or with S, the following estimates are obtained (Caires & Sterl, 2003;Stoffelen, 1998), Here, the σ 2 = var(ε) are the absolute uncertainties squared, angular brackets hÁi denote statistical averaging, and only days with data for all three SWE estimates are used. Relative uncertainties, δ, will be calculated by dividing the absolute uncertainty by the mean of the When the CRNS-based SWE estimate is used for y, the true SWE should be regarded as having a large footprint, as both the CRNS and snowmelt model have large footprints. Any lack of representivity in the depth-based estimate due to its smaller footprint will then be included in its error term. Conversely, when the buried-CRNS is used for y, the true SWE should be regarded as having a small footprint.

| Attribution of CRNS and buried-CRNS uncertainties
The triple collocation uncertainties in the CRNS and buried-CRNS estimates (Section 5.3) can be attributed to two main sources: the statistical uncertainty in the counts, σ N = ffiffiffiffiffiffiffiffiffiffiffi ffi N=24 p (where the factor 24 accounts for daily averaging), and the uncertainty in the estimated count rate over snow-free ground, σ Nθ . These combine to give the total uncertainty as For the CRNS, Equation (1) gives while for the buried-CRNS, where C = − 157 mm is the gradient ∂SWE/∂N * estimated using the 0-30 mm portion of the attenuation curve given in Howat et al. (2018). In Section 6.1, uncertainties given by Equation (14) will be equated with those from the triple collocation analysis so the causes and effects of the uncertainties can be discussed.

| ANALYSIS OF SWE ESTIMATES
6.1 | Comparison of SWE estimates and uncertainties Figure 9 shows the various SWE estimates for the example snow events at Easter Bush and Moor House. For the CRNS and buried-CRNS SWE estimates, Method 3 is used to estimate the snow-free count rate (see Section 3.1 and Section 4, respectively). Using Method 1 or Method 2 produces only limited differences (e.g., see snow-free count rates in Figures 2 and 6).
For both events shown in Figure 9, all four SWE estimates display quantitatively similar behaviour. The biggest exception to this is an apparently faster melting of the depth-based SWE estimate compared to all other SWE estimates, clearly visible in Figure 9a for the Easter Bush event after 4 March. This appears to be caused by the snowpack consolidating and becoming increasingly dense as it ages. This effect is not accounted for in the depth-based SWE estimate and led to the exclusion of this estimate for days at the end of snow events in which the depth has decreased by more than 20% from its maximum, as shown by the large hollow circles in Figure 9. This simple criteria does not work perfectly for complex snow events such as that shown in Figure 9b where the depth-based SWE estimates for 8-9 February have suffered considerably larger melting than all other estimates, and should be regarded with suspicion.
Scatter plots comparing daily CRNS/buried-CRNS SWE estimates to depth-based and snowmelt model estimates show a considerable degree of scatter as reflected by correlation coefficients of less than 1 (see Figure 10). A similar degree of scatter is also found comparing depth-based with snowmelt model estimates (see Figure 11a) or comparing CRNS with buried-CRNS estimates (see Figure 11b). Some of the scatter is caused by differences in the footprints for the various To understand the level of absolute uncertainty involved with each SWE estimate, the method of triple collocation (Section 5.3) was used to produce Table 2 which includes results for all three SWE estimation methods for the CRNS and buried-CRNS. Caution should be employed when interpreting this table. In particular, the requirement to have data for all three SWE estimates in order to use triple collocation, limits the analysis to just 162 snow days from 63 separate snow events with the availability of the depth-based estimate being the main limitation (compared to 874 days with CRNS-based SWE estimates across all sites, or 489 days with buried-CRNS SWE estimates across the snow sites). The uncertainty caused by this sample size was assessed using the bootstrap approach with 10,000 event-wise resamplings (Efron, 1979;Wilks, 2020) and is represented using a '±' for certain columns in Table 2. Additionally, the triple collocation analysis rarely uses data from near the end of snow events as depth-based SWE estimates are usually unavailable here. For the CRNS or buried-CRNS, this may result in an underestimation of the importance of correcting the snow-free count rate for changes in the soil moisture (Methods 2 and 3). Similarly, the uncertainty estimate quoted for the snowmelt model is primarily an estimation of the accuracy of the input precipitation.
Summarizing Table 2, triple collocation uncertainties for both the CRNS and buried-CRNS SWE estimates are σ y ≲ 4 mm, while for the depth-based and snowmelt model estimates the uncertainties are σ x ≈ 4 mm and σ z ≈ 5 mm, respectively. Compared to the depth-based estimate, there is a significant over-estimation of SWE for the buried-CRNS (a y ≈ 1.7), while a much smaller bias (a z ≈ 1.2) is found for the snowmelt model. For the CRNS estimate the bias is negligible. However, note that the use of Equation (6) to calculate N wat amounts to a partial calibration of the CRNS estimate to the depth-based estimate.
Also, note that if the buried-CRNS or snowmelt model SWE estimate were corrected for its bias, then according to the definitions in Equation (12), its uncertainty in Table 1 would be reduced by the same factor.
The uncertainty divided by the mean SWE recorded by each estimate for days included in the triple collocation analysis (the relative uncertainty) is also given in Table 2. According to this measure, the best performance is recorded by the buried-CRNS (δ y ≈ 0.3), followed by the CRNS (δ y ≈ 0.4), and then the depth-based and snowmelt model estimates (δ x ≈ 0.45 and δ z ≈ 0.5, respectively). Nevertheless, even for the CRNS and buried-CRNS the relative uncertainties are high. This appears to be a consequence of the low mean SWE depths, rather than high absolute errors. For the CRNS and buried-CRNS in particular, the absolute uncertainty is expected to increase only weakly with SWE depth for shallow SWE due to the almost linear neutron response for SWE ≲ Λ. As a result, the relative uncertainty would be expected to decrease almost inversely with increasing SWE (while SWE ≲ Λ).
For both the CRNS and buried-CRNS estimates, the choice of method for estimating the snow-free count rate has only a limited effect on the triple collocation uncertainties (see Table 2). In both cases, σ y indicates that Method 1 appears to be most accurate, followed by Method 3 and then Method 2. However in all cases the differences between σ y for each method fall within the overlap of their bootstrap sampling uncertainties. Additionally, the triple collocation analysis does not use data from near the end of the snow event, where an influx of meltwater into the soil increases the importance of correcting for the soil moisture (Methods 2 and 3).
For the CRNS and buried-CRNS the triple collocation uncertainty can be apportioned between statistical uncertainty due to the finite count rate, σ N , and the uncertainty in the estimated snow-free count rate, σ Nθ , using Equation (14). As σ N = ffiffiffiffiffiffiffiffiffiffiffi ffi N=24 p is known (the factor 24 accounts for daily averaging), σ Nθ can be estimated so that uncertainty predicted by Equation (14) reproduces that given in Table 2 when averaged over snow days included in the triple collocation analysis. For the CRNS (Method 3) this gives σ Nθ = 12:4 AE 2:4 where the uncertainties in this value reflect the bootstrap uncertainties in Table 2. For comparison, the average statistical uncertainty is σ N = 7.4. This suggests that uncertainty in the snow-free count rate accounts for the largest part of the combined uncertainty: an average of σ y = 2.7 ± 0.5 mm (carrying through the uncertainty in σ Nθ ), compared to σ y = 2.0 mm for σ N only. Also note that the CRNS uncertainty has a strong increasing dependence on the soil moisture (and to a lesser extent the SWE depth), which occurs because both N and N θ approach N wat in the denominators of Equation (15)  an average statistical uncertainty of σ N = 5.6. Again, the largest portion of the uncertainty is ascribed to uncertainty in the snow-free counts (an average of σ y = 3.6 ± 1.3 mm compared to σ y = 1.2 mm for σ N only). However, compared to the CRNS, the buried-CRNS uncertainty displays a much weaker sensitivity to the soil moisture because of the weaker dependencies seen in Equation (16).
For the CRNS and buried-CRNS SWE, the typical size of the discontinuity at the end of the snow event (which occurs when the albedo shows snow to be absent but the CRNS/buried-CRNS SWE estimates may be non-zero), provides a measure of the uncertainty at the end of the event. The mean discontinuity is given in This could be caused by an increasing uncertainty in the snow-free count rate, which is based on the count rate from before the start of the event. A very large increase in the uncertainties would, however, result in significant decreases in the correlation coefficients calculated using data from the end of the events (Figures 10 and 11). Additionally, and in contrast to the triple collocation uncertainties, the mean discontinuities are found to be marginally smaller when the estimate of the snow-free count rate is corrected for variations in soil moisture (Methods 2 and 3). Taken together, this suggests Method 3 may be (marginally) the more accurate.
The discontinuity in the CRNS or buried-CRNS SWE estimate can also be used to test the utility of Equation (14) in identifying snow events with particularly reliable/unreliable SWE estimates. For the F I G U R E 1 0 Comparison of daily CRNS (top row) and buried-CRNS (bottom row) Method 3 SWE estimates with depth-based (left column) and snowmelt model (right column) estimates. Dark/light blue points mark days from the 'start'/'end' of snow events at snow sites (depth not/has decreased by more than 20%), and pink points mark days at the non-snow sites (top right only). Error bars for the CRNS/buried-CRNS use Equation (14), error bars for the depth-based and snowmelt model estimates use Table 2, a green line displays the relative biases (Table 2), and dashes show the 1 to 1 line. Correlation coefficients are listed for various subsets of data at snow sites ('start, snow sites', 'end, snow sites' or 'all, snow sites'), or using all data for all sites ('all, all sites') CRNS (Method 3), selecting events with either lower or higher than median uncertainty at the end of the event (as calculated using Equation (14)) results in the mean discontinuity decreasing or increasing to 3.57 mm or 5.94 mm, respectively. In contrast, a much weaker effect is seen for the buried-CRNS, probably due to the weaker soil moisture dependence given in Equation (16). In this case, selecting events with either lower or higher than median uncertainty, results in only an increase/decrease of the mean discontinuity to 4.19 mm or 4.67 mm, respectively.

| Distribution of SWE across the UK
The total number of recorded snow days (as detected using albedo), by day of year across all sites and all winters, is shown in Figure 12.
Periods of snow affecting large numbers of sites can be discerned as peaks on the histogram. This is especially the case for winter 2017/18, indicated by the pink portion of the histogram bars, and accounting for 55% of all recorded snow days over the five winters. The two example snow events presented in Figure 9 are  T A B L E 2 Triple collocation bias parameters (a y , a z ), absolute uncertainties (σ x , σ y , σ z ), and relative uncertainties (δ x , δ y , δ z ) for SWE estimations based on the depth, x, the CRNS/ buried-CRNS, y, or the snowmelt model, z Note: The mean remaining CRNS/buried-CRNS SWE, Δy, at the end of the snow events also provides a measure of the uncertainty. Sampling uncertainties for a y , σ y , δ y , and Δy are standard deviations of 10 000 event-wise resamplings of the original data within the bootstrap approach. Figure 13 shows the annual distribution of accumulated snowfall across the UK, approximated by summing increases in the daily SWE for the various estimation methods. Note that the buried-CRNS and depth-based estimates are only available for the seven snow-sites (named on Figure 1); also that omission of accumulated snowfall due to excessive missing or removed SWE data is particularly restrictive for the depth-based estimate. This is signified using a cross (×) on Figure 13 if this occurs for greater than 30% of snow days. Increases in the depth-based SWE from the ends of snow events are included here, in contrast to the method outlined in Section 5.1 and used elsewhere in this paper.
The expansion of the COSMOS-UK network is evident in Figure 13, with sites added over the five winters and only one removed (the Redmere site in East England). For the three-winter average CRNS and snowmelt model estimates, greater snowfall accumulations generally occur over Scotland, Northern England and Wales, especially for those sites located at relatively higher elevations (compare to Figure 1). Also recall that no manual quality-control checks were applied to either the CRNS or buried-CRNS SWE estimates in order to assess their accuracy in an operationally relevant setting. This contrasts with the depth-based and snowmelt model estimates which did undergo quality-control checks, including manual comparison to phenocam images, in order to improve their accuracy for validation purposes.
Relative to the other estimates, a bias has been found in the buried-CRNS SWE estimates which could be corrected, although it should be recalled that there is no guarantee that the other estimates are bias-free.
Uncertainties in the SWE estimates from both the CRNS and buried-CRNS sensors are primarily attributed to variation in the underlying soil moisture and, to a lesser extent, statistical uncertainty in the 24 h average neutron count rate. For the CRNS SWE estimate in particular, the uncertainty in SWE is found to increase strongly with soil moisture, which allows the identification of snow events for which the CRNS-based estimate is either particularly reliable or particularly unreliable. In addition, there may be systematic errors in the method.
For example, the CRNS method for estimating SWE from the neutron counts is based on the relation from Desilets (2017) (see Equation (1)) which may suffer from inaccuracy in either the attenuation length, Λ = 48 mm, or the count rate over infinite depth of water, N wat . A measurement campaign in which SWE depths are assessed at distances representative of the CRNS footprint for several sites and for several snow periods would help expose and correct potential biases.
Additionally, the exponential dependence assumed in Equation (1) may be overly simplistic, which may become more apparent at larger SWE depths. One consequence is that Equation (1) predicts that the reduction in count rate will saturate for SWE depths greater than a few times Λ, while Schattan et al. (2017) reported no complete saturation for SWE depths up to 600 mm. However, this is not a major concern for the relatively shallow snow depths measured here (the maximum CRNS SWE estimate is 36.9 mm).
The CRNS approach to estimating SWE has several attractive qualities for hydrological monitoring, process understanding and modelling, and for verification of remote sensing snow products. The involving not only the water content of the soil but that of the overlying snow using the methods developed herein.