Sensing Area‐Average Snow Water Equivalent with Cosmic‐Ray Neutrons: The Influence of Fractional Snow Cover

Cosmic‐ray neutron sensing (CRNS) is a promising non‐invasive technique to estimate snow water equivalent (SWE) over large areas. In contrast to preliminary studies focusing on shallow snow conditions (SWE < 130 mm), more recently the method was shown experimentally to be sensitive also to deeper snowpacks providing the basis for its use at mountain experimental sites. However, hysteretic neutron response has been observed for complex snow cover including patchy snow‐free areas. In the present study we aimed to understand and support the experimental findings using a comprehensive neutron modeling approach. Several simulations have been set up in order to disentangle the effect on the signal of different land surface characteristics and to reproduce multiple observations during periods of snow melt and accumulation. To represent the actual land surface heterogeneity and the complex snow cover, the model used data from terrestrial laser scanning. The results show that the model was able to accurately reproduce the CRNS signal and particularly the hysteresis effect during accumulation and melting periods. Moreover, the sensor footprint was found to be anisotropic and affected by the spatial distribution of liquid water and snow as well as by the topography of the nearby mountains. Under fully snow‐covered conditions the CRNS is able to accurately estimate SWE without prior knowledge about snow density profiles or other spatial anomalies. These results provide new insights into the characteristics of the detected neutron signal in complex terrain and support the use of CRNS for long‐term snow monitoring in high elevated mountain environments.


