Page 1 of 72

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

1

ON THE INFORMATION CONTENT OF COSMIC-RAY NEUTRON DATA IN THE

2

INVERSE ESTIMATION OF SOIL HYDRAULIC PROPERTIES

3
Giuseppe Brunettia, Jiří Šimůnekb, Heye Bogenac, Roland Baatzc, Johan Alexander

4

Huismanc, Helen Dahlkea, and Harry Vereeckenc*

5
6
7
8
9
10
11

a

University of California, Davis, Department of Land, Air and, Water Resources, CA 95616, USAb
University of California, Riverside, Department of Environmental Sciences, CA 92521, USAc Agrosphere Institute, Forschungszentrum Jülich GmbH, Germany*Corresponding Author

12

CORE IDEAS

13•

Cosmic-ray neutron data are used to inversely estimate soil hydraulic properties

14•

The forward neutron operator COSMIC is coupled with the vadose zone model

15
16

HYDRUS-1D•

Bayesian analyses confirm the information content of cosmic-ray neutron data

17
18

ABSTRACT

19

Observations of soil moisture content from remote sensing platforms can be used in20

conjunction with hydrological models to inversely estimate soil hydraulic properties (SHPs). In

21

recent years, cosmic-ray neutron sensing (CRNS) has proven to be a reliable method for the22

estimation of area-average soil moisture at field scales. However, their use in the inverse estimation

23of the effective SHPs is largely unexplored. Thus, the main objective of this study is to assess the24

information content of aboveground fast-neutron counts to estimate SHPs using both a synthetic

25

modeling study and actual experimental data from the Rollesbroich catchment in Germany. For this,26

the forward neutron operator COSMIC was externally coupled with the hydrological model27

HYDRUS-1D. The coupled model was combined with the Affine Invariant Ensemble Sampler to28

calculate the posterior distributions of effective soil hydraulic parameters as well as the model1

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

29

predictive uncertainty for different synthetic and experimental scenarios. Measured water contents30

at different depths were used to assess estimated SHPs. The analysis of both synthetic and actual

31

CRNS data from homogenous and heterogeneous soil profiles, respectively, led to confident

32

estimations of the shape parameters α and n, while higher uncertainty was observed for the33

saturated hydraulic conductivity. Furthermore, results demonstrated that neutron data are less34

influenced by local sources of uncertainty compared to near-surface point measurements. The

35

simultaneous use of CRNS and water content data further reduced the overall uncertainty, opening36

up new perspectives for the combination of CRNS with other remote sensing techniques for the37

inverse estimation of the effective SHPs.

38
39

Keywords: Cosmic ray neutron probe, soil hydraulic properties, COSMIC, HYDRUS, Bayesian40

inference

2

Page 2 of 72

Page 3 of 72

41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

Table of Contents
ON THE INFORMATION CONTENT OF COSMIC-RAY NEUTRON DATA IN THE INVERSE
ESTIMATION OF SOIL HYDRAULIC PROPERTIES ................................................................. 1
CORE IDEAS ................................................................................................................................. 1
ABSTRACT ................................................................................................................................... 1
INTRODUCTION .......................................................................................................................... 4
MATERIALS AND METHODS..................................................................................................... 8
Theory ...................................................................................................................................... 8
Subsurface Water Flow and Root Water Uptake .................................................................. 8
COsmic-ray Soil Moisture Interaction Code (COSMIC) .................................................... 10
Model Coupling Strategy ................................................................................................... 11
Uncertainty Analysis ............................................................................................................... 12
Bayesian Inference ............................................................................................................ 12
MCMC Diagnostic ............................................................................................................ 15
Likelihood Function .......................................................................................................... 17
Rollesbroich grassland experimental catchment ....................................................................... 17
Numerical Domain and Boundary Conditions.......................................................................... 19
Inverse estimation with synthetic data ..................................................................................... 20
Model Set-up ..................................................................................................................... 20
Synthetic Modeling Scenarios ........................................................................................... 22
Inverse estimation with actual data .......................................................................................... 24
Model Set-up ..................................................................................................................... 24
Experimental Modeling Scenarios ..................................................................................... 24
RESULTS AND DISCUSSION .................................................................................................... 25
Inverse estimation with synthetic data ..................................................................................... 25
MCMC Diagnostic ............................................................................................................ 25
Uncertainty analysis .......................................................................................................... 26
Inverse estimation with actual data .......................................................................................... 33
MCMC Diagnostic ............................................................................................................ 33
Uncertainty analysis .......................................................................................................... 34
CONCLUSIONS AND SUMMARY ............................................................................................ 46
ACKNOWLEDGMENTS ............................................................................................................. 48
REFERENCES ............................................................................................................................. 48

75

3

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

76
77

INTRODUCTION
Soil hydraulic properties (SHPs) define the relationship between volumetric water content

78(θ), pressure head (h), and hydraulic conductivity (K), thus regulating the movement of water in the79

soil and influencing the water-energy cycle at the land surface. In particular, the near-surface soil

80

moisture has a profound influence on the partitioning of precipitation into surface runoff,81

evapotranspiration, and infiltration, and thus on the surface energy balance. In this view, there is82

strong evidence that land surface processes play a key role in climate models (Pitman, 2003). An

83

accurate description of hydrological processes at the land surface is thus of crucial importance.

84

Land Surface Models (LSMs) incorporating a mechanistic description of infiltration require85

precise estimation of large-scale SHPs. Typically, the measurement scale for the characterization of

86

SHPs is on the order of 10 cm, with a sample spacing of 100 m or larger (Hopmans et al., 2002). In

87

this perspective, the upscaling of point measurements is time-consuming and introduces significant

88

uncertainty in land surface hydrological modeling (Vereecken et al., 2007). To overcome this89

problem, pedotransfer functions (PTFs) are often used (Van Looy et al., 2017) to translate available

90

soil textural information into SHPs. However, uncertainty in PTF inputs often propagates into the91

poor accuracy of SHPs (Deng et al., 2009).

92

A valid alternative is the inverse estimation of effective soil hydraulic properties from93

transient data (Dane and Hruska, 1983; Kool and Parker, 1988; Šimůnek and Van Genuchten, 1996;

94

Iden and Durner, 2007; Ines and Mohanty, 2008; Vrugt et al., 2008a). For example, Gutmann and95

Small (2007) used the land surface model Noah (Ek et al., 2003) to compare the results obtained

96

using inversely estimated and PTF-based SHPs. They found a significant reduction in errors

97

between measured and simulated latent heat fluxes when inversely estimated SHPs were used.

98

Moreover, the use of emerging global optimization algorithms and statistically rigorous methods

99

based on Bayesian inference now allow obtaining both a reliable estimate of SHPs and a proper100

assessment of their uncertainty (e.g., Huisman et al., 2010; Scharnagl et al., 2011; Wöhling and4

Page 4 of 72

Page 5 of 72

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

101

Vrugt, 2011), which is directly related to the information content of the measured data used in the102

inverse estimation framework.

103

Traditionally, point measurements of pressure head and water content have been used for the104

inverse estimation of SHPs. However, during the last decade, soil moisture remote sensing

105

techniques have gained popularity (Bogena et al., 2015; Mohanty et al., 2017). Despite their106

restriction to near-surface soil moisture measurements, several studies have demonstrated that107

remote sensing data can be used to successfully estimate effective SHPs at the field scale (Mohanty,

1082013). For example, Santanello et al. (2007) inversely estimated soil hydraulic parameters of a109

semiarid watershed using the Noah Land Surface Model by combining soil moisture evaluated110

using passive microwave remote sensing by NASA’s push-broom microwave radiometer and active

111

microwave RADARSAT-1 imagery with a parameter estimation algorithm and then compared the112

estimated and PTF-based soil hydraulic parameters. In another study, Steenpass et al. (2010)

113

assimilated infrared-measured soil surface temperatures and time domain reflectometry–based soil

114

water contents to estimate SHPs for both synthetic and field scale conditions.

115

Among the ground-based remote sensing techniques, the Cosmic-Ray Neutron Sensing

116(CRNS) has shown good promise to capture soil moisture at relevant scales. Briefly, fast neutrons

117

are produced by nuclear interactions between the incoming cosmic rays and elements of the Earth’s

118

atmosphere (Zreda et al., 2008). When they reach the soil surface, they penetrate to a certain depth

119

and are scattered back into the atmosphere. Since fast neutrons are mainly moderated by hydrogen,120

the fast neutron intensity at a near-ground-level is negatively correlated with the near-surface soil

121

moisture. For this reason, soil moisture can be inferred from the fast neutron intensity measured by

122

Cosmic-Ray Neutron Probes (CRNPs) (Zreda et al., 2008). The CRNP measures integral soil

123

moisture in a circular footprint centered on a detector. This radius varies between 130 and 240 m

124

depending on the site conditions, which is significantly larger than the typical spatial correlation

125

length of soil moisture patterns (Western et al., 2004). The penetration depth of CRNP

126

measurements varies between 15 cm for wet soils to 55 cm for dry soils (Schrön et al., 2017). Due
5Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

127

to this large sample volume, the information content of CRNP data could potentially improve the128

estimation of effective SHPs for land-surface models because the scale mismatch between sample

129

volume and model resolution is much smaller.

130

The applicability of the CRNP method has been demonstrated for various land cover types,

131

soil types, and climates, and CRNPs have been installed in several national monitoring networks

132(Bogena et al., 2015). The first network was established in the USA by the University of Arizona133

and has already deployed more than 60 CRNPs at various locations (Zreda et al., 2012). Similar

134

networks have been established in Australia (Hawdon et al., 2014), the United Kingdom (Evans et135al., 2016), and Germany (Baatz et al., 2014). To date, almost 200 stationary CRNPs have been136

installed worldwide (Andreasen et al., 2017). The required conversion of neutron intensity into soil

137

moisture and the correction for important influencing factors (e.g., biomass, snow, and litter layer

138

water content) are described in detail in the literature (e.g., Andreasen et al., 2017; Desilets et al.,1392010). Few recent studies have focused on the use of CRNP data for the estimation of SHPs with140

promising results. Both Finkenbiner et al. (2018) and Gibson and Franz (2018) combined CRNP

141

data with empirical orthogonal functions to isolate spatial patterns of SHPs in selected agricultural

142

fields. A statistical analysis confirmed that the use of CRNP data can improve the estimate of SHPs

143

at the field scale.

144

Several types of neutron source models are available to study cosmic-ray neutron

145

interactions near the soil surface. One of the most widely used models is the Monte Carlo N-Particle

146

eXtended (MCNPX) neutron transport code (Pelowitz, 2011). Given specified chemistry of the147

atmosphere and soil (e.g., soil mineralogy, soil moisture), MCNPX simulates nuclear collisions of148

individual, randomly generated, incoming cosmic rays and their interaction byproducts through the149

atmosphere and in the soil using libraries of nuclear properties. In this way, MCNPX allows to150

accurately simulate the number of fast neutrons that enter a defined detector volume above the151

ground. In an attempt to simultaneously describe the neutrons’ dynamics and water flow in the152

unsaturated zone, Franz et al. (2012) successfully coupled MCNPX with the hydrological model
6

Page 6 of 72

Page 7 of 72

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

153

HYDRUS-1D (Šimůnek et al., 2016) for calculating the effective sensor depth of a CRNP while154

accounting for three sources of hydrogen and their vertical variability. However, the high155

computational cost of this Monte Carlo based model makes it unsuited for Bayesian analysis

156

frameworks, which require thousands of model executions. Recently, the Ultra Rapid Adaptable

157

Neutron-Only Simulation (URANOS) model (Köhli et al., 2015) was developed to reduce the158

computational cost of CRNP simulations and to provide a neutron transport code tailored to159

environmental applications. By neglecting less relevant physical processes, URANOS is more160

computationally efficient than other available Monte Carlo codes while maintaining a similar

161

accuracy in the simulation of neutron fluxes (Köhli et al., 2015). However, the computational gain

162

is still not sufficient to justify its use in Bayesian analysis.

163

The forward neutron operator COSMIC has taken a major step toward the efficient164

integration of CRNS data in models (Shuttleworth et al., 2013). It provides a simplified description

165

of neutron interactions, and the use of a physically-based analytical approach reduces the166

computational cost by four orders of magnitude compared to MCNPX. COSMIC calculates the fast167

neutron flux from the soil moisture profile and can be easily coupled with hydrological models.

168

Recently, Baatz et al. (2017) were able to assimilate data from a network of CRNPs installed in the169

River Rur catchment into the Common Land Model (CLM4.5) in order to update soil moisture and170

SHP of the catchment. In another study, the COSMIC and Noah models were coupled to assimilate

171

CRNP data from the Santa Rita Range field site (Shuttleworth et al., 2013). In particular, the Noah-172

COSMIC model was used in conjunction with the global optimization algorithm AMALGAM

173(Vrugt and Robinson, 2007) to recalibrate the Noah model by minimizing the errors between

174

simulated and observed neutron fluxes. Their results demonstrated that the assimilation of CRNP

175

