1 Evaluating hydrodynamic parameters accounting for water retention hysteresis
2 in a large sand column using surface GPR

3 Emmanuel Léger<sup>a,c,*</sup>, Albane Saintenoy<sup>a</sup>, Yves Coquet<sup>b</sup>, Piotr Tucholka<sup>a</sup>, Hermann Zeyen<sup>a</sup>

<sup>a</sup>Université Paris-Saclay, CNRS, GEOPS, 91405, Orsay, France

<sup>b</sup>Université Paris-Saclay, INRAE, AgroParisTech, UMR ECOSYS, 78850 Thiverval-Grignon, France

<sup>c</sup>Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley, CA 94720, USA

7 Abstract

This study presents laboratory experiments using surface Ground Penetrating Radar measurements to exhibit and monitor hysteresis of the water retention function of Fontainebleau sand. A commercial impulsionnal GPR system monitored the volumetric water content changes in a large sand column subject to different hydraulic heads applied at its bottom during drainage/wetting cycles. Coupled hydrodynamic and electromagnetic modeling was used to simulate radargrams whereas hydrodynamic modeling coupled with 1D optical ray tracing was used to estimate the hydrodynamic parameters of the sand from GPR reflection two-way travel times. Statistical and uncertainty analysis were performed on numerical and experimental data. The range of optimized parameters obtained from experimental data were compared to those obtained with classical laboratory methods such as the hanging water column and the constant head permeameter. The range of parameters retrieved using GPR monitoring was consistent for the hydraulic conductivity at saturation and the van Genuchten parameters $\alpha_d$, $\alpha_w$. The difference between the GPR method and classical methods for the retrieved $\lambda$ and $n$ van Genuchten's parameters are believed to be caused by the dynamic character of the proposed method.

8 **Keywords:** Ground Penetrating Radar, Hydrodynamic properties of soils, Hysteresis, Vadose Zone

9 1. Introduction

10 Soil moisture is a key state variable for a wide variety of hydrological processes acting over a range of 11 spatial and temporal scales. For example, soil moisture is one of the factors influencing the partitioning of 12 rain fall into infiltration and runoff and the partitioning of net radiation into sensible and latent heat [1]. In 13 addition, it controls the subsurface drainage of water and thereby the leaching of chemicals to the groundwater 14 as well as nutrient and nitrogen [2] cycles. It is therefore evident that a thorough understanding of the soil 15 moisture dynamics is of primary importance for soil hydrological research and water budget prediction under 16 climate forcing [3].

17 The movement of water and solutes in variably saturated soils is to a large extent controlled by the 18 retention and the permeability properties of the soil matrix. The Soil-Water Retention (SWR) property is 19 expressed by the relationship between the soil water potential ($h$, expressed in cm) and the soil water content 20 ($\theta$ expressed in $cm^3.cm^{-3}$). Permeability to water is expressed by the Unsaturated Hydraulic Conductivity

<footnote>*Corresponding author</footnote>

<footnote>Email address: eleger@lbl.gov (Emmanuel Léger)</footnote>
(UHC) function, $K(h)$. These soil water properties are commonly described using analytical expressions such as van Genuchten's [4], or Brooks and Corey's [5]. This classical definition of the hydrodynamic functions assumes that $h$ is a unique function of $\theta$.

However, as first pointed out by Haines [6], this assumption is not valid because of hysteresis phenomena. At a given water head, the water content differs depending on whether the soil gets water (imbibition) or loses water (drainage). Phenomena resulting in SWR hysteresis are multiple: the “ink-bottle effect” due to irregular shapes of pores, differences between ascending and descending contact angles at an interface [7], entrapped air bubbles [8] or Haines jumps and temperature [9].

Historically, the basis of hysteresis theory, first observed in magnetism and adsorption on solids, was set in the works of Neel [10], Everett and Whitton [11]. Everett and Whitton [11] developed in parallel the independent domain theory, which was later suggested to be applicable to soil capillary hysteresis by Collis-George [12]. The assumption behind the independent domain theory of soil-water hysteresis is that the pore space can be divided into discrete pores, each of which drains and fills independently of the state of other pores. Subsequently Poulovassilis [13] accounted for capillary hysteresis in terms of independent domains with experiments on glass bead materials. He observed hysteresis in the water retention function depending on wetting or drying processes and found good agreement between measurements and expected values. Talsma [14] pursued this approach with success on sands. However, disagreement between theory and experiments was also reported for glass bead material [15], for sand [16], for sandy loam [17] and for silty loam and clay loam [18]. These studies suggest that the independent domain theory applies only to data obtained during static-equilibrium or steady state flow conditions but fails to depict hysteresis properly in case of unsteady flow [18].

Following Topp's conclusion, Mualem [19] improved and developed a new formalism based on the similarity hypothesis. This hypothesis is based on the assumption that the pore water distribution function can be described as a product of two independent functions, the pore opening radii distribution and the pore body radii distribution. This new formalism resulted in a succession of papers [19, 20] in which Mualem's formalism of hysteresis was improved.

In parallel to the conceptual models derived from physical representations of hysteresis, empirical models have been proposed based on fitting water retention curve data during wetting cycles. They are defined for specific types of soils such as the model applicable for sand of Kool and Parker [21], based on the scanning model of Scott et al. [22]. Several laboratory methods have been developed to obtain the data needed for empirical models of the SWR and UHC curves. Depending on soil type, different techniques were designed such as, the hanging water column [23], the evaporation methods [24], or the multi-step outflow method [25]. The hanging water column is one of the ways to obtain the SWR function parameters for the drainage and wetting phase whereas the multi-step outflow experiment allows to obtain the UHC parameters as well as the SWR function parameters. These methods, which we refer to as “classical”, are robust and efficient. However they need several days, are facility-demanding and cannot be run *in situ* on the field. They give the main
57
58

drying and main wetting curves delimiting the water content range one can measure at one soil water head $h$.

59

60

61

62

63

64

Estimating the SWR and UHC functions at the field scale is challenging leading to parameter estimation [26] and simplifying assumptions such as homogeneity, isotropy and non hysteretic behaviour. However, different studies demonstrated that the hysteretic behaviour of soils undergoing dynamic wetting and drying processes has to be taken into account for runoff predictions at catchment scale [27] or slope stability studies [28]. To overcome the difficult task of getting both drainage and wetting SWR, Lamorski and al. [29] used a machine learning approach to estimate the main wetting curve from the main drying curve.

65

66

67

68

Recently, another technique was tested on pyroclastic soils [30], with a laboratory experiment to define the main wetting and drying SWR curves. The authors developed a novel empirical technique from in situ monitoring in order to estimate the main wetting and drying curves. They redefined the van Genuchten $n$ parameters allowing more flexibility to design the inner scanning curves.

69

70

71

72

73

74

75

Advances in electronics in the past thirty years have significantly reduced the costs and improved the acquisition rates of geophysical methods [31]. This lead to the development of the field of hydrogeophysics [32]. These methods provide physical properties, such as electrical conductivity, nuclear magnetic resonance relaxation time constants and dielectric permittivity. These parameters are highly sensitive to soil composition, structure, density and water content. Among them, electrical resistivity measurements [33], nuclear magnetic resonance [34] and Ground Penetrating Radar (GPR) [31] are commonly used. Their non-invasive nature and their sensitivity to fine water content variations suit them for the vadose zone characterization.

76

77

78

79

80

81

