A universal calibration function for determination of soil moisture with cosmic-ray neutrons

A cosmic-ray soil moisture probe is usually calibrated locally using soil samples collected within its support volume. But such calibration may be difficult or impractical, for example when soil contains stones, in presence of bedrock outcrops, in urban environments, or when the probe is used as a rover. Here we use the neutron transport code MCNPx with observed soil chemistries and pore water distribution to derive a universal calibration function that can be used in such environments. Reasonable estimates of pore water content can be made from neutron intensity measurements and by using measurements of the other hydrogen pools (water vapor, soil lattice water, soil organic carbon, and biomass). Comparisons with independent soil moisture measurements at one cosmic-ray probe site and, separately, at 35 sites, show that the universal calibration function explains more than 79 % of the total variability within each dataset, permitting accurate isolation of the soil moisture signal from the measured neutron intensity signal. In addition the framework allows for any of the other hydrogen pools to be separated from the neutron intensity measurements, which may be useful for estimating changes in biomass, biomass water, or exchangeable water in complex environments.


Introduction
Understanding the exchange of water between the land surface and atmosphere is critical for accurate initialization of general circulation models (Koster et al., 2004;Wang et al., 2006), understanding energy and water fluxes , and thus making short-term weather predictions. However, accurate and exhaustive soil moisture datasets in space and time are difficult to obtain (Robinson et al., 2008), hindering progress in fundamental understanding of the land-atmosphere coupling (Jung et al., 2010;Seneviratne et al., 2010). The recently developed cosmicray soil moisture method (Zreda et al., 2008) and probe by Hydroinnova LLC, Albuquerque, NM, USA, allow for near surface soil moisture measurements (∼ 12 to 70 cm) at intermediate horizontal scales (∼ 35 ha) (Desilets et al., 2010). Fifty probes have been deployed around the continental USA as part of the COsmic-ray Soil Moisture Observing System (COSMOS) (Zreda et al., 2012; data available at http://cosmos.hwr.arizona.edu/), and other networks are being installed elsewhere.
Previous work has concentrated on the relation between the neutron intensity in air above the land surface and soil moisture. However, that intensity is influenced not only by soil moisture, but also by hydrogen in other reservoirs. Desilets et al. (2010) presented a theoretical calibration function that required at least one independent estimate of areaaverage soil moisture to define the free parameter, N 0 . Analysis of 35 different COSMOS site calibration datasets indicates that N 0 varies significantly from site to site and in time within the same site (Table S1) due to presence of other timevarying sources of hydrogen, such as fast-growing maize (Hornbuckle et al., 2012).
While it is possible to collect soil calibration datasets at some sites, at others it may be impossible (e.g., where rock outcrops are present, in urban environments, in inaccessible areas, etc.), impractical (e.g., for large-scale mobile surveys; Desilets et al., 2010), or difficult to obtain a representative area-average soil water content (e.g., for sites that contain significant amounts of large cobbles or stones). However, water is still exchanged between the land surface and atmosphere at these sites, and observations of area-average moisture are necessary to understand the transfer of mass, momentum, and energy in these systems. Such observations may be possible with the aid of a universal calibration function. In this work we aim to show that measurements of neutron intensity and other more easily measurable hydrogen pools (water vapor, soil lattice water, soil organic carbon, and above-ground biomass) are adequate to provide reasonable estimates of soil moisture through a universal calibration function.
Universal calibration functions and algorithms have been developed in the past, with some of them being transformative in terrestrial hydrology. For example, Knyazikhin et al. (1998) utilized a radiative transfer model in six different vegetation classes to derive a global leaf area index and fraction of absorbed photosynthetically active radiation from satellite-derived spectroradiometer measurements. Topp et al. (1980) found a single relationship between the measured dielectric constant and volumetric water content across a wide range of soil types using coaxial transmission lines.
In this work, we develop a universal calibration function for the cosmic-ray neutron probe. It accounts for several sources of time-varying hydrogen signals that may be present in the probe's support volume in order to expand the potential use of the probe to hitherto difficult sites and novel applications. We first develop the function using the neutron transport code MCNPx (Pelowitz, 2005), for two different cases: uniform variations of pore water in 50 different soil chemistries, and vertical variations in pore water in four different soil textures using a numerical solution to the 1-D Richards equation. We then test the validity of the function by using observed neutron data from 35 COSMOS sites (where we have full calibration datasets) with a wide range of conditions.

Relationship between low-energy neutron intensity and hydrogen
The intensity of low-energy neutrons in air above the land surface is controlled mainly by the number of atoms of hydrogen in the combined soil-air system (Zreda et al., 2008). Hydrogen at the surface is present in three forms: static (soil mineral structure, mostly constant in time), quasi-static (vegetation and soil organic carbon, possibly varying in time) and transient (water vapor, pore or ponded water, snow, all changing in time). The support volume of a cosmic-ray probe is governed by the average travel distance of neutrons in air near the land surface (Desilets et al., 2010). Because of the additional sources of hydrogen, it is difficult to isolate a single transient signal (here, soil moisture) from the convoluted detected signal (here, neutron intensity) without additional information, constraints or simplifications. The influence of additional sources of hydrogen is evident from the wide range of N 0 values computed for COSMOS sites using the Desilets et al. (2010) calibration function: where θ is the volumetric pore water content (m 3 m −3 ), N is the neutron counting rate/flux normalized to a reference atmospheric pressure and solar activity level and N 0 is the counting rate/flux over dry soil under the same reference conditions. Figure 1 illustrates an example of Eq.
(1) for one value of N 0 and includes uncertainty bounds for two count rate uncertainties of 50 and 100 counts per hour, cph. We note that the uncertainty in the pore water content grows with increasing water content given the nonlinearity in Eq. (1). Using 45 calibration datasets from 35 different COSMOS sites, we find N 0 varies significantly (mean = 2632 cph, st.dev. = 433 cph, min = 1892 cph, max = 3394 cph; Table S1). This parameter should be approximately constant when accounting for all hydrogen sources; therefore the large spread indicates that not all sources of hydrogen are taken into account.

A practical framework
Given the deficiency of Eq.
(1) to account for additional or time-varying hydrogen signals, we present a simplified but general framework to account for all hydrogen sources present. We assume that a monotonic relationship exists between observed neutrons and the amount of hydrogen present in the support volume. In order to isolate the different transient signals, we employ neutron intensity correction factors and assign average properties within the support volume. From neutron transport simulations, Zreda et al. (2008) found that 86 % of the neutron signal occurs within a 335-m radius and is nearly independent of soil moisture. They also found that the vertical extent of the neutron signal depends on soil moisture, ranging from 12 cm in wet soils (0.40 m 3 m −3 ) to 70 cm in dry ones (0 m 3 m −3 ).
In the air, the support volume is approximately a hemisphere with a radius of 335 m. The neutron signal is normalized to the same reference pressure, geomagnetic latitude, and the incident high-energy neutron intensity as summarized in Zreda et al. (2012) and implemented in the COS-MOS project (Level 2 data available at http://cosmos.hwr. arizona.edu/). In addition, the neutron signal is corrected for variations in atmospheric water vapor (Franz et al., 2012b;Rosolem et al., 2013): where CWV is the scaling factor for temporal changes in cosmic-ray intensity as a function of changes in atmospheric water vapor (N · CWV), ρ 0 v (g m −3 ) is the absolute humidity at the surface, and ρ ref v (g m −3 ) is the absolute humidity at the surface at a reference condition (here we use dry air, Estimates of absolute humidity can be made with surface measurements of air temperature, air pressure, and relative humidity following Rosolem et al. (2013) (see Table S1). We note that the support radius is a weak function of water vapor and that it may be reduced by ∼ 10 % for fully saturated air as compared to dry air (Rosolem et al., 2013). Here we correct the neutron intensity to dry air conditions and do not consider the small changes in the support radius.
In the subsurface, the support volume is a cylinder with a fixed radius of 335 m and a depth that varies with pore water content. In order to calculate the measurement depth and depth-weighted average of properties within the volume, we modified the framework outlined in Franz et al. (2012a). For uniform distributions of pore water, lattice water, soil organic carbon, and bulk density, Franz et al. (2012a) found the effective depth z * (cm) as where 5.8 (cm) represents the 86 % cumulative sensitivity depth of low-energy neutrons in liquid water, 0.0829 is controlled by the nuclear cross sections of SiO 2 , ρ bd is the dry bulk density of soil (g cm −3 ), ρ w is the density of liquid water assumed to be 1 (g cm −3 ), τ is the weight fraction of lattice water in the mineral grains and bound water, defined as the amount of water released at 1000 • C preceded by drying at 105 • C (g of water per g of dry minerals, also known as lattice water), and SOC is the soil organic carbon (g of water per g of dry minerals, estimated from stoichiometry using measurements of total soil carbon and CO 2 , SOC = TC − 12 44 CO 2 ). For this work, measurements of lattice water, total soil carbon, and CO 2 measurements were made on a ∼ 100 g composite sample collected at the study site and analyzed at Actlabs Inc. of Ontario, Canada. The 100 g composite sample was made by collecting ∼ 1 g of soil from 108 individual samples at 18 locations with the footprint (samples collected at 6 depths, 0-5, 5-10, 10-15, 15-20, 20-25, 25-30 cm, and 18 horizontal locations, bearings of 0 • , 60 • , 120 • , 180 • , 240 • , 300 • , and radii of 25, 75, and 200 m). As a slight modification to the original equation (Franz et al., 2012a), we have added the effects of soil organic carbon (Table S1) on the effective depth of the sensor by treating it similarly to lattice water following Zreda et al. (2012), as it may contain a significant amount of hydrogen molecules. Given the complexity of soil organic carbon molecules, here we assume that the weight percent of organic carbon is the same as lattice water based on the stoichiometry of lignin (C 9 H 10 O 2 ) and equivalent amount of hydrogen. We note that the effects of SOC and lattice water on the effective depth will be most pronounced for low soil water contents.
With the estimates of sensor support volume and average properties within that volume, the total mass in the system and the total number of moles of each element can be calculated (Table S1). Following the neutron correction factor for variations in atmospheric water vapor Eq. (2), we assume the atmosphere is composed of only nitrogen (79 % by weight) and oxygen (21 % by weight). We estimate wet above-ground biomass, AGB (kg m −2 ), from US Forest Service maps in the continental USA (http://webmap.ornl.gov/ biomass/biomass.html) and assume vegetation is only composed of water (60 % by weight) and cellulose (C 6 H 10 O 5 , 40 % by weight). We assume the subsurface is composed of solid grains (pure quartz, SiO 2 , plus lattice water and SOC lattice water equivalent) and pore water. With the estimates of volume, mass, and chemical composition, hydrogen molar fraction, hmf (mol mol −1 ), is where H i is the sum of hydrogen moles from lattice water H τ , soil organic carbon lattice water equivalent H SOC , pore water H θ , and vegetation H AGB inside the support volume, and A i is the sum of all moles from air NO, soil SiO 2 , lattice water H 2 O τ , soil organic carbon lattice water equivalent H 2 O SOC , pore water H 2 O θ , and above-ground biomass C 6 H 10 O 5 + H 2 O AGB inside the support volume. Table S1 presents a summary of all calculations from the 45 calibration datasets at 35 different COSMOS locations (some locations have more than one calibration dataset).

Modeled neutron intensity for variations in soil mineral chemistry
To simulate the transport of cosmic-ray particles through the atmosphere and shallow subsurface, we used the MCNPx model (Pelowitz, 2005), a general purpose Monte Carlo code that tracks the individual life history of a particle and subsequent particles as it interacts with matter. At 1 to 2 m above the surface, the fast neutron flux, N (energy range 10-100 eV), is tallied; it corresponds to approximately the same energy neutrons that are measured by the cosmic-ray neutron detector moderated or "fast" channel (Desilets, 2011). Desilets (2011) notes that the moderated channel measures a range of neutron energies with the median value near 10 eV, where detection frequencies drop off as 1/e for higher-energy neutrons. We note that MCNPx simulations of neutron flux of each binned energy level between 1 and 10 5 eV indicate nearly identical response functions to soil moisture changes, thus justifying the use of the energy range used here. The neutron transport simulations used soil mineral chemistries from 50 different soil types, 49 COSMOS sites mostly across the continental USA ranging from sand to clays, and the hypothetical case of SiO 2 (Table S2). Using these soils we simulated 12 different uniformly distributed pore water contents (Fig. 2). Despite large variations in soil chemistry and soil texture, the framework presented in Sect. 2.2 explained 99 % of the variation (Table 1) between relative neutron flux and hydrogen molar fraction within the sensor support volume using a two-term exponential function: N N s = 4.486 exp(−48.1·hmf)+4.195 exp(−6.181·hmf). (5) As opposed to dry soil, which was used in previous work (Desilets et al., 2010), neutron flux is normalized to the case of an infinitely deep layer of water beneath the sensor (i.e., no soil) because of site to site differences in lattice water and soil organic carbon. The modeled neutron intensity over SiO 2 is 8.5 times higher than that over water, which is similar to the theoretical value reported in Fig. 1 of (Hendrick and Edge, 1966). To compare modeled neutron fluxes to observed neutron counts, the parameter N s , which represents the cosmicray neutron count rate normalized to that over water, must be specified. This parameter may be sensor dependent, but can be easily specified by measurements over a large water body (> ∼ 500 m on all sides and deeper than 1 m) and by following the standard correction factors summarized in Zreda et al. (2012). Most notably, the results in Fig. 2 illustrate that all 49 observed soil chemistries collapse to the SiO 2 case when using hmf, where we can account for the hydrogen in the lattice water and soil organic carbon. This result justifies the use of Eq. (4) and its simple representation of soil chemistry consisting of only SiO 2 , lattice water, and organic carbon. With respect to future work using stationary or mobile cosmic-ray neutron probes, we recommend site-specific estimates of lattice water and soil organic carbon for most accurate results.  (Table S2). The horizontal axis in (b) is the ratio of the sum of all hydrogen moles that are present in the support volume (vegetation, pore water, soil mineral water, and soil organic carbon) and the sum of all moles from all elements that are present in air, vegetation, soil, and pore water. Note the combined modeled chemistry and pore water profile cases were used for the best curve fitting analysis (see Table 1 for individual data fits).

Modeled neutron intensity for variations in pore water
To model variations in neutron intensity due to variable pore water content profiles, we used 1-D solutions of the Richards equation as input for the MCNPx modeling. Full details of the experiment are reported in Franz et al. (2012a), who simulated one complete infiltration and drying cycle in four different soil textures: sand, sandy loam, silt, and silty clay loam. Fitting the same two-term exponential function Eq. (5) to the model results, the explained fraction of variance decreased from 99.5 to 96.5 % and the RMSE increased from 63.9 to 75.1 cph ( Fig. 3 and Table 1). The increased variation is due to the slight hysteresis that exists in the neutron intensity because of the vertical averaging of the sharp wetting front that exists during infiltration. The hysteresis is most pronounced in the coarser soils where sharp wetting fronts are strongest as indicated by the RMSE for different soil types in Table 1. We note that all uncertainties are expressed in terms of cph, given the nonlinearity in converting cph to pore water content (Eq. 1). However, Fig. 1 (5), N/N s = a · exp(b · hmf) + c · exp(d · hmf). b Assuming N s = 1000 cph to convert between counts count −1 to cph. c Coefficients used in all figures and observed fitting. used to estimate the equivalent uncertainty in terms of pore water content for different count rates and uncertainties.

Neutron intensity observations at Santa Rita experimental range
To further validate the universal calibration function, we compare observed neutron data with estimates of hydrogen molar fraction for five different volumetric calibration datasets (Table S1, COSMOS site number 11) and for continuous measurements from a distributed sensor network in Southern Arizona over a six-month period (details of the distributed sensor network observations are presented in Franz et al., 2012b). Continuous measurements were taken at depths of 10, 20, 30, 50, and 70 cm, in the same spatial pattern as the volumetric calibration samples. Using the five volumetric sample data points, we estimate an N s value of 1037 cph using Eq. (5). Then, using the independently distributed sensor network measurements, we compute an RMSE of 110.1 cph and R 2 of 0.854 (Fig. 4, Table 1). However, we note that there is a slight bias at the wet end for the independent distributed sensor network measurements. Because the shallowest soil moisture sensor was placed at 10 cm, the dynamics in the top 5 cm are likely not well captured, leading to a bias when shallow wetting fronts exist shortly after precipitation (Franz et al., 2012b). Fig. 3. MCNPx modeled neutron flux versus hydrogen molar fraction for four different soil textures undergoing one wetting and drying cycle. Each soil moisture profile was generated using a numerical solution to the 1-D Richards equation for a top boundary condition of a 2.54 cm rain event that lasted 24 h followed by a 2 mm day −1 potential evapotranspiration for 9 days, a free drainage lower boundary condition, initial condition set to field capacity, and a vertical resolution of 2 cm. Note the combined modeled chemistry and pore water profile cases were used for the best curve fitting analysis (see Table 1 for individual data fits).  Table S1 for calibration datasets.
The problem is less pronounced for the calibration datasets in which the volumetric soil water content was determined using the gravimetric method.  Table S1 for calibration datasets. Note the 95 % confidence intervals are shown by the dotted line; data inside the dotted line red ellipse show dry sites with large amounts of above-ground biomass (site numbers 43, 51, 52, 53, see Table S1), and data inside the solid line blue ellipse show sites with hydrogen molar fraction values above 0.23 (site numbers 30, 31, 35).

Neutron intensity observations at multiple sites
The analysis of 45 calibration datasets from 35 different cosmic-ray neutron probe locations shows that Eq. (5) reasonably describes the data with an RMSE of 199 cph and R 2 of 0.791 (Fig. 5a), over a wide range of soil moisture, soil bulk density, soil texture classes, lattice water, soil organic carbon, vegetation, and water vapor conditions. As a direct comparison with Eq. (1), we find the best fit between total soil water θ + ρ bd ρ w (τ + SOC) and the observed neutrons counts leads to an RMSE of 310.2 cph and a R 2 of 0.493 (Fig. 5b). By including the differences in AGB between sites, Eq. (5) reduces RMSE from 310.2 to 199 cph and increases R 2 from 0.493 to 0.791.
Looking at the residuals between Eq. (5) and observed values, we find two sets of data that account for a large portion of the remaining uncertainty. The first subset (data inside solid line blue ellipse) comprises three sites where the observed hydrogen molar fraction is greater than the liquid water case at ∼ 0.23, as the sites have a large amount of aboveground biomass and were relatively wet on the day of the calibration. Above this value of 0.23, the count rate becomes flat, as the neutron probe is no longer sensitive to hydrogen being gained or lost to the system. The second subset of data (data inside dotted line red ellipse) is from four sites where large forests grow in dry sandy soils supported by shallow water tables (southeastern USA and northern Michigan); it will be discussed in the next section.

Remaining uncertainties
Numerical and observational results show that the overall uncertainty using the cosmic-ray method to detect time-varying hydrogen signals is small for a range of expected conditions (Table 1). We found from an inter-comparison analysis of 45 calibration datasets that Eq. (5) does poorest for sites that have a large amount of above-ground biomass in low soil moisture environments. In this work, the vegetation is presented as an equivalent layer of cellulose and water. However, as described in Sect. 2.1, the distribution of hydrogen above the surface may need to be explicitly accounted for in the life history of a neutron. Given the strong relationship for the five volumetric sample sets at Santa Rita Experimental Range (Fig. 3) (Franz et al., 2012b) and other local calibration functions (Desilets et al., 2010;Rivera Villarreyes et al., 2011), site-specific calibrations will implicitly include vegetation effects on observed neutron counts. But when comparing various sites with different geometries, there exists a systematic uncertainty that we are not properly accounting for with the equivalent layer assumption. Moreover, we do not include the effects of biomass below the ground surface, which may be significant but is difficult to quantify. Future observational and theoretical work should aim to understand the effect on neutron intensity of biomass above and below ground surface. In terms of calibration datasets, accurate spatial estimates of volumetric water content may be difficult to obtain because of a large uncertainty in the determination of soil bulk density (Table S1) (Dane and Topp, 2002). Additional experimental work on accurate determinations of bulk density at intermediate spatial scales, and better quantification of above-and below-ground biomass and soil organic matter will help reduce the overall uncertainty in the cosmicray moisture measurements.

Outlook
With the presented framework, reasonably accurate estimates of soil moisture at difficult sampling sites or many points when using mobile surveys are now possible using cosmicray neutron probes. The proposed method requires ancillary data on location (for the computation of the incident cosmicray intensity), surface pressure, air temperature, and relative humidity to perform the necessary neutron intensity corrections. In addition, estimates of soil bulk density and total biomass are needed, but spatially contiguous data are readily available from various sources. The most challenging aspect is a spatial map of lattice water and soil organic carbon. Here we provide data from 50 sites (Table S1), but additional measurements along the mobile survey need to be made, or, alternatively, relationships between bulk density and more readily available properties (Greacen, 1981) need to be established.
In addition to computing pore water from the neutron intensity measurements, it is possible to isolate any one of the hydrogen signals, given measurements of the other pools (Fig. 6). Figure 6 illustrates the expected neutron count rates solving Eqs. (3)-(5) for different combinations of pore water, lattice water, and above-ground biomass. This method and framework has potential applications to mapping total biomass or changes in biomass, which could be used as validation datasets for remote-sensing products (Lu et al., 2012;Weiss et al., 2007) or to help understand carbon cycle dynamics. Furthermore, it may be possible to isolate the fraction of the total neutron signal that corresponds to an exchange of water from the surface to the atmosphere. It is the exchange of all water from the surface that will influence the transfer of mass, momentum, and energy, and this total exchangeable water is the critical climatic variable in interactions between land surface and the atmosphere and in ecosystem processes.