data could effectively improve the accuracy of the model predictions. To date, the only available

176

studies on the use of CRNP data in the inverse estimation of soil hydraulic properties are those of177

Rivera Villarreyes et al. (2014) and Baatz et al. (2017). In the former, integral soil moisture178

estimates obtained for the CRNP support volume were transformed into the soil water storage for7

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

179

direct calibration of the effective soil hydraulic parameters using the HYDRUS-1D model. Despite180

the overall promising results of their study, the numerical approach based on the use of soil water

181

storage data led to a poor approximation of the volumetric water content at different depths. In

182

Baatz et al. (2017), soil water contents estimated from CRNS data were assimilated with the local

183

ensemble transform Kalman filter in the CLM 4.5 model to update soil texture and organic matter

184

used in the PTFs for the estimation of the SHPs. Results of that study demonstrated the potential of

185

CRNS networks to improve subsurface parameterization in the regional land surface model.

186

However, the use of PTFs and the indirect assimilation of CRNS data represent potential sources of187

additional uncertainty in the data assimilation framework.

188

Thus, the main aim of this study is to provide a comprehensive theoretical and experimental

189

assessment of the information content of CRNP data for the inverse estimation of effective SHP at

190

the field scale. The problem is addressed in the following way. First, a combined mechanistic

191

description of variably-saturated flow and soil-neutron interactions is provided by coupling the192

hydrological model HYDRUS-1D and the neutron transport code COSMIC. Before proceeding193

with an analysis of actual CRNP data from an experimental facility in Germany, an analysis of

194

synthetic data is performed by combining the coupled HYDRUS-COSMIC model with a Bayesian

195

inference algorithm. The synthetic example simplifies the experimental case study by describing the196

soil-neutrons interactions in a homogeneous soil profile and provides the first proof of concept for197

the estimation of SHP from CRNP data. Finally, the actual data from a heterogeneous soil profile is

198

used to evaluate the potential of CRNP data for the inverse estimation of SHPs in real conditions. A

199

comparison between inverse estimation using neutron fluxes and classical point measurements of200

water content completes and further enriches the analysis.

201
202

MATERIALS AND METHODS

203

Theory204

Subsurface Water Flow and Root Water Uptake
8Page 8 of 72

Page 9 of 72

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

205

The finite-element hydrological model HYDRUS-1D (Šimůnek et al., 2016) numerically206

solves the one-dimensional Richards equation describing variably-saturated water flow in the soil

207

profile:∂θ∂  ∂h 
=  K (h) + 1 − S∂t ∂z 
 ∂z


(1)

208

where θ is the volumetric water content [L3L-3], h is the soil water pressure head [L], K(h) is the209

unsaturated hydraulic conductivity [LT-1], t is time [T], z is the soil depth [L] and S is a sink term

210

[L3L-3T-1], defined as a volume of water removed from a unit volume of soil per unit of time due to211plant water uptake. A solution of the Richards equation requires knowledge of the SHPs, which212

should accurately describe the hydraulic characteristics of the porous medium and its heterogeneity.

213

It is well known that agriculture land management can significantly influence the soil structure and214

corresponding SHPs in the near surface. Under such circumstances, the shape of the retention curve215

may be constant with depth, but the soil hydraulic conductivity and porosity are expected to change

216

significantly. Thus, in the present study, the vertical variability of the soil is described by combining217

the unimodal van Genuchten-Mualem (VGM) (van Genuchten, 1980) functions with the scaling

218

factors proposed by Vogel et al. (1991). In particular, three independent scaling factors are used: φθ219(-), φh (-), and φK (-) for the water content, the pressure head, and the saturated hydraulic220

conductivity, respectively. These three scaling parameters are used to define a linear model of the221

actual spatial variability in the soil hydraulic properties as follows:[θ (h) = θ r + ϕθ θ * (h* ) − θ r**
 h = ϕh h
 K ( h) = ϕ K * ( h * )
K


]
(2)

222

where θ*(h*) and K*(h*) are the reference soil hydraulic functions, which are ideally those223

unaffected by surface management activities. This theoretical approach guarantees high modeling224

flexibility since it can describe both homogeneous and heterogeneous soil profiles using a single set

9Page 10 of 72

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

225

of reference soil hydraulic functions and only few more parameters (i.e., the scaling factors). If the226

soil is homogeneous, the scaling factors are set to 1.

227

Feddes et al. (1978) defined S as:
S (h ) = a (h ) ⋅ S p

(3)

228

where Sp is the potential root water uptake rate and a(h) is a dimensionless water stress response229

function that depends on the soil pressure head h and can vary between 0 and 1. Parameters of the230

stress response function for a majority of agricultural crops can be found in various databases (e.g.,231

Taylor and Ashcroft, 1972; Wesseling et al., 1991). Potential root water uptake Sp is calculated from232

the potential transpiration rate Tp and the root distribution function. Beer’s equation is used to233

partition the reference evapotranspiration calculated with the Penman-Monteith equation (Allen et234al., 1998) into potential transpiration and soil evaporation fluxes (e.g., Ritchie, 1972) using Leaf

235

Area Index (LAI) information. For a detailed explanation of the partitioning of evapotranspiration,

236

please refer to Sutanto et al. (2012).

237
238

COsmic-ray Soil Moisture Interaction Code (COSMIC)

239

The physically-based COsmic-ray Soil Moisture Interaction Code (COSMIC) (Shuttleworth

240

et al., 2013) is used to calculate the aboveground neutron intensity for a given soil moisture profile.

241

COSMIC includes simple descriptions of (a) degradation of the incoming high-energy neutron flux

242

with soil depth, (b) creation of fast neutrons at each depth in the soil, and (c) scattering of the243

resulting fast neutrons before they reach the soil surface. We refer to Shuttleworth et al. (2013) for a244

detailed description of the COSMIC model.

245

COSMIC requires several site-independent and site-specific time-invariable parameters. In

246

particular, the parameters L1 = 162.0 g/cm2, L2 = 129.1 g/cm2, and L4 = 3.16 g/cm2 were found to be247

site-independent (Shuttleworth et al., 2013), while the parameters L3 (g/cm2) and α (cm3/g) were248

found to be dependent on the soil bulk density ρb according to:

10

Page 11 of 72

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

L3 = −31.65 + 99.29ρb(4)

α = 0.404 − 0.101ρb

(5)

249

The remaining two parameters, the lattice water content θl [L3L-3] and the number of high-energy

250

neutrons at the soil surface NCOSMIC [-], are also site-dependent and must be measured and251

calibrated, respectively.

252

A FORTRAN version of the COSMIC code is freely available online