Few studies dealt with geophysical measurements as data to be fitted to quantify the hysteresis of water retention and hydraulic conductivity functions [35, 36, 37, 38, 39]. Knight [37], followed by several others (i.e. [40]), found that the resistivity measured during imbibition in sandstone samples was consistently less than the one measured at the same saturation during drainage. She attributed this difference to the presence of electrical conduction at the air/water interface and concluded that geophysical properties depended on saturation history.

82

83

84

85

86

87

88

GPR data of Lai *et al.* [36] exhibited dielectric hysteresis for both soil and asphalt, depending on drying or wetting processes. More recently, Weihnacht and Boerner [35] performed a GPR monitoring of the soil water content in a large rectangular Plexiglas® tank filled with homogeneous sand. They used transmission radar data, one antenna on each side of the tank, at different depths to obtain the volumetric water content during a moving water table experiment. They found different sand water retention functions for drainage and imbibition. Although no hydrodynamic modeling was used, coupling GPR measurements with a multi-step outflow experiment was already an advancement.

89

90

91

92

Dagenbach et al. [38] and Klenk et al. [39] used a GPR test site to monitor the water table during fluctuating water table experiments. Both studies were inspired by the precursory works of Tsolfias et al. [41] and Endres et al. [42] interested in capillary fringe observation with GPR. Using a stationary GPR antenna at the soil surface, Dagenbach et al. [38] recorded capillary fringe reflections during drainage and
imbibition and compared them to those obtained numerically with three different water retention models, namely Brooks and Corey [43] and the constrained and unconstrained van Genuchten [4] models. Their measurements of reflection, amplitude and phase are similar to those of Saintenoy and Hopmans [44]. Klenk et al. [39] performed radargram profiles at different experimental times and studied the wavelet associated with peculiar reflections during drainage/relaxation/imbibition cycles. Their comparison between numerical modeling and experiments allowed them to obtain information on the relative hydraulic properties of the materials. All of these studies paved the way for hydrodynamic parameter estimation using geophysical methods and particularly GPR, but did not get to the quantification of these parameters.

Our study was multi-purposes: i) to demonstrate the hysteresis of the SWR function at the macroscopic (or so-called Darcy) scale using a surface commercial time domain GPR, and ii) to obtain the Mualem-van Genuchten (M-vG) [45, 4] parameters that characterize drainage and imbibition water retention as accurately as the classical laboratory methods do but with a much shorter experimental time. We used the empirical model of Kool and Parker [21] applied to sand for the hysteresis formalism.

In our study, we used commercial transmitter-receiver antennae (Ramac Malå 1600 MHz) placed at the top of a column of sand subjected to controlled water table variations. We monitored the time of reflection from the water column bottom during different drainage-imbibition cycles, which we then inverted to obtain the parameters. We derived experimental uncertainty from multiple inversions of independent noise-added synthetic data. After numerical tests, we inverted experimental GPR data acquired in a laboratory experiment to retrieve hydraulic properties in drainage or wetting condition of a large sand column. We account for dispersive frequency shifts in our procedure. Statistical analysis used bootstrapping. We compared our results with the M-vG parameters estimated from classical hanging water column and constant head permeameter experiments.

## 2. Background

### 2.1. Soil hydraulic functions

The soil hydraulic functions can be described by several mathematical expressions requiring a different number of parameters depending on the chosen model. For the $\theta(h)$ relationship, one of the largest groups of models represented in the literature is the power function model of the form:

$$ (1 - S_e^{-a})^b S_e = \alpha h, \qquad (1) $$

where $S_e = \frac{\theta - \theta_r}{\theta_s - \theta_r}$ is the effective saturation, $\theta_s$ is the saturated volumetric water content, $\theta_r$ the residual volumetric water content, and $\alpha$ (>0), a and b, are fitting parameters. Many models based on this power form have been proposed to fit SWR data [43, 4].

In the current study, we used the equation from van Genuchten [4],

$$ S_e = (1 + (\alpha h)^n)^{-m}, \qquad (2) $$
where $\alpha$ is proportional to the inverse of the absolute value of the soil matric potential at the inflection point of the $\theta(h)$ curve and $n$ determines the slope of the curve at the inflection point. The two dimensionless fitting parameters, $n$ and $m$, are linked by

$$m = 1 - \frac{k}{n}, \qquad (3)$$

with $n > k$ [46]. For simplicity, $k$ takes integer values which correspond to different models giving closed form analytic expressions for the UHC functions. The case $k = 1$ corresponds to the model of Mualem [45], while $k = 2$ gives the conductivity model of Burdine [47]. Using $k = 1$, the UHC function is written as

$$K(\theta) = K_s S_e^{\lambda} \left[ 1 - \left[ 1 - S_e^{\frac{n}{n-1}} \right]^{\frac{n-1}{n}} \right]^2, \qquad (4)$$

with $K_s$ the saturated hydraulic conductivity and $\lambda$ a parameter accounting for pore tortuosity. The whole description of water flow in the unsaturated zone, considering no hysteresis, with equations (2) and (4) requires the determination of 6 parameters: $\theta_r, \theta_s, n, \alpha, K_s$ and $\lambda$.

The choice of van Genuchten's model for describing the SWR function was motivated by i) the good fitting to our hanging water column data and ii) hydrodynamic modeling reasons since it was the only one taking into account hysteresis in the Hydrus1D code that we used for hydrodynamic modeling (see section 2.5).

## 2.2. Hysteresis of the SWR function

A schematic representation of hysteresis for the SWR curve is presented in Figure 1. The boundary hysteresis loop consists of the main drying (red) and main wetting (green) curves. In this theoretical representation, the main drying curve is described by the van Genuchten parameter vector, $vG^d = \{\theta_r^d, \theta_s^d, n^d, \alpha^d\}$. The main wetting curve is described by $vG^w = \{\theta_r^w, \theta_s^w, n^w, \alpha^w\}$. Considering only the hysteresis on the SWR function, we are left with a total of eight unknown parameters to determine.

Based on the work of [48], we will assume

$$\theta_r^w = \theta_r^d = \theta_r, \qquad (5)$$

$$\theta_s^w = \theta_s^d = \theta_s, \qquad (6)$$

and

$$n^w = n^d = n. \qquad (7)$$

The two first conditions imply that we disregard changes in the volume of entrapped air during rewetting due to temperature differences [8]. The last condition imposes that the curvature of the two SWR functions are similar. Imposing these constraints results in some loss of flexibility in describing the hysteresis of the SWR function, since the main drying and main wetting curves differ only by their $\alpha$ values. However, we will see later in section 4.1, that this last assumption is representing well the hysteresis cycle measured for the Fontainebleau sand with the hanging water column method.
## 2.3. Electromagnetic wave velocity as a function of soil dielectric properties

Surface GPR consists of two antennae, being electric dipoles, positioned on the surface of the medium to be sounded. The transmitting antenna emits short pulses of spherical electromagnetic (EM) waves in response to an excitation current source. The receiving antenna converts the incoming EM fields to electrical signals. Following the work of Annan [49], the velocity of an electromagnetic wave propagating in a non-magnetic soil, with low electrical conductivity, is

