Published online February 1, 2019

Original ResearchCore Ideas• Cosmic-ray neutron data are used
to inversely estimate soil hydraulicproperties.
• The forward neutron operator COSMIC is coupled with the vadose zone
model HYDRUS-1D.
• Bayesian analyses confirm theinformation content of cosmic-ray
neutron data.

G. Brunetti, Dep. of Land, Air and Water
Resources, Univ. of California, Davis, CA
95616, and Institute of Hydraulics andRural Water Management, Univ. of Natural
Resources and Life Sciences (BOKU), Vienna,
1190, Austria; J. Šimůnek, Dep. of Environmental Sciences, Univ. of California, Riverside,
CA 92521; H. Bogena, R. Baatz, J.A. Huisman,and H. Vereecken, Agrosphere Institute,
Forschungszentrum Jülich GmbH, 52428
Jülich, Germany; H. Dahlke, Dep. of Land,
Air and Water Resources, Univ. of California,
Davis, CA 95616. *Corresponding author(h.vereecken@fz-juelich.de).
Received 26 June 2018.
Accepted 23 Sept. 2018.
Citation: Brunetti, G., J. Šimůnek, H. Bogena, R.
Baatz, J.A. Huisman, H. Dahlke, and H. Vereecken. 2019. On the information content ofcosmic-ray neutron data in the inverse estimation of soil hydraulic properties. Vadose Zone
J. 18:180123. doi:10.2136/vzj2018.06.0123© Soil Science Society of America.
This is an open access article distributed under
the CC BY-NC-ND license(http://creativecommons.org/licenses/by-ncnd/4.0/).

On the Information Content of CosmicRay Neutron Data in the Inverse
Estimation of Soil Hydraulic Properties
Giuseppe Brunetti, Jiří Šimůnek, Heye Bogena, Roland Baatz,
Johan Alexander Huisman, Helen Dahlke, and Harry Vereecken*
Observations of soil moisture content from remote sensing platforms can beused in conjunction with hydrological models to inversely estimate soil hydraulicproperties (SHPs). In recent years, cosmic-ray neutron sensing (CRNS) has provento be a reliable method for the estimation of area-average soil moisture at field
scales. However, its use in the inverse estimation of the effective SHPs is largely
unexplored. Thus, the main objective of this study was to assess the information content of aboveground fast-neutron counts to estimate SHPs using botha synthetic modeling study and actual experimental data from the Rollesbroichcatchment in Germany. For this, the forward neutron operator COSMIC was externally coupled with the hydrological model HYDRUS-1D. The coupled model was
combined with the Affine Invariant Ensemble Sampler to calculate the posterior
distributions of effective soil hydraulic parameters as well as the model-predictive
uncertainty for different synthetic and experimental scenarios. Measured watercontents at different depths were used to assess estimated SHPs. The analysis ofboth synthetic and actual CRNS data from homogenous and heterogeneous soil
profiles, respectively, led to confident estimations of the shape parameters a andn, while higher uncertainty was observed for the saturated hydraulic conductivity. Furthermore, results demonstrated that neutron data are less influenced by
local sources of uncertainty compared with near-surface point measurements.
The simultaneous use of CRNS and water content data further reduced the overall uncertainty, opening up new perspectives for the combination of CRNS withother remote sensing techniques for the inverse estimation of the effective SHPs.
Abbreviations: AIES, Affine Invariant Ensemble Sampler; CRNP, cosmic-ray neutron probe; CRNS, cosmicray neutron sensing; IAT, integrated autocorrelation time; LSM, land surface model; MCMC, Markov chain
Monte Carlo; MCNPX, Monte Carlo N-Particle eXtended; PTF, pedotransfer function; SHPs, soil hydraulicproperties; VGM, van Genuchten–Mualem.

Soil hydraulic properties (SHPs) define the relationship between volumetric watercontent (q), pressure head (h), and hydraulic conductivity (K), thus regulating the movement of water in the soil and influencing the water–energy cycle at the land surface. In
particular, the near-surface soil moisture has a profound influence on the partitioningof precipitation into surface runoff, evapotranspiration, and infiltration and thus on thesurface energy balance. In this view, there is strong evidence that land surface processes
play a key role in climate models (Pitman, 2003). An accurate description of hydrologicalprocesses at the land surface is thus of crucial importance.
Land surface models (LSMs) incorporating a mechanistic description of infiltration
require precise estimation of large-scale SHPs. Typically, the measurement scale for thecharacterization of SHPs is on the order of 10 cm, with a sample spacing of 100 m or larger(Hopmans et al., 2002). In this perspective, the upscaling of point measurements is time
consuming and introduces significant uncertainty into land surface hydrological modeling (Vereecken et al., 2007). To overcome this problem, pedotransfer functions (PTFs)
are often used (Van Looy et al., 2017) to translate available soil textural information into
SHPs. However, uncertainty in PTF inputs often propagates into the poor accuracy of
SHPs (Deng et al., 2009).
Vadose Zone Journal | Advancing Critical Zone ScienceA valid alternative is the inverse estimation of effective SHPs
from transient data (Dane and Hruska, 1983; Kool and Parker,
1988; Šimůnek and van Genuchten, 1996; Iden and Durner,
2007; Ines and Mohanty, 2008; Vrugt et al., 2008a). For example,
Gutmann and Small (2007) used the land surface model Noah (Eket al., 2003) to compare the results obtained using inversely estimated and PTF-based SHPs. They found a significant reduction
in errors between measured and simulated latent heat fluxes when
inversely estimated SHPs were used. Moreover, the use of emerging
global optimization algorithms and statistically rigorous methods
based on Bayesian inference now allow obtaining both a reliable
estimate of SHPs and a proper assessment of their uncertainty (e.g.,
Huisman et al., 2010; Scharnagl et al., 2011; Wöhling and Vrugt,2011), which is directly related to the information content of themeasured data used in the inverse estimation framework.
Traditionally, point measurements of pressure head andwater content have been used for the inverse estimation of SHPs.
However, during the last decade, soil moisture remote sensing
techniques have gained popularity (Bogena et al., 2015; Mohanty
et al., 2017). Despite their restriction to near-surface soil moisture measurements, several studies have demonstrated that remotesensing data can be used to successfully estimate effective SHPs
at the field scale (Mohanty, 2013). For example, Santanello et al.
(2007) inversely estimated soil hydraulic parameters of a semiarid
watershed using the Noah land surface model by combining soil
moisture evaluated using passive microwave remote sensing by
NASA’s push-broom microwave radiometer and active microwave
RADARSAT-1 imagery with a parameter estimation algorithm
and then compared the estimated and PTF-based SHPs. In
another study, Steenpass et al. (2010) assimilated infrared-measured soil surface temperatures and time domain reflectometry
based soil water contents to estimate SHPs for both synthetic andfield-scale conditions.
Among the ground-based remote sensing techniques, cosmicray neutron sensing (CRNS) has shown good promise to capturesoil moisture at relevant scales. Briefly, fast neutrons are producedby nuclear interactions between the incoming cosmic rays and elements of the Earth’s atmosphere (Zreda et al., 2008). When theyreach the soil surface, they penetrate to a certain depth and arescattered back into the atmosphere. Since fast neutrons are mainly
moderated by hydrogen, the fast neutron intensity at near groundlevel is negatively correlated with the near-surface soil moisture.
For this reason, soil moisture can be inferred from the fast neutron intensity measured by cosmic-ray neutron probes (CRNPs)(Zreda et al., 2008). The CRNP measures integral soil moisture
in a circular footprint centered on a detector. This radius varies
between 130 and 240 m depending on the site conditions, which
is significantly larger than the typical spatial correlation lengthof soil moisture patterns (Western et al., 2004). The penetration
depth of CRNP measurements varies between 15 cm for wet soils
to 55 cm for dry soils (Schrön et al., 2017). Due to this large sample
volume, the information content of CRNP data could potentially
improve the estimation of effective SHPs for LSMs because the

VZJ | Advancing Critical Zone Sciencescale mismatch between sample volume and model resolution ismuch smaller.
The applicability of the CRNP method has been demonstrated for various land cover types, soil types, and climates, and
CRNPs have been installed in several national monitoring networks (Bogena et al., 2015). The first network was established in
the United States by the University of Arizona and has alreadydeployed more than 60 CRNPs at various locations (Zreda etal., 2012). Similar networks have been established in Australia(Hawdon et al., 2014), the United Kingdom (Evans et al., 2016),
and Germany (Baatz et al., 2014). To date, almost 200 stationary
CRNPs have been installed worldwide (Andreasen et al., 2017).
The required conversion of neutron intensity into soil moistureand the correction for important influencing factors (e.g., biomass,snow, and litter layer water content) have been described in detailin the literature (e.g., Andreasen et al., 2017; Desilets et al., 2010).
Few recent studies have focused on the use of CRNP data for theestimation of SHPs with promising results. Both Finkenbiner et al.
(2018) and Gibson and Franz (2018) combined CRNP data with
empirical orthogonal functions to isolate spatial patterns of SHPs
in selected agricultural fields. A statistical analysis confirmed thatthe use of CRNP data can improve the estimate of SHPs at thefield scale.
Several types of neutron source models are available tostudy cosmic-ray neutron interactions near the soil surface. Oneof the most widely used models is the Monte Carlo N-Particle
eXtended (MCNPX) neutron transport code (Pelowitz, 2011).
Given the specified chemistry of the atmosphere and soil (e.g.,soil mineralogy, soil moisture), MCNPX simulates nuclear collisions of individual, randomly generated, incoming cosmic rays
and their interaction byproducts through the atmosphere and in
the soil using libraries of nuclear properties. In this way, MCNPX
allows accurate simulation of the number of fast neutrons thatenter a defined detector volume above the ground. In an attempt
to simultaneously describe the neutrons’ dynamics and waterflow in the unsaturated zone, Franz et al. (2012) successfully
coupled MCNPX with the hydrological model HYDRUS-1D(Šimůnek et al., 2016) for calculating the effective sensor depthof a CRNP while accounting for three sources of hydrogen andtheir vertical variability. However, the high computational costof this Monte Carlo based model makes it unsuited for Bayesiananalysis frameworks, which require thousands of model executions.
Recently, the Ultra Rapid Adaptable Neutron-Only Simulation(URANOS) model (Köhli et al., 2015) was developed to reducethe computational cost of CRNP simulations and to provide aneutron transport code tailored to environmental applications.
By neglecting less relevant physical processes, URANOS is morecomputationally efficient than other available Monte Carlo codes
while maintaining a similar accuracy in the simulation of neutron
fluxes (Köhli et al., 2015). However, the computational gain is still
not sufficient to justify its use in Bayesian analysis.
The forward neutron operator COSMIC has taken a major
step toward the efficient integration of CRNS data in models
p. 2 of 24

(Shuttleworth et al., 2013). It provides a simplified description ofneutron interactions, and the use of a physically based analytical
approach reduces the computational cost by four orders of magnitude compared with MCNPX. COSMIC calculates the fastneutron flux from the soil moisture profile and can be easily coupled with hydrological models. Recently, Baatz et al. (2017) wereable to assimilate data from a network of CRNPs installed in theRiver Rur catchment into the Common Land Model (CLM4.5)
to update soil moisture and SHPs of the catchment. In anotherstudy, the COSMIC and Noah models were coupled to assimilate
CRNP data from the Santa Rita Range field site (Shuttleworth etal., 2013). In particular, the Noah–COSMIC model was used inconjunction with the global optimization algorithm AMALGAM(Vrugt and Robinson, 2007) to recalibrate the Noah model byminimizing the errors between simulated and observed neutron
fluxes. Their results demonstrated that the assimilation of CRNP
data could effectively improve the accuracy of the model predictions. To date, the only available studies on the use of CRNP datain the inverse estimation of soil hydraulic properties are thoseof Rivera Villarreyes et al. (2014) and Baatz et al. (2017). In theformer, integral soil moisture estimates obtained for the CRNP
support volume were transformed into the soil water storage fordirect calibration of the effective soil hydraulic parameters using
the HYDRUS-1D model. Despite the overall promising resultsof their study, the numerical approach based on the use of soil
water storage data led to a poor approximation of the volumetricwater content at different depths. In Baatz et al. (2017), soil watercontents estimated from CRNS data were assimilated with thelocal ensemble transform Kalman filter in the CLM4.5 model to
update soil texture and organic matter used in the PTFs for theestimation of the SHPs. Results of that study demonstrated thepotential of CRNS networks to improve subsurface parameterization in the regional LSM. However, the use of PTFs and theindirect assimilation of CRNS data represent potential sources ofadditional uncertainty in the data assimilation framework.
Thus, the main aim of this study was to provide a comprehensive theoretical and experimental assessment of the information
content of CRNP data for the inverse estimation of effective
SHPs at the field scale. The problem was addressed in the following way. First, a combined mechanistic description of variably
saturated flow and soil-neutron interactions was provided bycoupling the hydrological model HYDRUS-1D and the neutrontransport code COSMIC. Before proceeding with an analysis ofactual CRNP data from an experimental facility in Germany, ananalysis of synthetic data was performed by combining the coupled
HYDRUS–COSMIC model with a Bayesian inference algorithm.
The synthetic example simplified the experimental case study bydescribing the soil-neutrons interactions in a homogeneous soil
profile and provided the first proof of concept for the estimationof SHPs from CRNP data. Finally, the actual data from a heterogeneous soil profile was used to evaluate the potential of CRNP
data for the inverse estimation of SHPs under real conditions. A
comparison between inverse estimation using neutron fluxes and

VZJ | Advancing Critical Zone Scienceclassical point measurements of water content completes and further enriches the analysis.

66Materials and Methods

TheorySubsurface Water Flow and Root Water Uptake
The finite-element hydrological model HYDRUS-1D(Šimůnek et al., 2016) numerically solves the one-dimensional
Richards equation describing variably saturated water flow in thesoil profile:æ ¶ h öù¶q ¶ é
= ê K ( h )ççç +1÷÷÷ú - Sè ¶zøúû¶ t ¶ z êë[1]

where q is the volumetric water content [L3 L−3], h is the soil waterpressure head [L], K(h) is the unsaturated hydraulic conductivity[L T−1], t is time [T], z is the soil depth [L], and S is a sink term
[L3 L−3 T−1], defined as a volume of water removed from a unit
volume of soil per unit of time due to plant water uptake. A solution of the Richards equation requires knowledge of the SHPs,
which should accurately describe the hydraulic characteristics ofthe porous medium and its heterogeneity. It is well known that
agricultural land management can significantly influence the soil
structure and corresponding SHPs in the near surface. Under suchcircumstances, the shape of the retention curve may be constantwith depth, but the soil hydraulic conductivity and porosity are
expected to change significantly. Thus, in the present study, thevertical variability of the soil is described by combining the unimodal van Genuchten–Mualem (VGM) (van Genuchten, 1980)
functions with the scaling factors proposed by Vogel et al. (1991).
In particular, three independent dimensionless scaling factors
are used: j q , j h, and j K for the water content, the pressure head,
and the saturated hydraulic conductivity, respectively. These threescaling parameters are used to define a linear model of the actual
spatial variability in the soil hydraulic properties as follows:ì q ( h ) = q r +j q é q * ( h * )-q r * ù
ï
ï
ë
û
ï
ïh = j h *
í
h
ï
ï
ï
K h = j K K *( h *)
ï
î ( )

[2]

where q*(h*) and K*(h*) are the reference soil hydraulic functions,which are ideally those unaffected by surface management activities. This theoretical approach guarantees high modeling flexibility
since it can describe both homogeneous and heterogeneous soil
profiles using a single set of reference soil hydraulic functions andonly a few more parameters (i.e., the scaling factors). If the soil ishomogeneous, the scaling factors are set to 1.
Feddes et al. (1978) defined S as
S ( h ) = a( h ) S p[3]

where Sp is the potential root water uptake rate and a(h) is a dimensionless water stress response function that depends on the soilpressure head h and can vary between 0 and 1. Parameters of thep. 3 of 24

stress response function for a majority of agricultural crops canbe found in various databases (e.g., Taylor and Ashcroft, 1972;
Wesseling et al., 1991). Potential root water uptake Sp is calculated
from the potential transpiration rate Tp and the root distributionfunction. Beer’s equation is used to partition the reference evapotranspiration calculated with the Penman–Monteith equation(Allen et al., 1998) into potential transpiration and soil evaporation fluxes (e.g., Ritchie, 1972) using leaf area index information.
For a detailed explanation of the partitioning of evapotranspiration, see Sutanto et al. (2012).

Cosmic-Ray Soil Moisture Interaction CodeThe physically based COsmic-ray Soil Moisture Interaction
Code (COSMIC) (Shuttleworth et al., 2013) is used to calculate
the aboveground neutron intensity for a given soil moisture profile. COSMIC includes simple descriptions of (i) degradation of
the incoming high-energy neutron flux with soil depth, (ii) creation of fast neutrons at each depth in the soil, and (iii) scatteringof the resulting fast neutrons before they reach the soil surface.
See Shuttleworth et al. (2013) for a detailed description of the
COSMIC model.
COSMIC requires several site-independent and site-specifictime-invariable parameters. In particular, the parameters L1 =162.0 g cm−2 , L2 = 129.1 g cm−2 , and L4 = 3.16 g cm−2 were foundto be site independent (Shuttleworth et al., 2013), while the parameters L3 (g cm−2) and a (cm3 g−1) were found to be dependent on
the soil bulk density r b according to

L3 =- 31.65 + 99.29r b[4]

a = 0.404 - 0.101r b[5]

The remaining two parameters, the lattice water content q l3[L L−3] and the number of high-energy neutrons at the soil surface NCOSMIC (dimensionless), are also site dependent and must
be measured and calibrated, respectively.
A FORTRAN version of the COSMIC code is freely available online (http://cosmos.hwr.arizona.edu). For this study, theFORTRAN code was rewritten in the Python language to enable
a more efficient model coupling with HYDRUS-1D. The opensource compiler Numba (Lam et al., 2015) was used to speed up
the Python code and to significantly reduce its execution time.
The Python version of the COSMIC code was verified against
the original FORTRAN code before further use.

Model Coupling Strategy
The COSMIC and HYDRUS models were coupled via a
Python interface to dynamically calculate subsurface water flow
by HYDRUS and neutron intensity by COSMIC. To this end,
several user-defined Python subroutines that interact with the
HYDRUS files were developed. The coupling of the two models
is accomplished in the following way:
1. Model initialization: HYDRUS is initialized by setting up
input data, the numerical domain, and boundary conditions.

VZJ | Advancing Critical Zone ScienceIn addition, the temporal resolution of the model output forthe soil moisture profile is defined, which determines the temporal resolution of the simulated neutron flux. Similarly, the
COSMIC parameters are defined, and the Python code is compiled using Numba.
2. HYDRUS execution and data processing: HYDRUS-1D isexternally executed directly from the Python code, and a checkis performed for numerical convergence and the final mass balance error. If the simulation is successful, the calculated soil
moisture profiles are read and stored in a matrix, in which rows
represent the vertical discretization of the model domain andcolumns represent the selected time steps for soil moisture profile information.
3. COSMIC execution: The matrix of soil moisture profiles is passed
to COSMIC. Two nested loops are used to iterate through depth
and time and to calculate a time series of simulated abovegroundneutron fluxes that can be compared with measured CRNP data.

Uncertainty Analysis
The main goal of the uncertainty analysis is to provide a statistical basis for the assessment of the information content of CRNP
data for the inverse estimation of soil hydraulic parameters. This isaccomplished by first using a global optimization algorithm for theidentification of the global optimum, which is subsequently used
as a starting point for the Bayesian uncertainty analysis based on a
Markov chain Monte Carlo (MCMC) algorithm to reduce its autocorrelation time, thus improving the computational efficiency of
the analysis. The global optimum is determined using the Particle
Swarm Optimization algorithm (Kennedy and Eberhart, 1995).

Bayesian Inference
Global optimization algorithms are generally used to find an
optimal set of parameters that maximize a likelihood function.
However, different sources of uncertainty can affect hydrological
data assimilation frameworks (i.e., model inadequacy, observations
and measurements errors, improper boundary conditions, etc.).
Therefore, there is no reason to assume that a single optimized
parameter set represents the only realistic model realization. It is
more reasonable to consider it as one of many acceptable solutionsof the problem. This idea underpins the equifinality concept first
introduced by Beven and Binley (1992) in the early 90s, whichemphasizes the fact that there are many acceptable (i.e., behavioral)
representations that cannot be excluded from the analysis and thatshould be considered when quantifying the uncertainty associatedwith model predictions.
Several methods have been proposed in the literature toassess the model predictive uncertainty (e.g., Beven and Binley,
1992; Kavetski et al., 2006; Vrugt et al., 2008a). However,during the last decade, the scientific community has focusedon the application of statistically rigorous Bayesian methods for the calibration of hydrological models. To explain theBayesian approach, let us consider a HYDRUS–COSMIC model
realization with a symbol Y, which is a function of q model
parameters, u = {u1, u2 ,…, uq}:p. 4 of 24

Y = f (u )[6]Bayes’ theorem combines the prior knowledge about theinvestigated process with the observed data to obtain the posterior distribution of the parameters. In a probabilistic context, thiscan be written as | u )´Pr ( u ) = Pr ( Y )´Pr ( u | Y)Pr ( Y
Likelihood ´ Prior = Evidence ´ Posterior[7]

L( u )´p( u )d u = Z ´ p ( u )d u
where Ỹ are the observed data, which modulate our prior beliefp(u)du into the posterior p(u)du. One of the main goals of Bayesiananalysis in hydrological model calibration is the estimation of theparameters’ posterior distribution. This is usually accomplished byapplying MCMC algorithms. One of the main problems faced in
the Bayesian analysis of hydrological models is their strong nonlinearity and their discontinuous nature, which makes it very difficult
to apply efficient MCMC algorithms such as Hamiltonian Monte
Carlo sampling (Betancourt, 2017). Although the classic Metropolis–
Hastings algorithm (Metropolis et al., 1953) is easy to implement, it
is unable to deal with highly correlated and multimodal posteriors
often encountered in the calibration of hydrological models. Thus, avariety of MCMC algorithms have been developed and successfully
applied in water-related problems in recent years.
In this study, we use the Affine Invariant Ensemble Sampler(AIES), which was first proposed by Goodman and Weare (2010).
The AIES algorithm is an ensemble MCMC algorithm that has
proven to be extremely effective when the dimensionality of the
problem is not very high (Huijser et al., 2015), as is the case inthis study. Each ensemble, W, consists of D walkers, which can beconsidered as a vector in the q-dimensional parameter space. The
core of AIES is its affine invariance property, which implies thatthe performance of the method is independent of the aspect ratio
in highly skewed distributions. This is a highly desirable property
when calibrating hydrological models because of the prevalence ofhighly correlated posterior distributions. In AIES, the ensemble
Markov chain is evolved using a stretch move step (Goodman andWeare, 2010), which significantly reduces the chain’s autocorrelation time. Briefly, the main idea is that the ensemble carries useful
information about the approximated distribution, and thus eachwalker is updated using a complementary ensemble. A parameter a,which is usually set to 2.0 (Foreman-Mackey et al., 2012), is used
to “stretch” the proposal. For a thorough description of the AIES,
see Goodman and Weare (2010), Foreman-Mackey et al. (2012),
and Huijser et al. (2015). It is worth noting that this study is thefirst one that uses the AIES algorithm in a hydrological context.
One of the main advantages of AIES is that it requires the
specification of only two parameters, i.e., the stretch factor, a, andthe number of walkers, D. There are no clear indications in the literature about the best values of D and a, and they typically dependon the investigated problem. Based on some preliminary tests, D
and a were set to 10q (q is the number of model parameters) and 2.0,

VZJ | Advancing Critical Zone Science

respectively. In this study, the Python package emcee (ForemanMackey et al., 2012) is coupled with the HYDRUS–COSMIC
model to carry out the Bayesian analysis. Emcee is a stable and welltested Python implementation of the AIES method that includes
several MCMC convergence diagnostic tools. The walkers are
initialized by random sampling from a multivariate normal distribution centered on the global optimum with a diagonal covariancematrix characterized by a small standard deviation of 0.01. Suchan initialization around the global optimum is expected to reduce
the autocorrelation time, thus increasing the overall computational
efficiency of the analysis.
The AIES algorithm calls the coupled HYDRUS–COSMIC
model at each iteration. However, the stochastic combination of various SHPs can lead to non-convergent model runs for HYDRUS-1D,
mainly due to the inability to appropriately adjust the HYDRUS-1D
settings (e.g., mesh size) during the Monte Carlo analysis. In thisstudy, non-convergent runs are identified and an unrealistically large
value is attributed to the sum of squared residuals. This procedure,
which is similar to what was reported by Brunetti et al. (2016) and
Wöhling and Vrugt (2011), can significantly slow down the convergence of the AIES because it introduces sudden and unrealistic
jumps in the response surface due to the artificially low value of thelikelihood. In future work, regularization of the response surface
could be considered through a local interpolation or a surrogate
model (e.g., Brunetti et al., 2017).
A non-informative uniform prior is used for all parametersin the uncertainty analysis. The only information provided bythis prior are the bounds of the feasible parameter space. Withinthese bounds, all parameter values have the same probability. The
bounds were set based on field measurements (e.g., particle size
distribution) and literature values.

Markov Chain Monte Carlo Diagnostic
Goodman and Weare (2010) suggested using the integrated
autocorrelation time (IAT), tint, to assess the accuracy of the AIES.
The IAT is the number of steps (i.e., ensemble moves) required
for the walker to produce an independent sample. This number
is expected to be relatively large if the walkers are initialized in a
distribution that is not equal to the stationary distribution. Sucha transient period is usually defined as burn-in and should be discarded from the analysis. In the present study, the burn-in period isassumed to be 2tint. It is worth noting that we avoid the use of thewidely applied Gelman–Rubin test because walkers in the AIES
are not formally independent.
The IAT determines statistical errors in the Monte Carloanalysis once a stationary distribution has been attained (Sokal,
1996). In particular, the statistical error e of the mean of the posterior depends on tint and is calculated ase=

2t intDm[8]where m is the number of steps. Obviously, a larger sample reduces
the errors but significantly increases the computational time. Thus,
p. 5 of 24

in this study, the admissible maximum statistical error was set to5%, which should guarantee a good trade-off between numericalaccuracy and computational effort.
The estimation of the IAT is notoriously difficult because it
is a two-point statistic. More specifically, the IAT should asymptotically reach its true value for a large number of samples. Thus,a preliminary IAT convergence analysis is performed to assess its
calculated value. This is accomplished in the following way:ʶʶ The AIES algorithm is preliminary run for 3000 steps to
obtain a significant statistical sample. This requires a few daysof computation time considering that the execution time of the
HYDRUS–COSMIC model is <1 s.
ʶʶ Each walker is divided into 20 parts, and the IATs are calculated
for increasing chain lengths. We adopt the numerical proceduredescribed by Sokal (1996). In particular, a mean normalized
autocorrelation function is computed by averaging the normalized autocorrelation functions calculated for each walker in theensemble. Then, the procedure described by Sokal (1996) isused to estimate the IAT for each soil hydraulic parameter. The
IAT is expected to exhibit a logarithmic growth.
ʶʶ The analysis is assumed to converge if
ìe < 5%
ï
ï max
í
max
ï
ï
î m > 10t int[9]

Otherwise more samples are generated by evolving the AIES. The
second condition in Eq. [9] is meant to reduce the MCMC error
and mitigate the effect of an imperfect estimation of the IAT.

Likelihood Function
If we assume that the error residuals are uncorrelated andnormally distributed with constant variance, s2 , the likelihood
function L(u) can be written askL( u ) = Õ
i =1

1

{

1
2
exp - s2 éë yi ( u )- yi ùû
2
2ps
2

}

[10]

where yi(u) and ỹi(u) are the ith model realization and its corresponding measured value, respectively. For algebraic simplicity, thelog-likelihood l(u) is used:

k
k
1 k
2
 ( u ) =- ln ( 2p )- ln ( s2 )- 2 å éë yi ( u )- yi ùû22
2s i=1[11]

Equation [11] is used in the inverse modeling with synthetic
data because the error variance is known. In the inverse modeling with measured data, s2 is treated as a nuisance parameter andintegrated out from the inference equation (e.g., Kavetski et al.,2006). The new likelihood is then written asük ïì k2ï ( u ) =- ln ïí å éë yi ( u )- yi ùû ïýïþï
2 ïîï i=1

[12]Rollesbroich Grassland Experimental Catchment
The Rollesbroich grassland experimental catchment is locatedin western Germany in the Eifel mountain range and is part of

VZJ | Advancing Critical Zone Science

the TERENO observatory Eifel/Lower Rhine Valley (Bogena etal., 2012, 2018). It covers an area of about 40 ha with altitudes
ranging from 474 to 518 m asl (the average slope is 1.63°). The
mean annual air temperature and precipitation are 7.7°C and1033 mm, respectively. The vegetation is dominated by perennial
ryegrass (Lolium perenne L.) and smooth meadow grass (Poa pratensis L.), and the predominant soils are Cambisols in the southernpart and Stagnosols in the northern part of the catchment. Allcomponents of the water balance (e.g., precipitation, evapotranspiration, runoff) are continuously monitored using state-of-the-art
instrumentation. A detailed description of the instrumentation
at the Rollesbroich experimental catchment and available hydrological datasets was given by Qu et al. (2016). A SoilNet wirelessnetwork (Bogena et al., 2010) with more than 500 soil moisturesensors (ring oscillator probes, Bogena et al., 2017) was installedin May 2011 to measure the soil water content at depths of 5, 20,
and 50 cm (Qu et al., 2014) at 87 locations. In particular, two
sensors were installed in parallel at each depth with a distance of~10 cm to increase the sensing volume and to allow examinationof inconsistencies in sensor reading. In May 2011 and May 2012,
two CRNPs (Type CRS1000, HydroInnova LLC) were installedin the southern and northern parts of the catchment. The CRNP
calibration and a comparison with in situ soil water content dataobtained with SoilNet were presented in Baatz et al. (2014). This
study uses data only from the CRNP located in the southern partof the catchment.
A 6-mo-long dataset (1 Apr. 2012–30 Sept. 2012) was
selected for further analysis (Fig. 1). In particular, daily precipitation and reference evapotranspiration were used as input data in
HYDRUS for both synthetic and experimental scenarios, while
daily measured water contents at three depths (i.e., z = −5, −20,
and −50 cm) and the aboveground neutron fluxes were used in
the inverse estimation of the soil hydraulic parameters in theexperimental scenarios.

Numerical Domain and Boundary Conditions
The 200-cm-deep soil profile was discretized in HYDRUS-1D
into 100 finite elements, refined at the top to accommodate
pressure head gradients induced by the atmospheric conditions. An
atmospheric boundary condition is applied at the soil surface using(i) precipitation and potential evaporation fluxes, (ii) a prescribedzero pressure head (i.e., full saturation) during ponding, and(iii) equilibrium between the soil surface pressure head and theatmospheric water vapor pressure when the atmospheric evaporative
demand cannot be met. The effect of groundwater was not simulated
in this study and thus a free-drainage boundary condition was usedat the bottom of the model domain (z = −200 cm).
The initial pressure head was assumed to be constant andequal to −100 cm in the entire model domain. A 3-mo-long spinup period was used to limit the influence of the initial conditionon the simulation results. Three observation points were set at
depths of −5, −20, and −50 cm to retrieve simulated volumetricwater contents.
p. 6 of 24

Fig. 1. (A) Precipitation (bars) and reference evapotranspiration ET0 (dashed line) used as input data in HYDRUS, and (B) measured neutron fluxesand (C,D,E) water contents at three depths in the Rollesbroich catchment between April and September 2012.

Inverse Estimation with Synthetic Data
Several studies have demonstrated how data assimilation
frameworks can be significantly affected by various sources of
uncertainty at the field scale that can lead to biased conclusions if
not properly accounted for in the analysis (Mantovan and Todini,
2006; Mockler et al., 2016; Montanari, 2007). Because the main
aim of this study was to provide an objective evaluation of theinformation content of the CNRP data for the inverse estimation of SHPs at the field scale, it was important to first analyze
synthetic modeling scenarios that are not affected by uncertainty
in the model structure and measurement errors before analyzing
actual experimental data. Here, the synthetic modeling scenario
is a simplification of the experimental case study and is focused

VZJ | Advancing Critical Zone Science

on soil-neutron interactions in a homogeneous sandy soil. The
strongly unsaturated conditions induced by the use of a sandy soil
will increase the aboveground neutron intensity and thus provide
different information than the experimental case study characterized by a silty soil.

Model Setup
Variably saturated water flow is simulated in a 200-cm-deep
homogeneous loamy sand soil profile covered by grass. The VGM
parameters reported by Carsel and Parrish (1988) for loamy sandwere used to generate the synthetic dataset. The scaling factors
were fixed to 1.0 to obtain a homogeneous soil profile. The rooting
depth and density were assumed to be 15 cm and 1.0, respectively,
p. 7 of 24

while a constant leaf area index of 2.0 (Qu et al., 2016) was used
to partition reference evapotranspiration. The Feddes parameters
for grass were chosen according to Taylor and Ashcroft (1972) andWesseling et al. (1991). The COSMIC parameters were set according to Baatz et al. (2014) (Table 1).
The synthetic dataset was generated in the following way.
First, a forward numerical simulation by the HYDRUS-1Dmodel was used to generate a time series of the volumetric watercontent at z = −5 cm and daily soil moisture profiles. The latter
were passed to the COSMIC model for the simulation of theaboveground neutron fluxes. Next, synthetic noise was added tothe original time series. The selection of s values has a significant influence on parameter uncertainty and thus they should beproperly set to reflect the accuracy of real measured data. Indeed,
an unbalanced choice of such values can lead to biased conclusions on the information content of the investigated data. These
values were thus set in this study based on a preliminary analysis of the noise in selected measured datasets. Briefly, a spectral
analysis based on the Fast Fourier Transform was used to separate the deterministic (signal) and stochastic (noise) componentsof the measured neutron fluxes and volumetric water contents.
The Shapiro–Wilk test was then used to assess the normality of
the noise, and the standard deviations were calculated. Based on
this analysis, random noise with zero mean and standard deviations equal to 30.0 N h−1 and 0.01 cm 3 cm−3 for the neutron
flux and water content time series, respectively, was used. The

Table 1. The COSMIC parameters used in the synthetic and experimental scenarios (Baatz et al., 2014).
ql

NCOSMIC L1

L2

L4

1.42

0.037

189

161.98129.14

3.16

Experimental 1.09

0.068

213

161.98129.14

3.16

Scenarior bd†
g cm−3Synthetic† r bd, soil bulk density; q l, lattice water content; NCOSMIC , number of highenergy neutrons at the soil surface; L1, L2 , and L3, empirical parameters.

input, original, and synthetic data for the selected time period
are shown in Fig. 2.

Synthetic Modeling Scenarios
Three synthetic modeling scenarios are considered for the
inverse estimation of soil hydraulic parameters (Table 2). Scenarios
S1 and S2 examine the information content of the temporal evolution of aboveground neutron flux and water content, respectively,
at one observation point. The comparison of the parameter uncertainty obtained by these two scenarios will clarify the potential of
CRNP data to estimate SHPs compared with traditional point
measurements. Scenario S3, which combines the two synthetic
datasets, is used to investigate the effect of the joint inversion ofwater content and neutron flux data. The parameter bounds usedin this analysis are reported in Table 3.

Fig. 2. (A) Precipitation (bars) and reference evapotranspiration ET0 (dashed line) from the Rollesbroich catchment used as input data in HYDRUS,
and simulated (solid lines) and synthetic (with added noise, diamonds) (B) neutron fluxes and (C) water contents at depth z = −5 cm between April
and September 2012.

VZJ | Advancing Critical Zone Sciencep. 8 of 24

Inverse Estimation with Actual Data
Model Setup
Measured particle size distributions reported by Qu et al.
(2016) for the Rollesbroich catchment revealed negligible differences in the textural composition of the soil at different depths.
On the other hand, bulk densities increased with depth from 0.94
to 1.52 g cm−3, and porosities decreased with depth from 0.65 and0.43 cm3 cm−3. This vertical soil variability can be attributed toprior agricultural land management (i.e., the plow layer is 20–30
cm thick), which reduced the compaction of the surface layer. The
measured soil properties highlight the existence of two different
soil horizons in the top layer (i.e., z > −30 cm). Indeed, the porosity and the bulk density for z > −10 cm are 20% higher and 25%
lower, respectively, than for the underlying topsoil. Thus, threesoil materials are assumed in the model: Soil Material 1 extends
to z = −10 cm; Soil Material 2 ranges from −10 to −30 cm; andSoil Material 3 extends below −30 cm. Since Soil Material 3 is notaffected by agricultural land management, the scaling factors were
set to 1.0 for this material.
The pedotransfer model ROSETTA (Schaap et al., 2001)
was used to identify the plausible ranges of the soil hydraulicparameters and the scaling factors. On the basis of this analysis, thewater content scaling factors j q were fixed to 1.4 and 1.1 for Soil
Materials 1 and 2, respectively, thus reducing the dimensionality of
the inverse problem. Because the residual water content exhibited
low variability, it was fixed to 0.07, which is a common value forsilty soils. The tortuosity and pore-connectivity parameter L(dimensionless) in the VGM equation was set to 0.5 in all numericalsimulations. The aforementioned model simplifications reduced
the number of unknown parameters to eight. The parameterbounds used in the analysis are reported in Table 3.

Experimental Modeling Scenarios
Three experimental data analysis scenarios are considered(Table 2). Scenarios E1 and E2, in analogy to Scenarios S1 andS2, consider only measured neutron flux and near-surface watercontent (i.e., z = −5 cm), respectively. The comparison of these
two scenarios will clarify the information content of the CRNP
data to estimate SHPs in layered soils relative to near-surface point
measurements. In addition, Scenario E3 uses measured watercontents at three different depths (i.e., z = −5, −20, and −50 cm).
This scenario should provide more accurate estimates of SHPs
because of the higher amount of measurements in multiple soil
layers and thus serves as a benchmark for the other two scenarios.

Table 2. Synthetic modeling and experimental data analysis scenarios.
Scenario
Data

Synthetic
Experimental(homogeneous soil) (heterogeneous soil)

Aboveground neutron flux time series

S1

E1

Near-surface water content time series(z = −5 cm)

S2

E2

Three water content time series (z = −5, –−20, and −50 cm)

E3

Near-surface water content and neutron S3
flux time series–

that a reliable estimation has been obtained for all three scenarios.
It can be seen that all three scenarios exhibit similar autocorrelation times. In particular, the dynamics of the IAT indicate thatapproximately 100 ensemble moves (i.e., steps) are sufficient for
the AIES to forget about its initial position and to start samplingfrom a stable posterior distribution. The saturated water content
exhibits the largest IAT in Scenarios S1 and S2 (Fig. 3A and 3B),
while in Scenario S3 (Fig. 3C) it is slightly larger for the shapeparameter n. The short autocorrelation times lead to a maximum
statistical error of the Monte Carlo analyses (Eq. [8]) of approximately 3.5% after 3000 steps. This is well below the prescribed
tolerance of 5%, thus indicating a good accuracy in the estimationof the posterior distributions.

Uncertainty Analysis
The joint and marginal posterior distributions for each soil
hydraulic parameter for Scenario S1 (aboveground neutron fluxes)
are shown in Fig. 4. The inverse problem is markedly unimodal,
with strong parameter interactions. More specifically, the joint
Table 3. The bounds of the van Genuchten–Mualem parameters andtheir true values used in the synthetic and experimental scenarios.
Synthetic

Experimental

Parameter†

Lower
bound

Upper
bound

True value(Carsel andParrish, 1988)

Lower
bound

Upper
bound

q r, cm3 cm−3

0.03

0.08

0.06

0.07

0.07

q s, cm3 cm−3

0.35

0.45

0.41

0.30

0.45

log10 a, cm−1

−1.50

−0.50

−0.91

−2.50

−0.85

n

1.50

3.00

2.28

1.10

2.00

log10Ks , cm d−1

1.00

3.00

2.54

0.50

3.00

j h1

1.00

1.00

1.00

1.00

5.00

j h2

1.00

1.00

1.00

1.00

5.00

Inverse Estimation with Synthetic Data

jK1

1.00

1.00

1.00

1.00

20.00

Markov Chain Monte Carlo Diagnostic

jK2

1.00

1.00

1.00

1.00

5.00

The development of the IAT in the MCMC analyses for thethree synthetic scenarios is reported in Fig. 3. The asymptotic
behavior of the IAT with increasing algorithmic steps suggests† q r, residual water content; q s, saturated water content; a and n, shape parameters; Ks, saturated hydraulic conductivity, j h1 and j K1, and j h2 and j K2 ,pressure head and saturated hydraulic conductivity scaling factors for Soil
Materials 1 and 2, respectively.

66Results and Discussion

VZJ | Advancing Critical Zone Sciencep. 9 of 24

Fig. 3. The diagnostic plot of the Affine Invariant Ensemble Sampler (AIES) algorithm for synthetic modeling scenarios showing the integrated autocorrelation time tint against the algorithm step for the soil hydraulic parameters of residual soil water content (qr), saturated soil water content (qs), thevan Genuchten–Mualem parameters n and a, and saturated hydraulic conductivity (Ks), along with the maximum statistical error (emax) for Scenarios(A) S1, (B), S2, and (C) S3.

posterior distribution of q s and Ks reveals an expected positive
correlation, while the VGM shape parameters are negatively
correlated. The interdependence of a and n is exacerbated, as
indicated by their slight banana-shaped joint posterior distribution. Interestingly, the posterior for Ks is unimodal, even though
it spans a relatively broad range of values from 83 and 602 cm
d−1. On the other hand, the marginal posterior distribution forq s is flat, indicating that the neutron data do not contain enoughinformation content to update the uniform prior distribution.
This is caused mainly by the fact that the soil water content isrelatively low (0.05–0.20 m3 m−3) in these synthetic data due tothe use of a highly permeable loamy sand. Therefore, water content
never approaches the saturated water content of 0.41 and little
information is available to constrain q s. Different results would
probably have been obtained for other (finer) soil textures. In thiscase, independent measurements of soil porosity would be needed
to better constrain q s. The strong interaction effect shown in the
joint posterior q s–Ks suggests the impossibility of fixing the saturated water content independently from the saturated hydraulicconductivity. As expected in strongly unsaturated conditions, q r is
precisely estimated. Similarly to q s, different types of Ks measurements are needed to reduce the uncertainty. Similar findings werereported by Gibson and Franz (2018), who used a mobile CRNP
to quantify spatial patterns of SHPs across a range of spatial scales.
In that study, the researchers underscored how an inability to map
near-saturated soil water contents increased the uncertainty in theestimation of Ks.
Overall, the marginal posterior distributions indicate a wellposed inverse problem characterized by acceptable parameter
identifiability and moderate uncertainty. This is confirmed by
the Bayesian credible intervals reported in Table 4. In particular,
the shape parameters a and n are generally well predicted, with thelatter exhibiting a slightly narrower posterior distribution.

VZJ | Advancing Critical Zone ScienceAs shown in Fig. 5, the use of near-surface water content
instead of neutron flux (Scenario S2) leads to similar posterior
distributions. Again, the inverse problem is mostly unimodal andaffected by significant parameter correlations that are similar to
those in S1 (i.e., a positive correlation between q s and Ks, and a negative correlation between a and n). The residual water content q r
is generally uncorrelated, except for a slight interdependence with
n. As in Scenario S1, the near-surface water content measurements
do not provide enough information to estimate q s. The marginal
posterior distribution for Ks again is unimodal and spans values
between 100 and 1000 cm d−1, thus indicating a somewhat higher
uncertainty than in Scenario S1. This effect could be related tothe increased information content of the neutron flux, which isrepresentative of the soil moisture dynamics in the entire topsoil
and not just for a specific depth (i.e., 5 cm in this study).
Overall, a comparison of Fig. 4 and 5 suggests a minor decreasein parameter uncertainty for Scenario S2 compared with Scenario S1(Table 4). This is particularly evident for the VGM shape parameters,
Table 4. Original values and 95% Bayesian credible intervals of the soil
hydraulic parameters for the three synthetic scenarios. The 2.5 and 97.5%
quantiles are the bounds of the credible intervals.
ScenarioParameter†

2.5%

97.5% 2.5%

97.5% 2.5%

True value(Carsel and
97.5% Parrish, 1988)

q r, cm3 cm−3

0.05

0.06

0.05

0.06

0.05

0.06

0.057

q s, cm3 cm−3

0.35

0.44

0.35

0.44

0.35

0.44

0.41

log10 a, cm−1

−1.36 −0.66 −1.19 −0.7

−1.19 −0.76 −0.91

n

2.08

2.88

1.97

2.63

2.06

2.66

2.28

log10Ks, cm d−1

1.92

2.78

2.06

2.93

2.00

2.86

2.54

S1

S2

S3

† q r, residual water content; q s, saturated water content; a and n, shape parameters;
Ks, saturated hydraulic conductivity.

p. 10 of 24

Fig. 4. Joint (below diagonal) and marginal (diagonal) posterior distributions of the soil hydraulic parameters of residual soil water content (qr), saturated soil water content (qs), the van Genuchten–Mualem parameters n and a, and saturated hydraulic conductivity (Ks) for Scenario S1. The red lines
indicate the original values of the parameters for the loamy sand, while the isolines indicate 12, 39, 68, and 84% probability levels. The dashed lines
indicate the Bayesian credible intervals. The axis ranges correspond to the parameter bounds reported in Table 3.

which have narrower Bayesian credible intervals and a good agreement between the most probable and true values. However, it must
be noted that these differences depend on the selected level of synthetic noise used to generate the synthetic datasets. In this view, thesensor accuracy, as well as different sources of noise, plays a fundamental role in the analysis of experimental data. Similar to Scenario
S1, the results suggest the need to use more informative prior distributions and to combine different data types to reduce the overall
parameter uncertainty. This is investigated in Scenario S3, where
both aboveground neutron fluxes and water contents are used.
The MCMC analysis for Scenario S3 (with both near-surface water contents and aboveground neutron fluxes) resulted instrongly unimodal and narrow posterior distributions (Fig. 6). The
correlations between parameters are similar to those for Scenarios
S1 and S2, and the relatively high uncertainty in the estimation

VZJ | Advancing Critical Zone Science

of q s remains. The q r is precisely estimated, while the uncertainty
for Ks persists. The overall reduction in the parameter uncertainty
for Scenario S3 is visually evident by comparing the joint posterior
distributions shown in Fig. 4, 5, and 6, and is further confirmed by
the comparison of the Bayesian credible intervals (Table 4). Indeed,
the VGM shape parameters have narrower credible intervals, suggesting that the simultaneous use of more and different data types
reduces the uncertainty in the estimated parameters. The resultsof Scenario S3 further confirm that the saturated water content
and hydraulic conductivity cannot be precisely inferred from data
representative of strongly unsaturated conditions and that independent measurements of soil porosity and hydraulic conductivity
may be needed in such a case.
As a result of this analysis, a number of conclusions can be
drawn:
p. 11 of 24

Fig. 5. Joint (below diagonal) and marginal (diagonal) posterior distributions of the soil hydraulic parameters of residual soil water content (qr), saturated soil water content (qs), the van Genuchten–Mualem parameters n and a, and saturated hydraulic conductivity (Ks) for Scenario S2. The red lines
indicate the original values of the parameters for the loamy sand, while the isolines indicate 12, 39, 68, and 84% probability levels. The dashed lines
indicate the Bayesian credible intervals. The axis ranges correspond to the parameter bounds reported in Table 3.

ʶʶ The three synthetic datasets provided sufficient information
content to obtain a satisfactory estimate of the SHPs. More specifically, the residual water content, the VGM shape parameters,
and the saturated hydraulic conductivity exhibited unimodal
posterior distributions, with the latter characterized by a higher
uncertainty. The saturated water content could not be inferredfrom the dataset (for loamy sand), suggesting that independentmeasurements of the soil porosity are needed to better constrain
the inverse problem for this soil texture.
ʶʶ The use of the aboveground neutron flux data led to a satisfactory estimate of SHPs. In particular, the use of neutron flux dataslightly reduced the uncertainty in the estimation of Ks. This
provides the first proof of concept for the use of CRNP data for
the inverse estimation of effective SHPs at the field scale.
ʶʶ The simultaneous use of water content and neutron flux data
reduced the parameter uncertainty. This suggests that the

VZJ | Advancing Critical Zone Science

combination of CRNP data with sparse point measurementsin the near-surface environment could reduce parameter
uncertainty, in particular for the VGM shape parameters.
However, the number of point measurements should be
sufficient to characterize the soil heterogeneity within the
CRNP footprint.

Inverse Estimation with Actual Data
Markov Chain Monte Carlo Diagnostic
The IAT values for the MCMC analyses for the three experimental scenarios are shown in Fig. 7. The maximum statistical error
ranges between 3 and 5%, which is under the prescribed tolerance
and indicates a reliable estimation of the posterior distributions.
As expected, the increased dimensionality of the inverse problem
resulted in longer autocorrelation times and more steps before convergence to the true posterior distribution. This behavior agrees
p. 12 of 24

Fig. 6. Joint (below diagonal) and marginal (diagonal) posterior distributions of the soil hydraulic parameters of saturated soil water content (qs), thevan Genuchten–Mualem parameters n and a, saturated hydraulic conductivity (Ks) for Scenario S3. The red lines indicate the original values of theparameters for the loamy sand, while the isolines indicate 12, 39, 68, and 84% probability levels. The dashed lines indicate the Bayesian credible intervals. The axis ranges correspond to the parameter bounds reported in Table 3.

with the properties of the AIES algorithm in high dimensions
reported by Huijser et al. (2015).
The maximum IAT for Scenario E2 (Fig. 7B) is approximately
three times larger than the IAT obtained for Scenarios E1 andE3, thus indicating convergence issues. This is related to multimodal posterior distributions, discussed below. Poor convergenceis particularly evident for the scaling factor j K1, which is characterized by an IAT of 1600 steps. Under such circumstances, theensemble samplers are known to perform poorly, mainly due totheir inability to traverse low-probability valleys between peaks(Matthews et al., 2017). It is worth noting that the performanceof the AIES for Scenario E2 was insensitive to the value of thestretch factor a, which was then set to the suggested value of 2.0.
Future work should consider using multimodal nested samplingtechniques (Feroz et al., 2009), which were found to be effective for

VZJ | Advancing Critical Zone Sciencesuch types of inverse problems. On the other hand, the relatively
large autocorrelation time for the shape parameter n in Scenario E1(Fig. 7A) is partially related to the position of the high-likelihood
region, which is close to the boundary of the numerical domain.
Indeed, bounded prior distributions can lead to convergence issues
when the posterior develops near the boundary.

Uncertainty Analysis
Because Scenario E3 serves as a benchmark for the other two
scenarios, this will be discussed first. The joint and marginal posterior distributions of the SHPs for this scenario are shown in Fig.
8 and have a markedly unimodal pattern characterized by a leptokurtic behavior. Consequently, the Bayesian credible intervals(Table 5) for the SHPs are extremely narrow, thus indicating their
good identifiability. The Bayesian analysis reveals, as expected, ap. 13 of 24

Fig. 7. The diagnostic plot of the Affine Invariant Ensemble Sampler (AIES) algorithm for the experimental data analysis scenario showing the integrated autocorrelation time tint against the algorithm step for the soil hydraulic parameters of residual soil water content (qr), saturated soil watercontent (qs), the van Genuchten–Mualem parameters n and a, saturated hydraulic conductivity (Ks), and scaling factors for hydraulic head and hydraulic conductivity (jh1 and jh2, and jK1 and jK2 for Soil Materials 1 and 2, respectively), along with the maximum statistical error (emax) for Scenarios(A) E1, (B) E2, and (C) E3.

well-posed inverse problem characterized by low uncertainty
mainly stemming from the simultaneous use of measured watercontents at different depths, which better constrains the inverse
estimation problem. Overall, the results agree well with thosereported by Qu et al. (2016), where the ROSETTA pedotransfer function was used to predict the spatial distribution of SHPs
using the soil textural information for the top 5 cm of soil. Indeed,
the saturated water content for Soil Material 1 ranges between
0.50 and 0.52 cm3 cm−3, which is in line with the values of theporosity reported by Qu et al. (2016), ranging between 0.46 and0.56 cm3 cm−3 for the southern part of the Rollesbroich catchment.
Results slightly deviate for the VGM shape parameters a and n,which are estimated to be higher and lower, respectively, than thoseof Qu et al. (2016), who reported a and n values ranging between−2.40 and −2.30 cm−1 and between 1.63 and 2.23, respectively.
Higher uncertainty is found for the saturated hydraulic conductivity of Soil Material 1, which ranges between 10 and 40 cm d−1.
It is worth noting that estimates of SHPs using PTFs are usuallyassociated with significant uncertainty. Thus, it is preferable to use
inversely estimated SHPs when the model inversion is well posed
and the associated uncertainty is properly assessed.
The effect of the parameter uncertainty on the simulated
water contents and neutron fluxes is shown in Fig. 9, which reports
the model predictive uncertainty (gray lines) obtained by random
sampling of 500 solutions from the posterior distribution. The
model can accurately describe the soil moisture dynamics in thethree observation depths with limited uncertainty. The qualityof the fitting is significantly better than that reported by Qu et al.
(2014), who inversely estimated the SHPs from soil water contents
measured with a sensor network at 41 locations in the Rollesbroichcatchment. This error reduction confirms that the use of the scaling factors is a valid approach for this dataset because it results
in a better description of water flow than the assumption of

VZJ | Advancing Critical Zone Science

homogeneous soil used by Qu et al. (2014). On the other hand,the comparison between simulated and measured abovegroundneutron fluxes (Fig. 9A) reveals an appreciable deviation. Morespecifically, the coupled model tends to overestimate the variabilityof neutron fluxes. This is particularly evident during August and
September 2012. The Rollesbroich site is managed grassland with arelatively low biomass. Except for an occasional snow cover during
the winter period, there are no additional transient hydrogen poolsthat have a strong influence on the measured neutron signal. As
a consequence, such behavior can be partially explained by thesimplifying assumptions of the COSMIC model, by an imperfectneutron flux correction (e.g., for variations due to fluctuations in
the atmospheric density and humidity or incoming neutrons), andby an inadequate calibration of the COSMIC model parameters.
In such circumstances, two approaches can be recommended:
1. Two-Step Sequential Calibration: First, the HYDRUS modelis preliminarily calibrated against measured data from different locations and depths in the measurement footprint of
the cosmic-ray probe. Next, the calibrated HYDRUS model
is coupled with the COSMIC model for the inverse estimation of the COSMIC parameters from the measured neutron
fluxes. It must be emphasized that the possibility of using
simulated high-resolution soil moisture profiles instead of afew point measurements at different soil depths can increase
the accuracy of the COSMIC model predictions. In thisperspective, it may be worth coupling the three-dimensionalmodel (e.g., HYDRUS 2D/3D), with the COSMIC modelto account for the horizontal and vertical soil heterogeneity
at the field scale (»300 m). The recent work of Gibson andFranz (2018) suggested that sparse point measurements of the
SHPs can significantly improve the description of spatial patterns of SHPs.
2. Joint Calibration: The soil hydraulic and COSMIC parameters
are jointly estimated using the coupled HYDRUS–COSMIC
p. 14 of 24

Fig. 8. Joint (below diagonal) and marginal (diagonal) posterior distributions of the soil hydraulic parameters of saturated soil water content (qs), thevan Genuchten–Mualem parameters n and a, saturated hydraulic conductivity (Ks), and scaling factors for hydraulic head and hydraulic conductivity(jh1 and jh2, and jK1 and jK2 for Soil Materials 1 and 2, respectively) for Scenario E3. The isolines and the dashed lines indicate the 12, 39, 68, and
84% probability levels and the Bayesian credible intervals, respectively. The axis ranges correspond to the parameter bounds reported in Table 3.

model using both measured neutron fluxes and volumetricwater contents. However, this will increase the dimensionalityof the inverse problem, which can lead to convergence issues in
the MCMC analysis.
The joint and marginal posterior distributions of the soil
hydraulic parameters for Scenario E1 that relies on measured
aboveground neutron flux are shown in Fig. 10. It is evident thatthe measured neutron fluxes do not provide sufficient information to update the prior distributions of j h2 , j K1, and j K2 ,which remain almost flat after the MCMC analysis. Indeed, their
Bayesian credible intervals (Table 5) are only marginally narrower
than the prior bounds listed in Table 3. A similar but less pronounced behavior is visible for Ks and j h1, although the marginal
posterior distributions are right skewed, thus indicating a slightinfluence of the measured data on these parameters. It is worth
noting that the MCMC analysis reveals a highly skewed bivariate
posterior distribution for the pair of parameters a–Ks. In such circumstances, the affine invariant properties of the AIES algorithm
significantly mitigate convergence issues usually encountered in
traditional MCMC algorithms.

VZJ | Advancing Critical Zone ScienceInterestingly, the use of measured neutron fluxes leads to wellconstrained estimates of VGM shape parameters, thus confirming
the findings of the synthetic modeling scenarios. The confident
estimation of q s stems from the use of near-saturated water contents, which better constrain the inverse problem. This effect ismore pronounced for n and q s, while a presents a more platykurtic
distribution. In particular, the Bayesian credible intervals for q s
and n range between 0.33 and 0.38 cm3 cm−3 and between 1.10
and 1.47, respectively, which are slightly larger than those obtainedin Scenario E3. On the other hand, a ranges between 0.004 and0.04 cm−1, thus indicating a tendency to overestimate the soil
water retention capacity. Indeed, very low values of a are usuallyassociated with a high retention capacity and delayed desaturationof the porous medium. In such circumstances, simulated watercontent fluctuations are limited and propagated in the calculation of the fast neutron fluxes. Thus, from a numerical point ofview, the Bayesian inference process adjusted the SHPs to compensate for other plausible sources of uncertainty (e.g., COSMIC
model assumptions, COSMIC parameters). Therefore, the resultsfor Scenario E1 further confirm the findings for Scenario E3 andp. 15 of 24

suggest the use of a joint Bayesian calibration of the HYDRUS–COSMIC model.
This aspect is further investigated in Fig. 11, which shows
a comparison between measured (black circles) and modeled
depth-weighted soil water contents. In particular, measured data
are compared with soil water contents calculated from the CRNP
using the standard calibration method (Baatz et al., 2014) (a grayline in Fig. 11) and predicted using the model with the median
SHPs from the posterior distribution for Scenario E1 (a black line
in Fig. 11), respectively. Depth-weighted soil water contents are
obtained using weights of 0.852, 0.146, and 0.002 at z = −5, −20,
and −50 cm, respectively, as reported by Baatz et al. (2014). Overall,both methods provide a good description of the near-surface soil
water contents as indicated by low RMSE values. An appreciable
overestimation of the water content can be observed for bothmethods between August and September 2012, with the magnitude being higher for the COSMIC-predicted water contents (ablack line in Fig. 11), which generally exhibit a lower variability
than water contents calculated using the standard calibration
method. Such difference confirms potential approximations in
the COSMIC modeling assumptions and an imperfect calibration of the COSMIC parameters, although deviation from themeasured data is limited.
Figure 12 shows the model predictive uncertainty (gray lines)
in the simulated neutron fluxes (Fig. 12A) and volumetric watercontents at three different depths. It is evident that both the discrepancy between measurement and model and the uncertainty
grow with depth. The model predictive uncertainty only partiallyencompasses the measured data, and the water contents are considerably underestimated at z = −50 cm. This behavior is expected
because the measured neutron fluxes are mostly informative ofnear-surface water dynamics. Furthermore, the results indicatethat the variability of water content is underestimated. This isparticularly visible between June and September 2012 at z = −5
Table 5. 95% Bayesian credible intervals of the soil hydraulic parameters
for the three experimental scenarios. The 2.5 and 97.5% quantiles are
the bounds of the credible intervals.
Scenario
E1

E2

E3

Parameter†

2.5%

97.5%

2.5%

97.5%

2.5%

97.5%

q s, cm3 cm−3

0.33

0.38

0.35

0.39

0.36

0.37

log10 a, cm−1

−2.40

−1.41

−1.87

−0.93

−1.67

−1.49

n

1.10

1.47

1.14

1.24

1.15

1.18

log10Ks, cm d−1

0.50

2.82

0.50

1.44

0.97

1.18

j h1

1.00

4.32

1.00

2.33

1.00

1.05

j h2

1.27

5.00

1.00

4.35

1.80

2.21

jK1

1.00

18.94

1.00

11.51

1.07

2.63

jK2

1.00

4.78

1.00

4.65

1.00

1.60

† q s, saturated water content; a and n, shape parameters; Ks, saturated hydraulic conductivity, j h1 and j K1, and j h2 and j K2 , pressure head and saturated
hydraulic conductivity scaling factors for Soil Materials 1 and 2, respectively.

VZJ | Advancing Critical Zone Science

and −20 cm, and it can be partially explained by the previously
described underestimation of the shape parameter a.
Overall, this analysis suggests that CRNP data can provide
reasonable estimates of the effective SHPs at the field scale, which
can be extremely valuable in land surface hydrological modeling.
In particular, the VGM shape parameters and the saturated watercontent are successfully inferred from the measured neutron flux,
thus leading to a coherent simulation of the soil moisture dynamics in the near surface, which is crucial for a proper description of
the surface energy balance in LSMs (Iwema et al., 2017). The large
uncertainty observed for Ks can significantly affect simulations ofoverland flow and the partitioning of precipitation into infiltration and surface runoff, and thus other types of measurements areneeded to better constrain this parameter. Nevertheless, the large
footprint of CRNPs represents an important advantage comparedwith classical near-surface point measurements, the information
content of which can be affected by various sources of uncertainty.
This aspect is better investigated in the final experimental data
analysis scenario.
Figure 13 shows the joint and marginal posterior distributions for Scenario E2 that relies on near-surface water content
measurements only. A multimodal posterior distribution withtwo peaks is observed for both a and Ks . The examination of
the joint posterior distribution of the pair a–Ks reveals the existence of two regions of high probability, one centered around thelocation [log10 a = −1.7, log10Ks = 0.5] and another one around[log10 a = −1.2, log10Ks = 1.3]. While the first region of high probability is in line with the values obtained for Scenarios E1 and E3
and thus probably reflects the true SHPs, the latter first region ofhigh probability significantly overestimates the shape parameter a.
Such predictions could be representative of another plausible poresystem induced by the plowing of the top layer. These sources of
uncertainty propagate in the credible interval (Table 5), which now
ranges between 0.01 and 0.11 cm−1, thus indicating poor identifiability of the a parameter. In this view, it must be emphasized thataboveground neutron fluxes are less influenced by local sources of
uncertainty than near-surface point measurements due to theirfootprint and penetration depth.
On the other hand, q s and n are well constrained and in linewith the values obtained in the previous two scenarios. A significant uncertainty is noticeable for the scaling factor j K1, whichranges between 1.0 and 11.51, despite being characterized by anapparently strong unimodal posterior distribution. In particular, aregion of high probability is evident around 1.0, followed by a flat
valley of low probability. Such posterior distributions are mainly
responsible for the large integrated autocorrelation time observedin Fig. 7. Nevertheless, a comparison of the posterior distributionsof the scaling factors of Soil Materials 1 and 2 confirms that themeasured water contents at z = −5 cm are mainly informative of
the soil hydraulic properties of the surface layer. Indeed, the posterior distributions of j h2 and j K2 are very similar to the priordistributions, whereas the posterior distributions clearly are muchsmaller for j h1 and j K1.
p. 16 of 24

VZJ | Advancing Critical Zone Science

p. 17 of 24

Fig. 9. A comparison between the model predictions (gray lines) obtained by random sampling of 500 solutions from the posterior distribution (Scenario E3) and measured (black circles) (A) neutron
fluxes (Nflux) and (B,C,D) volumetric water contents (q) at three different depths. The red lines indicate the model prediction with the median soil hydraulic parameters from the posterior distribution.

Fig. 10. Joint (below diagonal) and marginal (diagonal) posterior distributions of the soil hydraulic parameters of saturated soil water content (qs), thevan Genuchten–Mualem parameters n and a, saturated hydraulic conductivity (Ks), and scaling factors for hydraulic head and hydraulic conductivity(jh1 and jh2, and jK1 and jK2 for Soil Materials 1 and 2, respectively) for Scenario E1. The isolines and the dashed lines indicate the 12, 39, 68, and
84% probability levels and the Bayesian credible intervals, respectively. The axis ranges correspond to the parameter bounds reported in Table 3.

The propagation of uncertainty in the model predictionsfor Scenario E2 is shown in Fig. 14. It is again evident that theunderestimation of the soil retention capacity, caused by the lowinferred values of the shape parameter a, results in lower valuesof the water contents at z = −20 and −50 cm. The reduced retention capacity propagates in the simulated aboveground neutron
fluxes, which are significantly overestimated. However, the magnitude of the uncertainty is smaller than that obtained for thewater contents at deeper soil horizons, which is expected because

neutron fluxes are mostly influenced by the soil moisture dynamics in the near surface.
As a result of these analyses, a number of conclusions can be
drawn:
ʶʶ The results for Scenario E3 indicate that it is not possible tosimultaneously fit the measured neutron flux and soil watercontent data. This may be related to uncorrected effects ofaboveground biomass on the measured neutron data or inadequate modeling assumptions. In particular, the deviation

Fig. 11. A comparison between measured(black circles) and modeled (solid lines)
depth-weighted soil water contents (q). The
gray and black lines indicate soil water contents calculated from the cosmic-ray neutronprobes (CRNP) using the standard calibration method (Baatz et al., 2014) and predictedusing the model with the median soil hydraulic parameters from the posterior distribution
from Scenario E1, respectively.

VZJ | Advancing Critical Zone Sciencep. 18 of 24

VZJ | Advancing Critical Zone Sciencep. 19 of 24

Fig. 12. A comparison between the model predictions (gray lines) obtained by random sampling of 500 solutions from the posterior distribution (Scenario E1) and (A) measured (black circles) neutron
fluxes (Nfluxes) and (B,C,D) volumetric water contents (q) at three different depths. The red lines indicate the model prediction with the median soil hydraulic parameters from the posterior distribution.

Fig. 13. Joint (below diagonal) and marginal (diagonal) posterior distributions of the soil hydraulic parameters of saturated soil water content (qs), thevan Genuchten–Mualem parameters n and a, saturated hydraulic conductivity (Ks), and scaling factors for hydraulic head and hydraulic conductivity(jh1 and jh2, and jK1 and jK2 for Soil Materials 1 and 2, respectively) for Scenario E2. The isolines and the dashed lines indicate the 12, 39, 68, and
84% probability levels and the Bayesian credible intervals, respectively. The axis ranges correspond to the parameter bounds reported in Table 3.

between measured neutron fluxes and those simulated by
the calibrated HYDRUS model suggest potential inaccuracies in the COSMIC parameters. This could potentially be
addressed with a two-step sequential or a joint calibrationof the coupled HYDRUS–COSMIC model.
ʶʶ The use of measured aboveground neutron f luxes led toa reasonably confident estimation of the VGM shapeparameters and the saturated water content. From thispoint of view, the possibility of inversely estimating theeffective soil hydraulic parameters at the field scale fromthe CRNP data is confirmed. The combination of neutron
data with other independent types of measurements (e.g.,soil porosity, saturated conductivity) can further reduce
the parameter uncertainty.
ʶʶ A comparison of inversion results obtained for near-surface
volumetric water contents or neutron fluxes revealed thatthe latter are less influenced by typical sources of uncertainty
associated with land surface management activities (e.g., soil
tillage), thus leading to more reliable and representative estimates of soil hydraulic parameters.

VZJ | Advancing Critical Zone Science

66Conclusions and Summary
The main goal of this study was to assess the information
content of cosmic-ray neutron data for the inverse estimation of
the effective soil hydraulic parameters. An extensive theoreticaland experimental assessment based on the combination of the coupled HYDRUS–COSMIC model with a Bayesian data analysis
framework revealed that neutron fluxes could be used effectively toestimate the SHPs. In particular, the analysis of both synthetic andactual data from homogeneous and heterogeneous soil profiles led
to a confident estimation of the VGM shape parameters. Moreover,
the analysis demonstrated that the residual and saturated watercontents can be successfully inferred from measured neutron
fluxes if the soil approximately reaches a dry or wet state, respectively. On the other hand, high uncertainty was observed for thesaturated hydraulic conductivity, which suggests that the joint consideration of CRNP data with independent measurements of thisparameter may be required. The simultaneous use of both neutron
flux and water content data resulted in an appreciable reductionin uncertainty, in particular for the VGM shape parameters andp. 20 of 24

VZJ | Advancing Critical Zone Science

p. 21 of 24

Fig. 14. A comparison between the model predictions (gray lines) obtained by random sampling of 500 solutions from the posterior distribution (Scenario E2) and (A) measured (black circles) neutron
fluxes (Nfluxes) and (B,C,D) volumetric water contents (q) at three different depths. The red lines indicate the model prediction with the median soil hydraulic parameters from the posterior distribution.

for the saturated hydraulic conductivity. This result opens new
perspectives to estimate soil hydraulic properties from the combination of CRNP and other remote sensing techniques. However,to have a precise assessment of the benefits of such joint inversion
for LSMs, it is worth investigating numerically the propagationof the uncertainty in values of SHPs using LSMs. In this perspective, a model-based Bayesian experimental design could be used to
optimize the location and the number of cosmic-ray probes andtraditional sensors to minimize the effect of the uncertainty on
the land surface processes analyzed.
Furthermore, the comparison between the results obtainedwith the use of traditional near-surface point measurements and
CRNP data indicates that the latter data are less influenced by thetypical sources of uncertainty induced by land surface management
activities (e.g., plowing) and lead to a more representative estimation of the soil hydraulic parameters, which is fundamental when
dealing with land surface hydrological modeling. In particular, theshape parameter a, which characterizes the soil retention capacity, was significantly overestimated when inversions were solely
based on near-surface water content measurements, leading to abimodal posterior distribution representative of another plausible
pore system induced by former land management activities andmacropores. On the other hand, the CRNP penetration depth andfootprint significantly mitigate these issues. Nevertheless, to what
extent this difference in the estimation of the SHPs is significant
for practical purposes depends on the type of analysis performedby the modeler (e.g., surface runoff, infiltration, solute fluxes) andcan be clarified only by propagating the estimated uncertainty. In
this perspective, further studies are encouraged to better assess theinformation content gain of CRNP data in hydrological modeling.
In addition, the analysis presented here suggests that the
coupled HYDRUS–COSMIC model can be used to detectand correct an imperfect calibration of the COSMIC parameters. Avery et al. (2016) provided extensive results about thecorrelation between lattice water and the soil clay fraction. The
relationship reported in this study can be used to better constrain
COSMIC parameters. Moreover, future research should consider the coupling of the COSMIC model with HYDRUS-3D tobetter account for horizontal and vertical soil heterogeneity at thefield scale. In this perspective, the use of CRNP rover data could
improve the description of the spatial patterns of SHPs (Gibson
and Franz, 2018). A computationally cheap alternative would be
to use the kriging technique (Brunetti et al., 2017) to interpolate
the HYDRUS-1D predictions.

Acknowledgments

We thank the TERENO (Terrestrial Environmental Observatories) project fundedby the Helmholtz Association and the Transregional Collaborative Research Center 32 “Patterns in Soil–Vegetation–Atmosphere Systems” funded by Deutsche
Forschungsgemeinschaft (DFG) for providing data for this study.

References
Allen, R.G., L.S. Pereira, D. Raes, and M. Smith. 1998. Crop evapotranspiration.
Irrig. Drain. Pap. 56. FAO, Rome. doi:10.1016/j.eja.2010.12.001

VZJ | Advancing Critical Zone Science

Andreasen, M., K.H. Jensen, D. Desilets, T.E. Franz, M. Zreda, H.R. Bogena, and
M.C. Looms. 2017. Status and perspectives on the cosmic-ray neutron
method for soil moisture estimation and other environmental science
applications. Vadose Zone J. 16(8). doi:10.2136/vzj2017.04.0086
Avery, W.A., C. Finkenbiner, T.E. Franz, T. Wang, A.L. Nguy-Robertson,
A. Suyker, et al. 2016. Incorporation of globally available datasets
into the roving cosmic-ray neutron probe method for estimating
field-scale soil water content. Hydrol. Earth Syst. Sci. 20:3859–3872.
doi:10.5194/hess-20-3859-2016
Baatz, R., H.R. Bogena, H.J. Hendricks Franssen, J.A. Huisman, W. Qu, C.
Montzka, and H. Vereecken. 2014. Calibration of a catchment scalecosmic-ray probe network: A comparison of three parameterization
methods. J. Hydrol. 516:231–244. doi:10.1016/j.jhydrol.2014.02.026
Baatz, R., H.J.H. Franssen, X. Han, T. Hoar, H.R. Bogena, and H. Vereecken.
2017. Evaluation of a cosmic-ray neutron sensor network for improved
land surface model prediction. Hydrol. Earth Syst. Sci. 21:2509–2530.
doi:10.5194/hess-21-2509-2017
Betancourt, M. 2017. A conceptual introduction to Hamiltonian Monte
Carlo. https://arxiv.org/pdf/1701.02434.pdf.
Beven, K., and A. Binley. 1992. The future of distributed models: Modelcalibration and uncertainty prediction. Hydrol. Processes 6:279–298.
doi:10.1002/hyp.3360060305
Bogena, H.R., M. Herbst, J.A. Huisman, U. Rosenbaum, A. Weuthen, and H.
Vereecken. 2010. Potential of wireless sensor networks formeasuring soil water content variability. Vadose Zone J. 9:1002–1013.
doi:10.2136/vzj2009.0173
Bogena, H.R., J.A. Huisman, A. Güntner, C. Hübner, J. Kusche, F.
Jonard, et al. 2015. Emerging methods for noninvasive sensing of soil
moisture dynamics from field to catchment scale: A review. Wiley
Interdiscip. Rev. Water 2:635–647. doi:10.1002/wat2.1097
Bogena, H.R., J.A. Huisman, B. Schilling, A. Weuthen, and H. Vereecken.
2017. Effective calibration of low-cost soil water content sensors.
Sensors 17:208. doi:10.3390/s17010208
Bogena, H.R., R. Kunkel, S. Zacharias, T. Pütz, M. Schwank, O. Bens, et al.
2012. TERENO: Long-term monitoring network for terrestrial research.
Hydrol. Wasserbewirtsch. 56:138–143.
Bogena, H.R., C. Montzka, J.A. Huisman, A. Graf, M. Schmidt, M. Stockinger, etal. 2018. The TERENO-Rur Hydrological Observatory: A multi-scale multicompartment research platform for the advancement of hydrologicalscience. Vadose Zone J. 17:180055. doi:10.2136/vzj2018.03.0055
Brunetti, G., J. Šimůnek, and P. Piro. 2016. A comprehensive numericalanalysis of the hydraulic behavior of a permeable pavement. J. Hydrol.
540:1146–1161. doi:10.1016/j.jhydrol.2016.07.030
Brunetti, G., J. Šimůnek, M. Turco, and P. Piro. 2017. On the use of surrogatebased modeling for the numerical analysis of low impact development
techniques. J. Hydrol. 548:263–277. doi:10.1016/j.jhydrol.2017.03.013
Carsel, R.F., and R.S. Parrish. 1988. Developing joint probability
distributions of soil water retention characteristics. Water Resour. Res.
24:755–769 doi:10.1029/WR024i005p00755
Dane, J.H., and S. Hruska. 1983. In-situ determination of soil hydraulicproperties during drainage. Soil Sci. Soc. Am. J. 47:619–624.
doi:10.2136/sssaj1983.03615995004700040001x
Deng, H., M. Ye, M.G. Schaap, and R. Khaleel. 2009. Quantification of
uncertainty in pedotransfer function-based parameter estimation
for unsaturated flow modeling. Water Resour. Res. 45:W04409.
doi:10.1029/2008WR007477
Desilets, D., M. Zreda, and T.P.A. Ferré. 2010. Nature’s neutron probe: Landsurface hydrology at an elusive scale with cosmic rays. Water Resour.
Res. 46:W11505. doi:10.1029/2009WR008726
Ek, M.B., K.E. Mitchell, Y. Lin, E. Rogers, P. Grunmann, V. Koren, et al. 2003.
Implementation of Noah land surface model advances in the National
Centers for Environmental Prediction operational mesoscale Eta
model. J. Geophys. Res. 108:8851. doi:10.1029/2002JD003296
Evans, J.G., H.C. Ward, J.R. Blake, E.J. Hewitt, R. Morrison, M. Fry, et al. 2016.
Soil water content in southern England derived from a cosmic-ray soil
moisture observing system: COSMOS-UK. Hydrol. Process. 30:4987–
4999. doi:10.1002/hyp.10929
Feddes, R.A., P.J. Kowalik, and H. Zaradny. 1978. Simulation of field waterp. 22 of 24

use and crop yield. PUDOC, Wageningen, the Netherlands.
Feroz, F., M.P. Hobson, and M. Bridges. 2009. MultiNest: An efficient and robust
Bayesian inference tool for cosmology and particle physics. Mon. Not. R.
Astron. Soc. 398:1601–1614. doi:10.1111/j.1365-2966.2009.14548.x
Finkenbiner, C.E., T.E. Franz, J. Gibson, D.M. Heeren, and J. Luck. 2018.
Integration of hydrogeophysical datasets and empirical orthogonal
functions for improved irrigation water management. Precis. Agric. (inpress). doi:10.1007/s11119-018-9582-5
Foreman-Mackey, D., D.W. Hogg, D. Lang, and J. Goodman. 2012.
emcee: The MCMC hammer. Publ. Astron. Soc. Pac. 125:306–312.
doi:10.1086/670067
Franz, T.E., M. Zreda, T.P.A. Ferre, R. Rosolem, C. Zweck, S. Stillman, X.
Zeng, and W.J. Shuttleworth. 2012. Measurement depth of the cosmic
ray soil moisture probe affected by hydrogen from various sources.
Water Resour. Res. 48:W08515. doi:10.1029/2012WR011871
Gibson, J., and T.E. Franz. 2018. Spatial prediction of near surface soil water
retention functions using hydrogeophysics and empirical orthogonal
functions. J. Hydrol. 561:372–383. doi:10.1016/j.jhydrol.2018.03.046
Goodman, J., and J. Weare. 2010. Ensemble samplers with affine invariance.
Comm. App. Math. Comp. Sci. 5:65–80. doi:10.2140/camcos.2010.5.65
Gutmann, E.D., and E.E. Small. 2007. A comparison of land surface
model soil hydraulic properties estimated by inverse modelingand pedotransfer functions. Water Resour. Res. 43:W05418.
doi:10.1029/2006WR005135
Hawdon, A., D. McJannet, and J. Wallace. 2014. Calibration andcorrection procedures for cosmic-ray neutron soil moistureprobes located across Australia. Water Resour. Res. 50:5029–5043.
doi:10.1002/2013WR015138
Hopmans, J.W., D.R. Nielsen, and K.L. Bristow. 2002. How useful are smallscale soil hydraulic property measurements for large-scale vadosezone modeling? In P.A.C. Raats et al., editor, Environmental mechanics:
Water, mass and energy transfer in the biosphere. Geophys.
Monogr. Ser. 129. Am. Geophys. Union, Washington, DC. p. 247–258.
doi:10.1029/129GM20
Huijser, D., J. Goodman, and B.J. Brewer. 2015. Properties of the affine
invariant ensemble sampler in high dimensions. arXiv:1509.02230.
Huisman, J.A., J. Rings, J.A. Vrugt, J. Sorg, and H. Vereecken. 2010.
Hydraulic properties of a model dike from coupled Bayesian andmulti-criteria hydrogeophysical inversion. J. Hydrol. 380:62–73.
doi:10.1016/j.jhydrol.2009.10.023
Iden, S.C., and W. Durner. 2007. Free-form estimation of the unsaturated
soil hydraulic properties by inverse modeling using global optimization.
Water Resour. Res. 43:W07451. doi:10.1029/2006WR005845
Ines, A.V.M., and B.P. Mohanty. 2008. Near-surface soil moistureassimilation for quantifying effective soil hydraulic properties
under different hydroclimatic conditions. Vadose Zone. J. 7:39–52.
doi:10.2136/vzj2007.0048
Iwema, J., R. Rosolem, M. Rahman, E. Blyth, and T. Wagener. 2017. Landsurface model performance using cosmic-ray and point-scale soil
moisture measurements for calibration. Hydrol. Earth Syst. Sci.
21:2843–2861. doi:10.5194/hess-21-2843-2017
Kavetski, D., G. Kuczera, and S.W. Franks. 2006. Bayesian analysis of input
uncertainty in hydrological modeling: 1. Theory. Water Resour. Res.
42:W03407. doi:10.1029/2005WR004368
Kennedy, J., and R. Eberhart. 1995. Particle swarm optimization. Eng.
Technol. 1942–1948. doi:10.1109/ICNN.1995.488968
Köhli, M., M. Schrön, M. Zreda, U. Schmidt, P. Dietrich, and S. Zacharias.
2015. Footprint characteristics revised for field-scale soil moisturemonitoring with cosmic-ray neutrons. Water Resour. Res. 51:5772–
5790 . doi:10.1002/2015WR017169
Kool, J.B., and J.C. Parker. 1988. Analysis of the inverse problem
for transient unsaturated flow. Water Resour. Res. 24:817–830.
doi:10.1029/WR024i006p00817
Lam, S.K., A. Pitrou, and S. Seibert. 2015. Numba: A LLVM-based Python JIT
compiler. In: Proceedings of the 2nd Workshop on the LLVM Compiler
Infrastructure in HPC, Austin, TX. 15–20 Nov. 2015. Assoc. Comput.
Mach., New York. Article no. 7. doi:10.1145/2833157.2833162
Mantovan, P., and E. Todini. 2006. Hydrological forecasting uncertainty

VZJ | Advancing Critical Zone Science

assessment: Incoherence of the GLUE methodology. J. Hydrol.
330:368–381. doi:10.1016/j.jhydrol.2006.04.046
Matthews, C., J. Weare, A. Kravtsov, and E. Jennings. 2017. Umbrellasampling: A powerful method to sample tails of distributions. Mon. Not.
R. Astron. Soc. 480:4069–4079. doi:10.1093/mnras/sty2140
Metropolis, N., A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller, and E. Teller.
1953. Equation of state calculations by fast computing machines. J.
Chem. Phys. 21:1087–1092. doi:10.1063/1.1699114
Mockler, E.M., K.P. Chun, G. Sapriza-Azuri, M. Bruen, and H.S. Wheater. 2016.
Assessing the relative importance of parameter and forcing uncertainty
and their interactions in conceptual hydrological model simulations.
Adv. Water Resour. 97:299–313. doi:10.1016/j.advwatres.2016.10.008
Mohanty, B.P. 2013. Soil hydraulic property estimation using remotesensing: A review. Vadose Zone J. 12(4). doi:10.2136/vzj2013.06.0100
Mohanty, B.P., M.H. Cosh, V. Lakshmi, and C. Montzka. 2017. Soil
moisture remote sensing: State-of-the-science. Vadose Zone J. 16(1).
doi:10.2136/vzj2016.10.0105
Montanari, A. 2007. What do we mean by ‘uncertainty’? The need fora consistent wording about uncertainty assessment in hydrology.
Hydrol. Process. 21:841–845. doi:10.1002/hyp.6623
Pelowitz, D.B. 2011. MCNPX user’s manual version 2.7.0. LA-CP-11-00438.
Los Alamos Natl. Lab., Los Alamos, NM.
Pitman, A.J. 2003. The evolution of, and revolution in, land surface
schemes designed for climate models. Int. J. Climatol. 23:479–510.
doi:10.1002/joc.893
Qu, W., H.R. Bogena, J.A. Huisman, G. Martinez, Y.A. Pachepsky, and H.
Vereecken. 2014. Effects of soil hydraulic properties on the spatial
variability of soil water content: Evidence from sensor network data andinverse modeling. Vadose Zone J. 13(12). doi:10.2136/vzj2014.07.0099
Qu, W., H.R. Bogena, J.A. Huisman, M. Schmidt, R. Kunkel, A. Weuthen,et al. 2016. The integrated water balance and soil data set of theRollesbroich Hydrological Observatory. Earth Syst. Sci. Data 8:517–529.
doi:10.5194/essd-8-517-2016
Ritchie, J.T. 1972. Model for predicting evaporation from a row
crop with incomplete cover. Water Resour. Res. 8:1204–1213.
doi:10.1029/WR008i005p01204
Rivera Villarreyes, C.A., G. Baroni, and S.E. Oswald. 2014. Inverse modellingof cosmic-ray soil moisture for field-scale soil hydraulic parameters.
Eur. J. Soil Sci. 65:876–886. doi:10.1111/ejss.12162
Santanello, J.A., C.D. Peters-Lidard, M.E. Garcia, D.M. Mocko, M.A. Tischler,
M.S. Moran, and D.P. Thoma. 2007. Using remotely-sensed estimatesof soil moisture to infer soil texture and hydraulic properties
across a semi-arid watershed. Remote Sens. Environ. 110:79–97.
doi:10.1016/j.rse.2007.02.007
Schaap, M.G., F.J. Leij, and M.Th. van Genuchten. 2001. Rosetta: A
computer program for estimating soil hydraulic parameterswith hierarchical pedotransfer functions. J. Hydrol. 251:163–176.
doi:10.1016/S0022-1694(01)00466-8
Scharnagl, B., J.A. Vrugt, H. Vereecken, and M. Herbst. 2011. Inverse
modelling of in situ soil water dynamics: Investigating the effect ofdifferent prior distributions of the soil hydraulic parameters. Hydrol.
Earth Syst. Sci. 15:3043–3059. doi:10.5194/hess-15-3043-2011
Schrön, M., M. Köhli, L. Scheiffele, J. Iwema, H.R. Bogena, L. Lv, et al. 2017.
Improving calibration and validation of cosmic-ray neutron sensors
in the light of spatial sensitivity. Hydrol. Earth Syst. Sci. 21:5009–5030.
doi:10.5194/hess-21-5009-2017
Shuttleworth, J., R. Rosolem, M. Zreda, and T. Franz. 2013. The COsmic-ray
Soil Moisture Interaction Code (COSMIC) for use in data assimilation.
Hydrol. Earth Syst. Sci. 17:3205–3217. doi:10.5194/hess-17-3205-2013Šimůnek, J., and M.Th. van Genuchten. 1996. Estimating unsaturated soil
hydraulic properties from tension disc infiltrometer data by numericalinversion. Water Resour. Res. 32:2683–2696. doi:10.1029/96WR01525
Šimůnek, J., M.Th. van Genuchten, and M. Šejna. 2016. Recent
developments and applications of the HYDRUS computer software
packages. Vadose Zone J. 15(7). doi:10.2136/vzj2016.04.0033
Sokal, A.D. 1996. Monte Carlo methods in statistical mechanics:
Foundations and new algorithms. NATO ASI Ser., Ser. B 361:131–192.
doi:10.1007/978-1-4899-0319-8_6

p. 23 of 24

Steenpass, C., J. Vanderborght, M. Herbst, J. Šimůnek, and H. Vereecken.
2010. Estimating soil hydraulic properties from infrared measurementsof soil surface temperatures and TDR data. Vadose Zone J. 9:910–924.
doi:10.2136/vzj2009.0176

Vrugt, J.A., P.H. Stauffer, T. Wöhling, B.A. Robinson, and V.V. Vesselinov.
2008a. Inverse modeling of subsurface flow and transport properties:
A review with new developments. Vadose Zone J. 7:843–864.
doi:10.2136/vzj2007.0078

Sutanto, S.J., J. Wenninger, A.M.J. Coenders-Gerrits, and S. Uhlenbrook.
2012. Partitioning of evaporation into transpiration, soil evaporationand interception: A comparison between isotope measurements
and a HYDRUS-1D model. Hydrol. Earth Syst. Sci. 16:2605–2616.
doi:10.5194/hess-16-2605-2012

Vrugt, J.A., C.J.F. ter Braak, M.P. Clark, J.M. Hyman, and B.A. Robinson.
2008b. Treatment of input uncertainty in hydrologic modeling: Doinghydrology backward with Markov chain Monte Carlo simulation. Water
Resour. Res. 44: W00B09. doi:10.1029/2007WR006720

Taylor, S.A., and G.L. Ashcroft. 1972. Physical edaphology: The physics ofirrigated and nonirrigated soils. W.H. Freeman, San Francisco.
van Genuchten, M.Th. 1980. A closed-form equation for predicting thehydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44:892–
898. doi:10.2136/sssaj1980.03615995004400050002x
Van Looy, K., J. Bouma, M. Herbst, J. Koestel, B. Minasny, U. Mishra, et al.
2017. Pedotransfer functions in Earth system science: Challenges andperspectives. Rev. Geophys. 55:1199–1256. doi:10.1002/2017RG000581
Vereecken, H., R. Kasteel, J. Vanderborght, and T. Harter. 2007. Upscaling
hydraulic properties and soil water flow processes in heterogeneous
soils. Vadose Zone J. 6:1–28. doi:10.2136/vzj2006.0055
Vogel, T., M. Cislerova, and J.W. Hopmans. 1991. Porous media withlinearly variable hydraulic properties. Water Resour. Res. 27:2735–2741.
doi:10.1029/91WR01676
Vrugt, J.A., and B.A. Robinson. 2007. Improved evolutionary optimizationfrom genetically adaptive multimethod search. Proc. Natl. Acad. Sci.
104:708–711. doi:10.1073/pnas.0610471104

VZJ | Advancing Critical Zone Science

Wesseling, J., J. Elbers, P. Kabat, and B. Van den Broek. 1991. SWATRE:
Instructions for input. Winand Staring Ctr., Wageningen, the Netherlands.
Western, A.W., S.-L. Zhou, R.B. Grayson, T.A. McMahon, G. Blöschl, and D.J.
Wilson. 2004. Spatial correlation of soil moisture in small catchments
and its relationship to dominant spatial hydrological processes. J.
Hydrol. 286:113–134. doi:10.1016/j.jhydrol.2003.09.014
Wöhling, T., and J.A. Vrugt. 2011. Multiresponse multilayer vadose zone
model calibration using Markov chain Monte Carlo simulationand field water retention data. Water Resour. Res. 47:W04510.
doi:10.1029/2010WR009265
Zreda, M., D. Desilets, T.P.A. Ferré, and R.L. Scott. 2008. Measuring
soil moisture content non-invasively at intermediate spatial
scale using cosmic-ray neutrons. Geophys. Res. Lett. 35:L21402
doi:10.1029/2008GL035655
Zreda, M., W.J. Shuttleworth, X. Zeng, C. Zweck, D. Desilets, T. Franz, and R.
Rosolem. 2012. COSMOS: The cosmic-ray soil moisture observing system.
Hydrol. Earth Syst. Sci. 16:4079–4099. doi:10.5194/hess-16-4079-2012

p. 24 of 24