253(http://cosmos.hwr.arizona.edu). For this study, the FORTRAN code was rewritten in the Python254

language to enable a more efficient model coupling with HYDRUS-1D. The open-source compiler

255

Numba (Lam et al., 2015) was used to speed-up the Python code and to significantly reduce its

256

execution time. The Python version of the COSMIC code was verified against the original

257

FORTRAN code before further use.

258
259

Model Coupling Strategy

260

The COSMIC and HYDRUS models were coupled via a Python interface to dynamically261

calculate subsurface water flow by HYDRUS and neutron intensity by COSMIC. To this end,

262

several user-defined Python subroutines that interact with the HYDRUS files were developed. The

263

coupling of the two models is accomplished in the following way:264

1. Model Initialization: HYDRUS is initialized by setting up input data, the numerical265

domain, and boundary conditions. In addition, the temporal resolution of the model

266

output for the soil moisture profile is defined, which determines the temporal resolution

267

of the simulated neutron flux. Similarly, the COSMIC parameters are defined, and the268

Python code is compiled using Numba.

269

2. HYDRUS execution and data processing: HYDRUS-1D is externally executed directly270from the Python code, and a check is performed for numerical convergence and the final

271

mass balance error. If the simulation is successful, the calculated soil moisture profiles
11Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

272

are read and stored in a matrix, in which rows represent the vertical discretization of the273

model domain and columns represent the selected time steps for soil moisture profile274

information.

275

3. COSMIC execution: The matrix of soil moisture profiles is passed to COSMIC. Two

276

nested loops are used to iterate through depth and time and to calculate a time series of277

simulated aboveground neutron fluxes that can be compared with measured CRNP data.

278
279

Uncertainty Analysis

280

The main goal of the uncertainty analysis is to provide a statistical basis for the assessment281

of the information content of CRNP data for the inverse estimation of soil hydraulic parameters.

282

This is accomplished by first using a global optimization algorithm for the identification of the283

global optimum, which is subsequently used as a starting point for the Bayesian uncertainty analysis

284

based on a Markov Chain Monte Carlo algorithm to reduce its autocorrelation time, thus improving

285

the computational efficiency of the analysis. The global optimum is determined using the Particle

286

Swarm Optimization algorithm (Kennedy and Eberhart, 1995).

287
288

Bayesian Inference289

Global optimization algorithms are generally used to find an optimal set of parameters that290

maximize a likelihood function. However, different sources of uncertainty can affect hydrological291

data assimilation frameworks (i.e., model inadequacy, observations and measurements errors,292

improper boundary conditions, etc.). Therefore, there is no reason to assume that a single optimized

293

parameter set represents the only realistic model realization. It is more reasonable to consider it as294

one of many acceptable solutions of the problem. This idea underpins the equifinality concept first295

introduced by Beven and Binley (1992) in the early 90s, which emphasizes the fact that there are296

many acceptable (i.e., behavioral) representations that cannot be excluded from the analysis and297

that should be considered when quantifying the uncertainty associated with model predictions.
12Page 12 of 72

Page 13 of 72

298

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

Several methods have been proposed in the literature to assess the model predictive

299

uncertainty (e.g., Beven and Binley, 1992; Kavetski et al., 2006; Vrugt et al., 2008). However,300

during the last decade, the scientific community has focused on the application of statistically301

rigorous Bayesian methods for the calibration of hydrological models. To explain the Bayesian302

approach, let us consider a HYDRUS-COSMIC model realization with a symbol Y, which is a303function of q model parameters, u={u1, u2,…, uq}:

Y = f (u)

(6)

304

Bayes’ theorem combines the prior knowledge about the investigated process with the observed

305

data to obtain the posterior distribution of the parameters. In a probabilistic context, this can be

306

written as:~~~
Pr(Y | u) × Pr(u) = Pr(Y) × Pr(u | Y)
Likelihood× Prior = Evidence× Posterior
L(u) × π (u)du = Z × p(u)du

(7)

307

where Ỹ are the observed data, which modulate our prior belief π(u)du into the posterior p(u)du.

308

One of the main goals of Bayesian analysis in hydrological model calibration is the estimation of309

the parameters’ posterior distribution. This is usually accomplished by applying Markov Chain

310

Monte Carlo (MCMC) algorithms. One of the main problems faced in the Bayesian analysis of311

hydrological models is their strong nonlinearity and their discontinuous nature, which makes it very

312

difficult to apply efficient MCMC algorithms such as Hamiltonian Monte Carlo sampling

313(Betancourt, 2017). Although the classic Metropolis-Hastings algorithm (Metropolis et al., 1953) is314

easy to implement, it is unable to deal with highly correlated and multimodal posteriors often

315

encountered in the calibration of hydrological models. Thus, a variety of MCMC algorithms have316

been developed and successfully applied in water-related problems in recent years.

317

In this study, we use of the Affine Invariant Ensemble Sampler (AIES), which was first318

proposed by Goodman and Weare (2010). The AIES algorithm is an ensemble MCMC algorithm,

319

which has proven to be extremely effective when the dimensionality of the problem is not very high13

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

320(Huijser et al., 2015), as is the case in this study. Each ensemble, Ω, consists of D walkers, which

321

can be considered as a vector in the q-dimensional parameter space. The core of AIES is its affine

322

invariance property, which implies that the performance of the method is independent of the aspect323

ratio in highly skewed distributions. This is a highly desirable property when calibrating324hydrological models because of the prevalence of highly correlated posterior distributions. In AIES,325

the ensemble Markov Chain is evolved using a stretch move step (Goodman and Weare, 2010),

326

which significantly reduces the chain’s autocorrelation time. Briefly, the main idea is that the327

ensemble carries useful information about the approximated distribution, and thus each walker is328

updated using a complementary ensemble. A parameter a, which is usually set to 2.0 (Foreman-329

Mackey et al., 2012), is used to “stretch” the proposal. For a thorough description of the AIES,

330

please refer to Goodman and Weare (2010), Foreman-Mackey et al. (2012), and Huijser et al.

331(2015). It is worth noting that this study is the first one that uses the AIES algorithm in a332hydrological context.

333

One of the main advantages of AIES is that it requires the specification of only two334

parameters, i.e., the stretch factor, a, and the number of walkers, D. There are no clear indications in335

the literature about the best values of D and a, and they typically depend on the investigated336

problem. Based on some preliminary tests, D and a were set to 10q (q is the number of model

337

parameters) and 2.0, respectively. In this study, the Python package emcee (Foreman-Mackey et al.,3382012) is coupled with the HYDRUS-COSMIC model to carry out the Bayesian analysis. Emcee is a339

stable and well-tested Python implementation of the AIES method, which includes several MCMC

340

convergence diagnostic tools. The walkers are initialized by random sampling from a multivariate

341

normal distribution centered on the global optimum with a diagonal covariance matrix characterized342

by a small standard deviation of 0.01. Such an initialization around the global optimum is expected343

to reduce the autocorrelation time, thus increasing the overall computational efficiency of the344

analysis.

14

Page 14 of 72

Page 15 of 72

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

345

The AIES algorithm calls the coupled HYDRUS-COSMIC model at each iteration.

346

However, the stochastic combination of various soil hydraulic parameters can lead to non-347

convergent model runs for HYDRUS-1D, mainly due to the inability to appropriately adjust the348

HYDRUS-1D settings (e.g., mesh size) during the Monte Carlo analysis. In this study, non-349

convergent runs are identified and an unrealistically large value is attributed to the sum of squared

350

residuals. This procedure, which is similar to what was reported in Brunetti et al. (2016b) and351

Wöhling and Vrugt (2011), can significantly slow down the convergence of the AIES since it

352

introduces sudden and unrealistic jumps in the response surface due to the artificially low value of353

the likelihood. In future work, regularization of the response surface could be considered through a354

local interpolation or a surrogate model (e.g., Brunetti et al., 2017).

355

A non-informative uniform prior is used for all parameters in the uncertainty analysis. The

356

only information provided by this prior are the bounds of the feasible parameter space. Within these357bounds, all parameter values have the same probability. The bounds were set based on field358

measurements (e.g., particle size distribution) and literature values (see Table 3).

359
360

MCMC Diagnostic

361

Goodman and Weare (2010) suggested using the Integrated Autocorrelation Time, τint (IAT),

362

to assess the accuracy of the AIES. The IAT is the number of steps (i.e., ensemble moves) required

363

for the walker to produce an independent sample. This number is expected to be relatively large if364

the walkers are initialized in a distribution that is not equal to the stationary distribution. Such a

365transient period is usually defined as burn-in and should be discarded from the analysis. In the366

present study, the burn-in period is assumed to be 2τint. It is worth noting that we avoid the use of367

the widely applied Gelman-Rubin test since walkers in the AIES are not formally independent.

368

The IAT determines statistical errors in the Monte Carlo analysis once a stationary369

distribution has been attained (Sokal, 1996). In particular, the statistical error ε of the mean of the370

posterior depends on τint and is calculated as:15

Page 16 of 72

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

2τ int
Dm

ε=

(8)

371

where m is the number of steps (-). Obviously, a larger sample reduces the errors but significantly

372

increases the computational time. Thus, in this study, the admissible maximum statistical error was

373

set to 5%, which should guarantee a good trade-off between numerical accuracy and computational374

effort.

375

The estimation of the IAT is notoriously difficult since it is a two-point statistics. More376

specifically, the IAT should asymptotically reach its true value for a large number of samples. Thus,377

a preliminary IAT convergence analysis is performed to assess its calculated value. This is378

accomplished in the following way:379-

The AIES algorithm is preliminary run for 3,000 steps to obtain a significant statistical

380

sample. This requires a few days of computation time considering that the execution time of381

the HYDRUS-COSMIC model is less than 1 second.

382-

Each walker is divided into 20 parts, and the IATs are calculated for increasing chain

383

lengths. We adopt the numerical procedure described by Sokal (1996). In particular, a mean

384

normalized autocorrelation function is computed by averaging the normalized

385

autocorrelation functions calculated for each walker in the ensemble. Then, the procedure386

described by Sokal (1996) is used to estimate the IAT for each soil hydraulic parameter. The

387

IAT is expected to exhibit a logarithmic growth.

388-

The analysis is assumed to converge if:

 ε max < 5%max
m > 10τ int

(9)

389

Otherwise more samples are generated by evolving the AIES. The second condition in eq. 9390

is meant to reduce the MCMC error and mitigate the effect of an imperfect estimation of the

391

IAT.

392
16

Page 17 of 72

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

393

Likelihood Function394

If we assume that the error residuals are uncorrelated and normally distributed with constant395

variance, σ2, the likelihood function L(u) can be written as:k

L(u) = ∏
i =1

 1
2
exp − σ 2 ( yi (u ) − ~yi ) 
 2

2πσ
1

(10)

2

396

where yi(u) and ỹi(u) are the ith model realization and its corresponding measured value,397

respectively. For algebraic simplicity, the log-likelihood ℓ(u) is used:

( )

1 k
k
k
2
l(u) = − ln(2π ) − ln σ 2 − 2 ∑ ( yi (u) − ~yi )222σ i=1(11)398

Eq. 11 is used in the inverse modeling with synthetic data since the error variance is known. In the399

inverse modeling with measured data, σ2 is treated as a nuisance parameter and integrated out from400

the inference equation (e.g., Kavetski et al., 2006). The new likelihood is then written as:k k2l(u) = − ln ∑ ( yi (u ) − ~yi ) 
2  i=1


(12)

401

Rollesbroich grassland experimental catchment

402

The Rollesbroich grassland experimental catchment is located in Western Germany in the403

Eifel mountain range and is part of the TERENO observatory Eifel/Lower Rhine Valley (Bogena et404al., 2012, 2018). It covers an area of about 40 ha with altitudes ranging from 474 to 518 m a.s.l. (the405average slope is 1.63°). The mean annual air temperature and precipitation are 7.7 °C and 1033 mm,

406

respectively. The vegetation is dominated by perennial ryegrass (Lolium perenne) and smooth407

meadow grass (Poa pratensis) and the predominant soils are Cambisols in the southern part and408

Stagnosols in the northern part of the catchment. All components of the water balance (e.g.,409

precipitation, evapotranspiration, runoff) are continuously monitored using state-of-the-art

410

instrumentation. A detailed description of the instrumentation at the Rollesbroich experimental

411

catchment and available hydrological datasets is given in Qu et al. (2016). A SoilNet wireless17

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

412

network (Bogena et al., 2010) with more than 500 soil moisture sensors (ring oscillator probes,

413

Bogena et al., 2017) was installed in May 2011 to measure the soil water content at depths of 5, 20414

and 50 cm (Qu et al., 2014) at 87 locations. In particular, two sensors were installed in parallel at415

each depth with a distance of ⁓10 cm to increase the sensing volume and to allow examination of416

inconsistencies in sensor reading. In May 2011 and May 2012, two CRNPs (type CRS1000,

417

HydroInnova LLC) were installed in the southern and northern parts of the catchment. The CRNP

418

calibration and a comparison with in-situ soil water content data obtained with SoilNet are419

presented in Baatz et al. (2014). This study uses data only from the CRNP located in the southern420part of the catchment.

421

A six-month long dataset (1 April 2012 to 30 September 2012) is selected for further

422

analysis (Figure 1). In particular, daily precipitation and reference evapotranspiration are used as

423

input data in HYDRUS for both synthetic and experimental scenarios, while daily measured water424

contents in three depths (i.e., z=-5, -20, and -50 cm) and the aboveground neutron fluxes are used in425

the inverse estimation of the soil hydraulic parameters in the experimental scenarios.

426

18

Page 18 of 72

Page 19 of 72

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

427
428
429
430

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

431
432

Numerical Domain and Boundary Conditions

433

The 200 cm deep soil profile was discretized in HYDRUS-1D into 100 finite elements,

434

refined at the top to accommodate pressure head gradients induced by the atmospheric conditions.

435

An atmospheric boundary condition is applied at the soil surface using: (a) precipitation and436

potential evaporation fluxes, (b) a prescribed zero pressure head (i.e., full saturation) during437

ponding, and (c) equilibrium between the soil surface pressure head and the atmospheric water

19

Page 20 of 72

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

438

vapor pressure when the atmospheric evaporative demand cannot be met. The effect of groundwater439

is not simulated in this study and thus a ‘Free Drainage’ BC is used at the bottom of the model

440

domain (z=-200 cm).

441

The initial pressure head was assumed to be constant and equal to -100 cm in the entire

442

model domain. A three-month long spin-up period was used to limit the influence of the initial443

condition on the simulation results. Three observation points were set at depths of -5, -20, and -50

444

cm to retrieve simulated volumetric water contents.

445
446

Inverse estimation with synthetic data

447

Several studies have demonstrated how data assimilation frameworks can be significantly

448

affected by various sources of uncertainty at the field scale that can lead to biased conclusions if not449

properly accounted for in the analysis (Mantovan and Todini, 2006; Mockler et al., 2016;

450

Montanari, 2007). Since the main aim of this study is to provide an objective evaluation of the

451

information content of the CNRP data for the inverse estimation of soil hydraulic parameters at the452

field scale, it is important to first analyze synthetic modeling scenarios that are not affected by453

uncertainty in the model structure and measurement errors before analyzing actual experimental

454

data. Here, the synthetic modeling scenario is a simplification of the experimental case study and is455

focused on soil-neutron interactions in homogeneous sandy soil. The strongly unsaturated

456

conditions induced by the use of sandy soil will increase the aboveground neutron intensity and thus

457

provides different information compared to the experimental case study characterized by a silty soil.

458
459
460
461
462

Model Set-up

Table 1. The COSMIC parameters used in the synthetic and experimental scenarios (Baatz et al.,2014).
Scenarioρbd (g/cm3)

θl (-)NCOSMIC (-)

20

L1 (-)

L2 (-)

L4 (-)

Page 21 of 72

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

Synthetic

1.42

0.037

189

161.98 129.14

3.16

Experimental

1.09

0.068

213

161.98 129.14

3.16

463
464

Variably-saturated water flow is simulated in a 200 cm deep homogeneous loamy sand soil

465

profile covered by grass. The VGM parameters reported in Carsel and Parrish (1988) for loamy

466

sand are used to generate the synthetic dataset. The scaling factors are fixed to 1.0 to obtain a

467

homogeneous soil profile. The rooting depth and density are assumed to be 15 cm and 1.0,468

respectively, while a constant LAI equal to 2.0 (Qu et al., 2016) is used to partition reference469

evapotranspiration. The Feddes parameters for grass are chosen according to Taylor and Ashcroft

470(1972) and Wesseling et al. (1991). The COSMIC parameters are set according to Baatz et al.

471(2014) (472

Table 1).

473

The synthetic dataset is generated in the following way. First, a forward numerical474

simulation by the HYDRUS-1D model is used to generate a time series of the volumetric water475

content at z = -5 cm and daily soil moisture profiles. The latter are passed to the COSMIC model for476

the simulation of the aboveground neutron fluxes. Next, synthetic noise is added to the original time477

series. The selection of σ values has a significant influence on parameter uncertainty and thus they478

should be properly set to reflect the accuracy of real measured data. Indeed, an unbalanced choice479

of such values can lead to biased conclusions on the information content of the investigated data.

480

These values are thus set in this study based on a preliminary analysis of the noise in selected481

measured datasets. Briefly, a spectral analysis based on the Fast Fourier Transform was used to482

separate the deterministic (signal) and stochastic (noise) components of the measured neutron fluxes

483

and volumetric water contents. The Shapiro-Wilk test was then used to assess the normality of the484

noise, and the standard deviations were calculated. Based on this analysis, random noise with zero485

mean and standard deviations equal to 30.0 N/hour and 0.01 cm3/cm3 for the neutron flux and water

21

Page 22 of 72

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

486

content time series, respectively. The input, original, and synthetic data for the selected time period487

are shown in Figure 2.

488

489
490
491
492
493

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

494
495
496

Synthetic Modeling Scenarios
Table 2. Synthetic modeling and experimental data analysis scenarios.
Scenario ID
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

22

Page 23 of 72

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

Near-surface water content and neutron flux time series

S3

-

497
498

Three synthetic modeling scenarios are considered for the inverse estimation of soil

499

hydraulic parameters (Table 2). Scenarios S1 and S2 examine the information content of the500

temporal evolution of aboveground neutron flux and water content in one observation point,

501

respectively. The comparison of the parameter uncertainty obtained by these two scenarios will502

clarify the potential of CRNP data to estimate SHP compared to traditional point measurements.

503

Scenario S3, which combines the two synthetic datasets, is used to investigate the effect of the joint

504

inversion of water content and neutron flux data. The parameter bounds used in this analysis are505

reported in Table 3.

506
507
508

Table 3. The bounds of the VGM parameters and their true values used in the synthetic andexperimental scenarios.
Synthetic

Experimental
True Value

Parameter

Lower

Upper

Lower

Upper

Bound

Bound

(Carsel and Parrish,
Bound

Bound
1988)

θr (cm3/cm3)

0.03

0.08

0.06

0.07

0.07

θs (cm3/cm3)

0.35

0.45

0.41

0.30

0.45

log10α (1/cm)

-1.50

-0.50

-0.91

-2.50

-0.85

n (-)

1.50

3.00

2.28

1.10

2.00

1.00

3.00

2.54

0.50

3.00

φh1 (-)

1.00

1.00

1.00

1.00

5.00

φh2 (-)

1.00

1.00

1.00

1.00

5.00

φk1 (-)

1.00

1.00

1.00

1.00

20.00

φk2 (-)

1.00

1.00

1.00

1.00

5.00

log10Ks(cm/day)23

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

509
510

Inverse estimation with actual data

511

Model Set-up

512

Measured particle size distributions (PSD) reported in Qu et al. (2016) for the Rollesbroich

513

catchment revealed negligible differences in the textural composition of the soil at different depths.

514

On the other hand, bulk densities increased with depth from 0.94 to 1.52 g/cm3, and porosities

515

decreased with depth from 0.65 and 0.43 (cm3/cm3). This vertical soil variability is attributed to516

prior agriculture land management (i.e., the plow layer is 20 to 30 cm thick), which reduced the517

compaction of the surface layer. The measured soil properties highlight the existence of two

518

different soil horizons in the top layer (i.e., z>-30 cm). Indeed, the porosity and the bulk density for519

z>-10 cm are 20% higher and 25% lower than for the underlining topsoil, respectively. Thus, three520

soil materials are assumed in the model: soil material 1 extends until z=-10 cm; soil material 2521

ranges from -10 to -30 cm; and soil material 3 extends below -30 cm. Since material 3 is not522

affected by agriculture land management, the scaling factors are set to 1.0 for this material.

523

The pedotransfer model ROSETTA (Schaap et al., 2001) is used to identify the plausible524

ranges of the soil hydraulic parameters and the scaling factors. On the basis of this analysis, the525water content scaling factors φθ are fixed to 1.4 and 1.1 for layers 1 and 2, respectively, thus526

reducing the dimensionality of the inverse problem. Since the residual water content exhibited low

527

variability, it was fixed to 0.07, which is a common value for silty soils. The tortuosity and pore-528

connectivity parameter L [-] in the VGM equation is set to 0.5 in all numerical simulations. The

529

aforementioned model simplifications reduced the number of unknown parameters to 8. The

530

parameter bounds used in the analysis are reported in Table 3.

531
532

Experimental Modeling Scenarios

533

Three experimental data analysis scenarios are considered (Table 2). Scenarios E1 and E2,534

in analogy to scenarios S1 and S2, consider only measured neutron flux and near-surface water
24

Page 24 of 72

Page 25 of 72

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

535

content (i.e., z = -5 cm), respectively. The comparison of these two scenarios will clarify the536

information content of the CRNP data to estimate SHP in layered soils relative to near-surface point537

measurements. In addition, Scenario E3 uses measured water contents at three different depths (i.e.,538

z = -5, -20, and -50 cm). This scenario should provide more accurate estimates of SHP because of539

the higher amount of measurements in multiple soil layers, and thus serves as a benchmark for the540

other two scenarios.

541
542

RESULTS AND DISCUSSION

543

Inverse estimation with synthetic data

544

MCMC Diagnostic

545
546
547
548

Figure 3. The diagnostic plot of the AIES algorithm for synthetic modeling scenarios (scenarios S1,
S2, and S3 in A, B, and C, respectively) showing the integrated autocorrelation time against thealgorithm step for each soil hydraulic parameter.

549
550

The development of the IAT in the MCMC analyses for the three synthetic scenarios is551

reported in Figure 3. The asymptotic behavior of the IAT with increasing algorithmic steps suggests

552

that a reliable estimation has been obtained for all three scenarios. It can be seen that all three553

scenarios exhibit similar autocorrelation times. In particular, the dynamics of the IAT indicate that554

approximately 100 ensemble moves (i.e., steps) are sufficient for the AIES to forget about its initial

555

position and to start sampling from a stable posterior distribution. The saturated water content
25Page 26 of 72

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

556

exhibits the largest IAT in scenarios S1 and S2 (subplots A and B in Figure 3), while in Scenario S3

557(subplot C in Figure 3) it is slightly larger for the shape parameter n. The short autocorrelation

558

times lead to a maximum statistical error of the Monte Carlo analyses (eq. 8) of approximately 3.5

559% after 3,000 steps. This is well below the prescribed tolerance of 5%, thus indicating a good560

accuracy in the estimation of the posterior distributions.

561
562
563
564

Uncertainty analysis
Table 4. Original values and 95% Bayesian Credible Intervals of the soil hydraulic parameters forthe three synthetic scenarios. The 2.5 and 97.5% quantiles are the bounds of the credible intervals.
Scenario
True Values
Parameter

S1

S2

S3
(Carsel and Parrish, 1988)

2.5%

97.5%

2.5%

97.5% 2.5% 97.5%

θr (-)

0.05

0.06

0.05

0.06

0.05

0.06

0.057

θs (-)

0.35

0.44

0.35

0.44

0.35

0.44

0.41

log10α (1/cm)

-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/day)