$$v = \frac{c}{\sqrt{\varepsilon'}}$$ (8)

where $\varepsilon'$ denotes the real part of the relative dielectric permittivity and $c = 0.3$ m/ns is the velocity of EM waves in air.

## 2.4. Petrophysical relationships

The bulk soil has a frequency-dependent relative dielectric permittivity $\varepsilon_b$, expressed by a complex number with the loss factor represented by the imaginary part, defined as:

$$\varepsilon_b = \varepsilon_b'(f) - i \left[ \varepsilon_b''(f) + \frac{\sigma_b}{2\pi f \varepsilon_0} \right],$$ (9)

where $f$ is the frequency, $\varepsilon_b'(f)$ the real part of soil permittivity, $\varepsilon_b''(f)$ the imaginary part, including the relaxation loss, and $\sigma_b$ the bulk conductivity leading to conductivity losses. The permittivity of free space is denoted $\varepsilon_0$.

The dielectric permittivity depends on the water content. We used the volumetric CRIM model [50], where the relative dielectric permittivity of the porous geological material, $\varepsilon_b$, is a function of the material porosity, its state of saturation and the respective permittivity of each of its individual constituents. We carried out our experiments over a Fontainebleau sand which is a tri-phasic medium composed of water, air and silica. Then,

$$\varepsilon_b^\gamma = \theta\varepsilon_w^\gamma + (1-\phi)\varepsilon_s^\gamma + (\phi-\theta)\varepsilon_a^\gamma,$$ (10)

where $\varepsilon_w$, $\varepsilon_s$ and $\varepsilon_a$ are respectively the relative dielectric permittivity of water, silica and air, $\phi$ the porosity and $\theta$ the volumetric water content. $\gamma$ is an empirical coefficient that depends on soil structure, set to 0.5 in this study. The relative dielectric permittivity of air is equal to 1. We measured $\varepsilon_s = 2.23$ using a time-domain reflectometer in laboratory experiments. We assumed that the solid phase was non conductive and therefore its dielectric permittivity was independent of frequency. The permittivity of water $\varepsilon_w$ is a function of frequency and conductivity, and is expressed in time domain by Debye's equation:

$$\varepsilon_w = \varepsilon_\infty + \frac{\varepsilon_{st} - \varepsilon_\infty}{1 + i\omega\tau},$$ (11)

in which $\varepsilon_\infty = 4.9$ is the high-frequency limit, and $\varepsilon_{st} = 80.088$ the static value, of the real relative dielectric permittivity at 20°C [51], $\tau$ is the relaxation time of water at 20°C [52]. We further assumed that there was no loss due to bound water, or at least it was negligible in our Fontainebleau sand.
## 2.5. Hydrodynamical modeling

We considered one-dimensional vertical soil water flow described by the 1D-Richards' equation expressed in terms of water content as

$$ \frac{\partial \theta}{\partial t} = \frac{\partial K(\theta)}{\partial z} + \frac{\partial}{\partial z} \left[ K(\theta) \frac{dh}{d\theta} \frac{\partial \theta}{\partial z} \right], \qquad (12) $$

where $K(\theta)$, and $\theta(h)$ are described by the M-vG model, presented in equations (2) and (4).

We used the Hydrus-1D code [53] to solve this equation, with an atmospheric (Neumann) boundary condition at the top. We changed the bottom boundary conditions with time (through a Dirichlet boundary condition) using the evolution of the applied hydraulic head we used in our experiments. Our domain was divided into 5-mm thick horizontal layers, thin compared to the shortest wavelength of the electromagnetic waves propagating through it.

## 2.6. Electromagnetic modeling

Among the numerous techniques available for simulating GPR data [54], we chose the open source software GprMax2D [55]. It solves Maxwell's equations in two dimensions with Finite Difference Time Domain (FDTD) techniques. All synthetic radargrams presented hereafter have been simulated using GprMax2D (See, e.g., Figure 4).

However, to avoid time-consuming FDTD modeling in our inversion algorithm, we computed the Two Way Travel-time (TWT) of the wave reflected at the bottom of the tank assuming 1D propagation through a stack of 5 mm thick layers. The electromagnetic wave inside each layer was computed with equation (8) taking the dielectric permittivity value defined by equations (10) and the water content distribution modeled by HYDRUS-1D. Further details on this procedure are presented in [56].

# 3. Numerical experiments

## 3.1. Set-up

We considered a homogeneous cylindrical medium (71 cm height and 60 cm diameter), with a porosity set to 0.39 and hydrodynamic properties in the range of what we measured with the water hanging column and the UMS Ksat System on Fontainebleau sand (section 4.1). These parameters are summarized as $p_{in}$ in Table 1. A transmitter (emitting a Ricker wavelet centered on 1000 MHz) and a receiver were positioned at the sand column surface. The draft of the apparatus is presented in Figure 2. The tube on the right indicates the different water table levels set during the experiment. The two triangles symbolize the GPR transmitter and receiver. A very thin layer of Perfect Electrical Conductor (PEC) was positioned at the bottom of the sand column. At the initial state, the water level was set at the top of the column, corresponding to $h_0 = 71$ cm.

Different hydraulic heads were applied at the bottom of the soil cylinder. The hydraulic heads were chosen in agreement with experimental data, in order to later compare numerical radargrams with experimental ones.
Two types of numerical experiments were conducted: i) a multi-step (subscript MS) experiment where we applied successively $h_1 = 41.5$ cm, $h_2 = 31.5$ cm, $h_3 = 22.5$ cm, $h_4 = 0$ cm, $h_5 = h_3$, $h_6 = h_2$, $h_7 = h_1$ and $h_8 = h_0$ at experimental times $t_{MS} = [100; 1600; 2700; 4000; 5800; 6900; 7900; 9200; 10500]$ seconds, and ii) a one-step (subscript OS) experiment where the heads applied to the bottom of the cylindrical medium were $h_0 = 71$ cm, $h_1 = 0$ cm and $h_2 = h_0$ at experimental times $t_{OS} = [700; 6600; 8100]$.

A volumetric water content profile was computed at each experimental time step using Hydrus-1D with the M-vG parameters $p_{in}$ presented in Table 1. As presented in the part 2.2, the hysteresis effect was fully represented by the parameter $\alpha$, taking $\alpha_w = 2\alpha_d$ according to Kool and Parker [21] and Simunek [48]. The two water retention boundary curves (wetting and drying curves) corresponding to the parameters $p_{in}$ are shown on Figure 1. The GPR antenna set at the top was used to monitor the water content of the sand, using the TWT of the reflection on the PEC surface at the bottom of the tank.

## 3.2. Forward modeling

### 3.2.1. 1D velocity modeling

The water content profiles generated at each experimental time step were converted to permittivity profiles using relations (10). The TWT of the PEC reflection was computed using 1D velocity modeling. Figure 3 displays the TWTs obtained when varying the water table levels with (a) the multi-step protocol, and the (b) one-step protocol (b). In both cases, we present the TWTs obtained from the simulated water content profile accounting for hysteresis of the water retention function (red circles) and without hysteresis (blue curves). As expected, since the simulated water content profiles during the drying and wetting cycle were different for the same applied heads, the TWTs are different as well. These differences are more emphasized in the multi-step protocol (Figure 3-a) compared with the one-step (Figure 3-b). In this last case, the difference, between hysteresis and non-hysteresis, resides in only 6 points in the wetting phase. For this reason, we decided to continue with only the multi-step experimental protocol.

### 3.2.2. FDTD modeling highlighting frequency dispersion