Introduction
In regions relying on snow-fed mountain rivers, snow water equivalent (SWE) is a fundamental environmental variable for the management of water resources (Clark et al., 2011;Sturm et al., 2017;Viviroli et al., 2007). Mountain snowpacks are characterized by a complex vertical structure (Sturm et al., 1995) as well as high spatial heterogeneity due to the combined effects of vegetation, solar radiation, and wind-induced snow redistribution Sturm et al., 1995;Webster et al., 2015). In particular, the small-scale variability of SWE can be significant while the corresponding spatial patterns of SWE distribution are often unstable over space and time (Jost et al., 2007). Since conventional techniques for SWE estimation exhibit a small measurement footprint (Goodison et al., 1987;Kinar & Pomeroy, 2015), the complex spatial snow patterns can considerably hamper the determination of representative SWE (Grünewald & Lehning, 2015). Therefore, SWE monitoring in mountain regions is challenging (Dozier et al., 2016), and the best results could be obtained by a combination of snow-hydrological modeling, remote sensing, and in situ data (Sturm, 2015).
In contrast to single SWE measurements that are affected by small-scale variability, data aggregated over distances of up to 400 m have been found to correlate well with terrain parameters (Grünewald et al., 2013;Helfricht et al., 2014;Jost et al., 2007). This is of high relevance for the calibration and validation of snow-hydrological models, as most modeling approaches dealing with snow redistribution are designed to capture processes at an intermediate scale of hundreds of meters, neglecting the sub-grid variability (Freudiger et al., 2017). Hence, differences between single point measurements and modeled SWE values can be substantial, even when the mean value of a grid cell is captured accurately by the model A promising technique to close the scale gap between in situ measurements, hydrological models, and remote sensing is cosmic-ray neutron sensing (CRNS). It was originally introduced for monitoring area-averaged soil moisture dynamics that are insensitive to small-scale heterogeneity (Zreda et al., 2008. Today, CRNS comprises an established method for intermediate scale soil moisture monitoring (Mohanty et al., 2017). Applications range from agricultural fields (Baroni & Oswald, 2015;Franz et al., 2015;Rivera Villarreyes et al., 2011), over measuring soil moisture in forested areas (Baatz et al., 2015;Heidbüchel et al., 2016;Nguyen et al., 2017), to the determination of the water content of frozen soils (Zhu et al., 2016). Operational CRNS networks for measuring soil moisture exist in the United States , the United Kingdom (Evans et al., 2016), and Australia .
In addition to soil moisture, CRNS is also sensitive to other landscape elements containing hydrogen-rich molecules like vegetation or snow (Desilets et al., 2010). While approaches to separate the signals of soil moisture and biomass are emerging (Baroni & Oswald, 2015;Baatz et al., 2015;Jakobi et al., 2018), the knowledge about using CRNS for measuring SWE is still limited. For high-energy neutrons it is known that detectors, such as neutron monitors, are not only influenced by snow on the roof (Korotkov et al., 2011(Korotkov et al., , 2013 but also by the snow in the surrounding area (Cheminet et al., 2014;Eroshenko et al., 2010). Neutron spectrography with Bonner spheres confirmed that, below 20 MeV, the cosmic-ray induced neutron flux is highly influenced by the amount of SWE in the surrounding area (Rühm et al., 2012). The sensitive energy range of CRNS detectors lies approximately between 1 eV and 100 keV , and thus, the signal should be even more sensitive to the hydrogen in the water molecules of the snowpack. First applications of monitoring SWE using CRNS in the United States (Desilets et al., 2010;Rasmussen et al., 2012), Northern Germany (Rivera Villarreyes et al., 2011), China (Tian et al., 2016), and Canada (Sigouin & Si, 2016) confirmed the sensitivity of CRNS detectors to snow, but measurements were limited to shallow, rather uniform snowpacks. And initially, a measurement limit of around 100 mm of SWE was presumed (Desilets et al., 2010;Desilets, 2014). In contrast, a recent study in the Austrian Alps showed that the CRNS signal could be related quantitatively to SWE of several hundreds of mm, and there is no evidence for a complete signal saturation even for SWE values of up to 600 mm (Schattan et al., 2017). A major finding of this study was that the relationship between neutron counts and SWE differs considerably between the accumulation and the melting season, respectively (Desilets et al., 2010;Schattan et al., 2017).
Simulating the response of the neutron flux to different boundary conditions can give important insights into the characteristics of the CRNS signal. In principle, it is well-known since decades that the signal of a neutron detector is influenced by the hydrogen content of its surrounding area (Hendrick & Edge, 1966). But only Monte Carlo neutron transport simulations were able to prove that a cosmic-ray neutron detector could be used for assessing soil moisture at an intermediate spatial scale of several hectares (Zreda et al., 2008) and to quantify the footprint of the signal with regard to atmospheric pressure and detector height (Desilets & Zreda, 2013). More recently, neutron transport modeling has proven to be a valuable tool for understanding the CRNS response to more complex environments and to spatial heterogeneity in the footprint. For instance, it was used to reproduce field measurements in a forest site (Andreasen et al., 2016(Andreasen et al., , 2017, to quantify variations in footprint size caused by changes in moisture content in soil and atmosphere (Köhli et al., 2015), to understand the effects of land surface heterogeneity in the CRNS footprint (Köhli et al., 2015;, and even to generate spatial maps of the neutron response in a whole landscape . These neutron transport simulations have shown that the sensitivity of the sensor is non-linearly dependent on the distance between the surface and the detector and is affected by other hydrogen sources like atmospheric moisture or vegetation (Köhli et al., 2015;Schrön et al., 2017). Accounting for this spatial sensitivity during the calibration and validation of CRNS substantially improves its performance (Cai et al., 2018;Heidbüchel et al., 2016;Schattan et al., 2017;Schrön et al., 2017). While CRNS has been generally considered insensitive to local anomalies, complex spatial structures resulting in large local differences in signal attenuation, like a heterogeneous mix of soil, roads and buildings in urban areas, can also influence the signal . Measurements of soil moisture within a forest stand could only be reproduced by neutron transport modeling, when considering the spatial structure of the forest (Andreasen et al., 2017). In mobile applications,  found that the dry road itself can introduce a substantial bias to CRNS-derived soil moisture compared to off-road measurements on agricultural land. Thus, the effect of a spatial heterogeneity on the CRNS signal constitutes a possible explanation for the empirically observed increase in neutron count rates during snow melt.
In the present work, we aim to advance the understanding of the interaction between epithermal neutron and land surface based on Monte Carlo neutron simulation and to support the use of CRNS approach for snow estimation in complex snow pack conditions. In particular, a rigorous analysis using scenario-based neutron modeling and long empirical data sets is conducted to assess the application of CRNS for monitoring intermediate scale alpine water resources. The main scientific goal is to evaluate possible influencing factors and to enhance the validity of its interpretation. Therefore, step-wise simulation experiments have been conducted using the URANOS (Ultra RApid Neutron-Only Simulation) model (Köhli et al., 2015) to assess the effect of alpine topography, the heterogeneity of the snow-free land surface, snow density profiles, and spatially non-uniform snow patterns including fractional snow coverage. In particular, the simulations aim at addressing the following research questions: 1. Can neutron simulations reproduce field measurements of a mountain snowpack with up to 600 mm of SWE? 2. Does the vertical and horizontal heterogeneity of the snowpack influence the CRNS signal? 3. How large is the potential effect of the underlying soil moisture conditions? 4. Can the increase in neutron count rates causing a hysteresis of the neutron response between melting and accumulation periods be explained by the presence of snow-free patches?

Site Description
The domain used in the neutron transport simulations consists of an area of 1 km × 1 km centered at the Weisssee Snow Research Site in the Austrian Alps. This site was chosen as it represents an alpine environment with spatially and temporally heterogeneous mountain snowpacks. The Weisssee station (Figure 1b) is located in a side valley of the upper Kaunertal in level area with sandy deposits at 2,464 m a.s.l. and is surrounded by mountain peaks (see Figure 1e). The southern part of the valley is occupied by the Weissseeferner glacier with a ski resort and opens in northward direction. During little ice age, the glacier terminus was within the model domain (see Figure 1a).
The site shows a typical alpine snowpack evolution with snow accumulation starting in October to November, peak accumulation in terms of SWE around April, and depleting snowpack in summer (Schattan et al., 2017). The variability of the height of snow (HS) is affected by preferential deposition of snow, wind-induced snow redistribution, and differences in solar irradiance. As a result, the mean values in individual radial sectors of a 250 m radius around the weather station differ from mean SWE in the same area by around 0.1 m in winter and up to 0.3 m in spring, with differing spatial patterns particularly when the winter seasons of 2014/2015 and 2015/2016 are compared (Fey et al., 2019a). In contrast, the variability of snow density in the area around the weather station is considerably lower (Schattan et al., 2017). Though strongly affecting the representativeness of the point-scale measurements at the snow research site, the inter-and intra-annual differences in the snow accumulation patterns allow for analyzing the CRNS response to contrasting natural snow distributions.
The inner 500m × 500m area of the model domain, fully covering the empirically derived footprint of 230 m (Schattan et al., 2017), is shown in a 3D plot in Figure 2c. Within this area, the elevation ranges between 2,417 in the north and 2,556 m a.s.l in the southeast. The mean elevation of 2,464 m a.s.l. matches exactly with the elevation of the snow research site. The inner domain is characterized by a moderately steep mountain environment with mean slopes of about 20 • , comparable to the neutron simulation experiments conducted by Köhli et al. (2015). The outer part of the model domain was additionally included as it is considerably steeper and possibly might affect the incoming neutron flux. Furthermore, this larger domain allows for simulating neutrons originating in larger distances from the detector. Due to the steeper terrain in the outer part of the model domain, the elevation range is considerably larger for the entire area (2,349 to 2,675 m a.s.l., Figure 2a) than for the inner model domain.
The natural landscape is dominated by glacial deposits and alpine vegetation cover. In close vicinity of the station the terrain is relatively flat but characterized by moraine material accumulated by the Weissseeferner glacier during little ice age (see Figure 1). Soil pockets with a high skeleton fraction coexist with a notable number of ferrous boulders ( Figure 1c). Vegetation cover is generally low and changes with soil stability, soil type. and soil water availability from bare ground to alpine meadows. These differences are clearly visible in Figure 2b. A Sentinel-2A scene (European Space Agency (ESA), 2019) during the peak of the vegetation season on 3 August 2015 was used to calculate the normalized difference vegetation index (NDVI). For that, the bands 4 (central wavelength 664.5 nm, bandwidth 38 nm) and 8 (central wavelength 835.1 nm, bandwidth 145 nm) were used in their native spatial resolution of 10m × 10m. Though being primarily used for characterizing vegetation, NDVI is also suitable for detecting fractional vegetation cover as the spectral differences between bare soil or rock and vegetation cover are very high (Huete & Jackson, 1988;Ormsby et al., 1987). Areas with bare ground, as indicated by low NDVI values, are mostly found in steep slopes with unconsolidated moraine material in the eastern part of the scene but also in the rocky slopes south of lake Weisssee and along the road. Furthermore, plant productivity and thus NDVI are physiologically related to water availability. Therefore, in areas with comparable plant composition NDVI can be used to map spatial patterns of soil moisture (Engstrom et al., 2008;Grzywna et al., 2018). Highest NDVI values are found in areas with more stable soils covered by alpine meadows, predominantly in the moderately steep slopes in the northwestern part of the scene and in the concave area upslope of the lateral moraine in the northeastern ridge of the scene. In the center, intermediate NDVI values are associated with a mixture of boulders and meadow-covered soil pockets. Although the majority of the areas with low NDVI consist of surfaces with no or shallow soils, the effects of topography on the spatial distribution of soil moisture might not be fully reflected in these areas. Further landscape features include surface water bodies and artificial landscape elements. In addition to the natural lake Weisssee, surface water bodies exist along two creeks flowing northward into the Fagge River. Artificial structures are formed by a cable car station and the paved Kaunertal glacier road.

Available Field Observations 2.2.1. Cosmic-Ray Neutron Sensing
Dedicated to measure the CRNS response to high snow accumulations, the Weisssee Snow Research Site was equipped with a 3 He gas filled CRS-1000 neutron detector (Hydroinnova LLC, USA), similar to the standard probes of the COSMOS network . Continuous CRNS measurements exist from March 2014 to June 2014 and from October 2014 to June 2016 (Schattan et al., 2017). The raw data are available from Schattan et al. (2019). Based on historical snow depth measurements, the detector was installed 2.7 m above snow-free ground to ensure that the installation is always above the snowpack. With increasing snow accumulation the distance between the surface and the detector diminished to up to only 0.4 m during the highest snow accumulation in April 2015. The neutron detector is located at the weather station in the very center of the domain.
High density polyethylene (HDPE) is commonly used as a moderator to reduce the response of neutron detectors (typically 3 He or BF 3 ) to thermal neutrons and to extend the sensitivity toward epithermal and fast neutron energies (Hutcheson et al., 2017;Woolf et al., 2019). Similarly, the CRS-1000 detector is moderated by 2.5 cm of HDPE to improve its response to neutrons in the hydrogen-sensitive energy range yielding the highest detection efficiency between 10 2 and 10 4 eV (Desilets et al., 2010;Köhli et al., 2018). Still, a contribution from slower thermal neutrons of up to 30% is very likely (Andreasen et al., 2016;Hutcheson et al., 2017;Köhli et al., 2018;McJannet et al., 2014). Thus, for comparing measured neutron count data with simulation results, actual energy response of the detector as modeled by Köhli et al. (2018) is used.
The measured CRNS data are corrected for the effects of incoming neutron flux, barometric pressure, and atmospheric humidity according to the functions reported in the literature (Rosolem et al., 2014;. Data from the Jungfraujoch neutron monitor data have been used as the incoming cosmic-ray flux. Reference barometric pressure and atmospheric humidity are set to the mean local conditions of 750 mbar and 1 g m −3 , respectively.

Meteorological and Snow Related Data
Being part of the regional hydrometeorological network, the site is equipped with standard meteorological instruments and sensors for measuring snow state variables. Snow density is measured continuously using a SnowPackAnalyzer (Stähli et al., 2004). In addition, snow pits were dug close to the weather station to manually measure snow density. Snow density data are used for processing distributed SWE maps, while meteorological data are needed for processing CRNS data. The data for the period 1 October 2014 to 30 September 2018 is provided by Schöber et al. (2019).
Distributed snow depth maps based on terrestrial laser scanning (TLS) are available from Fey et al. (2019b). Within the period covered by CRNS data, a total of 17 TLS-based snow maps in the winter seasons 2014/2015 and 2015/2016 covering the center of the model domain exists (Fey et al., 2019a;Schattan et al., 2017). This data set includes a snow-free data acquisition, 10 snow maps with full snow cover, and seven snow maps during snow melt. The processing of the TLS data and the accuracy of the HS maps derived from difference surface models are described in Fey et al. (2019a). For details on the processing of the distributed HS maps into SWE maps, see Schattan et al. (2017). Although largely covering the inner 500 m × 500 m area of the model domain (marked as inner domain in Figure 2), data gaps still exist in the TLS-based SWE maps as exemplified for 28 April 2016 (peak of accumulation) and 7 June 2016 (melting period) in Figures 3a and 3c. In the inner domain this is mainly caused by obstacles in the line-of-flight. In the outer domain data gaps result predominantly from a limited viewshed, that is, the area in the line-of-sight from the position of the laserscanner, due to the mountainous topography. For individual dates, the viewshed was even smaller as high snow accumulations close to the scan positions hampered the TLS data acquisition (see Fey et al., 2019a).
A multi-linear regression based on proxy data is thus applied for every scene to produce gap-filled snow maps suitable as input for the neutron transport model. The regression includes a gap-free airborne laserscan acquisition of April 2010  to cover interannually persistent snow patterns. Preferential deposition and redistribution of snow due to gravity and wind transport is accounted for by including elevation, slope, and curvature of the terrain (Blöschl et al., 1991) and topographic openness (Hanzer et al., 2016) into the regression. As the detector is considerably more sensitive to the water pools in the area close to the detector, the domain was additionally subdivided into zones based on 16 sectors and 50 m distance bands around the neutron detector, see also section 2.4. If the regression coefficient is higher for a given

Soil Moisture Dynamics
The large number of boulders and the high skeleton fraction considerably hamper the measurement of soil moisture. Thus, the CRNS signal is also used to reconstruct the temporal dynamics of soil water storage. To calibrate the CRNS, in situ measurements were conducted during a dry-spell on 13 August 2015 with a portable TDR instrument (6050X1 TRASE System I, Soilmoisture Equipment Corp.). Within a 50 m radius around the detector, a radial spatial sampling scheme  was applied to account for the spatial sensitivity of the CRNS. The TDR measurements yielded an average soil volumetric moisture of 10%. Based on that, the parameter N 0 of the equation proposed by Desilets et al. (2010) was calibrated to local conditions. The absolute soil moisture values might be systematically biased toward wetter conditions as measurements were taken in the soil only. Still the temporal soil moisture dynamics of the site can be illustrated. The CRNS-based soil moisture is mostly below 20% and ranges between 8% and 25% during the snow-free period (Figure 4). Nominally higher values are caused by the presence of snow on the ground resulting in a contribution of SWE to the signal. Overall, this indicates a rather low water holding capacity of the soils in the area. Throughout the paper, volumetric soil moisture values refer to the overall percentage in respect to the total soil volume.

Neutron Transport Modeling
The Monte Carlo tool URANOS, freely online available at http://ufz.de/uranos/ (Köhli et al., 2015), is designed specifically for modeling neutrons interacting with the environment in the frame of CRNS. The standard calculation routine features a ray-casting algorithm for a single neutron propagation and a voxel engine. The physics model follows the implementation declared by the ENDF (Evaluated Nuclear Data File, see also Brown et al., 2018) database standard and as was described by OpenMC (Romano & Forget, 2013). It features the treatment of elastic collisions in the thermal and epithermal regime, as well as inelastic collisions, absorption, and emission processes such as evaporation. Cross sections, energy distributions, and angular distributions were taken from the databases ENDF/B-VII.1 (Chadwick et al., 2011) and JENDL/HE-2007 (Shibata et al., 2011). For energies above 20 MeV an effective high-energy model is used which emulates the shower cone from empirical data, resulting in an attenuation length in water as reported in the literature (Korotkov et al., 2011(Korotkov et al., , 2013. Generally, instead of propagating particle showers in atmospheric cascades, URANOS reduces the computational effort and makes use of the analytically defined cosmic-ray neutron spectrum by Sato (Sato & Niita, 2006;Sato, 2015Sato, , 2016, which has been calculated by PARMA (Sato, 2016). Specifically important is the voxel engine, which extrudes the whole domain to a stack of 3D pixels, each of which can contain materials like soil, rock, or snow with variable densities. This allows to directly transfer the snow distribution measured by the 3D laser scanner to an input file for the geometry definition. In this work the URANOS model version 0.99 was used.

Simulation Setups
A step-wise simulation experiment using the URANOS model is set up to better understand the interaction between the signal and land surface conditions and to evaluate influencing factors possibly affecting the application of CRNS in an alpine environment. Similar to a sensitivity analysis, the factor modified between two scenarios can be tested in its influence on the functional relationship between neutron counts and snow water based on the magnitude of its impact in the simulations. In general, the simulations reflect the mean local conditions in terms of cutoff-rigidity (4.5 GeV) and barometric pressure (750 mbar). The atmospheric humidity is fixed at a value of 1 g m −3 . Taking some of the environmental conditions as constant not only reduces the number of influencing factors but also allows for a direct comparison with the measured neutron count data. This is due to the fact that the empirical data are corrected for variations in barometric pressure and atmospheric humidity with the same values as for the reference conditions. From this basic configuration being identical in all simulations, nine variants of the URANOS model setup are deducted ( Table 1). The setups basically differ in the assumptions regarding the soil and the snow layers and in the consideration of landscape features and topography. A soil moisture of 15% is assumed to represent mean conditions at this specific site according to the measurements described above. Nevertheless, on that, also higher and lower values are tested to reflect the uncertainty in determining absolute values of soil moisture and for assessing the sensitivity to soil moisture changes.
In a first set, setups 1 to 4 aim at analyzing the possible effect of alpine topography and variations of the snow-free land surface. These modeling scenarios thus have no snow layer. The lower the differences between these scenarios, the less the signal is affected by topographic shielding due to the surrounding mountain slopes and heterogeneity of the land surface. The first setup represents entirely homogeneous soil conditions with uniformly distributed soil moisture and without the consideration of landscape elements other than soil. In the second setup, the same uniform soil moisture distribution is used but including a number of landscape features described above (i.e., open water, the paved road, and the building) in the soil layer. This is implemented by defining material codes other than soil for the respective pixels of the soil layer. Additionally, soil moisture is distributed according to the relative variations in NDVI in the third setup. CRNS-derived soil moisture represents the water content of the surface layer including soil skeleton and rock fractions. Therefore, the variations of NDVI are used to represent spatial patterns of soil moisture and the presence of stones and boulders. Being otherwise identical to the third setup, only the fourth setup is explicitly accounting for topographical effects. For that, the 3D geometry of the soil layer is represented in 2 m elevation steps.
Setups 5 to 8, that is, set two, now with a snow layer, are aimed to test the overall effects of vertical and horizontal snow heterogeneity in fully snow-covered conditions. In setup 5 snow density is uniformly 100 kg m −3 . In contrast, a vertical snow density profile conceptualizing measured data from the peak  (300 and 500 kg m −3 , respectively), and a depth-hoar influenced layer of 0.3 m thickness with 300 kg m −3 density. While in setup 7 landscape features are included, but snow is distributed homogeneously; in setup 8 spatially non-correlated Gaussian noise with a standard deviation of 30 mm was added to the snow layer. Gaussian noise was selected to reflect that, as other variables in nature, SWE is distributed randomly around a mean value. The standard deviation constitutes a first guess based on the measured variability of SH during the early accumulation season. Snow-free pixels are not allowed and covered instead with small amounts of SWE. As this constraint affects the average SWE in the domain, SWE is calculated from the modified snow layer for comparison with setup 7.
Based on the knowledge gained from the first two sets, the additional effect of partly snow-free areas is evaluated against measured neutron count data. This constitutes the most complex setup 9 having measured snow distributions implemented. As in setup 3, the soil moisture is distributed according to NDVI values from the Sentinel-2A scene including the additional landscape features. The snow layer consists of gap-filled TLS data (see Figure 3), resulting in one realization for each of the 17 dates with TLS data. The detection layer representing the CRNS sensor is adapted according to the measured snow depth at the weather station, as neutron intensity decreases with elevation above the surface while thermal intensity decreasing faster than epithermal (Andreasen et al., 2017).
Finally, the changes in footprint are evaluated, as this is of high practical relevance for comparing CRNS data to remote sensing products or hydrological models. The CRNS footprint is defined as the area within which 86% of the detected neutrons originate (Desilets & Zreda, 2013;Köhli et al., 2015). It is relevant for relating SWE data derived from neutron sensing to other measured or simulated snow data. The spatial sensitivity of CRNS probes with respect to sensing distance is highly non-linear. With a weighting function exhibiting roughly a double-exponential shape, one can identify a local contribution from the near-field and a far-field contribution resulting from a long-range neutron transport. For selected simulation runs, the origins of the detected neutrons are tracked and analyzed to illustrate the effects of topography and contrasting snowpack conditions. The footprint is calculated based on the modeling domain of 1 km × 1 km, excluding all distances larger than 500 m to avoid artifacts due to the non-circular domain geometry. Based on sensitivity calculations with the footprint functions presented by Schrön et al. (2017), the exclusive consideration of neutrons originating within the simulated domain (500 m radius) is expected to bias the footprint result by around 8% as compared to using the full range of the function (up to 600 m in radius). Furthermore, the geometry of the simulated detector is defined as a spherical entity in the center of the model domain with a model-based response function. It differs from a real detector due to its larger extent as the main objective of the simulations is to investigate the influence of spatial heterogeneity on neutron flux intensities rather than footprint characteristics. Thus, in absolute numbers the resulting footprints are biased toward lower values but still give insights into differences related to spatial heterogeneity of the domain and the presence of snow.

Heterogeneity of the Snow-Free Landscape
The simulations of the different setups of the snow-free landscape are shown in Figure 5. For a better comparability, the data are scaled to the respective scenario minimum. The 3D scenario shows lower total neutron counts being caused by topographic shielding (Dunne et al., 1999) of parts of the neutron flux. Topographic shielding can also occur in the small scale due to obstacles like large boulders (Balco, 2014). However, focusing on CRNS with a footprint of several hectares, only the large-scale effects caused by elevation changes were considered. In addition to a decreased total neutron flux, also the attenuation length is altered (Dunne et al., 1999).
All in all, the response of the different scenarios to changing soil moisture conditions is similar. The largest differences are found between the first and the fourth setups, especially in very dry conditions (0% to 5% soil moisture). In entirely dry conditions, setup 1 is around 7.5% higher and setup 4 is 4.5% lower than setup 3. Setup 2, in contrast, is almost identical (+0.5%). Accordingly, the dynamics between dry and wet conditions are highest in the setup without landscape features and lowest in the 3D setup. The consideration of the lake and the road has thus a larger effect on the neutron response to changing soil moisture conditions than the topography. The magnitude of the difference is, however, site-specific and is most likely influenced by the large water body of the lake. Compared to the mean value for each soil moisture condition, the differences of the soil layer scenarios are around ±1% in the range between 10% and 25% soil moisture. Using setup 2 as a reference, the RMSD for the entire range from 0% to 25% soil moisture is 8.5% for setup 1, 5.0% for setup 3, and 5.5% for setup 4. In the range between 5% and 25% soil moisture the RMSD is reduced to 3.8%, 5.4%, and 1.9%, respectively. This means that the effect of the surrounding topography and the land surface heterogeneity on the functional relationships between neutron counts and soil moisture and/or snow is negligible for the complete range of soil moisture values observed at the site. Similarly, the impact of the uncertainty associated with the NDVI based method for distributing soil moisture values is likely very small.

Vertical and Horizontal Heterogeneity of the Snowpack
A fundamental question for using CRNS in hydrological applications is whether the technique is sensitive to local anomalies. The influence of a vertical and horizontal heterogeneity of the snowpack is illustrated based on setups 5 and 6 ( Figure 6a) and setups 7 and 8 ( Figure 6b). The snow density assumptions, as described in setups 5 and 6, have no effect on the simulated neutron response. The RMSD between the two realizations is 2.15% for SWE values between 30 and 250 mm. That is, neither the bulk snow density nor the typical horizontal profile of snow density alter the CRNS signal. Similar results have been found by Dunai et al. (2014) in an empirical experiment using neutron detectors covered by HDPE plates and various air gaps. Similarly, disturbances of the vertical SWE distribution (Figure 6b) would not alter the CRNS response under the assumption of full snow cover (RMSD in the range of SWE values between 30 and 250 mm: 0.13%). The difference for small SWE values is caused by the altered mean SWE value in setup 8 (as snow-free areas are not allowed, see above). Considering the effective SWE, both setups reflect the same functional relationship between SWE and neutron numbers.

TLS-Based Snow Cover Distribution Including Patchy Snow Cover
In addition to the theoretical snow distributions presented above, the results from setup 9 allow for assessing differences in the neutron flux due to natural snow patterns. Figure 7 shows maps of the modeled neutron fluxes for different situations. The data are normalized to the domain average for the highest snow accumulation on 9 April 2015 to illustrate relative changes of the neutron flux in comparison to a situation with the lowest neutron count rates (Figure 7b).
On that date (9 April 2015) the values are very similar throughout the domain. The only exception is the cleared road with around 50% higher values. The contrast between neutron fluxes in the air over snow and ground is very high (Delunel et al., 2014;Korotkov et al., 2011;Paquet et al., 2007). On average, the neutron flux doubled in snow-free conditions with local differences arising from the landscape features and the heterogeneous soil moisture distribution (Figure 7a). While the neutron flux is decreased over the lake, the highest values are found in areas with dry soil conditions such as the lateral moraine. Due to low snow accumulations and high solar irradiance the lateral moraine on the northeastern hillslope is one of the first areas to become snow free. On 28 April 2016 small snow-free patches on the lateral moraine locally increase the neutron flux. The affected area is not limited to the snow-free areas (Figure 3b) but expands to neighboring snow-covered regions (Figure 7c). This effect is considerably more pronounced during snow melt (Figure 7d). The widespread absence of snow on the southwestern and northeastern hillslopes results in areas with higher neutron fluxes. Moreover, with decreasing amounts of SWE the snow-covered part is more heterogeneous.
An important aspect common to all modeling scenarios is the influence of the surrounding area on the local neutron flux even for larger geographic features. This is fundamentally caused by the size and the shape of the footprint (Desilets & Zreda, 2013;Köhli et al., 2018). The neutron flux over linear features such as the road or the creeks is clearly affected by the presence or absence of snow. In snow-free conditions, the creeks lower the neutron flux locally but to a lesser degree than in the case of, for example, the lake or snow cover. For linear elements this effect is very local (Köhli et al., 2015) but can introduce biases if the CRNS instrument is located directly on a linear feature like a road . In the bare area the snow-free neutron counts are around 250% higher. The creek reduces this value to 200%. Despite its notable size, the neutron flux over the lake is around 150% in snow-free conditions, while in snow-covered conditions it is around 100%. The same features are observable for the snow-free areas. The road shows neutron fluxes of 150% in winter, while this value changes to 200% in the absence of snow. The lateral moraine illustrates both the increase of the neutron flux over small snow-covered patches (as on 28 April 2016) and the attenuation of the neutron flux due to the presence of snow in adjacent areas. A reduced neutron flux was also measured during a field campaign in an area with mixed conditions of snow cover and bare ground in northern Canada (Woolf et al., 2019). Furthermore, the magnitude of neutron flux change depends on the size of the snow-free area, similar to the effect of road width . The neutron flux over that slope increases from around 150% in late April 2016 to 200% in early June 2016 and is around 250% under entirely snow-free conditions. Figure 8 shows URANOS simulation results for all 17 TLS scenes in comparison to the measured data. To assess the effects of the assumption of a soil moisture dependency, results with 10%, 15%, and 20% soil moisture are shown (Figures 8b-8d). The plot depicts the neutron flux as simulated by the virtual detector in orange (according to the actual CRNS energy response modeled by Köhli et al., 2018) and additionally shows values for the epithermal energy range (in blue) and values based on the weighted mean of 30% thermal and 70% epithermal neutrons as proposed by McJannet et al. (2014) (in red). Comparing modeled neutron counts with measured data shows that RMSE values are lowest for the simulations with 10% soil moisture being 7.6% in the epithermal energy range, 5.6% for the virtual detector, and 6.9% for the weighted mean of 30% thermal and 70% epithermal neutrons. With 15% soil moisture, the RMSE values rise to 6.3%, 12.2%, and 11.7%, respectively. Even higher RMSE values are found when using 20% soil moisture, RMSE (13.6%,17.1%,and 15.9%). Confirming prior studies (Andreasen et al., 2016;McJannet et al., 2014), little differences lie in the URANOS results using the modeled detector sensitivity and the assumption of a thermal contribution of 30%. However, the data of the epithermal flux without thermal contribution are considerably higher under dry conditions. Neglecting the thermal contribution to the CRNS signal would lead to an overestimation of the snow-induced signal reduction.
The measured data (Figure 8a) feature a different functional relationship for the accumulation and melting seasons, respectively (Schattan et al., 2017). Similar effects were also reported from a site in the United States (Desilets et al., 2010). This behavior can be reproduced well by the TLS-based simulations. Vertical and horizontal heterogeneity of snow cover is shown above to have no influence on the functional relationship between neutron flux and SWE at all. The differences hence result from the presence of snow-free patches during snow melt. In consequence, the empirical correction function for the melting season (Schattan et al., 2017) is likely only valid for sites with similar development of snow-free patches and should not be transferred to other sites without modifications.
The simulations reproduce the empirical neutron count data well up to 400 mm of SWE. Above this threshold further reduction of the neutron flux found in the empirical data is not represented by the simulations. Similar simulations by Eroshenko et al. (2010) showed a decrease even up to 800 mm, predominantly in the thermal energy region, by modeling the response of the neutron flux in the presence of water on the ground using the FLUKA-2006 model, a multi-particle Monte Carlo code simulating the interactions of different particles including neutrons. This underlines that above 400 mm of SWE uncertainties regarding the slope of the function relating SWE and neutron counts are considerably larger than below this threshold. Similarly, empirical findings have shown that due to the steep slope of the function relating neutron counts to SWE there is a higher uncertainty in the CRNS-based SWE signal for SWE values above 400 mm (Schattan et al., 2017).
For the entire range of SWE values used for model application (22 to 548 mm of SWE) further sources of uncertainty include the assumption regarding the energy response of the detector and the soil moisture conditions. The dynamics of the response of the neutron flux to additional hydrogen in snow are highly sensitive to the assumed soil moisture conditions. Thus, the largest differences are found below 200 mm of SWE. At our site the best agreement with measured data is found for the simulations with 10% soil moisture. With 20% soil moisture the neutron flux is lower than measured, resulting in a steeper slope of the function relating neutron counts to SWE. This indicates even lower soil moisture conditions than reconstructed from measured CRNS data. Considering the difficulties in measuring absolute values of soil moisture at the site, essentially caused by the large amounts of boulders and stones, this is not surprising. As the neutron response is highly non-linear, another factor could be the presence of ferrous boulders (see Figure 1) contributing to higher neutron evaporation rates. Furthermore, the differences due to the soil moisture assumptions are most notable for low SWE values. This confirms empirical results showing that while shallow snowpacks are influenced by dynamics of the underlying soil moisture content, this effect is smaller for larger SWE values (Sigouin & Si, 2016;Schattan et al., 2017).

Footprint Anisotropy
The previous analysis has shown that snow-free patches in the environment can significantly contribute to the CRNS signal. In contrast to snow-covered areas, neutrons are able to leave the bare ground unhindered, leading to a higher fraction of detected neutrons from the direction of the snow-free patches. Following this argument, the angular sensitivity of the sensor cannot be expected to be isotropic in a highly heterogeneous environment. To understand from which direction and distance the detected neutrons come from, we have analyzed the angular footprint distribution, defined as the radius within which 86% of detected neutrons originated for 30 individual sectors with a 12 • opening angle. Figure 9 illustrates the angular footprint distribution for the snow-free simulations in setups 3 (2D) and 4 (3D) and simulations based on TLS data on a snowy and a patchy day (setup 9). The average footprint of the CRNS signal under snow-free conditions is 206 m for the 2D model and 213 m for the 3D variant. While the 3% difference is in the range of the model uncertainties (Köhli et al., 2015), the angular footprint distributions are significantly different. This supports the assumption of an influence of geometric effects (see also Dunne et al., 1999).
In the 3D simulation, more neutrons originated from the uphill terrain in the northwest and south, leading to a stronger influence of the low contribution from the lake. We suppose that the CRNS detector in the valley is more exposed to neutrons from mountainous slopes because they can reach it more directly. The path between origin and detection goes far above the ground, leading to less intermediate collisions with the soil and thus allows for the neutron to have traveled larger distances when reaching the sensor. This explanation is supported by the simulations of the fully covered and the patchy snow days shown in Figures 9c and 9d. Here, high footprint anisotropy is evident while the angular distribution is mainly skewed toward the road and high mountain slopes, where snowpack is either thin, patchy, or even absent. The average footprint for the fully covered scenario in April is 152 m, and the patchy snowpack in June yields an average footprint of 182 m. Although the absolute numbers are probably rather slightly underestimated (see section 2.4), they confirm the first approximations by Schattan et al. (2017), and they are also in line with the theory of the footprint dependency on water content found by Köhli et al. (2015).

Conclusions and Outlook
The present study represents the so far most complete analysis of how spatially heterogeneous neutron density distributions can be modeled and how the model can be related to measurements to understand the signal of CRNS in complex alpine terrain. The URANOS neutron transport model was used to simulate the neutron response to changing SWE patterns at the Weisssee Snow Research Site. The model was able to reproduce the CRNS observations for complex snow patterns based on 3D laser scanning data. In particular, the observed hysteresis effect during snow melt and accumulation found by Schattan et al. (2017) has been confirmed by neutron simulations.
The results demonstrate that under fully covered snow conditions the CRNS signal is only influenced by the total SWE while being insensitive to vertical and horizontal heterogeneity, such as different snow density profiles or local anomalies. In contrast, rocks and other snow-free areas have a substantial effect on the neutron count rate, particularly during snow melt. Snow-free patches can increase the neutron flux, and its intensity is governed by the size of the depleted area. Since this effect is probably site-specific, data on other locations should be corrected during the melting season based on the individual site-specific patterns and with the help of neutron simulations.
To achieve a good performance on reproducing observation results, it has been shown that it is important to model landscape features such as buildings, roads, and lakes and to account for thermal neutrons particularly under dry conditions. Furthermore, uncertainties in the neutron modeling have been identified with regard to the assumption of the soil moisture content. This underlines that in particular for shallow snowpacks, soil moisture dynamics can significantly alter the neutron response. Non-invasive techniques for measuring local soil moisture or an accurate estimation of the volumetric fractions of soil and rocks would be needed to improve the assessment of soil moisture dynamics. With a better in situ data base and a vertical and horizontal weighting procedure as proposed by Schrön et al. (2017), it would be feasible to refine the CRNS calibration for soil moisture measurements.
The footprint of the CRNS signal was found to be anisotropic, mainly due to the complex topography and the contribution of the nearby lake and road. The average footprint for SWE measurements at our site lies between 150 m for fully snow-covered scenarios and 210 m for the complete absence of snow. The work shows that further research is needed to fully understand the footprint sensitivity to terrain features, the detector geometry, and the different neutron energies involved. From a practical point of view, the results show that the CRNS signal is well understood, reliable, and applicable even in complex environments.