1.92

2.78

2.06

2.93

2.00

2.86

2.54

565
566

The joint and marginal posterior distributions for each soil hydraulic parameter for Scenario

567

S1 (above-ground neutron fluxes) are shown in Figure 4. The inverse problem is markedly

568

unimodal with strong parameter interactions. More specifically, the joint posterior distribution of θs

569

and Ks reveals an expected positive correlation, while the VGM shape parameters are negatively

570

correlated. The interdependence of α and n is exacerbated, as indicated by their slight banana-571

shaped joint posterior distribution. Interestingly, the posterior for Ks is unimodal, even though it

572

spans a relatively broad range of values from 83 and 602 cm/day. On the other hand, the marginal

573

posterior distribution for θs is flat, indicating that the neutron data do not contain enough26

Page 27 of 72

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

574

information content to update the uniform prior distribution. This is caused mainly by the fact that575

soil water content is relatively low (between 0.05 and 0.20 m3/m-3) in this synthetic data due to the576

use of a highly permeable loamy sand. Therefore, water content never approaches the saturated

577

water content of 0.41 and little information is available to constrain θs. Different results would

578

likely have been obtained for other (finer) soil textures. In this case, independent measurements of579

soil porosity would be needed to better constrain θs. The strong interaction effect shown in the joint

580

posterior θs-Ks suggests the impossibility to fix the saturated water content independently from the581

saturated hydraulic conductivity. As expected in strongly unsaturated conditions, θr is precisely

582

estimated. Similarly to θs, different types of Ks measurements are needed to reduce the uncertainty.

583

Similar findings were reported in Gibson and Franz (2018), who used a mobile CRNP to quantify

584

spatial patterns of SHPs across a range of spatial scales. In that study, the authors underlined how

585

inability to map near-saturated soil water contents increased the uncertainty in the estimation of Ks.

586

Overall, the marginal posterior distributions indicate a well-posed inverse problem

587

characterized by acceptable parameter identifiability and moderate uncertainty. This is confirmed

588

by the Bayesian credible intervals reported in Table 4. In particular, the shape parameters α and n

589

are generally well predicted, with the latter exhibiting a slightly narrower posterior distribution.

27

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

590
591
592
593
594
595
596

Figure 4. Joint (below diagonal) and marginal (diagonal) posterior distributions of soil hydraulicparameters 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.

597

As shown in Figure 5, the use of near-surface water content instead of neutron flux

598

(Scenario S2) leads to similar posterior distributions. Again, the inverse problem is mostly599

unimodal and affected by significant parameter correlations that are similar as in S1 (i.e., a positive600

correlation between θs and Ks, and a negative correlation between of α and n). The residual water
28

Page 28 of 72

Page 29 of 72

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

601

content θr is generally uncorrelated, except for a slight interdependence with n. As in Scenario S1,602

the near-surface water content measurements do not provide enough information to estimate θs. The

603

marginal posterior distribution for Ks again is unimodal, and spans values between 100 and 1,000604

cm/day, thus indicating a somewhat higher uncertainty compared to Scenario S1. This effect could

605

be related to the increased information content of the neutron flux, which is representative of the606

soil moisture dynamics in the entire top soil and not just for a specific depth (i.e., 5 cm in this607

study).

608

Overall, a comparison of Figure 4 and Figure 5 suggests a minor decrease in parameter

609

uncertainty for Scenario S2 compared to Scenario S1 (Table 4). This is particularly evident for the610

VGM shape parameters, which have narrower Bayesian credible intervals and a good agreement611

between the most probable and true values. However, it must be noted that these differences depend612

on the selected level of synthetic noise used to generate the synthetic datasets. In this view, the613

sensor accuracy, as well as different sources of noise, play a fundamental role in the analysis of614

experimental data. Similar to Scenario S1, the results suggest the need to use more informative prior

615

distributions and to combine different data types to reduce the overall parameter uncertainty. This is616

investigated in Scenario S3, where both aboveground neutron fluxes and water contents are used.

29

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

617
618
619
620
621
622

Figure 5. Joint (below diagonal) and marginal (diagonal) posterior distributions of the soil hydraulicparameters 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.

623
624

The MCMC analysis for Scenario S3 (with both near-surface water contents and above

625

ground neutron fluxes) resulted in strongly unimodal and narrow posterior distributions (Figure 6).

626

The correlations between parameters are similar to Scenarios S1 and S2, and the relatively high627

uncertainty in the estimation of θs remains. The θr is precisely estimated, while the uncertainty for
30

Page 30 of 72

Page 31 of 72

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

628

Ks persists. The overall reduction in the parameter uncertainty for Scenario S3 is visually evident by

629

comparing the joint posterior distributions shown in Figure 4, Figure 5, and Figure 6, and is further

630

confirmed by the comparison of the Bayesian credible intervals (Table 4). Indeed, the VGM shape631

parameters have narrower credible intervals, suggesting that the simultaneous use of more and632

different data types reduces the uncertainty in the estimated parameters. The results of Scenario S3

633

further confirm that the saturated water content and hydraulic conductivity cannot be precisely

634

inferred from data representative of strongly unsaturated conditions and that independent635

measurements of soil porosity and hydraulic conductivity may be needed in such a case.

636

As a result of this analysis, a number of conclusions can be drawn:

637-

The three synthetic datasets provided sufficient information content to obtain a638

satisfactory estimate of the SHP. More specifically, the residual water content, the VGM

639

shape parameters, and the saturated hydraulic conductivity exhibited unimodal posterior

640

distributions, with the latter characterized by a higher uncertainty. The saturated water641

content could not be inferred from the dataset (for loamy sand), suggesting that642

independent measurements of the soil porosity are needed to better constrain the inverse

643

problem for this soil texture.

644-

The use of the aboveground neutron flux data led to a satisfactory estimate of SHPs. In

645

particular, the use of neutron flux data slightly reduced the uncertainty in the estimation

646

of Ks. This provides the first proof of concept for the use of CRNP data for the inverse

647

estimation of effective SHPs at the field scale.

648-

The simultaneous use of water content and neutron flux data reduced the parameter

649

uncertainty. This suggests that the combination of CRNP data with sparse point

650

measurements in the near-surface environment could reduce parameter uncertainty, in651

particular for the VGM shape parameters. However, the number of point measurements

652

should be sufficient to characterize the soil heterogeneity within the CRNP footprint.

653
31

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

654
655
656
657
658
659
660

Figure 6. Joint (below diagonal) and marginal (diagonal) posterior distributions of the soil hydraulicparameters for Scenario S3. 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.

661

32

Page 32 of 72

Page 33 of 72

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

662

Inverse estimation with actual data

663

MCMC Diagnostic

664
665
666
667

Figure 7. The diagnostic plot of the AIES algorithm for the experimental data analysis scenarios(scenarios E1, E2, and E3 in A, B, and C, respectively) showing the integrated autocorrelation timeagainst the algorithm step for each soil hydraulic parameter.

668
669

The IAT for the MCMC analyses for the three experimental scenarios is shown in Figure 7.

670

The maximum statistical error ranges between 3 and 5%, which is under the prescribed tolerance671

and indicates a reliable estimation of the posterior distributions. As expected, the increased672

dimensionality of the inverse problem resulted in longer autocorrelation times and more steps

673

before convergence to the true posterior distribution. This behavior agrees with the properties of the

674

AIES algorithm in high-dimensions reported in Huijser et al. (2015).

675

The maximum IAT for Scenario E2 (subplot B in Figure 7) is approximately three times

676

larger than the IAT obtained for scenarios E1 and E3, thus indicating convergence issues, This is677

related to multimodal posterior distributions that will be discussed later. Poor convergence is678

particularly evident for the scaling factor φk1, which is characterized by an IAT of 1,600 steps.

679

Under such circumstances, the ensemble samplers are known to perform poorly, mainly due to their680

inability to traverse low-probability valleys between peaks (Matthews et al., 2017). It is worth

681

noting that the performance of the AIES for Scenario E2 was insensitive to the value of the stretch

682

factor a, which was then set to the suggested value of 2.0. Future work should consider using
33

Page 34 of 72

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

683

multimodal Nested Sampling techniques (Feroz et al., 2009), which were found to be effective for684

such types of inverse problems. On the other hand, the relatively large autocorrelation time for the685

shape parameter n in Scenario E1 (subplot A in Figure 7) is partially related to the position of the686high likelihood region, which is close to the boundary of the numerical domain. Indeed, bounded