The water content distribution outputs from Hydrus-1D were used as inputs for the GprMax2D code using equations (10) and (11). The simulated GPR monitoring of the dynamic water level variation is presented in Figure 4 for the multi-step case. The reflection on the PEC layer at the base of the sand is arriving between 15 and 22 ns. On this figure we see that the width of the reflected pulse is changing during the experiment. This is because of relaxation losses due to the presence of water, and to the dielectric permittivity gradient created by the water retention properties of the medium the electromagnetic wave is traveling through. The dominant wave period associated with the bottom tank reflection found between 15 and 22 ns evolves with the applied head. The dispersion of the wave happening along the experiment was evaluated by determining the changes in the frequency spectrum computed with an instantaneous Fast Fourier Transform (FFT) applied over a selected time window containing the reflected wavelet. This is represented in Figure 5, where the maximum of the instantaneous frequency associated with this reflection is presented for each experimental
time step. On this figure, we display the picking of the maximum instantaneous frequency as circles and the TWT of the bottom tank reflection picked from the radargrams presented in Figure 4 (red line). The drier the column, the higher is the frequency content, and, similarly, the wetter the column, the more present are the relaxation effects, the steeper is the dielectric permittivity gradient, and the smaller is the frequency content of the bottom tank reflection.

### 3.3. Uncertainty analysis

#### 3.3.1. Data uncertainty quantification

As presented above, even if we are taking into account the frequency shift in our 1D-velocity modeling through the computation of the dielectric permittivity at a given frequency (equations 10 and 11), the 1D-velocity modeling does not represent all the electromagnetic phenomena taking place, and we observe a small discrepancy between picked arrival times of the reflection present in the radargram of Figure 4 and our 1D velocity algorithm results. We evaluated this numerical discrepancy, by using the RMSE, to 0.09 ns.

An other source of uncertainty in our data is the one potentially associated with the applied head. For determining the TWT uncertainties associated with the applied hydraulic head uncertainties, we considered the RMSE between the TWT resulting from the hydraulic head presented in section 3.1 for the multi-step experiment, and other TWT generated with the hydraulic heads decreased by 1 cm. The RMSE between these TWT is 0.11 ns.

At last, we considered the picking uncertainties by considering that time picking can not be better than one time step, which was 0.058 ns for experiments to be presented in section 4.3. Even if the time-step was much smaller for the numerical experiment we chose to take the experimental time step in order to better represent the uncertainty.

Taking these three sources of uncertainties in our multi-step numerical experiments, we estimated our overall uncertainty associated with each TWT data point, using a quadratic summation, to $\sigma_{MS} = 0.13$ ns.

#### 3.4. Inverting numerical data perturbed by Gaussian noise

We used the TWT of the bottom tank reflection as presented in Figure 3 as data to be fitted to get the set of M-vG parameter values $p = \{\theta_r, \alpha_d, \alpha_w, n, K_s, \lambda\}$. We did not invert the volumetric water content at saturation $\theta_s$, fixed at 0.39 cm³/cm³.

The chosen optimization procedure was the Shuffled Complex Evolution (SCE-UA) algorithm [57] minimizing the objective function expressed as the RMSE between the modeled and observed TWT of the bottom tank reflection to be fitted. This RMSE was given by

$$RMSE = \sqrt{\frac{1}{N} \sum_{j=1}^{N} [\tau_{xp}(t_j) - \tau_{mod}(t_j, p)]^2}, \quad (13)$$

where $\tau_{xp}$ is the radargram-measured TWT generated with the hydrodynamic parameters $p_{in}$ at experimental time $t_j$, and $\tau_{mod}$ is the modeled TWT reflection at time $t_j$ subject to soil hydraulic parameter vector $p$, and $N$ is the number of points used for comparison, i.e. the number of considered experimental times $t_j$.
We sampled the parameters in ranges defined according to the values found in the literature. The range for $\theta_r$ was [0.01; 0.1] cm$^3$/cm$^3$. $K_s$ interval was [0.05; 3] cm.min$^{-1}$, $n$ was allowed to get values between 2 and 10, and $\alpha_W$ and $\alpha_d$ ranges were both set to [0.01; 0.1] cm$^{-1}$. We set the $\lambda$ range to [0; 10], considering this parameter as a fitting one. The historical paper of [45] proposed $\lambda$ values in the range of [-1; 2.5] and concludes on good fitting for most soils with $\lambda = 0.5$. This value is the most used in the literature. The studies of Schaap et al [26] were among the first to develop a wider range for $\lambda$ parameter [-1; 1]. We based our $\lambda$ range on the study of [58], on the database of European soils, where $\lambda$ could reach values higher than unity.

For each generated parameter set $p$, we computed the TWT of the bottom tank reflection using the 1D velocity model. This inversion modeling procedure is presented as a flow diagram in Figure 6, where the objective function is the RMSE. We used 20 complexes, 50,000 iterations for each loop and a constraint on 30,000 iterations to stop if no better model was found during the Complex search.

Since the inversion algorithm, SCE-UA, is not following the Bayesian rule, the statistical analysis performed on one single inversion cannot be considered as mathematically correct [59]. We decided to perform a statistical analysis on models obtained by independent SCE-UA optimization procedures carried out on different noise-added generated TWT data. Examples of the perturbed TWT are presented in Figure 5 as purple points. These data correspond to the original simulated ones, in blue, to which we added small time delays sampled from a Gaussian distribution centered on the original value, with standard deviation $\sigma_{MS} = 0.13$ ns.

We performed 80 inversions on 80 different Gaussian noise-added data. We selected the models with associated RMSE smaller than or equal to $2\sigma_{MS} = 0.26$ ns. Possible double models were taken out in order to keep the same weight for all models, and not influence the histograms. This resulted in 968,291 models. The statistical distributions of each parameter of all the selected models are presented in Figure 8. The value for each parameter set as input is represented by the vertical dashed orange bar and can be found in Table 1. These histograms show that the input parameters are included in every distribution.

Figure 9 shows objective function plots for four couples of parameters ($K_s$, $\alpha_w$), ($\alpha_w$, $n$), ($\alpha_w$, $\alpha_d$), ($\lambda$, $n$). These objective function plots show that the regions with lowest RMSE (dark blue areas) always contain the input parameter, except for $\lambda$.

The case of the $\lambda$ parameter and its link with the other parameters require another analysis. On Figures 5-a) we displayed 3 TWT curves of the bottom tank reflection. The set of parameters taken for the original blue line data is the one set as input (Table 1). The two other lines represent the result of an increase or a decrease of the $n$ value by 2. Increasing $n$ results in increasing the TWT variations and vice-versa. Figure 5-b) shows as well the effect of an increase or a decrease of the parameter $\lambda$. At the opposite of $n$, increasing $\lambda$ by 2 results in decreasing the TWT variations. From these figures, we see that different couples of ($n$, $\lambda$) will fit as well the noisy data (purple points). A too small $n$ value associated with a too small $\lambda$ value will result in the same curve as too high values for both $n$ and $\lambda$. This is illustrated by the objective function
312 plot Figure 9-d).

This numerical analysis showed the attainable accuracy of our modeling for inverting most of the hydro- 313 dynamic parameters as well as some non-unicity due to inter-correlations between parameters such as $n$ and 314 $\lambda$. However, we see that working under dynamical conditions, without waiting for static equilibrium at each 315 step, did not prevent us from retrieving most of the hydrodynamical parameters even if we applied a quite 316 simple uncertainty analysis. For instance, we applied the same uncertainty to every point of the TWT data 317 not considering the influence of the frequency shift uncertainties being stronger when the soil is becoming 318 drier. In addition, an uncertainty of $\pm 1$ cm on the hydraulic head is a pessimistic estimation. 319