687

prior distributions can lead to convergence issues when the posterior develops near the boundary.

688
689
690
691

Uncertainty analysis
Table 5. 95% Bayesian Credible Intervals of the soil hydraulic parameters for the threeexperimental scenarios. The 2.5 and 97.5% quantiles are the bounds of the credible intervals.
Scenario
E1

Parameter

E2

E3

2.5%

97.5%

2.5%

97.5%

2.5%

97.5%

θs (-)

0.33

0.38

0.35

0.39

0.36

0.37

log10α (1/cm)

-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/day)

0.50

2.82

0.50

1.44

0.97

1.18

φh1 (-)

1.00

4.32

1.00

2.33

1.00

1.05

φh2 (-)

1.27

5.00

1.00

4.35

1.80

2.21

φk1 (-)

1.00

18.94

1.00

11.51

1.07

2.63

φk2 (-)

1.00

4.78

1.00

4.65

1.00

1.60

692
693

Since Scenario E3 serves as a benchmark for the other two scenarios, this will be discussed694

first. The joint and marginal posterior distributions of the SHP for this scenario are shown in Figure695

8 and have a markedly unimodal pattern characterized by a leptokurtic behavior. Consequently, the696

Bayesian credible intervals (Table 5) for the soil hydraulic parameters are extremely narrow, thus

697

indicating their good identifiability. The Bayesian analysis reveals, as expected, a well-posed

698

inverse problem characterized by low uncertainty mainly stemming from the simultaneous use of34

Page 35 of 72

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

699

measured water contents at different depths, which better constrain the inverse estimation problem.

700

Overall, the results agree well with those reported in Qu et al. (2016), where the ROSETTA

701

pedotransfer function was used to predict the spatial distribution of soil hydraulic parameters using702

the soil textural information for the top 5-cm of soil. Indeed, the saturated water content for the soil

703

material 1 ranges between 0.50 and 0.52 cm3/cm3, which is in line with the values of the porosity704

reported in Qu et al. (2016) ranging between 0.46 and 0.56 cm3/cm3 for the southern part of the705

Rollesbroich catchment. Results slightly deviate for the VGM shape parameters α and n, which are

706

estimated to be higher and lower, respectively, compared to Qu et al. (2016), who reported α and n

707

values ranging between -2.40 and -2.30 1/cm and 1.63 and 2.23, respectively. Higher uncertainty is708

found for the saturated hydraulic conductivity of soil material 1, which ranges between 10 and 40

709

cm/day. It is worth noting that SHP estimates using pedotransfer functions are usually associated710

with significant uncertainty. Thus, it is preferable to use inversely estimated soil hydraulic711

parameters when the model inversion is well-posed and the associated uncertainty is properly712

assessed.

35

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

713
714
715
716
717
718

Figure 8. Joint (below diagonal) and marginal (diagonal) posterior distributions of the soil hydraulicparameters 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.

719

The effect of the parameter uncertainty on the simulated water contents and neutron fluxes is

720

shown in Figure 9, which reports the model predictive uncertainty (grey lines) obtained by random721

sampling of 500 solutions from the posterior distribution. The model can accurately describe the722soil moisture dynamics in the three observation depths with limited uncertainty. The quality of the723

fitting is significantly better than that reported in Qu et al. (2014), who inversely estimated the724

SHPs from soil water contents measured with a sensor network at 41 locations in the Rollesbroich

725

catchment. This error reduction confirms that the use of the scaling factors is a valid approach for726

this dataset since it results in a better description of water flow compared to the assumption of727

homogeneous soil used in Qu et al. (2014). On the other hand, the comparison between simulated

728

and measured aboveground neutron fluxes (Figure 9A) reveals an appreciable deviation. More

36

Page 36 of 72

Page 37 of 72

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

729

specifically, the coupled model tends to overestimate the variability of neutron fluxes. This is730

particularly evident during August and September 2012. The Rollesbroich site is managed grassland

731

with a relatively low biomass. Except for an occasional snow cover during the winter period, there732

are no additional transient hydrogen pools that have a strong influence on the measured neutron733

signal. As a consequence, such behavior can be partially explained by the simplifying assumptions

734

of the COSMIC model, by an imperfect neutron flux correction (e.g., for variations due to735

fluctuations in the atmospheric density and humidity or incoming neutrons), and by an inadequate736

calibration of the COSMIC model parameters. In such circumstances, two approaches can be737

recommended.

738

1. Two-Step Sequential Calibration: First, the HYDRUS model is preliminarily calibrated

739

against measured data from different locations and depths in the measurement footprint740

of the cosmic-ray probe. Next, the calibrated HYDRUS model is coupled with the741

COSMIC model for the inverse estimation of the COSMIC parameters from the742

measured neutron fluxes. It must be emphasized that the possibility to use simulated

743

high-resolution soil moisture profiles instead of few point measurements at different soil

744

depths can increase the accuracy of the COSMIC model predictions. In this perspective,745