## 4. Laboratory Measurements

### 4.1. Characteristics of the sand

The van Genuchten water retention parameters were obtained by fitting the data obtained with three 322 Hanging Water Column (HWC) experiments performed on Fontainebleau sand core (250 cm³) samples packed 323 at $1.67 \text{ g.cm}^{-3}$ bulk density. We considered uncertainties of 1 cm on $h$ and $0.03 \text{ cm}^3/\text{cm}^3$ on $\theta$. The van 324 Genuchten parameters obtained are presented in Table 2, under the label “Lab. meas.”. The uncertainties 325 on the van Genuchten parameters retrieved with HWC were determined by generating 10,000 models fitting 326 10,000 WRC generated from the noised WRC, with the uncertainties on $\theta$ and $h$ presented above. 327

Figure 10 shows the retention curve data acquired from the HWC experiments. They were fitted using the 328 van Genuchten model. As expected for this type of soil, $\alpha_w$ was close to $2\alpha_d$. The HWC results agree with 329 the hypothesis that the $n$ parameter could be considered as independent from drying or wetting processes, 330 since the curvatures of the water retention curves are very similar for both drying and wetting branches. 331

The saturated hydraulic conductivity was determined on several 250 cm³ soil core samples packed at 332 $1.6 \text{ g/cm}^3$, using the UMS Ksat System (https://www.metergroup.com/) under constant head. The range of 333 values obtained are presented in Table 2. The $\lambda$ parameter was not determined, but its value was initially 334 expected to be close to 0.5 [45] without having any other measurements. 335

Granulometric determinations of the Fontainebleau sand used for the experiments showed that 13% of its 336 grains had a diameter larger than 0.25 mm and 98% larger than 0.125 mm. Its average particle density was 337 $2.8 \pm 0.1 \text{ g.cm}^{-3}$. The sand was packed to a porosity estimated to be 0.39. 338

### 4.2. Experimental Set-up

The experimental apparatus consisted of two cylindrical tanks (80 cm height and 60 cm diameter), one 340 tank, $T_S$, filled with the Fontainebleau sand and the other one, $T_w$, filled with water at a level corresponding 341 to the hydraulic head. $T_s$ was filled with 321 kg of sand, compacted manually in order to obtain a bulk 342 density of $\rho_b = 1.7 \text{ g.cm}^{-3}$. The sand layer was set on top of a 5 cm gravel layer in order to drain and 343 wet homogeneously the sand. To obtain a strong reflection at the bottom of the sand we put between the 344 gravel layer and the sand a thin aluminum sheet, finely meshed, to let the water go through. Three pictures 345
of the experimental set-up are presented in Figure 11. Figure 11-a) shows $T_w$, where we see the holes used as overflow to maintain the water table during drainage and where we see the green water supply pipe. $T_s$ and $T_w$ bottoms were linked through four independent plastic tubes with taps (Figure 11-b). Figure 11-c) displays $T_s$ with the antenna set at the top of the sand, maintained in position by a load. The white thin layer is a geo-textile enclosing the sand. Piezometers were connected at the bottom of each tank to indicate the applied head.

### 4.3. *Experimental measurements*

Several experiments were carried out using this experimental setup. For each experiment, the whole soil profile was initially water saturated during hours in order to be on the main drainage curve during the first cycle. Both tanks were hydraulically separated by closing the taps on the linking tubes. The water level was set to the desired height (for example 42 cm) in $T_w$ using an overflow hole. Then the hydraulic contact between both tanks was established again, opening the taps. The water drained out of $T_s$ was evacuated through the overflow hole in $T_w$. This procedure was repeated for all steps during draining and wetting processes. The radargram measured during one of the multi-step experiment is shown in Figure 12-a). The horizontal axis is the experimental time in seconds, and the vertical axis is the TWT in nanoseconds.

During the drainage and imbibition cycle, the direct wave signal evolves through time, due to the evolution of the media properties surrounding the antennae, which affects the antennae coupling. The reflection on the aluminum sheet at the base of the sand is arriving around 28 ns, other reflections are tank side effects. As expected, the TWT of this reflection evolves with the water table level. In the multi-step experiment, the heads were changed at [800;1500;2000;2900] seconds for drainage and [4400;5200;5800;6800] seconds for re-wetting. On Figure 12-b, we display the instantaneous maximum frequency of the considered reflection. On this figure we display as well the TWT of this reflection, relative to its initial TWT, when saturated. We see that the behavior of this reflection is comparable to the one we simulated (Figures 3, 4 and 5). Despite the fact that we used a Malå Ramac antenna with a nominal frequency of 1.6 GHz, we clearly see that the center frequency is far from 1.6 GHz, due to antenna coupling, Debye dispersion and the dielectric permittivity gradient inside the sand column.

### 4.4. *Inversion of Experimental Data*

#### 4.4.1. *Frequency shift correction*

The frequency shift observed in the experimental data of Figure 12-b) has to be taken into account in order to correct for the higher frequency picking, since all our TWT are relative to the saturated case (low frequency). Using the following procedure, we corrected the TWTs due to the frequency variations during the experiment.

We considered that the wavelet associated with the bottom tank reflection could be assimilated as a time derived Ricker wavelet. We then measured the time delay between the maximum peak arrival of a wavelet computed with the lower frequency of 460 MHz (observed when the sand column is fully saturated, see
Figure 12), and the maximum peak arrival of a wavelet computed with the frequency associated with the experimental time at which we want to apply a correction. This principle is represented in Figure 13. As expected, when the frequency difference between the initial measurement with the current one is larger the maximum time difference is larger. The maximum time delay to be corrected is 0.7 ns. It is too large to be taken into account as uncertainty on the data to be inverted. Instead, all numerical TWT were corrected with the computed time delays before inversion.

### 4.4.2. Data inversion and uncertainty analysis with bootstrapping

Using the same inversion loop as the one presented in the numerical case summarized in Figure 6, we inverted the experimental TWT data, re-sampled every 20 seconds, to obtain the hydrodynamic parameters of the considered soil. The water content at saturation was fixed to the value of porosity, $\theta_s = 0.39 \text{ cm}^3/\text{cm}^3$. The parameters minimizing the objective function (Eq. 13) between the simulated and the experimental data are presented in Table 2.

We will present the uncertainty analysis only on one experiment, the multi-step experiment labeled number 1.

The uncertainty analysis was based on running multiple SCE-UA inversions on down-sampled data. This approach is similar to bootstrapping without replacement. For each experimental case, we chose randomly 50 sets of 80 points (40 points during the drainage and 40 points during the rewetting). We performed inversions on each bootstrapped data subset. We selected the models with RMSE values up to 10% larger than the RMSE of the best fitting model. Possible double models were taken out. This resulted in 231,192 parameter sets. The distributions of all model parameters are presented as histograms in Figure 14.

### 4.5. Discussions on the experimental results

The best fitting parameter sets obtained from two different laboratory experiments, summarized in Table 2 and in the histograms issued from the bootstrapping analysis show differences with the classical results that are challenging to explain. The main differences between results concern the parameters $\lambda$ and $n$. Figures 14 show that the factor between $\alpha_d$ and $\alpha_w$ is more than the expected factor two. The histograms of $\theta_r$, $n$ and $\lambda$ optimized parameters do not match the HWC parameters distribution determined in the laboratory. Many reasons have to be considered while comparing GPR-derived parameters and classical laboratory parameters (HWC). Three main reasons, which we consider to be the most influencing, are presented below.