it may be worth coupling the three-dimensional (3D) model (e.g., HYDRUS (2D/3D),

746with the COSMIC model in order to account for the horizontal and vertical soil

747

heterogeneity at the field scale (≈300 m). The recent work of Gibson and Franz (2018)748

suggests that sparse point measurements of the SHPs can significantly improve the749

description of spatial patterns of SHPs.

750

2. Joint Calibration: The soil hydraulic and COSMIC parameters are jointly estimated751

using the coupled HYDRUS-COSMIC model using both measured neutron fluxes and752

volumetric water contents. However, this will increase the dimensionality of the inverse

753

problem, which can lead to convergence issues in the MCMC analysis.

754
37

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

755
756
757
758
759

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

760
761

The joint and marginal posterior distributions of the soil hydraulic parameters for Scenario

762

E1 that relies on measured aboveground neutron flux are shown in Figure 10. It is evident that the763

measured neutron fluxes do not provide sufficient information to update the prior distributions of764

φh2, φk1, and φk2, which remain almost flat after the MCMC analysis. Indeed, their Bayesian credible

765

intervals (Table 5) are only marginally narrower than the prior bounds listed in Table 3. A similar

766

but less pronounced behavior is visible for Ks and φh1, although the marginal posterior distributions

767

are right-skewed, thus indicating a slight influence of the measured data on these parameters. It is768

worth noting that the MCMC analysis reveals a highly skewed bivariate posterior distribution for769

the pair of parameters α-Ks. In such circumstances, the affine invariant properties of the AIES

770

algorithm significantly mitigate convergence issues usually encountered in traditional MCMC

771

algorithms.

38

Page 38 of 72

Page 39 of 72

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

772

Interestingly, the use of measured neutron fluxes leads to well-constrained estimates of773

VGM shape parameters, thus confirming the findings of the synthetic modeling scenarios. The

774

confident estimation of θs stems from the use of near-saturated water contents, which better

775

constrain the inverse problem. This effect is more pronounced for n and θs, while α presents a more

776

platykurtic distribution. In particular, the Bayesian credible intervals for θs and n range between

777

0.33 and 0.38 cm3/cm3 and 1.10 and 1.47, respectively, which are slightly larger than those obtained778

in Scenario E3. On the other hand, α ranges between 0.004 and 0.04 1/cm, thus indicating a779

tendency to overestimate soil water retention capacity. Indeed, very low values of α are usually780

associated with a high retention capacity and delayed desaturation of the porous medium. In such

781

circumstances, simulated water content fluctuations are limited and propagated in the calculation of782

the fast neutron fluxes. Thus, from a numerical point of view, the Bayesian inference process783

adjusted the soil hydraulic parameters to compensate for other plausible sources of uncertainty (e.g.,

784

COSMIC model assumptions, COSMIC parameters). Therefore, the results for Scenario E1 further

785

confirm the findings for Scenario E3 and suggest the use of a joint Bayesian calibration of the786

HYDRUS-COSMIC model.

787
788

39

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

789
790
791
792
793

Figure 10. Joint (below diagonal) and marginal (diagonal) posterior distributions of the soil
hydraulic parameters 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 correspondto the parameter bounds reported in Table 3.

794
795

This aspect is further investigated in Figure 11 which shows a comparison between

796

measured (black circles in Figure 11) and modeled depth-weighted soil water contents. In

797

particular, measured data are compared with soil water contents calculated from the CRNP using798

the standard calibration method (Baatz et al., 2014) (a grey line in Figure 11), and predicted using799

the model with the median soil hydraulic parameters from the posterior distribution for scenario E1

800(a black line in Figure 11), respectively. Depth-weighted soil water contents are obtained using

801

weights of 0.852, 0.146 and 0.002 at z=-5, -20 and -50 cm, respectively, as reported in Baatz et al.

802(2014). Overall, both methods provide a good description of the near-surface soil water contents as

803

indicated by low RMSE values. An appreciable overestimation of the water content is observed for804

both methods between August and September 2012, with the magnitude being higher for the40

Page 40 of 72

Page 41 of 72

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

805

COSMIC predicted water contents (a black line in Figure 11), which generally exhibit a lower

806

variability, compared to water contents calculated using the standard calibration method. Such807

difference confirms potential approximations in the COSMIC modeling assumptions and an

808

imperfect calibration of the COSMIC parameters, although the deviation from measured data is809

limited.

810
811
812
813
814
815

Figure 11. A comparison between measured (black circles) and modeled (solid lines) depthweighted soil water contents. The grey and black lines indicate soil water contents calculated fromthe CRNP using the standard calibration method (Baatz et al., 2014), and predicted using the modelwith the median soil hydraulic parameters from the posterior distribution from scenario E1,
respectively.

816
817

Figure 12 shows the model predictive uncertainty (grey lines) in the simulated neutron

818

fluxes (subplot A) and volumetric water contents at three different depths. It is evident that both the819

discrepancy between measurement and model and the uncertainty grow with depth. The model

820

predictive uncertainty only partially encompasses measured data, and the water contents are821

considerably underestimated at z=-50 cm. This behavior is expected since the measured neutron

822

fluxes are mostly informative of near-surface water dynamics. Furthermore, the results indicate that823

the variability of water content is underestimated. This is particularly visible between 2012-06 and824

2012-09 at z=-5 and -20 cm, and it can be partially explained by the previously described

825

underestimation of the shape parameter α.

41

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

826

Overall, this analysis suggests that CRNP data can provide reasonable estimates of the827

effective soil hydraulic parameters at the field scale, which can be extremely valuable in land828

surface hydrological modeling. In particular, the VGM shape parameters and the saturated water829

content are successfully inferred from the measured neutron flux, thus leading to a coherent

830

simulation of the soil moisture dynamics in the near surface, which is crucial for a proper

831

description of the surface energy balance in LSMs (Iwema et al., 2017). The large uncertainty

832

observed for Ks can significantly affect simulations of overland flow and the partitioning of

833

precipitation into infiltration and surface runoff, and thus other types of measurements are needed to834

better constrain this parameter. Nevertheless, the large footprint of CRNPs represents an important

835

advantage against classical near-surface point measurements, the information content of which can

836

be affected by various sources of uncertainty. This aspect is better investigated in the final

837

experimental data analysis scenario.

838
839
840
841
842

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

843
42

Page 42 of 72

Page 43 of 72

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

844

Figure 13 shows the joint and marginal posterior distributions for Scenario E2 that relies on

845

near-surface water content measurements only. A multi-modal posterior distribution with two peaks

846

is observed for both α and Ks. The examination of the joint posterior distribution of the pair α-Ks847

reveals the existence of two regions of high probability, one centered around the location [log10α=-848

1.7, log10Ks=0.5] and another one around [log10α=-1.2, log10Ks=1.3]. While the first region of high849

probability is in line with the values obtained for Scenarios E1 and E3 and thus likely reflects the850

true SHP, the latter first region of high probability significantly overestimates the shape parameter851

α. Such predictions could be representative of another plausible pore system induced by the plowing852

of the top layer. These sources of uncertainty propagate in the credible interval (Table 5), which853

now ranges between 0.01 and 0.11 1/cm, thus indicating poor identifiability of the α parameter. In

854

this view, it must be emphasized that aboveground neutron fluxes are less influenced by local855

sources of uncertainty compared to near-surface point measurements due to their footprint and856

penetration depth.

857

On the other hand, θs and n are well constrained and in line with the values obtained in the858previous two scenarios. A significant uncertainty is noticeable for the scaling factor φk1, which859

ranges between 1.0 and 11.51, despite being characterized by an apparently strong unimodal

860

posterior distribution. In particular, a region of high-probability is evident around 1.0, followed by a861

flat valley of low probability. Such posterior distributions are mainly responsible for the large

862

integrated autocorrelation time observed in Figure 7. Nevertheless, a comparison of the posterior

863

distributions of the scaling factors of soil materials 1 and 2 confirms that the measured water864

contents at z=-5 cm are mainly informative of the soil hydraulic properties of the surface layer.

865

Indeed, the posterior distributions of φh2 and φk2 are very similar to the prior distributions, whereas866

the posterior distributions clearly are much smaller for φh1 and φk1.

867

43

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

868
869
870
871
872

Figure 13. Joint (below diagonal) and marginal (diagonal) posterior distributions of the soil
hydraulic parameters 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 correspondto the parameter bounds reported in Table 3.

873
874

The propagation of uncertainty in the model predictions for Scenario E2 is shown in Figure875

14. It is again evident that the underestimation of the soil retention capacity, caused by the low

876

inferred values of the shape parameter α, results in lower values of the water contents at z=-20 and -877

50 cm. The reduced retention capacity propagates in the simulated aboveground neutron fluxes,

878which are significantly overestimated. However, the magnitude of the uncertainty is smaller than879

that obtained for the water contents at deeper soil horizons, which is expected since neutron fluxes

880

are mostly influenced by the soil moisture dynamics in the near-surface.

881

As a result of these analyses, a number of conclusions can be drawn:

882

-

883

The results for Scenario E3 indicate that it is not possible to simultaneously fit themeasured neutron flux and soil water content data. This may be related to uncorrected
44

Page 44 of 72

Page 45 of 72

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

884

effects of aboveground biomass on the measured neutron data, or inadequate modeling885

assumptions. In particular, the deviation between measured neutron fluxes and those886

simulated by the calibrated HYDRUS model suggest potential inaccuracies in the887

COSMIC parameters. This could potentially be addressed with a two-step sequential or a888

joint calibration of the coupled HYDRUS-COSMIC model.

889-

The use of measured aboveground neutron fluxes led to a reasonably confident

890

estimation of the VGM shape parameters and the saturated water content. From this891

point of view, the possibility of inversely estimating the effective soil hydraulic892

parameters at the field scale from the CRNP data is confirmed. The combination of893

neutron data with other independent types of measurements (e.g., soil porosity, saturated

894

conductivity) can further reduce the parameter uncertainty.

895-

A comparison of inversion results obtained for near-surface volumetric water contents or

896

neutron fluxes revealed that the latter are less influenced by typical sources of897

uncertainty associated with land surface management activities (e.g., soil tillage), thus

898

leading to more reliable and representative estimates of soil hydraulic parameters.

899
900

45

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

901
902
903
904
905

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

906
907
908

CONCLUSIONS AND SUMMARY
The main goal of this study was to assess the information content of cosmic-ray neutron data909

for the inverse estimation of the effective soil hydraulic parameters. An extensive theoretical and910

experimental assessment based on the combination of the coupled HYDRUS-COSMIC model with911

a Bayesian data analysis framework revealed that neutron fluxes could be used effectively to912

estimate the soil hydraulic parameters. In particular, the analysis of both synthetic and actual data

913

from homogeneous and heterogeneous soil profiles led to a confident estimation of the VGM shape914

parameters. Moreover, the analysis demonstrated that the residual and saturated water contents can

915

be successfully inferred from measured neutron fluxes if the soil approximately reaches a dry or wet

916

state, respectively. On the other hand, high uncertainty was observed for the saturated hydraulic917

conductivity, which suggests that the joint consideration of CRNP data with independent918

measurements of this parameter may be required. The simultaneous use of both neutron flux and46

Page 46 of 72

Page 47 of 72

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

919

water content data resulted in an appreciable reduction of uncertainty, in particular for the VGM

920

shape parameters and for the saturated hydraulic conductivity. This result opens new perspectives to921

estimate soil hydraulic properties from the combination of CRNP and other remote sensing

922

techniques. However, in order to have a precise assessment of the benefits of such joint inversion

923

for LSMs, it is worth investigating numerically the propagation of the uncertainty in values of soil

924

hydraulic parameters using LSMs. In this perspective, a model-based Bayesian experimental design925

could be used to optimize the location and the number of cosmic-ray probes and traditional sensors

926

to minimize the effect of the uncertainty on the land surface processes analyzed.

927

Furthermore, the comparison between the results obtained with the use of traditional near-928

surface point measurements and CRNP data indicates that the latter data are less influenced by the929

typical sources of uncertainty induced by land surface management activities (e.g., ploughing) and

930

lead to a more representative estimation of the soil hydraulic parameters, which is fundamental

931

when dealing with land surface hydrological modeling. In particular, the shape parameter α, which932

characterizes the soil retention capacity, was significantly overestimated when inversions were933

solely based on near-surface water content measurements, leading to a bimodal posterior

934

distribution representative of another plausible pore system induced by former land management

935

activities and macropores. On the other hand, the CRNP penetration depth and footprint

936

significantly mitigate these issues. Nevertheless, to what extent this difference in the estimation of937

the SHPs is significant for practical purposes depends on the type of analysis carried out by the938

modeler (e.g., surface runoff, infiltration, solute fluxes), and can be clarified only by propagating939

the estimated uncertainty. In this perspective, further studies are encouraged in order to better assess940

the information content gain of CRNP data in hydrological modeling.

941

In addition, the analysis presented in this manuscript suggests that the coupled HYDRUS-

942

COSMIC model can be used to detect and correct an imperfect calibration of the COSMIC

943

parameters. Avery et al. (2016) provided extensive results about the correlation between lattice944

water and the soil clay fraction. The relationship reported in this study can be used to better
47

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

945

constrain COSMIC parameters. Moreover, future research should consider the coupling of the946

COSMIC model with HYDRUS-3D to better account for the horizontal and vertical soil

947

heterogeneity at the field scale. In this perspective, the use of CRNP rover data could improve the948

description of the spatial patterns of SHPs (Gibson and Franz, 2018). A computationally cheap

949

alternative would be to use the kriging technique (Brunetti et al., 2017) to interpolate the HYDRUS-

950

1D predictions.

951
952
953

ACKNOWLEDGMENTS
The authors thank the TERENO (Terrestrial Environmental Observatories) project funded954

by the Helmholtz-Association and the Transregional Collaborative Research Center 32 “Patterns in955

Soil-Vegetation-Atmosphere Systems” funded by Deutsche Forschungsgemeinschaft (DFG) for

956

providing data for this study.

957
958

REFERENCES

959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981

Andreasen, M., Jensen, K.H., Desilets, D., Franz, T.E., Zreda, M., Bogena, H.R., Looms, M.C.,
2017. Status and Perspectives on the Cosmic-Ray Neutron Method for Soil Moisture
Estimation and Other Environmental Science Applications. Vadose Zo. J. 16, 0.
https://doi.org/10.2136/vzj2017.04.0086
Avery, W.A., Finkenbiner, C., Franz, T.E., Wang, T., Nguy-Robertson, A.L., Suyker, A.,
Arkebauer, T., Muñoz-Arriola, F., 2016. Incorporation of globally available datasets into theroving cosmic-ray neutron probe method for estimating field-scale soil water content. Hydrol.
Earth Syst. Sci. 20, 3859–3872. https://doi.org/10.5194/hess-20-3859-2016
Baatz, R., Bogena, H.R., Hendricks Franssen, H.J., Huisman, J.A., Qu, W., Montzka, C.,
Vereecken, H., 2014. Calibration of a catchment scale cosmic-ray probe network: A
comparison of three parameterization methods. J. Hydrol. 516, 231–244.
https://doi.org/10.1016/j.jhydrol.2014.02.026
Betancourt, M., 2017. A Conceptual Introduction to Hamiltonian Monte Carlo.
Beven, K., Binley, A., 1992. The future of distributed models: Model calibration and uncertainty
prediction. Hydrol. Process. 6, 279–298. https://doi.org/10.1002/hyp.3360060305
Bogena, H.R., Herbst, M., Huisman, J.A., Rosenbaum, U., Weuthen, A., Vereecken, H., 2010.
Potential of Wireless Sensor Networks for Measuring Soil Water Content Variability. Vadose
Zo. J. 9, 1002. https://doi.org/10.2136/vzj2009.0173
Bogena, H.R., Huisman, J.A., Güntner, A., Hübner, C., Kusche, J., Jonard, F., Vey, S., Vereecken,
H., 2015. Emerging methods for noninvasive sensing of soil moisture dynamics from field tocatchment scale: a review. Wiley Interdiscip. Rev. Water 2, 635–647.
https://doi.org/10.1002/wat2.1097
Bogena, H.R., Huisman, J.A., Schilling, B., Weuthen, A., Vereecken, H., 2017. Effective
48

Page 48 of 72

Page 49 of 72

982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

calibration of low-cost soil water content sensors. Sensors (Switzerland) 17.
https://doi.org/10.3390/s17010208
Bogena, H.R., Kunkel, R., Zacharias, S., Pütz, T., Schwank, M., Bens, O., Borg, E., Brauer, A.,
Dietrich, P., Hajnsek, I., Kunstmann, H., Munch, J.C., Papen, H., Priesack, E., Schmid, H.P.,
Wollschläger, U., Vereecken, H., 2012. TERENO – Long-term monitoring network forterrestrial research. Hydrol. und Wasserbewirtschaftung 56, 138–143.
Bogena, H.R., Montzka, C., Huisman, J.A., Graf, A., Schmidt, M., Stockinger, M., von Hebel, C.,
Hendricks-Franssen, H.J., van der Kruk, J., Tappe, W., Lücke, A., Baatz, R., Bol, R., Groh, J.,
Pütz, T., Jakobi, J., Kunkel, R., Sorg, J., Vereecken, H., 2018. The TERENO - The Rurhydrological observatory: A multi-scale multi-compartment research platform for theadvancement of hydrological science. Vadose Zo. J. https://doi.org/10.2136/vzj2018.03.0055
Brunetti, G., Šimůnek, J., Piro, P., 2016. A comprehensive numerical analysis of the hydraulicbehavior of a permeable pavement. J. Hydrol. 540, 1146–1161.
https://doi.org/10.1016/j.jhydrol.2016.07.030
Brunetti, G., Šimůnek, J., Turco, M., Piro, P., 2017. On the use of surrogate-based modeling for thenumerical analysis of Low Impact Development techniques. J. Hydrol. 548, 263–277.
https://doi.org/10.1016/j.jhydrol.2017.03.013
Carsel, R.F., Parrish, R.S., 1988. Developing joint probability distributions of soil water retention
characteristics. Water Resour. Res. 24, 755–769. https://doi.org/10.1029/WR024i005p00755
Dane, J.H., Hruska, S., 1983. In-Situ Determination of Soil Hydraulic Properties during Drainage.
Soil Sci. Soc. Am. J. https://doi.org/10.2136/sssaj1983.03615995004700040001x
Deng, H., Ye, M., Schaap, M.G., Khaleel, R., 2009. Quantification of uncertainty in pedotransfer
function-based parameter estimation for unsaturated flow modeling. Water Resour. Res. 45.
https://doi.org/10.1029/2008WR007477
Desilets, D., Zreda, M., Ferré, T.P.A., 2010. Nature’s neutron probe: Land surface hydrology at an
elusive scale with cosmic rays. Water Resour. Res. 46.
https://doi.org/10.1029/2009WR008726
Ek, M.B., Mitchell, K.E., Lin, Y., Rogers, E., Grunmann, P., Koren, V., Gayno, G., Tarpley, J.D.,
2003. Implementation of Noah land surface model advances in the National Centers forEnvironmental Prediction operational mesoscale Eta model. J. Geophys. Res. 108, 8851.
https://doi.org/10.1029/2002JD003296
Evans, J.G., Ward, H.C., Blake, J.R., Hewitt, E.J., Morrison, R., Fry, M., Ball, L.A., Doughty, L.C.,
Libre, J.W., Hitt, O.E., Rylett, D., Ellis, R.J., Warwick, A.C., Brooks, M., Parkes, M.A.,
Wright, G.M.H., Singer, A.C., Boorman, D.B., Jenkins, A., 2016. Soil water content insouthern England derived from a cosmic-ray soil moisture observing system – COSMOS-UK.
Hydrol. Process. 30, 4987–4999. https://doi.org/10.1002/hyp.10929
Feddes, R.A., Kowalik, P.J., Zaradny, H., 1978. Simulation of field water use and crop yield.
PUDOC, Wageningen.
Feroz, F., Hobson, M.P., Bridges, M., 2009. MultiNest: An efficient and robust Bayesian inference
tool for cosmology and particle physics. Mon. Not. R. Astron. Soc. 398, 1601–1614.
https://doi.org/10.1111/j.1365-2966.2009.14548.x
Finkenbiner, C.E., Franz, T.E., Gibson, J., Heeren, D.M., Luck, J., 2018. Integration ofhydrogeophysical datasets and empirical orthogonal functions for improved irrigation water
management. Precis. Agric. 1–23. https://doi.org/10.1007/s11119-018-9582-5
Foreman-Mackey, D., Hogg, D.W., Lang, D., Goodman, J., 2012. emcee: The MCMC Hammer.
https://doi.org/10.1086/670067
Franz, T.E., Zreda, M., Ferre, T.P.A., Rosolem, R., Zweck, C., Stillman, S., Zeng, X., Shuttleworth,
W.J., 2012. Measurement depth of the cosmic ray soil moisture probe affected by hydrogenfrom various sources. Water Resour. Res. https://doi.org/10.1029/2012WR011871
Gibson, J., Franz, T.E., 2018. Spatial prediction of near surface soil water retention functions usinghydrogeophysics and empirical orthogonal functions. J. Hydrol. 561, 372–383.
49

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083

https://doi.org/10.1016/J.JHYDROL.2018.03.046
Goodman, J., Weare, J., 2010. Ensemble samplers with affine invariance. Commun. Appl. Math.
Comput. Sci. 5, 65–80. https://doi.org/10.2140/camcos.2010.5.65
Gutmann, E.D., Small, E.E., 2007. A comparison of land surface model soil hydraulic properties
estimated by inverse modeling and pedotransfer functions. Water Resour. Res. 43.
https://doi.org/10.1029/2006WR005135
Hawdon, A., McJannet, D., Wallace, J., 2014. Calibration and correction procedures for cosmic-ray
neutron soil moisture probes located across Australia. Water Resour. Res. 50, 5029–5043.
https://doi.org/10.1002/2013WR015138
Hopmans, J.W., Nielsen, D.R., Bristow, K.L., 2002. How useful are small-scale soil hydraulicproperty measurements for large-scale vadose zone modeling? American Geophysical Union(AGU), pp. 247–258. https://doi.org/10.1029/129GM20
Huijser, D., Goodman, J., Brewer, B.J., 2015. Properties of the Affine Invariant Ensemble Sampler
in high dimensions.
Huisman, J.A., Rings, J., Vrugt, J.A., Sorg, J., Vereecken, H., 2010. Hydraulic properties of amodel dike from coupled Bayesian and multi-criteria hydrogeophysical inversion. J. Hydrol.
380, 62–73. https://doi.org/10.1016/J.JHYDROL.2009.10.023
Iden, S.C., Durner, W., 2007. Free-form estimation of the unsaturated soil hydraulic properties by
inverse modeling using global optimization. Water Resour. Res. 43.
https://doi.org/10.1029/2006WR005845
Ines, A.V.M., Mohanty, B.P., 2008. Near-Surface Soil Moisture Assimilation for Quantifying
Effective Soil Hydraulic Properties under Different Hydroclimatic Conditions. Vadose Zo. J.
7, 39. https://doi.org/10.2136/vzj2007.0048
Iwema, J., Rosolem, R., Rahman, M., Blyth, E., Wagener, T., 2017. Land surface model
performance using cosmic-ray and point-scale soil moisture measurements for calibration.
Hydrol. Earth Syst. Sci. 21, 2843–2861. https://doi.org/10.5194/hess-21-2843-2017
Kavetski, D., Kuczera, G., Franks, S.W., 2006. Bayesian analysis of input uncertainty inhydrological modeling: 1. Theory. Water Resour. Res. 42.
https://doi.org/10.1029/2005WR004368
Kennedy, J., Eberhart, R., 1995. Particle Swarm Optimization. Eng. Technol. 1942–1948.
Köhli, M., Schrön, M., Zreda, M., Schmidt, U., Dietrich, P., Zacharias, S., 2015. Footprint
characteristics revised for field-scale soil moisture monitoring with cosmic-ray neutrons.
Water Resour. Res. 51, 5772–5790. https://doi.org/10.1002/2015WR017169
Kool, J.B., Parker, J.C., 1988. Analysis of the inverse problem for transient unsaturated flow. Water
Resour. Res. https://doi.org/10.1029/WR024i006p00817
Lam, S.K., Pitrou, A., Seibert, S., 2015. Numba: a LLVM-based Python JIT compiler. Proc. Second
Work. LLVM Compil. Infrastruct. HPC - LLVM ’15 1–6.
https://doi.org/10.1145/2833157.2833162
Mantovan, P., Todini, E., 2006. Hydrological forecasting uncertainty assessment: Incoherence of
the GLUE methodology. J. Hydrol. 330, 368–381.
https://doi.org/10.1016/J.JHYDROL.2006.04.046
Matthews, C., Weare, J., Kravtsov, A., Jennings, E., 2017. Umbrella sampling: a powerful method
to sample tails of distributions.
Metropolis, N., Rosenbluth, A.W., Rosenbluth, M.N., Teller, A.H., Teller, E., 1953. Equation ofstate calculations by fast computing machines. J. Chem. Phys. 21, 1087–1092.
https://doi.org/10.1063/1.1699114
Mockler, E.M., Chun, K.P., Sapriza-Azuri, G., Bruen, M., Wheater, H.S., 2016. Assessing therelative importance of parameter and forcing uncertainty and their interactions in conceptualhydrological model simulations. Adv. Water Resour. 97, 299–313.
https://doi.org/10.1016/J.ADVWATRES.2016.10.008
Mohanty, B.P., 2013. Soil Hydraulic Property Estimation Using Remote Sensing: A Review.
50

Page 50 of 72

Page 51 of 72

1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

Vadose Zo. J. 12, 0. https://doi.org/10.2136/vzj2013.06.0100
Mohanty, B.P., Cosh, M.H., Lakshmi, V., Montzka, C., 2017. Soil Moisture Remote Sensing: Stateof-the-Science. Vadose Zo. J. 16, 0. https://doi.org/10.2136/vzj2016.10.0105
Montanari, A., 2007. What do we mean by ‘uncertainty’? The need for a consistent wording about
uncertainty assessment in hydrology. Hydrol. Process. 21, 841–845.
https://doi.org/10.1002/hyp.6623
Pelowitz, D.B., 2011. MCNPX User’s Manual Version 2.7.0.
Pitman, A.J., 2003. The evolution of, and revolution in, land surface schemes designed for climate
models. Int. J. Climatol. 23, 479–510. https://doi.org/10.1002/joc.893
Qu, W., Bogena, H.R., Huisman, J.A., Martinez, G., Pachepsky, Y.A., Vereecken, H., 2014. Effectsof Soil Hydraulic Properties on the Spatial Variability of Soil Water Content: Evidence from
Sensor Network Data and Inverse Modeling. Vadose Zo. J. 13, 0.
https://doi.org/10.2136/vzj2014.07.0099
Qu, W., Bogena, H.R., Huisman, J.A., Schmidt, M., Kunkel, R., Weuthen, A., Schiedung, H.,
Schilling, B., Sorg, J., Vereecken, H., 2016. The integrated water balance and soil data set of
the Rollesbroich hydrological observatory. Earth Syst. Sci. Data 8, 517–529.
https://doi.org/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. https://doi.org/10.1029/WR008i005p01204
Rivera Villarreyes, C.A., Baroni, G., Oswald, S.E., 2014. Inverse modelling of cosmic-ray soil
moisture for field-scale soil hydraulic parameters. Eur. J. Soil Sci. 65, 876–886.
https://doi.org/10.1111/ejss.12162
Santanello, J.A., Peters-Lidard, C.D., Garcia, M.E., Mocko, D.M., Tischler, M.A., Moran, M.S.,
Thoma, D.P., 2007. Using remotely-sensed estimates of soil moisture to infer soil texture andhydraulic properties across a semi-arid watershed. Remote Sens. Environ. 110, 79–97.
https://doi.org/10.1016/j.rse.2007.02.007
Schaap, M.G., Leij, F.J., Van Genuchten, M.T., 2001. Rosetta: A computer program for estimating
soil hydraulic parameters with hierarchical pedotransfer functions. J. Hydrol. 251, 163–176.
https://doi.org/10.1016/S0022-1694(01)00466-8
Scharnagl, B., Vrugt, J.A., Vereecken, H., Herbst, M., 2011. Inverse modelling of in situ soil water
dynamics: investigating the effect of different prior distributions of the soil hydraulicparameters. Hydrol. Earth Syst. Sci. 15, 3043–3059. https://doi.org/10.5194/hess-15-30432011
Schrön, M., Köhli, M., Scheiffele, L., Iwema, J., Bogena, H.R., Lv, L., Martini, E., Baroni, G.,
Rosolem, R., Weimar, J., Mai, J., Cuntz, M., Rebmann, C., Oswald, S.E., Dietrich, P.,
Schmidt, U., Zacharias, S., 2017. Improving calibration and validation of cosmic-ray neutron
sensors in the light of spatial sensitivity. Hydrol. Earth Syst. Sci. 21, 5009–5030.
https://doi.org/10.5194/hess-21-5009-2017
Shuttleworth, J., Rosolem, R., Zreda, M., Franz, T., 2013. The COsmic-ray Soil MoistureInteraction Code (COSMIC) for use in data assimilation. Hydrol. Earth Syst. Sci. 17, 3205–
3217. https://doi.org/10.5194/hess-17-3205-2013Šimůnek, J., Van Genuchten, M.T., 1996. Estimating unsaturated soil hydraulic properties fromtension disc infiltrometer data by numerical inversion. Water Resour. Res.
https://doi.org/10.1029/96WR01525
Šimůnek, J., van Genuchten, M.T., Šejna, M., 2016. Recent Developments and Applications of the
HYDRUS Computer Software Packages. Vadose Zo. J. 15, 25.
https://doi.org/10.2136/vzj2016.04.0033
Sokal, A.D., 1996. Monte Carlo Methods in Statistical Mechanics : Foundations and New
Algorithms. Funct. Integr. 131--192. https://doi.org/10.1007/978-1-4899-0319-8_6
Steenpass, C., Vanderborght, J., Herbst, M., Šimůnek, J., Vereecken, H., 2010. Estimating Soil
Hydraulic Properties from Infrared Measurements of Soil Surface Temperatures and TDR
51

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180

Data. Vadose Zo. J. 9, 910. https://doi.org/10.2136/vzj2009.0176
Sutanto, S.J., Wenninger, J., Coenders-Gerrits, A.M.J., Uhlenbrook, S., 2012. Partitioning of
evaporation into transpiration, soil evaporation and interception: A comparison between
isotope measurements and a HYDRUS-1D model. Hydrol. Earth Syst. Sci. 16, 2605–2616.
https://doi.org/10.5194/hess-16-2605-2012
Taylor, S.A., Ashcroft, G.L., 1972. Physical Edaphology: The Physics of Irrigated and NonirrigatedSoils. W.H. Freeman, San Francisco.
van Genuchten, M.T., 1980. A Closed-form Equation for Predicting the Hydraulic Conductivity of
Unsaturated Soils1. Soil Sci. Soc. Am. J. 44, 892–898.
https://doi.org/10.2136/sssaj1980.03615995004400050002x
Van Looy, K., Bouma, J., Herbst, M., Koestel, J., Minasny, B., Mishra, U., Montzka, C., Nemes,A., Pachepsky, Y.A., Padarian, J., Schaap, M.G., Tóth, B., Verhoef, A., Vanderborght, J., vander Ploeg, M.J., Weihermüller, L., Zacharias, S., Zhang, Y., Vereecken, H., 2017. Pedotransfer
Functions in Earth System Science: Challenges and Perspectives. Rev. Geophys.
https://doi.org/10.1002/2017RG000581
Vereecken, H., Kasteel, R., Vanderborght, J., Harter, T., 2007. Upscaling Hydraulic Properties andSoil Water Flow Processes in Heterogeneous Soils. Vadose Zo. J. 6, 1.
https://doi.org/10.2136/vzj2006.0055
Vogel, T., Cislerova, M., Hopmans, J.W., 1991. Porous Media With Linearly Variable Hydraulic
Properties. Water Resour. Res. 27, 2735–2741. https://doi.org/10.1029/91WR01676
Vrugt, J.A., Robinson, B.A., 2007. Improved evolutionary optimization from genetically adaptive
multimethod search. Proc. Natl. Acad. Sci. 104, 708–711.
https://doi.org/10.1073/pnas.0610471104
Vrugt, J.A., Stauffer, P.H., Wöhling, T., Robinson, B.A., Vesselinov, V. V., 2008a. Inverse
Modeling of Subsurface Flow and Transport Properties: A Review with New Developments.
Vadose Zo. J. 7, 843. https://doi.org/10.2136/vzj2007.0078
Vrugt, J.A., ter Braak, C.J.F., Clark, M.P., Hyman, J.M., Robinson, B.A., 2008b. Treatment of input
uncertainty in hydrologic modeling: Doing hydrology backward with Markov chain Monte
Carlo simulation. Water Resour. Res. 44. https://doi.org/10.1029/2007WR006720
Wesseling, J., Elbers, J., Kabat, P., Broek, B. Van den, 1991. SWATRE: instructions for input.
Intern. Note, Winand Star. Cent.
Western, A.W., Zhou, S.-L., Grayson, R.B., McMahon, T.A., Blöschl, G., Wilson, D.J., 2004.
Spatial correlation of soil moisture in small catchments and its relationship to dominant spatialhydrological processes. J. Hydrol. 286, 113–134.
https://doi.org/10.1016/J.JHYDROL.2003.09.014
Wöhling, T., Vrugt, J.A., 2011. Multiresponse multilayer vadose zone model calibration using
Markov chain Monte Carlo simulation and field water retention data. Water Resour. Res. 47.
https://doi.org/10.1029/2010WR009265
Zreda, M., Desilets, D., Ferré, T.P.A., Scott, R.L., 2008. Measuring soil moisture content noninvasively at intermediate spatial scale using cosmic-ray neutrons. Geophys. Res. Lett. 35,
L21402. https://doi.org/10.1029/2008GL035655
Zreda, M., Shuttleworth, W.J., Zeng, X., Zweck, C., Desilets, D., Franz, T., Rosolem, R., 2012.
COSMOS: The cosmic-ray soil moisture observing system. Hydrol. Earth Syst. Sci. 16, 4079–
4099. https://doi.org/10.5194/hess-16-4079-2012

52

Page 52 of 72

Page 53 of 72

1181

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

Figure 1. Precipitation (bars) and reference evapotranspiration ET0 (dashed line) used as

1182

input data in HYDRUS (A). Measured neutron fluxes (B) and water contents at three depths (C, D,

1183E) in the Rollesbroich catchment between April and September 2012.

1184

Figure 2. Precipitation (bars) and reference evapotranspiration ET0 (dashed line) from the1185

Rollesbroich catchment used as input data in HYDRUS (A). Simulated (solid lines) and synthetic1186(with added noise, diamonds) neutron fluxes (B) and water contents at z= -5 cm (C) between April

1187

and September 2012.

1188

Figure 3. The diagnostic plot of the AIES algorithm for synthetic modeling scenarios

1189(scenarios S1, S2, and S3 in A, B, and C, respectively) showing the integrated autocorrelation time

1190

against the algorithm step for each soil hydraulic parameter.

1191

Figure 4. Joint (below diagonal) and marginal (diagonal) posterior distributions of soil

1192

hydraulic parameters for Scenario S1. The red lines indicate the original values of the parameters

1193

for the loamy sand, while the isolines indicate 12, 39, 68, and 84% probability levels. The dashed1194

lines indicate the Bayesian credible intervals. The axis ranges correspond to the parameter bounds

1195

reported in Table 3.

1196

Figure 5. Joint (below diagonal) and marginal (diagonal) posterior distributions of the soil

1197

hydraulic parameters for Scenario S2. The red lines indicate the original values of the parameters

1198

for the loamy sand, while the isolines indicate 12, 39, 68, and 84% probability levels. The dashed

1199

lines indicate the Bayesian credible intervals. The axis ranges correspond to the parameter bounds

1200

reported in Table 3.

1201

Figure 6. Joint (below diagonal) and marginal (diagonal) posterior distributions of the soil

1202

hydraulic parameters for Scenario S3. The red lines indicate the original values of the parameters

1203

for the loamy sand, while the isolines indicate 12, 39, 68, and 84% probability levels. The dashed1204

lines indicate the Bayesian credible intervals. The axis ranges correspond to the parameter bounds

1205

reported in Table 3.

53

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

1206

Figure 7. The diagnostic plot of the AIES algorithm for the experimental data analysis

1207

scenarios (scenarios E1, E2, and E3 in A, B, and C, respectively) showing the integrated

1208

autocorrelation time against the algorithm step for each soil hydraulic parameter.

1209

Figure 8. Joint (below diagonal) and marginal (diagonal) posterior distributions of the soil

1210

hydraulic parameters for Scenario E3. The isolines and the dashed lines indicate the 12, 39, 68, and1211

84% probability levels and the Bayesian credible intervals, respectively. The axis ranges correspond1212

to the parameter bounds reported in Table 3.

1213

Figure 9. A comparison between the model predictions (grey lines) obtained by random

1214

sampling of 500 solutions from the posterior distribution (Scenario E3), and measured (black

1215

circles) neutron fluxes (A) and volumetric water contents at three different depths (B, C, D). The

1216

red lines indicate the model prediction with the median soil hydraulic parameters from the posterior

1217

distribution.

1218

Figure 10. Joint (below diagonal) and marginal (diagonal) posterior distributions of the soil

1219

hydraulic parameters for Scenario E1. The isolines and the dashed lines indicate the 12, 39, 68, and

1220

84% probability levels and the Bayesian credible intervals, respectively. The axis ranges correspond1221

to the parameter bounds reported in Table 3.

1222

Figure 11. A comparison between measured (black circles) and modeled (solid lines) depth-1223

weighted soil water contents. The grey and black lines indicate soil water contents calculated from1224

the CRNP using the standard calibration method (Baatz et al., 2014), and predicted using the model

1225

with the median soil hydraulic parameters from the posterior distribution from scenario E1,

1226

respectively.

1227

Figure 12. A comparison between the model predictions (grey lines) obtained by random

1228

sampling of 500 solutions from the posterior distribution (Scenario E1), measured (black circles)

1229

neutron fluxes (A), and volumetric water contents at three different depths (B, C, D). The red lines1230

indicate the model prediction with the median soil hydraulic parameters from the posterior

1231

distribution.
54

Page 54 of 72

Page 55 of 72

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

1232

Figure 13. Joint (below diagonal) and marginal (diagonal) posterior distributions of the soil

1233

hydraulic parameters for Scenario E2. The isolines and the dashed lines indicate the 12, 39, 68, and1234

84% probability levels and the Bayesian credible intervals, respectively. The axis ranges correspond1235

to the parameter bounds reported in Table 3.

1236

Figure 14. A comparison between the model predictions (grey lines) obtained by random

1237

sampling of 500 solutions from the posterior distribution (Scenario E2), measured (black circles)

1238

neutron fluxes (A), and volumetric water contents at three different depths (B, C, D). The red lines1239

indicate the model prediction with the median soil hydraulic parameters from the posterior

1240

distribution.

1241
1242

55

Page 56 of 72

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

1243
1244
1245

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

1246
1247
1248

Scenarioρbd (g/cm3)

θl (-)NCOSMIC (-)

Synthetic

1.42

0.037

Experimental

1.09

0.068

L1 (-)

L2 (-)

L4 (-)

189

161.98 129.14

3.16

213

161.98 129.14

3.16

Table 7. Synthetic modeling and experimental data analysis scenarios.
Scenario ID
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 flux time series

S3

-

1249
1250
1251

Table 8. The bounds of the VGM parameters and their true values used in the synthetic andexperimental scenarios.
Synthetic

Experimental
True Value

Parameter

Lower

Upper

Lower

Upper

Bound

Bound

(Carsel and Parrish,
Bound

Bound
1988)

θr (cm3/cm3)

0.03

0.08

0.06

0.07

0.07

θs (cm3/cm3)

0.35

0.45

0.41

0.30

0.45

log10α (1/cm)

-1.50

-0.50

-0.91

-2.50

-0.85

n (-)

1.50

3.00

2.28

1.10

2.00

1.00

3.00

2.54

0.50

3.00

log10Ks(cm/day)56

Page 57 of 72

1252
1253
1254
1255

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

φh1 (-)

1.00

1.00

1.00

1.00

5.00

φh2 (-)

1.00

1.00

1.00

1.00

5.00

φk1 (-)

1.00

1.00

1.00

1.00

20.00

φk2 (-)

1.00

1.00

1.00

1.00

5.00

Table 9. Original values and 95% Bayesian Credible Intervals of the soil hydraulic parameters forthe three synthetic scenarios. The 2.5 and 97.5% quantiles are the bounds of the credible intervals.
Scenario
True Values
Parameter

S1

S2

S3
(Carsel and Parrish, 1988)

1256
1257
1258
1259

2.5%

97.5%

2.5%

97.5% 2.5% 97.5%

θr (-)

0.05

0.06

0.05

0.06

0.05

0.06

0.057

θs (-)

0.35

0.44

0.35

0.44

0.35

0.44

0.41

log10α (1/cm)

-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/day)