First, the parameters obtained in laboratory are derived from $250 \text{ cm}^3$ core samples. Since about 330 kg of Fontainebleau sand were used in the large tank, the compaction was definitely not as homogeneous as in the case of small $250 \text{ cm}^3$ core samples, despite the fact that we tried to reach the same bulk density. In addition, evaporation could have influenced the HWC measurements.

Second, the hydraulic head range covered by the HWC experiment is different than the one covered by the large-scale experiment. Indeed, looking back at the HWC data presented in Figure 10, we see that the first drainage curve reaches 140 cm of suction, whereas in the case of our dynamical experiment we are dealing
with at the most 70 cm. Then it is possible that during the dynamic experiment we are not on the main wetting curve, but rather on an inner cycle curve. However this effect should be taken into account by our hydrodynamic modeling.

Third, the HWC measurements indicate a high value of *n*, even higher than the upper value found with bootstrapped data inversion. The HWC experiment is giving data at hydrostatic equilibrium. Therefore they do not involve hydraulic conductivity and the *n* parameter, fitting these data, is independent of the λ value, which is different from our dynamic GPR experiment.

As a result of these three facts, differences between the retention curve parameters obtained by our GPR experimental data inversion and those obtained through classical laboratory experiments (HWC and Ksat) do not mean that the quantification from GPR data is wrong. At the contrary, we think that it is pointing out that classical laboratory experiments may not give retention curve fitting parameters that are usable for dynamic process models at the macroscopic-scale.

At the opposite, the hydraulic conductivity at saturation $K_s$ is consistently estimated in all our bootstrapped data inversions in the range [0.6;1.5] cm.min⁻¹. This interval is consistent with the constant head permeameter estimations, [0.96;1.5] cm.min⁻¹.

Estimated values of λ parameter are larger than 12 from our experimental GPR data or smaller than 0.25 for a second set of experiment (Table 1). The numerical example shows that this parameter is correlated to *n* values but also that these two parameters are not well constrained by our GPR TWT data. Taking into account the amplitudes of the reflections through a full waveform inversion might improve the estimation of this couple of parameters but it was beyond the scope of this paper.

## 5. Conclusion

We demonstrated in a numerical and experimental example that static GPR data could be used to quantify the M-vG parameters taking into account hysteresis in the SWR properties of sand. The major advantage of our experiments resides in making measurements during dynamical processes, without waiting for hydrostatic equilibrium. Moreover, using a GPR gives non invasive data integrating the whole volume of soil in between the position of the antennas and a reflector. It is complementary to *in situ* sensors such as Time-Domain Reflectometers. Our experiment could be eventually considered outdoors in the specific situation where a reflector exists in a sandy soil submitted to drainage/imbibition cycles.

In our study, the hysteresis of the water retention function is represented by the *α* parameter only. The exponent parameter λ appearing in the hydraulic conductivity function is not very well constrained by our data. Our results show that the empirical parameters of the retention curve of a sand estimated from classical laboratory experiments are not correctly fitting our hydrodynamic drainage-imbibition experiment, which we interpret as indication that static parameters do not necessarily describe correctly dynamic processes.

Further research could be developed on different type of soils, especially sands with different grain sizes, as well as hydrodynamic modeling representing the hysteresis effect with a different set of parameters (i.e. with
both $\alpha$ and $n$ van Genuchten) or using a model based on contact angle-dependent hysteresis [60]. From an
hydrogeophysical point of view, full waveform inversion and/or coupled geophysical inversion (SIP and GPR
for instance), could help to better constrain the link between the parameters used to model the hysteresis
effect.

## Acknowledgment

The authors thank Guillaume Giot and Isabelle Cousin from INRA Orléans for their measurements of
the hydraulic conductivity at saturation on sand samples. All data are available upon request from the
corresponding author.

This study is part of Kri-Terres project, supported by the French National Radioactive Waste Management
Agency (ANDRA) under the French "Investments for the Future" Program.

## References

[1] I. T. Stewart, D. R. Cayan, M. D. Dettinger, Changes toward earlier streamflow timing across western north america, Journal of climate 18 (8) (2005) 1136–1155.

[2] D. A. Wedin, D. Tilman, Species effects on nitrogen cycling: a test with perennial grasses, Oecologia 84 (4) (1990) 433–441.

[3] D. B. G. Collins, R. Bras, G. E. Tucker, Modeling the effects of vegetation-erosion coupling on landscape evolution, Journal of Geophysical Research: Earth Surface 109 (F3).

[4] M. T. van Genuchten, A closed form equation for predicting the hydraulic conductivity of unsaturated soils, Soil Sci. Soc. Am. J. 44 (5) (1980) 892–898.

[5] R. H. Brooks, A. T. Corey, Properties of porous media affecting fluid flow, Journal of the irrigation and drainage division 92 (2) (1966) 61–90.

[6] W. Haines, Studies in the physical properties of soils: 5 the hysteresis effect in capillary properties and the modes of moisture distribution associated therewith, J. Agr. Sci. 20 (1930) 91–116.

[7] T. W. Bauters, T. S. Steenhuis, J.-Y. Parlange, D. A. DiCarlo, Preferential flow in water-repellent sands, Soil Science Society of America Journal 62 (5) (1998) 1185–1190.

[8] J. Hopmans, J. Dane, Temperature dependence of soil hydraulic properties, Soil Sci. Soc. Am. J. 50 (1) (1986) 4–9.

[9] S. A. Grant, A. Salehzadeh, Calculation of temperature effects on wetting coefficients of porous solids and their capillary pressure functions, Water Resources Research 32 (2) (1996) 261–270.
[10] L. Néel, Théorie des lois d'aimantation de Lord Rayleigh,2, Cahiers Phys. 13 (1943) 18-30.

[11] D. Everett, A general approach to hysteresis, 4, Trans. Faraday Soc. 51 (1955) 1551-1557.

[12] N. Collis-George, Hysteresis in moisture content-suction relationships in soils, Proc. Nat. Acad. Sci. 24A (1955) 80-85.

[13] A. Poulovassilis, Hysteresis of pore water, an application of the concept of independent domains, Soil Sci. 93 (1962) 405-412.

[14] T. Talsma, Hysteresis in two sands and independent domain model, Water Resour. Res. 6 (3) (1970) 964-970.

[15] G. Topp, E. Miller, Hysteretic moisture characteristics and hydraulic conductivities for glass-bead media, Soil Sci. Soc. Am. J. 30 (2) (1966) 156-162.

[16] G. Vachaud, J. Thony, Hysteresis during infiltration and redistribution in a soil column at different initial water contents, Water Resour. Res. 7 (1971) 111-127.

[17] G. Topp, Soil water hysteresis measured in sandy loam compared with the hysteric domain model, Soil Sci. Sco. Am. J. 33 (1969) 645-651.

[18] G. Topp, Soil water hysteresis in silt loam and clay loam soils, Water Resour. Res. 7 (4) (1971) 914-920.

[19] Y. Mualem, Modified approach to capillary hysteresis based on a similarity hypothesis, Water Resour. Res. 9 (5) (1973) 1324-1331.

[20] Y. Mualem, A modified dependent-domain theory of hysteresis, Soil Sci. 137 (5) (1984) 283-291.

[21] J. Kool, J. Parker, Development and evaluation of closed-form expressions for hysteretic soil hydraulic properties, Water Resour. Res. 23 (1) (1987) 105-114.

[22] P. Scott, G. Farquhar, N. Kouwen, Hysteretic effects on net infiltration, Advances in infiltration, Vol. 11-83, Am. Soc. Agric. Eng., (1983), Ch. Advances in infiltration, pp. 163-170.

[23] J. H. Dane, J. W. Hopmans, Method of soil analysis, Part 4, Physical methods, Soil Science Society of America, Inc., Madison, WI., USA, (2002), Ch. Hanging water column, pp. 680-684.

[24] A. Peters, W. Durner, Simplified evaporation method for determining soil hydraulic properties, J. Hydrol. 356 (1) (2008) 147-162.

[25] S. Eching, J. Hopmans, Optimization of hydraulic functions from transient outflow and soil water pressure data, Soil Sci. Soc. Am. J. 57 (5) (1993) 1167-1175.

[26] M. Schaap, J. Feike, M. van Genuchten, ROSETTA: a computer program for estimating soil hydraulic parameters with hierarchical pedutransfert functions, Journal of Hydrology 251 (2001) 163-176.
[27] B. B. Mirus, Evaluating the importance of characterizing soil structure and horizons in parameterizing a hydrologic process model, Hydrological Processes 29 (21) (2015) 4611–4623.

[28] B. A. Ebel, K. Loague, R. I. Borja, The impacts of hysteresis on variably saturated hydrologic response and slope failure, Environmental Earth Sciences 61 (6) (2010) 1215–1225.

[29] K. Lamorski, J. Simunek, C. Slawinski, J. Lamorska, An estimation of the main wetting branch of the soil water retention curve based on its main drying branch using the machine learning method, Water Resources Research 53 (2) (2017) 1539–1552.

[30] G. Capparelli, G. Spolverino, An empirical approach for modeling hysteresis behavior of pyroclastic soils, Hydrology 7 (1) (2020) 14.

[31] J. Huisman, S. Hubbard, J. Redman, A. Annan, Measuring soil water content with Ground-Penetrating Radar., Vadose Zone J. 2 (4) (2003) 476–491.

[32] Y. Rubin, S. S. Hubbard, Hydrogeophysics, Vol. 50, Springer Science & Business Media, (2006).

[33] Q. Zhou, J. Shimada, A. Sato, Three-dimensional spatial and temporal monitoring of soil water content using electrical resistivity tomography, Water Resour. Res. 37 (2) (2001) 273–285.

[34] R. Knight, E. Grunewald, T. Irons, K. Dlubac, Y.-Q. Song, H. Bachman, B. Grau, D. Walsh, J. Abraham, J. Cannia, Field experiment provides ground truth for surface NMR measurement, Geophy. Res. Lett. 39 (3).

[35] B. Weihnacht, F. Boerner, Measurement of retention functions with hysteresis using ground-penetrating radar, Near Surface Geophysics 12 (4) (2014) 539–548.

[36] W. Lai, W. Tsang, H. Fang, D. Xiao, Experimental determination of bulk dielectric properties and porosity of porous asphalt and soils using GPR and a cyclic moisture variation technique, Geophysics 71 (4) (2006) K93–K102.

[37] R. Knight, Hysteresis in the electrical resistivity of partially saturated sandstones, Geophysics 56 (12) (1991) 2139–2147.

[38] A. Dagenbach, J. Buchner, P. Klenk, K. Roth, Identifying a parametrisation of the soil water retention curve from on-ground GPR measurements., Hydrol. Earth Syst. Sci. 17 (1) (2013) 611–618.

[39] P. Klenk, S. Jaumann, K. Roth, Quantitative high-resolution observations of soil water dynamics in a complicated architecture using time-lapse ground-penetrating radar, Hydrol. Earth Syst. Sci. 19 (2015) 1125–1139.

[40] A. Moss, X. Jing, J. Archer, Laboratory investigation of wettability and hysteresis effects on resistivity index and capillary pressure characteristics, Journal of Petroleum Science & Engineering 24 (1999) 231–242.
[41] G. Tsolfias, T. Halihan, J. S. Jr., Monitoring pumping test response in a fractured aquifer using ground-penetrating radar, Water Resour. Res. 37 (2001) 1221-1229.

[42] A. Endres, W. Clement, D. Rudolph, Ground penetraing radar imaging of an aquifer during a pumping test, Ground Water 38 (2000) 566-576.

[43] R. Brooks, A. Corey, Hydraulic properties of porous media, Hydrol. pap. 3, Colo. State Univ., Fort Collins, Colo. (1964).

[44] A. Saintenoy, J. Hopmans, Ground Penetrating Radar: Water table detection sensitivity to soil water retention properties, IEEE J. Sel. Topics Appl. Earth Observ. in Remote Sens. 4 (4) (2011) 748-753.

[45] Y. Mualem, A new model for predicting the hydraulic conductivity of unsaturated porous media., Water Resour. Res. 12 (3) (1976) 513-522.

[46] R. Haverkamp, F. Leij, C. Fuentes, A. Sciortino, P. Ross, Soil water retention: I. introduction of a shape index, Soil Sci. Soc. Am. J. 69 (6) (2005) 1881-1890.

[47] N. Burdine, Relative permeability calculations from pore-size distribution data, Journ. Petr. Tech. 5 (3) (1953) 71-77.

[48] J. Šimunek, R. Kodešová, M. Gribb, M. van Genuchten, Estimating hysteresis in the soil water retention function from cone permeameter experiments, Water Resour. Res. 35 (5) (1999) 1329-1345.

[49] A. Annan, Ground Penetrating Radar: Workshop notes, Tech. rep., Sensors and Software Inc., Ontario, Canada (1999).

[50] J. Birchak, L. Gardner, J. Hipp, J. Victor, High dielectric constant microwave probes for sensing soil moisture, Proceedings IEEE 35 (1) (1974) 85-94.

[51] L. Klein, C. Swift, An improved model for the dielectric constant of sea water at microwave frequencies, IEEE Trans. Geosci. Remote Sens. 25 (1) (1977) 104-111.

[52] A. Stogryn, Equations for calculating the dielectric constant of saline water, IEEE Trans. Microw. Theory Tech. MTT-19 (8) (1971) 733-736.

[53] J. Šimunek, M. T. van Genuchten, M. Sejna, Development and applications of the HYDRUS and STAN-MOD software packages, and related codes, Vadose Zone J. 7 (2) (2008) 587-600.

[54] J. Irving, R. Knight, Numerical modeling of Ground-Penetrating Radar in 2D using Matlab, Comput. Geosci. 32 (9) (2006) 1247-1258.

[55] C. Warren, A. Giannopoulos, I. Giannakis, gprMax: Open source software to simulate electromagnetic wave propagation for Ground Penetrating Radar, Computer Physics Communications 209 (2016) 163-170.
574 [56] E. Léger, A. Saintenoy, Y. Coquet, Hydrodynamic parameters of a sandy soil determined by Ground Penetrating Radar inside a single ring infiltrometer, Water Resour. Res. 50 (7) (2014) 5459–5474. doi:10.1002/2013WR014226.

577 [57] Q. Duan, S. Sorooshian, V. Gupta, Effective and efficient global optimization for conceptual rainfall-runoff models, Water Resour. Res. 28 (4) (1992) 1015–1031.