1.92

2.78

2.06

2.93

2.00

2.86

2.54

Table 10. 95% Bayesian Credible Intervals of the soil hydraulic parameters for the threeexperimental scenarios. The 2.5 and 97.5% quantiles are the bounds of the credible intervals.
Scenario
E1

Parameter

E2

E3

2.5%

97.5%

2.5%

97.5%

2.5%

97.5%

θs (-)

0.33

0.38

0.35

0.39

0.36

0.37

log10α (1/cm)

-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/day)

0.50

2.82

0.50

1.44

0.97

1.18

φh1 (-)

1.00

4.32

1.00

2.33

1.00

1.05

φh2 (-)

1.27

5.00

1.00

4.35

1.80

2.21

57

Page 58 of 72

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

φk1 (-)

1.00

18.94

1.00

11.51

1.07

2.63

φk2 (-)

1.00

4.78

1.00

4.65

1.00

1.60

1260
1261

58

Page 59 of 72

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

Figure 1. Precipitation (bars) and reference evapotranspiration ET0 (dashed line) used as input data in
HYDRUS (A). Measured neutron fluxes (B) and water contents at three depths (C, D, E) in the Rollesbroichcatchment between April and September 2012.
210x190mm (300 x 300 DPI)

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

Figure 2. Precipitation (bars) and reference evapotranspiration ET0 (dashed line) from the Rollesbroichcatchment used as input data in HYDRUS (A). Original (simulated, solid lines) and synthetic (with added
noise, diamonds) neutron fluxes (B) and water contents at z= -5 cm (C) between April and September
2012.
201x120mm (300 x 300 DPI)

Page 60 of 72

Page 61 of 72

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

Figure 3. The diagnostic plot of the AIES algorithm for synthetic modeling scenarios (scenarios S1, S2, andS3 in A, B, and C, respectively) showing the integrated autocorrelation time against the algorithm step foreach soil hydraulic parameter.
213x82mm (300 x 300 DPI)

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

Figure 4. Joint (below diagonal) and marginal (diagonal) posterior distributions of soil hydraulic parameters
for Scenario S1. The red lines indicate the original values of the parameters for the loamy sand, while theisolines indicate 12, 39, 68, and 84% probability levels. The dashed lines indicate the Bayesian credibleintervals. The axis ranges correspond to the parameter bounds reported in Table 3.
299x299mm (200 x 200 DPI)

Page 62 of 72

Page 63 of 72

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

Figure 5. Joint (below diagonal) and marginal (diagonal) posterior distributions of the soil hydraulicparameters 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 Bayesiancredible intervals. The axis ranges correspond to the parameter bounds reported in Table 3.
299x299mm (200 x 200 DPI)

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

Figure 6. Joint (below diagonal) and marginal (diagonal) posterior distributions of the soil hydraulicparameters for Scenario S3. 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 Bayesiancredible intervals. The axis ranges correspond to the parameter bounds reported in Table 3.
299x299mm (200 x 200 DPI)

Page 64 of 72

Page 65 of 72

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

Figure 7. The diagnostic plot of the AIES algorithm for the experimental data analysis scenarios (scenarios
E1, E2, and E3 in A, B, and C, respectively) showing the integrated autocorrelation time against thealgorithm step for each soil hydraulic parameter.
218x83mm (300 x 300 DPI)

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

Figure 8. Joint (below diagonal) and marginal (diagonal) posterior distributions of the soil hydraulicparameters 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.
459x309mm (200 x 200 DPI)

Page 66 of 72

Page 67 of 72

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

Figure 9. A comparison between the model predictions (grey lines) obtained by random sampling of 500
solutions from the posterior distribution (Scenario E3), and measured (black circles) neutron fluxes (A) andvolumetric water contents at three different depths (B, C, D). The red lines indicate the model prediction
with the median soil hydraulic parameters from the posterior distribution.
457x254mm (200 x 200 DPI)

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

Figure 10. Joint (below diagonal) and marginal (diagonal) posterior distributions of the soil hydraulicparameters 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.
459x309mm (200 x 200 DPI)

Page 68 of 72

Page 69 of 72

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

Figure 11. A comparison between measured (black circles) and modeled (solid lines) depth-weighted soil
water contents. The grey and black lines indicate soil water contents calculated from the CRNP using thestandard calibration method (Baatz et al., 2014), and predicted using the model with the median soil
hydraulic parameters from the posterior distribution from scenario E1, respectively.
171x66mm (300 x 300 DPI)

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

Figure 12. A comparison between the model predictions (grey lines) obtained by random sampling of 500
solutions from the posterior distribution (Scenario E1), measured (black circles) neutron fluxes (A), andvolumetric water contents at three different depths (B, C, D). The red lines indicate the model prediction
with the median soil hydraulic parameters from the posterior distribution.
457x254mm (200 x 200 DPI)

Page 70 of 72

Page 71 of 72

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

Figure 13. Joint (below diagonal) and marginal (diagonal) posterior distributions of the soil hydraulicparameters 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.
459x309mm (200 x 200 DPI)

Vadose Zone J. Accepted Paper, posted 09/29/2018. doi:10.2136/vzj2018.06.0123

Figure 14. A comparison between the model predictions (grey lines) obtained by random sampling of 500
solutions from the posterior distribution (Scenario E2), measured (black circles) neutron fluxes (A), andvolumetric water contents at three different depths (B, C, D). The red lines indicate the model prediction
with the median soil hydraulic parameters from the posterior distribution.
457x254mm (200 x 200 DPI)

Page 72 of 72