579 [58] J. Wösten, A. Lilly, A. Nemes, C. L. Bas, Development and use of a database of hydraulic properties of European soils, Geoderma 90 (1999) 169–185.

581 [59] J. A. Vrugt, H. V. Gupta, W. Bouten, S. Sorooshian, A shuffled complex evolution metropolis algorithm for optimization and uncertainty assessment of hydrologic model parameters, Water Resources Research 39 (8).

584 [60] A.-N. Zhou, A contact angle-dependent hysteresis model for soil-water retention behaviour, Comput. Geotech. 49 (2013) 36–42.
<table><caption>Table 1: Hydrodynamical M-vG parameters of the numerical experiments. $p_{in}$ is the set of parameters used as input to generate the radargrams presented in Figures 4-a) and 4-b). The set of best fitting parameters obtained from the SCE-UA global method is $p_{MS}^{op}$ in the multi-step case.</caption><thead><tr><th></th><th>θ<sub>r</sub><br>(cm<sup>3</sup>.cm<sup>-3</sup>)</th><th>θ<sub>s</sub><br>(cm<sup>3</sup>.cm<sup>-3</sup>)</th><th>α<sup>d</sup><br>(cm<sup>-1</sup>)</th><th>α<sup>w</sup><br>(cm<sup>-1</sup>)</th><th>n</th><th>K<sub>s</sub><br>(cm/min)</th><th>λ</th><th>RMSE<br>(ns)</th></tr></thead><tbody><tr><td>p<sub>in</sub></td><td>0.05</td><td>0.35</td><td>0.025</td><td>0.05</td><td>6</td><td>1</td><td>0.5</td><td>-</td></tr><tr><td colspan="9">multi-step - heads = [71; 42; 32; 23; 0; 23; 32; 42; 71] cm</td></tr><tr><td>p<sub>MS</sub><sup>op</sup></td><td>0.048</td><td>-</td><td>0.024</td><td>0.052</td><td>5.29</td><td>1.08</td><td>0.55</td><td>0.014</td></tr></tbody></table>

<table><caption>Table 2: Hydrodynamical M-vG parameter sets from classical laboratory measurements and best fitting inverted laboratory GPR data monitoring. * obtained from Ksat. ** according to [45].</caption><thead><tr><th></th><th>θ<sub>r</sub><br>(cm<sup>3</sup>/cm<sup>-3</sup>)</th><th>α<sup>d</sup><br>(cm<sup>-1</sup>)</th><th>α<sup>w</sup><br>(cm<sup>-1</sup>)</th><th>n</th><th>K<sub>s</sub><sup>*</sup><br>(cm/min)</th><th>λ<sup>**</sup></th><th>RMSE<br>(ns)</th></tr></thead><tbody><tr><td>Lab. meas.</td><td>0.03 ± 0.005</td><td>0.024 ± 0.002</td><td>0.043 ± 0.003</td><td>8.5 ± 0.5</td><td>[0.96; 1.5]</td><td>0.5</td><td>-</td></tr><tr><td colspan="8">multi-step - heads = [71; 41.5; 31.5; 22.5; 0; 22.5; 31.5; 41.5; 71] cm</td></tr><tr><td>1</td><td>0.052</td><td>0.033</td><td>0.094</td><td>6.6</td><td>0.82</td><td>10.4</td><td>0.14</td></tr><tr><td>2</td><td>0.04</td><td>0.023</td><td>0.047</td><td>4.4</td><td>0.54</td><td>0.25</td><td>0.09</td></tr></tbody></table>
<figure><img src="image_1.png" /><figcaption>Figure 1: Schematic representation of hysteresis for the soil water retention curve. The drainage curve is red, the wetting curve is green. Black arrows symbolize the direction of the cycle. The hysteresis is only determined by the $\alpha$ parameters, with $\alpha_w = 2\alpha_d$. See Table 1 for the full parameter set.</figcaption></figure>
<figure><img src="image_2.png" /><figcaption>Figure 2: Draft of the experiment. The tank is 80 cm high and filled with sand. The small diameter reservoir on the right represents the different heads applied at the bottom of the tank. GPR antennae are set on top of the sand reservoir.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Figure 3: TWT of the bottom tank reflection. a): The heads applied at the bottom boundary condition were 71 cm (saturated)-42-32-23-0 cm, then re-wetting following the same steps. b) The heads applied were 71 cm (saturated)-0 cm, then re-wetting to 71 cm. Blue curves represent the computed TWT using a hydrodynamic function without taking into account hysteresis effects, whereas the red circles are modeled using hysteresis. The parameters used for the hydrodynamical modeling are presented in Table 1.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 4: Radargrams generated using GprMax2D from the simulated water content profiles in the multi-step experiment. The parameters used as input for the hydrodynamic modeling are presented in Table 1.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 5: Maximum frequency of the bottom tank reflection for multi-step experiment. Circles represent the automatically picked maximum instantaneous frequency from the signal at the bottom of the tank. The red line are picked TWTs relative to the initial state.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 6: Flow diagram of the inversion procedure. "SCE-UA": Shuffled Complex Algorithm, "M-vG": Mualem-van Genuchten, "CRIM": Complex Refractive Index Method, "TWT": Two Way Travel-time, "RMSE": Root Mean Squared Error.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 7: TWT with input parameters (Table 1) in blue, red and black curves represent TWT computed with different *n* values in a), and λ values in b). Purple points are the perturbed input data adding noise sampled in a Gaussian distribution with standard deviation 0.13 ns.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 8: Histograms realized on models fitting simulated data presented on radargram 4. For each model, the RMSE is smaller than $2\sigma = 0.26$ ns). The orange vertical lines symbolize the values used as starting parameters.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 9: Objective function plots realized from noise-added simulated data inversion keeping all parameter sets giving a RMSE value smaller or equal than $2\sigma$. The color bars displayed are RMSE (ns) and are valid for all objective function plots in each column.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 10: Results of water hanging column experiments including a drying and wetting cycle. Samples were prepared with the same bulk density as the sand in the cylindrical tank experiment. Red triangles: drying data; green circles: wetting data. Plain curves represent the model fitting the data with the parameters *p*<sub>lab</sub> in Table 2.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 11: Pictures of the experiments. On the left is presented the tank filled with water, with its overflow connections. Middle picture shows the tubing linking both tanks. The picture on the right displays the sand tank with the GPR set at the top with a load on it.</figcaption></figure>
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>Figure 12: a) Radargram acquired along a multi-step experiment. The reflection coming from the bottom of the tank is arriving between 25 and 30 ns. Other reflections are coming from tank edges and the capillary fringe. b) Maximum frequency of the reflection at the bottom of the tank (circles), after moving window averaging (green plain line), TWT relative to its initial value (red line).</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 13: a) Time derived Ricker wavelet with two different maximum frequencies, b) computed time delays to be applied to measured TWT of the reflection coming from the bottom of the tank, in the multi-step case.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 14: Histograms of the models fitting the experimental multi-step data with a RMSE smaller than or equal to 0.154 ns. Orange histograms represent the laboratory data obtained with HWC by the methods described in section 4.1, orange dashed lines limit the range for $K_s$ obtained with KSAT measurement and $\lambda = 0.5$ from [45]. The pink lines represent the optimum parameters obtained from inversion of GPR measurements (see Table 2).</figcaption></figure>