# On the interest of ground penetrating radar data for the estimation of unsaturated soil parameters

Rohianuu Moua, Nolwenn Lesparre, Jean-François Girard, Benjamin Belfort, François Lehmann, Anis Younès

<figure><img src="image_1.png" /></figure>

## To cite this version:

Rohianuu Moua, Nolwenn Lesparre, Jean-François Girard, Benjamin Belfort, François Lehmann, et al.. On the interest of ground penetrating radar data for the estimation of unsaturated soil parameters. Hydrology and Earth System Sciences Discussions, 2022, pp.936. 10.5194/egusphere-2022-936. hal-04296404

HAL Id: hal-04296404

https://hal.science/hal-04296404v1

Submitted on 20 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_3.png" /></figure>

# On the interest of ground penetrating radar data for the estimation of unsaturated soil parameters

Moua Rohianuu<fnref n="1" />, Lesparre Nolwenn<fnref n="1" />, Girard Jean-François<fnref n="1" />, Belfort Benjamin<fnref n="1" />, Lehmann François<fnref n="1" />, and Younes Anis<fnref n="1" />

**Correspondence:** Moua Rohianuu (rmoua@unistra.fr)

## Abstract.

In this study, the interest of ground penetrating radar (GPR) time-lapse measurements for the estimation of hydrodynamic unsaturated soil parameters is investigated using synthetic infiltration experiments. Water movement and electromagnetic wave propagation in the unsaturated zone are modeled using a one-dimensional hydrogeophysical model. The GPR travel time data are analyzed for different reflectors: a moving reflector (the infiltration wetting front) and three fixed reflectors located at different depths in the soil. Global sensitivity analysis (GSA) is employed to assess the influence of the saturated hydraulic conductivity, the saturated and residual water contents, and the Mualem–van Genuchten shape parameters $\alpha$ and $n$ of the soil on the GPR travel time data of the reflectors. Statistical calibration of the soil parameters is then performed using the Markov chain Monte Carlo (MCMC) method. The impact of the type of reflector (moving or fixed) is then evaluated by analyzing the calibrated model parameters and their confidence intervals for different scenarios. GSA results show that the sensitivities of the moving and fixed reflectors data to the hydrodynamic soil parameters are different whereas the fixed reflectors have similar sensitivities. Results of parameter estimation show that the use of only data from the moving or fixed reflectors does not allow a good identification of all soil parameters. When both data are combined, all soil parameters can be well estimated with narrow confidence intervals.

**Keywords:** coupled hydrogeophysical model; time-lapse ground penetrating radar; unsaturated soil parameters; global sensitivity analysis; Bayesian parameter estimation; uncertainty quantification.

## 1 Introduction

The vadose zone is defined by the region between the ground surface and the groundwater table. Because of its location, it is at the center of the interactive atmospheric-surface-underground water system. Hence, understanding water flow in the vadose zone is crucial for hydrological modeling and forecasting that can be useful for water resources management, agricultural practices optimization, or geotechnical studies. The porous medium in the vadose zone is filled by both water and air phases. The air phase is considered infinitely mobile and remains at atmospheric pressure. The movement of water has a non-linear behavior and is characterized by two fundamental hydraulic relationships, namely the water retention and the hydraulic con-

<footnote n="1">¹Université de Strasbourg, CNRS, ENGEES, EOST, ITES UMR 7063, F-67000 Strasbourg, France</footnote>
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

25

conductivity functions. Various mathematical expressions can describe these functions in terms of dependent variables and fitting parameters. In this work, we use the Mualem-van Genuchten (Mualem 1976, van Genuchten 1980) hydraulic conductivity and water retention models. These models include the following unsaturated soil hydraulic parameters: the saturated hydraulic conductivity, the saturated and residual water contents, and the Mualem-van Genuchten shape parameters $\alpha$ and $n$.

30

Different approaches can be applied to estimate the unsaturated soil parameters. A typical and prevalent approach relies on laboratory measurements conducted on soil core samples. The typical measurement approaches can be of various types but common practices rely on hydraulic fluxes measurements (Vereecken et al., 2008). If this approach can be accurate at the column scale, it is prone to certain limitations when the objective is to deduce the soil parameter values at larger scales. Indeed, sample analysis through laboratory experiments is unlikely to provide parameter estimates at field conditions since the volume of the analyzed samples is often not representative of the field heterogeneity at the mesoscale (Scharnagl et al., 2011). In addition, the method is invasive and can be labor-intensive for deep or large scales investigations (Binley et al., 2015). Furthermore, the conservation of collected samples can be challenging because of issues of compaction and changes in porosity.

35

40

A more appropriate approach for large scales is based on the monitoring of underground water content. In fact, water content measurements can show significant variations because of their sensitivity to different hydrological processes. As a consequence, they are convenient for the estimation of soil parameters of the subsurface at the field scale by inverse modeling. Water content measuring techniques can be classified into two categories, based on whether they provide information on water content directly or not. The first group uses direct measurements conducted at the point scale with water content sensing techniques, for instance using thermal or electromagnetic sensors (e.g., capacitance or time domain reflectometry, Jones et al., 2005; Belfort et al., 2019). These techniques can yield data with great resolution at one location and give information on the dynamics at the field scale (Vereecken et al., 2008). In addition, measurements taken at various locations can help to describe the distribution of water content, and thus, allow a good characterization of the state of the soil. However, the installation of sensors is often laborious, time-consuming, and destructive (Huisman et al., 2003; Dal Bo et al., 2019). Furthermore, their reliability requires an accurate calibration (Robinson et al., 2008).

45

The second category uses indirect measurements of the water content, like remote sensing and hydrogeophysical methods with non-invasive devices. Remote sensing techniques use devices that are not in direct contact with the ground, such as unmanned aerial vehicles thermal infrared imagery (Zhang et al., 2019) or airborne ground penetrating radar (Edemsky et al., 2021). These methods provide the mapping of water content at a large scale and in locations where direct sensing measurements cannot be conducted. However, remote sensing methods exhibit a penetration depth of only a few centimeters and are often limited by the vegetation density (Vereecken et al., 2008; Robinson et al., 2008).

50

Common hydrogeophysical methods include electromagnetic induction (Doolittle and Brevik, 2014), direct current resistivity (de Jong et al., 2020), nuclear magnetic resonance (Costabel and Günther, 2014), and ground penetrating radar (GPR) (Huisman et al., 2003; Klotzsche et al., 2018) methods. These techniques supply indirect information on hydraulic properties or states, at various scales, from estimated geophysical properties. As mentioned by Binley et al. (2015), such conversion from

55
<figure><img src="image_6.png" /></figure>

geophysical to hydraulic properties or states requires the use of robust petrophysical relationships to provide reliable estimates of hydraulic parameters.

Nowadays, GPR is highly used in the field of hydrogeophysics. Different techniques have been reviewed and discussed by Huisman et al. (2003) and Klotzsche et al. (2018). GPR is highly sensitive to water content, and, as such, it can close the gap between the spatial scales covered by direct and remote sensing techniques. Furthermore, the temporal variability of the soil water content can be characterized from time-lapse measurements. For instance, GPR data can be collected during artificial experiments (e.g., infiltration, runoff, drainage, imbibition) that can provide interesting information on the flow characteristics (e.g., Saintenoy et al., 2008; Moysey, 2010; Scholer et al., 2011; Busch et al., 2013; Jonard et al., 2015; Jaumann and Roth, 2018; Léger et al., 2014; 2016; 2020). Note however that the hydraulic properties estimated from GPR data are subject to an inherent compromise between a deep investigation and a fine spatial resolution. For instance, the lowest frequencies (typically from 1 GHz down to 100 MHz) allow deeper penetrations (until a maximum depth between 1 m and 3 m in most organic media).

Saintenoy and Hopmans (2011) performed a sensitivity analysis to evaluate the importance of the soil parameters on the reflections caused by the water table. They found that the detectability of the GPR reflections is mainly dependent on the slope of the retention curve in the capillary zone. Moysey (2010) showed that the Mualem-van Genuchten parameter *n* is the most poorly constrained among all unsaturated soil parameters when considering the two-way travel time (TWT) from various sources of reflection for a laboratory infiltration experiment conducted in a sandbox. At the laboratory scale, Léger et al. (2020) have monitored imbibition-drainage experiments using a single-offset surface GPR. At the field scale, infiltration experiments have been monitored with a single-offset surface GPR (Léger et al., 2014; 2016) or off-ground GPR (Jadoon et al., 2012). Saito et al. (2018) used a more complex multi-offset and multi-channel surface GPR to directly monitor the wetting front progression. Mono-channel multi-offset technique is usually not suited for monitoring experiments with high temporal variability, as the offset must be adjusted between each measurement. The multi-channel technique has the advantage to be multi-offset and is, therefore, able to simultaneously determine the propagation speed and the depths of reflectors.

Time-lapse GPR monitoring of artificial infiltrations is one of the cheapest non-destructive methods that can be easily conducted to estimate the unsaturated soil parameters in field conditions. Léger et al. (2014) have demonstrated the relevance of such a methodology to evaluate the hydraulic parameters of sandy soil. The authors have investigated synthetic and field examples and showed that the inverted parameters were in agreement with the values obtained in the laboratory for soil samples and with disk infiltrometer measurements. However, in their study, Léger et al. (2014) didn't assess the reliability of the estimated values since the uncertainty associated with the calibrated parameters has not been evaluated. Furthermore, Léger et al. (2014) employed only the TWT data obtained from the GPR reflection on the wetting front for the calibration of the soil parameters. The present study extends the work of Léger et al. (2014) by:

- Considering different reflectors at different depths: a moving reflector which corresponds to the infiltration dynamic wetting front and two fixed reflectors located at different depths in the soil.
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

- Investigating the influence of all soil parameters (the saturated hydraulic conductivity, the saturated and residual water contents, and the Mualem-van Genuchten shape parameters $\alpha$ and $n$) on the GPR TWT data of the three reflectors using Global Sensitivity Analysis (GSA).

- Performing statistical calibration of soil parameters using the Markov Chain Monte Carlo (MCMC) method and evaluating the reliability of the estimated parameters by analyzing not only the calibrated model parameters but also their associated uncertainty.

- Evaluating the impact of the type of reflector (moving or fixed) by analyzing the calibrated model parameters and their confidence intervals for different scenarios.

The plan of the paper is as follows: Section 2 describes the test case as well as the mathematical and numerical hydrogeophysical models. Section 3 reports on the GSA results of the different TWT signals. Then, Section 4 discusses the results of soil parameter estimation with MCMC for different scenarios.

# 2 Test case description and numerical solution

## 2.1 Test case description

The test case considered in this work is a hypothetical one-dimensional experiment of water infiltration in a homogeneous sandy soil of 150 cm (Fig.1a). The infiltration is driven by a constant pressure head of 10 cm applied at the surface of the soil (i.e., a 10 cm water ponding condition is maintained at the top). The medium is initially at the hydrostatic equilibrium with a water table maintained at 100 cm below the soil surface (Fig.1b). The domain is initially formed by an unsaturated zone of 100 cm thick above a saturated zone of 50 cm thick. We assume the experiment to be monitored with a surface GPR. The propagating time (i.e., the TWT) of the GPR waves reflected by two types of reflectors are considered (Fig.1c): (i) the moving infiltration wetting front and (ii) two fixed reflectors corresponding to a local heterogeneity at two different depths. For instance, these can be small objects that are artificially buried (e.g., moisture sensing probes) or naturally embedded (such as small rocks) in the porous medium. The fixed reflectors are supposed to be small enough compared to the section of the infiltrated area, so they do not significantly perturb the vertical flow. The upper fixed reflector, R50, is located in the initially unsaturated zone at 50 cm depth. The reflector R120 is located in the saturated zone, under the water table, at a distance of 120 cm from the soil surface (Fig.1a). In the following, the time-lapse TWT signal for reflection caused by the infiltration wetting front is noted TWT$_f$ and that from the two immovable diffracting points R50 and R120, are respectively noted TWT$_{50}$ and TWT$_{120}$.
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 1.** Test case and experimental device illustration at an advanced time step (a). R50 and R120 are fixed reflectors considered in this experiment. TX and RX refer to the transmitter and receiver antennas of the GPR system. Effective saturation $S_e$ (b) and reflection coefficient $r$ (c) profiles with depth.</figcaption></figure>

## 2.2 The mathematical model

### 2.2.1 Unsaturated flow model

Water infiltration in unsaturated/saturated soils is governed by the one-dimensional Richards' equation (Richards, 1931):

$$ \frac{\partial \theta}{\partial t} = \frac{\partial}{\partial z} \left[ K(\theta) \left( \frac{\partial h}{\partial z} - 1 \right) \right] \quad (1) $$

where $h$ (cm) is the pressure head; $z$ is the depth (cm), taken positive in the downward direction; $t$ is the time (s), $\theta$ (cm³/cm³) is the actual water content, and $K(\theta)$ (cm/s) is the hydraulic conductivity which is a function of water content. The initial condition is a hydrostatic pressure distribution corresponding to a water table at 100 cm depth. The boundary condition at the top of the domain is a fixed Dirichlet condition of 10 cm maintained during the experiment. The boundary condition at the bottom is a piezometric head fixed at -100 cm which corresponds to the water table position (Fig.1).
<figure><img src="image_11.png" /></figure>

The interdependencies of the pressure head, conductivity, and water content are described using the standard models of Mualem (1976) and van Genuchten (1980):

$$130 \qquad S_e(h) = \frac{\theta(h) - \theta_r}{\theta_s - \theta_r} = \begin{cases} [1 + (\alpha|h|)^n]^{-m} & \text{if } h < 0 \\ 1 & \text{if } h \ge 0 \end{cases} \qquad (2)$$

$$135 \qquad K(h) = \begin{cases} K_s \times S_e(h)^{1/2} [1 - (1 - S_e(h)^{1/m})^m]^2 & \text{if } h < 0 \\ K_s & \text{if } h \ge 0 \end{cases} \qquad (3)$$

where $S_e(h)$ (-) is the effective saturation, $\theta_s$ and $\theta_r$ (cm³/cm³) are the saturated and residual water contents respectively, $K_s$ (cm/s) is the saturated conductivity, $m = 1 - 1/n$, $\alpha$ (1/cm) and $n$ (-) are the Mualem-van Genuchten shape parameters.

## 2.2.2 Petrophysical and Geophysical relationships

In GPR sounding, pulses of radiofrequency (MHz to GHz) electromagnetic waves are emitted from a transmitting antenna through the sounded medium. The electromagnetic response is then acquired with a receiving antenna. With a surface GPR, both antennas are installed at the surface of the soil Fig.1). To monitor the experiment of water infiltration with time-lapse GPR, the sounding system is set immobile above the infiltration zone in order to capture the time variation of the electromagnetic response due to the change of saturation.

To describe the dependency of the dielectric permittivity on the water content, we use the complex refractive index model (Birchak et al., 1974). This petrophysical relationship relates the dielectric constant $\epsilon$ (-) of a three-phase (water-solid-air) medium to its water content by:

$$145 \qquad \sqrt{\epsilon(z,t)} = \theta(z,t)\sqrt{\epsilon_w} + (1 - \phi)\sqrt{\epsilon_s} + [\phi - \theta(z,t)]\sqrt{\epsilon_a} \qquad (4)$$

where $\phi$ (-) is the porosity, $\epsilon_w = 80$, $\epsilon_s = 2.5$ and $\epsilon_a = 1$ are the relative dielectric constants of water, silica (sand) and air, respectively.

In a low electrical conductivity and non-magnetic soil, the electromagnetic waves propagate at a speed $V$ (cm/ns) (Annan, 2003):

$$150 \qquad V = \frac{c}{\sqrt{\epsilon}} \qquad (5)$$

where $c \approx 30$ cm/ns is the speed of electromagnetic waves in air, and $\epsilon$ (-) is the dielectric constant of the porous medium. Equations (4) and (5) evidence that GPR waves propagate at a much lower speed in wet conditions. Any source of reflection in the sounded soil produces a reflected wave that is recorded at a time corresponding to the duration of its propagation, from the transmitting antenna, down to the source of reflection, then back up to the receiving antenna, i.e., the TWT of the reflected
<figure><img src="image_12.png" /></figure>

wave.

We consider a one-dimensional scenario (the offset between the antennas is null) and discretize the domain into N cells i, centered at a depth z_i, with element boundaries at z_{i-1/2} and z_{i+1/2}. The TWT for the reflection occurring at the interface (i - 1/2) between the elements i - 1 and i can be expressed as the sum of the vertical TWT in each element above i:

$$160 \quad \text{TWT}(z_{i-1/2}) = 2 \sum_{j=1}^{i-1} \frac{|l_j|}{V_j} \qquad (6)$$

in which $|l_j|$ (cm) is the length of the element j above i and $V_j$ (cm/ns) is the GPR propagation speed in the element j.

A reflection occurs at the interface between two successive elements if the reflection coefficient is not zero. The reflection coefficient expresses the contrast of dielectric constant (due to the contrast of water content) at the interface between the two elements i - 1 and i. When the offset between transmitting and receiving antennas is null, the reflection coefficient at interface (i - 1/2) is defined by:

$$165 \quad r(z_{i-1/2}) = \frac{\epsilon(z_i) - \epsilon(z_{i-1})}{\epsilon(z_i) + \epsilon(z_{i-1})} \qquad (7)$$

where $\epsilon(z_i)$ is the dielectric constant of the element i.

For an 800 MHz antenna, the wavelength can typically vary from 6 cm in a wet medium to around 18 cm in a dry medium. The abrupt change in the reflection coefficient makes the wetting front easily detectable whatever the hydraulic parameters, contrarily to the water table which may be hidden due to the capillary fringe (Bano, 2006; Saintenoy and Hopmans, 2011).

## 2.3 The numerical model

The variation of the water content in the soil during the infiltration is computed using the WAMOS-1D code (Belfort et al., 2018). The model describes the water movement in the porous medium using Richards' equation (1), and the constitutive relationships between the pressure, the hydraulic conductivity, and the volumetric water content given by Eq. (2) and Eq. (3).

The domain of 150 cm depth is discretized with uniform elements of 1 cm thick with homogeneous properties. The WAMOS-1D code solves the system of Eqs. (1)-(3) and yields the vertical distribution of water content at each time step. This distribution is then converted into a vertical dielectric permittivity profile $\epsilon$ using the petrophysical relationship Eq. (4) and into a GPR wave propagation speed profile $V$ using Eq. (5). Then, the time-lapse TWT signals for the fixed objects, TWT$_{50}$ and TWT$_{120}$, are calculated at each time step using Eq. (6) (dashed and dotted curves in Fig.2).

The time-lapse signal TWT$_{f}$, induced by wave reflection on the wetting front because of the sharp water content variation at the front position is calculated in two steps. First, we search the wetting front position $z_{i-1/2}^*$, which corresponds to the interface position having the maximum reflection coefficient from Eq. (7) as illustrated in Fig.1. Then, the TWT signal of the wetting front is obtained using TWT$_{f}$ = TWT($z_{i-1/2}^*$) from Eq. (6) (solid curve Fig.2).

Note that TWT$_{50}$ and TWT$_{120}$ signals are induced by fixed objects, thus, these signals exist regardless of the position of the infiltration front. On the other hand, TWT$_{f}$ is induced by the infiltration wetting front whose position varies over time. Besides, contrarily to TWT$_{50}$, and TWT$_{120}$, the TWT$_{f}$ signal disappears when the wetting front reaches the water table. To
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 2.** Hydrogeophysical model responses for $K_s = 0.08$ cm/s, $\theta_s = 0.4$ (-), $\theta_r = 0.07$ (-), $\alpha = 0.145$ cm$^{-1}$, $n = 2.68$ (-). TWT$_f$ corresponds to the TWT signal for the wetting front, while TWT$_{50}$ and TWT$_{120}$ are the TWT signals for fixed objects located at 50 and 120 cm below the surface, respectively.</figcaption></figure>

avoid numerical issues when simulations are performed with different soil parameter sets, the value of TWT$_f$ when the water table is reached, is artificially maintained for the remaining time steps until the end of simulation time. The water table is assumed to be reached when the maximum reflection coefficient of Eq. (7) is under a threshold of $10^{-2}$. This reflects a fully saturated domain with an almost uniform water content distribution (solid curve Fig.2). An explanation of the computation of all TWT signals is summarized in Fig.3.

# 3 Global sensitivity analysis of TWT signals

## 3.1 GSA method

The GSA method evaluates how the outputs of a model are influenced by the variation of the input parameters (Mara and Tarantola, 2008). Among the various forms of GSA, a variance-based sensitivity analysis, allowing the calculation of Sobol sensitivity indices (Sobol', 2001) is employed. Such indices depict the contribution of the variation of any input variable $x$ to the total variance of an output variable $y$. In our case, the input variables are the unsaturated soil parameters ($K_s, \theta_s, \theta_r, \alpha, n$) and the output variables are the TWT signals (TWT$_f$, TWT$_{50}$, TWT$_{120}$).

Given a model with a set of $p$ independent random parameters $\mathbf{X} = \{x_1, x_2, ..., x_p\}$ that yields a random response $y(\mathbf{X})$, the two variance-based sensitivity measures, also called Sobol indices (Sobol', 2001) are:
<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>Figure 3. Summary of the working process of the forward hydrogeophysical model and how it is used to build the PCE surrogate model.</figcaption></figure>

- the first-order sensitivity index:

$$S_i = \frac{\mathrm{Var}[E[y(\mathbf{X})|x_i]]}{\mathrm{Var}[y(\mathbf{X})]} \in [0, 1] \qquad (8)$$

- the total sensitivity index:

$$ST_i = \frac{E[\mathrm{Var}[y(\mathbf{X})|x_{-i}]]}{\mathrm{Var}[y(\mathbf{X})]} \in [0, 1] \qquad (9)$$

where $x_{-i} = \mathbf{X} \setminus x_i$ is the set of all parameters except $x_i$, $E(\cdot)$ and $E(\cdot|\cdot)$ are respectively the expectation and the conditional expectation operators, $\mathrm{Var}(\cdot)$ and $\mathrm{Var}(\cdot|\cdot)$ are the variance and the conditional variance, respectively. The first-order index $S_i$ quantifies the contribution of the parameter $x_i$ alone to the total variance of $y(\mathbf{X})$, while $ST_i$ also includes all interactions of $x_i$ with the other parameters $x_{-i}$.

To perform a variance-based GSA, a practical approach (to save computational time) is to use Polynomial Chaos Expansion (PCE; Wiener, 1938). The PCE approach consists in developing any signal $y(\mathbf{X})$ as a set of orthonormal multivariate polynomials of a degree not exceeding $D$:

$$y(\mathbf{X}) = \sum_{|\beta| \le D} s_{\beta} \Psi_{\beta}(\mathbf{X}) \qquad (10)$$
<figure><img src="image_17.png" /></figure>

<table><caption>**Table 1.** Prior intervals of the unsaturated soil parameters for both GSA and Bayesian estimation.</caption><thead><tr><th></th><th>K<sub>s</sub> (cm/s)</th><th>θ<sub>s</sub> (cm<sup>3</sup>/cm<sup>3</sup>)</th><th>θ<sub>r</sub> (cm<sup>3</sup>/cm<sup>3</sup>)</th><th>α (1/cm)</th><th>n (-)</th></tr></thead><tbody><tr><td>[x<sub>min</sub> - x<sub>max</sub>]</td><td>[0.001 - 0.15]</td><td>[0.32 - 0.48]</td><td>[0.01 - 0.13]</td><td>[0.01 - 0.28]</td><td>[1.5 - 10]</td></tr></tbody></table>

<table><caption>**Table 2.** Variance of TWT<sub>f</sub>, TWT<sub>50</sub> and TWT<sub>120</sub> signals at t = 50 s, 150 s and 2000 s calculated with the PCE surrogate model and with the hydrogeophysical model.</caption><thead><tr><th></th><th>TWT<sub>f</sub></th><th>TWT<sub>50</sub></th><th>TWT<sub>120</sub></th></tr></thead><tbody><tr><td colspan="4">t = 50 s</td></tr><tr><td>Var_HYD_model</td><td>14</td><td>5</td><td>15.9</td></tr><tr><td>Var_PCE_model</td><td>13.4</td><td>5</td><td>15.8</td></tr><tr><td>Var_error</td><td>4.3%</td><td>0.9%</td><td>0.5%</td></tr><tr><td colspan="4">t = 150 s</td></tr><tr><td>Var_hyd_model</td><td>54.6</td><td>5.5</td><td>23.3</td></tr><tr><td>Var_PCE_model</td><td>53.9</td><td>5.4</td><td>23.2</td></tr><tr><td>Var_error</td><td>1.3%</td><td>1.4%</td><td>0.4%</td></tr><tr><td colspan="4">t = 2000 s</td></tr><tr><td>Var_hyd_model</td><td>28.8</td><td>1.5</td><td>9.9</td></tr><tr><td>Var_PCE_model</td><td>27.1</td><td>1.4</td><td>9.3</td></tr><tr><td>Var_error</td><td>5.7%</td><td>7.4%</td><td>5.2%</td></tr></tbody></table>

where $β = β_1, β_2, ..., β_p \in R^p$ is a p-th-dimensional index, $s_β$ are the PC coefficients, $Ψ_β$ are the generalized polynomial chaos of degree $|β| = \sum_{i=1}^p β_i$.

In this work, Legendre polynomials are used since uniform distributions are assumed for all uncertain parameters. Uniform distributions express the absence of prior information. This makes all parameter values in the given prior intervals equally likely. Large prior distribution intervals are considered for all unsaturated soil parameters (Table 1).

The number of coefficients for a full PCE representation is $P = (p+D)!/p!D!$. A training dataset of M realizations of the forward coupled hydrogeophysical model is used to build the PCE surrogate model of order D (Fajraoui et al, 2011; Shao et al., 2017; Younes et al., 2013). The coefficients of the PCE are obtained by searching the best fit (in the least square sense) of the PCE surrogate model to the hydrogeophysical model for the M realizations. To work with low-discrepancy sets, the M realizations correspond to sets of input parameters sampled from their prior probability distributions, using quasi-random Sobol sequences (Shao et al., 2017). We illustrate the principle of the construction of the PCE with our hydrogeophysical model in Fig.3.
<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>

225

A PCE is constructed at each time step for all model responses (TWT<sub>f</sub>, TWT<sub>50</sub>, and TWT<sub>120</sub>) since we deal with transient simulations. In this work, M = 2048 hydrogeophysical model realizations are employed to obtain PCEs of degrees D = 5 containing P = 252 coefficients. The obtained PCEs are sufficiently accurate as the variance of the TWT output signals is calculated with the surrogate PCE model and the forward hydrogeophysical model at three different times t = 50 s, 150 s, and 2000 s. The results of Table 2 show that the relative difference between the two variances is very small for all investigated times. Note that although the relative variance error for the TWT<sub>50</sub> at t = 2000 s is the largest (around 7%), it remains insignificant since the total variance of the signal at this time is negligible (less than 2 ns<sup>2</sup>). The variance of the forward hydrogeophysical model is therefore well reproduced by the PCE surrogate model which will be employed for the GSA of the TWT signals using the variance decomposition.

230

## 3.2 GSA results

235

The temporal distribution of the output variance of the three TWT signals (TWT<sub>f</sub>, TWT<sub>50</sub> and TWT<sub>120</sub>) are represented in Fig.4. For each TWT signal, the variance is represented by the black curve and the relative contributions of the uncertain parameters to the variance are represented by the shaded area. The blank region between the variance curves and the shaded area represents interactions between parameters.

240

TWT<sub>f</sub> has a different behavior from the TWT signals of fixed reflectors TWT<sub>50</sub> and TWT<sub>120</sub> (Fig.4). Although the TWT signals of fixed reflectors have different variance magnitudes, they exhibit similar behavior (Fig.4b and 4c). The variance of the TWT signal is five times more significant for TWT<sub>120</sub> than for TWT<sub>50</sub>. This is in agreement with the physics since the zone of the porous medium affecting the GPR wave is more important for the TWT<sub>120</sub> signal than for the TWT<sub>50</sub>. In addition, the period of influence of the unsaturated parameters (θ<sub>r</sub>, α, n) is also more important for TWT<sub>120</sub> than for TWT<sub>50</sub> since saturated conditions for the reflector R120 are reached much later than for R50. Since fixed reflectors exhibit similar behavior, in the following, we comment on the results of TWT<sub>f</sub> and TWT<sub>120</sub> signals.

245

### 3.2.1 GSA of the TWT<sub>f</sub> signal

250

TWT<sub>f</sub> variance is zero at the beginning of the infiltration (Fig.4a) which means that the TWT<sub>f</sub> signal is not affected by the initial conditions. Indeed, the infiltration wetting front and the TWT<sub>f</sub> signal start at zero for all parameter sets. Then, the variance of the signal increases until a maximum of 60 ns<sup>2</sup>, reached at around 3 min. After that, the variance decreases, but keeps a significant value of around 25 ns<sup>2</sup> (Fig.4a). Concerning parameter sensitivities, at the beginning, the TWT<sub>f</sub> signal is mainly affected by K<sub>s</sub>. The influence of this parameter decreases over time and reaches zero for long times when steady-state conditions (corresponding to a fully saturated soil) are reached. The parameter θ<sub>s</sub> has a moderate influence on the TWT<sub>f</sub> signal. Its influence is not observable at short times since unsaturated conditions occur. Overall, the most influential parameter on the TWT<sub>f</sub> signal is the van Genuchten parameter α. This parameter seems influential even for saturated conditions. Note that this numerical artifact is observed because the value of TWT<sub>f</sub> is artificially maintained when the infiltration wetting front reaches the water table, while physically the signal disappears. The effects of the parameters θ<sub>r</sub> and n are not observable (Fig.4a). The blank region between the variance curve and the shaded area in this figure is due to the interaction between the parameters.

255
<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 4.** Time distribution of the variance of TWT<sub>f</sub> (a), TWT<sub>50</sub> (b) and TWT<sub>120</sub> (c). The shaded area under the variance curve represents the partial marginal contributions of the uncertain parameters; the blank region between the shaded area and the variance curves represents the contribution of interactions between the parameters. The marginal effects shown in Fig. 6 are represented at three time steps *t*₁ = 1 min, *t*₂ = 5 min, and *t*₃ = 200 min, highlighted here (dotted black lines).</figcaption></figure>

To estimate this interaction, we plot the difference between the total ($ST_i$) and the first order ($S_i$) sensitivity indices for all parameters (Fig.5a). This difference reflects the interaction between the parameters over time. Interactions between parameters are negligible for all parameters ($ST_i \approx S_i$), except for $K_s$ and $\alpha$ (Fig.5a). Hence, the interaction between these two parameters affects the variance of the TWT<sub>f</sub> signal as represented by the blank region between the variance curve and the shaded area (Fig.4a).
<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 5.** Difference between the total ($ST_i$) and the first order ($S_i$) sensitivity indices for all parameters for the TWT$_f$ (a) and the TWT$_{120}$ (b) signals.</figcaption></figure>

To evaluate further the effect of the unsaturated soil parameters on the TWT$_f$, we plot the marginal effect of each parameter (Fig.6). The marginal effect can be easily derived from the PCE coefficients and reflects the effect of one parameter on the output signal. Fig.6 depicts the marginal effects of each hydraulic parameter, i.e. their influence on the TWT signals as a function of their value when they vary over the range of their prior distribution interval, while the other hydraulic parameters are kept fixed at their center value. This representation allows determining the regions of influence of the hydraulic parameters, given that the stronger the slope of the marginal effect curve, the higher the influence of the parameter. These marginal effects can vary over time, so we represent them at the three time steps ($t_1$ = 1 min, $t_2$ = 5 min, and $t_3$ = 200 min) highlighted with dotted vertical black lines in Fig.4. The oscillations are caused by numerical artifacts related to the degree of the polynomials used in the PCE model. From Fig.6a, it can be noticed that:

265

- $K_s$ is highly influential at the beginning of the experiment. At $t_1$ = 1 min, the TWT$_f$ signal varies almost linearly with $K_s$. Indeed, at the beginning of the experiment, when $K_s$ increases, the wetting front is more advanced, thus, the GPR wave propagates at a lower speed and the TWT$_f$ signal increases. At $t_2$ = 5 min, the TWT$_f$ signal is sensitive only for small $K_s$ values. Indeed, for high $K_s$ values, the soil is fully saturated and the perturbation of the high value of $K_s$ doesn't change the TWT$_f$ signal. At $t_3$ = 200 min, the soil is fully saturated for almost all $K_s$ values, thus, the TWT$_f$ signal becomes insensitive to $K_s$.

275

- $\theta_s$ has no influence at the first times ($t_1$ = 1 min) since unsaturated conditions occur. For long times, the TWT$_f$ signal is very sensitive to $\theta_s$ with an almost linear behavior. Indeed, when the soil is fully saturated, the dielectric permittivity and thus the TWT$_f$ signal is almost proportional to $\theta_s$.

280

- the sensitivity of TWT$_f$ to $\theta_r$ is moderate and can be observed only at the beginning of the experiment (unsaturated conditions) with an almost linear behavior observable at $t$ = 1 min and 5 min. The positive slope of the curve is consistent
<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /><figcaption>**Figure 6.** Marginal effects of the unsaturated soil parameters $K_s$, $\theta_s$, $\theta_r$, $\alpha$ and $n$ on the TWT$_f$ and TWT$_{120}$ signals at three different times $t_1 = 1$ min, $t_2 = 5$ min, and $t_3 = 200$ min, highlighted in Fig. 4.</figcaption></figure>

with the physics of the process (when $\theta_r$ increases, the speed of the electromagnetic wave decreases, and the TWT$_f$ signal increases).

- The van Genuchten parameter $\alpha$ is highly influential notably for long times ($t_3 = 200$ min). A small variation of the parameter $\alpha$ can induce a strong variation of the TWT$_f$ signal. Notably, the sensitivity of $\alpha$ is very high for $\alpha \le 0.05$ cm$^{-1}$.

- The sensitivity of TWT$_f$ to the parameter $n$ is almost zero (flat curves) at all times ($t = 1$ min, 5 min, and 200 min). The parameter $n$ has therefore a negligible effect on the TWT$_f$ signal and, as a consequence, it is expected to be poorly identifiable from the TWT$_f$ data.

### 3.2.2 GSA of the TWT$_{120}$ signal

The variance of the TWT$_{120}$ signal is nonzero at the beginning of the experiment which means that the TWT$_{120}$ signal is affected by the initial conditions (Fig. 4c). Indeed, at the very beginning, the pressure distribution is hydrostatic and the water content distribution in the column is obtained from Eq. 2 which depends on all soil parameters except $K_s$. Therefore, the speed of the GPR wave depends on the initial water content distribution which is dependent on the unsaturated soil parameters $\theta_s$,
<figure><img src="image_26.png" /></figure>

$\theta_r$, $\alpha$, and $n$. The most influential parameter at the beginning of the experiment is the parameter $\alpha$. Over time, the effect of this parameter reduces, whereas the effect of $\theta_s$ increases. For long times, $\theta_s$ becomes the only sensitive parameter. The parameter $K_s$ is also very sensitive. Its effect starts at zero, and increases until a maximum is reached at around 3 min, then it slowly decreases and becomes negligible after 100 min. As with the TWT$_f$ signal, interactions between parameters are moderate. The difference between the total and first-order Sobol indices is negligible for all parameters except after 1 min for the parameters $K_s$, $\alpha$ and $\theta_s$ (see Fig.5b). This interaction corresponds to the blank region, between the variance curve and the shaded area in Fig.4c. The marginal effects of the soil parameters on the TWT$_{120}$ signal are plotted in Fig.6b for $t = 1$ min, 5 min, and 200 min. The curves in this figure show that:

- As for the TWT$_f$ signal, $K_s$ is highly sensitive, especially for $t = 1$ min and 5 min.

- The saturated water content $\theta_s$ is very influential for all times. The TWT$_{120}$ varies almost linearly with $\theta_s$ even at the beginning ($t_1 = 1$ min), since the fixed reflector is located in the lower saturated region.

- As for the TWT$_f$ signal, $\theta_r$ is sensitive only at the beginning of the experiment (unsaturated conditions) with an almost linear behavior at $t = 1$ min and 5 min. When $\theta_r$ increases, the water content increases, and hence, the TWT$_{120}$ increases.

- The van Genuchten parameter $\alpha$ is highly sensitive. However, contrarily to the TWT$_f$ signal where $\alpha$ is highly sensitive at long times ($t_3 = 200$ min), the sensitivity of $\alpha$ for the TWT$_{120}$ signal is high at short times ($t_1 = 1$ min). For long times, the influence of $\alpha$ disappears since the soil becomes fully saturated. The negative slope of the curve of the TWT$_{120}$ signal as a function of $\alpha$ observed at the beginning of the experiment is consistent with the physics of the process. Indeed, when $\alpha$ increases, the capillary fringe thickness decreases, hence, the water content in the unsaturated zone decreases, and thus the TWT$_{120}$ signal decreases.

- Surprisingly, and contrarily to the TWT$_f$ signal which showed a flat curve for the marginal effect of the parameter $n$ for all parameter values and at all investigated times, the TWT$_{120}$ signal is sensitive to $n$ at the beginning of the experiment ($t_1 = 1$ min) with a high sensitivity for $n < 3.5$ and a moderate sensitivity (the curve has a small slope) for $n \ge 3.5$.

# 4 Bayesian soil parameter estimation from the TWT signals

In this section, we estimate the unsaturated soil parameters in a Bayesian framework using the Markov Chain Monte Carlo (MCMC) sampler (Vrugt and Bouten, 2002; Vrugt et al., 2008). The statistical calibration is performed for a GPR monitored infiltration experiment in order to address the following questions:

1. Can we obtain an appropriate estimation of all unsaturated soil parameters from TWT data?

2. What is the impact of the kind of TWT data (moving/fixed reflectors) and of the number of reflectors on the calibrated model parameters and their confidence intervals?

3. What is the optimal set of TWT measurements that yields good reliability of all unsaturated soil parameters?
<figure><img src="image_27.png" /></figure>

The MCMC method has been successfully employed in various inverse hydrological problems (e.g., Fajraoui et al., 2011; Younes et al., 2016; Younes et al., 2017; Younes et al., 2018). The method generates random sequences of parameter sets that asymptotically converge toward the target joint posterior distribution by searching the ensemble of possible parameter sets that satisfactorily fit the observations. The converged sets can then be used to assess the quality of the parameter estimation such as the optimal parameter values and the 95% Confidence Intervals (CIs) which allow for evaluating the reliability of the parameters via uncertainty quantification.

330

In the sequel, the MCMC method is performed with the DREAM(ZS) software (Laloy and Vrugt, 2012), which is an efficient MCMC sampler. The vector of unknowns is formed by the five unsaturated soil parameters ($K_s$, $\theta_s$, $\theta_r$, $\alpha$, $n$). A reference solution is generated by simulating the hydrogeophysical problem formed by the system of equations (1)-(6) using the following reference parameter values $K_s^* = 0.08$ cm/s, $\theta_s^* = 0.4$, $\theta_r^* = 0.07$, $\alpha^* = 0.145$ cm$^{-1}$, $n^* = 2.68$, as shown in Table 3. The TWT$_f$, TWT$_{50}$, and TWT$_{120}$ signals used as calibration data are deduced from the results of the simulation using the reference parameter values. These TWT signals are then independently corrupted using a normally distributed noise with a standard deviation $\sigma = 0.5$ ns.

335

The TWT$_f$, TWT$_{50}$ and TWT$_{120}$ calibration signals, illustrated before noise corruption in Fig. 2, increase almost linearly until reaching a plateau. For the TWT$_{50}$ signal, the plateau is reached when the infiltration front attains the R50 reflector and the value of the plateau corresponds to the time needed by the electromagnetic wave to make a round trip from the surface to a 50 cm depth of a full saturated porous medium. For the TWT$_{120}$ signal, the plateau signal is reached when the infiltration front attains the water table (the domain becomes fully saturated) and the value of the plateau corresponds to the time needed by the electromagnetic wave to make a round trip from the surface to a 120 cm deep of a fully saturated porous medium. For TWT$_f$, the plateau value is also reached when the infiltration front attains the water table and the value of the plateau corresponds to the time needed by the electromagnetic wave to make a round trip from the surface to the water table at 100 cm deep.

340

345

The reliability of the unsaturated soil parameters is assessed for 5 different scenarios of measurement sets. In the first scenario, only data of the wetting-front TWT$_f$ signal are used for the calibration. The second and third scenarios use only the TWT$_{50}$ and TWT$_{120}$ signal, respectively, obtained from reflection on the fixed reflector R50 and R120. The fourth scenario uses both data of TWT$_f$ and TWT$_{120}$ as fitting data. The last scenario investigates the benefit of adding a fixed reflector by using data of the TWT$_f$, TWT$_{50}$ and TWT$_{120}$ signals as conditioning information.

350

In the five scenarios, the MCMC sampler uses three parallel chains and a total number of 50000 runs. The last 25% of the runs that adequately fit the model onto observations are used to estimate the joint posterior distribution.

The MCMC results of the five studied scenarios are given in Table 3 which depicts, for each parameter, the mean estimated value, its posterior CI size, and the ratio of prior to posterior intervals. Note that the CI and the last indicator are calculated from the standard deviation by assuming a Gaussian posterior distribution. A small CI indicates an accurate estimation of the parameter. A significant difference between the prior and posterior intervals is a sign of the high sensitivity of the model to that parameter (Dusek et al., 2015).

355

Results of table 3 for scenario 1 using only data of the TWT$_f$ signal for the estimation of the unsaturated soil parameters show that:

360
<figure><img src="image_28.png" /></figure>

<table><caption>**Table 3.** First line: Reference values used to build the synthetic calibration data. Then for the different scenarios: estimated mean values (bold), size of the posterior confidence intervals (CIs) (between brackets), and ratio of prior to posterior intervals (italic).</caption><thead><tr><th></th><th>K<sub>s</sub> (cm/s)</th><th>θ<sub>s</sub> (-)</th><th>θ<sub>r</sub> (-)</th><th>α (1/cm)</th><th>n (-)</th></tr></thead><tbody><tr><td>X*</td><td>0.08</td><td>0.40</td><td>0.07</td><td>0.145</td><td>2.68</td></tr><tr><td>Scenario 1</td><td><strong>0.081</strong></td><td><strong>0.39</strong></td><td><strong>0.076</strong></td><td><strong>0.211</strong></td><td><strong>2.75</strong></td></tr><tr><td>TWT<sub>f</sub></td><td>(0.037)</td><td>(0.031)</td><td>(0.14)</td><td>(0.167)</td><td>(0.93)</td></tr><tr><td></td><td>4</td><td>5</td><td>1</td><td>2</td><td>9</td></tr><tr><td>Scenario 2</td><td><strong>0.074</strong></td><td><strong>0.4</strong></td><td><strong>0.081</strong></td><td><strong>0.173</strong></td><td><strong>5.79</strong></td></tr><tr><td>TWT<sub>50</sub></td><td>(0.023)</td><td>(0.008)</td><td>(0.061)</td><td>(0.269)</td><td>(9.99)</td></tr><tr><td></td><td>6</td><td>19</td><td>2</td><td>1</td><td>1</td></tr><tr><td>Scenario 3</td><td><strong>0.078</strong></td><td><strong>0.4</strong></td><td><strong>0.089</strong></td><td><strong>0.167</strong></td><td><strong>5.93</strong></td></tr><tr><td>TWT<sub>120</sub></td><td>(0.011)</td><td>(0.007)</td><td>(0.053)</td><td>(0.195)</td><td>(9.36)</td></tr><tr><td></td><td>13</td><td>24</td><td>2</td><td>1</td><td>1</td></tr><tr><td>Scenario 4</td><td><strong>0.08</strong></td><td><strong>0.4</strong></td><td><strong>0.074</strong></td><td><strong>0.151</strong></td><td><strong>2.72</strong></td></tr><tr><td>TWT<sub>f</sub></td><td>(0.003)</td><td>(0.004)</td><td>(0.015)</td><td>(0.029)</td><td>(0.5)</td></tr><tr><td>TWT<sub>120</sub></td><td>46</td><td>37</td><td>8</td><td>9</td><td>17</td></tr><tr><td>Scenario 5</td><td><strong>0.079</strong></td><td><strong>0.4</strong></td><td><strong>0.073</strong></td><td><strong>0.149</strong></td><td><strong>2.68</strong></td></tr><tr><td>TWT<sub>f</sub></td><td>(0.003)</td><td>(0.004)</td><td>(0.015)</td><td>(0.027)</td><td>(0.49)</td></tr><tr><td>TWT<sub>50</sub></td><td>49</td><td>44</td><td>8</td><td>10</td><td>17</td></tr><tr><td>TWT<sub>120</sub></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

- An accurate estimation of $K_s$, the most sensitive parameter (Fig.4a), is obtained with a CI of 0.037 cm/s and a variation interval reduced by 4.

- A fair estimate of the parameters $\theta_s$ with a standard deviation of 0.031 (-) and a reduction of the interval of variation by 5. This result is relatively surprising as this parameter did not show a strong influence on TWT<sub>f</sub> sensitivity (Fig.4a).

- The parameter $\theta_r$ is not well estimated. Indeed, although its mean estimated value is very close to its reference value, the associated uncertainty of 0.14 is large and the posterior interval is as large as the prior one, which indicates the low reliability of the estimation.

- A poor estimation of $\alpha$, while the sensitivity analysis showed it has a strong influence on TWT<sub>f</sub> (Fig.4a). Its CI is large, with a value of 0.167 cm<sup>-1</sup> and its posterior interval size is half the prior one.
<figure><img src="image_29.png" /></figure>

<figure><img src="image_30.png" /></figure>

370

- The TWT<sub>f</sub> signal yields a mean estimated value *n* = 2.75 ± 0.47 which is close to the reference value *n*<sup>*</sup> = 2.68. The parameter *n* is quite well identified since its posterior interval is 9 times smaller than the prior interval. This is relatively surprising since the sensitivity of *n* is negligible (Fig.4a and 6a5).

375

In summary, using only data of the TWT<sub>f</sub> signal as conditioning information for the hydrogeophysical model calibration yielded well mean estimated parameter values, close to the reference values for all unsaturated soil parameters. However, the examination of the associated uncertainties, showed that only K<sub>s</sub>, θ<sub>s</sub>, and *n* are correctly identified (with narrow posterior intervals with respect to the prior ones). This points out the importance of statistical calibration methods for highly nonlinear problems to investigate not only estimated parameter values but also the associated uncertainties.

The estimation of the unsaturated soil parameters for scenarios 2 and 3, using only data of the TWT<sub>50</sub> or TWT<sub>120</sub> signal for the calibration shows that:

380

- The parameters K<sub>s</sub> and θ<sub>s</sub>, which are the most sensitive parameters during most of the experiment (Fig.4b and 4c), are well identified with small CI size and strong reductions by at least 6 for K<sub>s</sub>, and 19 for θ<sub>s</sub>, of their intervals of variation. We note that the TWT<sub>120</sub> signal allows a much better estimate of both K<sub>s</sub> and θ<sub>s</sub> as their CIs are smaller than the ones estimated with TWT<sub>50</sub>. It is especially true for K<sub>s</sub> where there is almost a factor 2 between the reduction ratios.

385

- The soil parameters θ<sub>r</sub>, α, and *n*, although sensitive (Fig.4b and 4c), cannot be identified from the TWT<sub>50</sub> and TWT<sub>120</sub> signals since their posterior intervals are as large as, or at best two times smaller than their prior ones.

The results of scenario 4 which combines data of TWT<sub>f</sub> and TWT<sub>120</sub> signals show that:

390

- Both parameters K<sub>s</sub>, θ<sub>s</sub>, and *n* are very well identified, with very narrow posterior intervals showing a strong reduction by 46, 37, and 17 of their prior intervals, respectively.

- The parameters θ<sub>r</sub> and α are reasonably well estimated with mean values very close to their reference and intervals of variation reduced by 8 and 9, respectively.

Finally, the results of the last scenario which combines data of TWT<sub>f</sub>, TWT<sub>50</sub> and TWT<sub>120</sub> signals, show performances very similar to scenario 4. Additional information from TWT<sub>50</sub> helped to reduce slightly the posterior intervals of K<sub>s</sub>, θ<sub>s</sub>, and α that in that case show a reduction of 49, 44, and 10 times their prior intervals, respectively.

The results of MCMC for this last scenario are shown in Fig.7 where diagonal plots depict the inferred posterior parameter distributions and the off-diagonal scatterplots represent the pairwise correlations in the MCMC draws. Almost bell-shaped posterior distributions are obtained for all unsaturated soil parameters. Negligible correlations are observed between the parameters, except moderate correlations observed between K<sub>s</sub> and θ<sub>r</sub> (r=-0.78) and between *n* and θ<sub>r</sub> (r=0.64).

395

Note that the parameter *n* is relatively well estimated as the target reference value 2.68 is located in the high sensitivity region (*n* < 3.5) (Fig.6). In the case of a reference value located in the low sensitivity region (*n* ≥ 3.5), the calibration of the hydrogeophysical model using TWT<sub>f</sub> and TWT<sub>120</sub> signals yields a much poorer identification of the parameter *n*. For instance,

400
<figure><img src="image_31.png" /></figure>

<figure><img src="image_32.png" /><figcaption>**Figure 7.** MCMC solutions using TWT<sub>r</sub>, TWT<sub>50</sub> and TWT<sub>120</sub> signals for the calibration of the hydrogeophysical model. The diagonal plots represent the inferred posterior parameter distributions, showing the estimated mean value (dotted black line) and the target value (hard red line). The off-diagonal represents the pairwise correlations between parameters.</figcaption></figure>

using scenario 5 with a reference value $n^* = 4.25$, the estimated mean value is 4.84 with a posterior CI size of 3.6, which corresponds to a reduction of the interval of variation by only 2.

These results point to the high benefit, for the identification of the unsaturated soil parameters, of combining the GPR signal data of a fixed reflector, preferably located sufficiently deep in the soil, with the TWT signal of the moving infiltration wetting front. This combination allows good reliability of almost all soil parameters with very narrow posterior intervals in comparison with the prior ones. In particular, the van Genuchten parameter *n* is relatively well identified for investigated sandy soil where *n* < 3.5.

405
<figure><img src="image_33.png" /></figure>

<figure><img src="image_34.png" /></figure>

# 5 Conclusions

In this work, we investigated the interest in ground penetrating radar (GPR) time data for the estimation of hydrodynamic unsaturated soil parameters. To this aim, water infiltration into an initially unsaturated sandy soil has been simulated using a one-dimensional hydrogeophysical model. GPR time signals have been analyzed from the reflection of the electromagnetic wave on the moving wetting front and on two fixed reflectors located at different depths. GSA, based on PCE decomposition, has been used to assess the effect of the unsaturated soil parameters (saturated hydraulic conductivity, saturated and residual water contents, and Mualem-van Genuchten shape parameters $\alpha$ and $n$) on the different TWT signals. Statistical calibration of the unsaturated soil parameters has been performed with the MCMC sampler using corrupted synthetic observations to evaluate the reliability of the soil parameters from the TWT signals.

The results of GSA showed that the TWT$_f$ signal of the wetting front is different from that of the two fixed reflectors which had similar behavior. For the fixed reflectors, the magnitude of the variance (and therefore the sensitivity of the soil parameters) is more pronounced for deeper reflectors. The TWT$_f$ signal is highly sensitive to $K_s$ and $\alpha$ and moderately sensitive to $\theta_s$. A low sensitivity was observed for $\theta_r$, whereas the parameter $n$ was insensitive. The TWT$_{120}$ signal of the fixed reflector located at 120 cm depth is highly sensitive to $K_s$, $\theta_s$ and $\alpha$, and moderately sensitive to $\theta_r$. The van Genuchten parameter $n$ has a high sensitivity for $n < 3.5$ and a poor sensitivity for $n \ge 3.5$.

The reliability of the unsaturated soil parameters has been assessed for 5 different scenarios of TWT measurement sets. When only data of the TWT$_f$ signal are used as conditioning information for the model calibration, all estimated parameter values were very close to the reference values. However, analyzing the associated uncertainties showed that only $K_s$, $\theta_s$, and $n$ were correctly identified (with narrow posterior intervals). Further, using only data of the TWT$_{50}$ or TWT$_{120}$ signals for the calibration allows also only a good identification of $K_s$ and $\theta_s$ with a strong reduction of their intervals of variation. The best results, in terms of parameter reliability, are obtained with the combination of TWT$_f$ with fixed reflectors. In this case, the four parameters $K_s$, $\theta_s$, $\theta_r$, and $\alpha$ are very well identified with very narrow posterior interval. The van Genuchten parameter $n$ is estimated with a low uncertainty but its estimation degrades in the low sensitivity region $n \ge 3.5$. We note that the deeper reflectors provide more information as the inversion of its signal furnishes parameters with lower uncertainty. Then using two or three reflectors in addition to the wetting front signal doesn't reduce consequently the uncertainty of the parameters.

The results of this study highlight the high benefit of combining TWT signals of fixed and moving (infiltration wetting front) reflectors for very good identification of all the unsaturated soil parameters. It also points out the role of GSA to assess the influence of the parameters on the output signals and the necessity to perform statistical calibration to assess the reliability of model parameters by evaluating not only estimated parameter values but also their associated uncertainties.

# Acknowledgements

The doctoral position of the first author is co-funded by the Grand-Est Region and the University of Strasbourg. This research work is partly funded by the French National Research Agency through the Exciting project ANR-17-CE06-0012-01. Finally, computing time was provided by the HPC-UdS.
<figure><img src="image_35.png" /></figure>

<figure><img src="image_36.png" /></figure>

**Code availability.** The hydrological code WAMOS-1D, the coupled hydrogeophysical model, and the global sensitivity analysis tools can be provided upon request by the authors.

**Data availability.** All data presented in this paper is available upon request from the authors.

**Author contributions.** GJF and LN formulated the aims of the research. BB, LF, and YA provided the hydrological code. YA provided the resources for sensitivity analysis and Bayesian estimation studies. MR established the hydrogeophysical coupled model and ran all experiments under the supervision of GJF and LN, and validation of all co-authors. MR prepared the manuscript with reviewing and editing contributions from all co-authors.

**Competing interests.** The authors declare that they have no conflict of interest.

**Acknowledgements.** The doctoral position of the first author is co-funded by the Grand-Est Region and the University of Strasbourg. This research work is partly funded by the French National Research Agency through the Exciting project ANR-17-CE06-0012-01. The authors would like to acknowledge the High-Performance Computing Center of the University of Strasbourg for supporting this work by providing scientific support and access to computing resources. Part of the computing resources was funded by the Equipex Equip@Meso project (Programme Investissements d'Avenir) and the CPER Alsac calcul/Big Data.
<figure><img src="image_37.png" /></figure>

<figure><img src="image_38.png" /></figure>

References

<figure><img src="image_39.png" /></figure>

Annan, P.: Ground penetrating radar: Principles, procedures and applications, Sensors and Software Inc, 2003.

Bano, M.: Effects of the transition zone above a water table on the reflection of GPR waves, Geophysical Research Letters, 33, https://doi.org/https://doi.org/10.1029/2006GL026158, 2006.

Belfort, B., Toloni, I., Ackerer, P., Cotel, S., Viville, D., and Lehmann, F.: Vadose Zone Modeling in a Small Forested Catchment: Impact of Water Pressure Head Sampling Frequency on 1D-Model Calibration, Geosciences, 8, https://doi.org/10.3390/geosciences8020072, 2018.

<figure><img src="image_40.png" /></figure>

Belfort, B., Weill, S., Fahs, M., and Lehmann, F.: Laboratory Experiments of Drainage, Imbibition and Infiltration under Artificial Rainfall Characterized by Image Analysis Method and Numerical Simulations, Water, 11, https://doi.org/10.3390/w11112232, 2019.

Binley, A., Hubbard, S. S., Huisman, J. A., Revil, A., Robinson, D. A., Singha, K., and Slater, L. D.: The emergence of hydrogeophysics for improved understanding of subsurface processes over multiple scales, Water Resources Research, 51, 3837–3866, https://doi.org/https://doi.org/10.1002/2015WR017016, 2015.

<figure><img src="image_41.png" /></figure>

Birchak, J., Gardner, C., Hipp, J., and Victor, J.: High dielectric constant microwave probes for sensing soil moisture, Proceedings of the IEEE, 62, 93–98, https://doi.org/10.1109/PROC.1974.9388, 1974.

Busch, S., Weihermüller, L., Huisman, J. A., Steelman, C. M., Endres, A. L., Vereecken, H., and van der Kruk, J.: Coupled hydrogeophysical inversion of time-lapse surface GPR data to estimate hydraulic properties of a layered subsurface, Water Resources Research, 49, 8480–8494, https://doi.org/https://doi.org/10.1002/2013WR013992, 2013.

<figure><img src="image_42.png" /></figure>

Costabel, S. and Günther, T.: Noninvasive Estimation of Water Retention Parameters by Observing the Capillary Fringe with Magnetic Resonance Sounding, Vadose Zone Journal, 13, vzj2013.09.0163, https://doi.org/https://doi.org/10.2136/vzj2013.09.0163, 2014.

Dal Bo, I., Klotzsche, A., Schaller, M., Ehlers, T. A., Kaufmann, M. S., Fuentes Espoz, J. P., Vereecken, H., and van der Kruk, J.: Geophysical imaging of regolith in landscapes along a climate and vegetation gradient in the Chilean coastal cordillera, CATENA, 180, 146–159, https://doi.org/https://doi.org/10.1016/j.catena.2019.04.023, 2019.

<figure><img src="image_43.png" /></figure>

de Jong, S. M., Heijenk, R. A., Nijland, W., and van der Meijde, M.: Monitoring Soil Moisture Dynamics Using Electrical Resistivity Tomography under Homogeneous Field Conditions, Sensors, 20, https://doi.org/10.3390/s20185313, 2020.

Doolittle, J. A. and Brevik, E. C.: The use of electromagnetic induction techniques in soils studies, Geoderma, 223-225, 33–45, https://doi.org/https://doi.org/10.1016/j.geoderma.2014.01.027, 2014.

<figure><img src="image_44.png" /></figure>

Dusek, J., Dohnal, M., Snehota, M., Sobotkova, M., Ray, C., and Vogel, T.: Transport of bromide and pesticides through an undis-
turbed soil column: A modeling study with global optimization analysis, Journal of Contaminant Hydrology, 175-176, 1–16,
https://doi.org/https://doi.org/10.1016/j.jconhyd.2015.02.002, 2015.

Edemsky, D., Popov, A., Prokopovich, I., and Garbatsevich, V.: Airborne Ground Penetrating Radar, Field Test, Remote Sensing, 13, https://doi.org/10.3390/rs13040667, 2021.

<figure><img src="image_45.png" /></figure>

Fajraoui, N., Ramasomanana, F., Younes, A., Mara, T. A., Ackerer, P., and Guadagnini, A.: Use of global sensitivity analysis and polynomial chaos expansion for interpretation of nonreactive transport experiments in laboratory-scale porous media, Water Resources Research, 47, https://doi.org/https://doi.org/10.1029/2010WR009639, 2011.

Huisman, J. A., Hubbard, S. S., Redman, J. D., and Annan, A. P.: Measuring Soil Water Content with Ground Penetrating Radar: A Review, Vadose Zone Journal, 2, 476–491, https://doi.org/https://doi.org/10.2136/vzj2003.4760, 2003.
<figure><img src="image_46.png" /></figure>

Jadoon, K. Z., Weihermüller, L., Scharnagl, B., Kowalsky, M. B., Bechtold, M., Hubbard, S. S., Vereecken, H., and Lambot, S.: Estimation of Soil Hydraulic Parameters in the Field by Integrated Hydrogeophysical Inversion of Time-Lapse Ground-Penetrating Radar Data, Vadose Zone Journal, 11, vzj2011.0177, https://doi.org/10.2136/vzj2011.0177, 2012.

490

Jaumann, S. and Roth, K.: Soil hydraulic material properties and layered architecture from time-lapse GPR, Hydrology and Earth System Sciences, 22, 2551–2573, https://doi.org/10.5194/hess-22-2551-2018, 2018.

Jonard, F., Weihermüller, L., Schwank, M., Jadoon, K. Z., Vereecken, H., and Lambot, S.: Estimation of Hydraulic Properties of a Sandy Soil Using Ground-Based Active and Passive Microwave Remote Sensing, IEEE Transactions on Geoscience and Remote Sensing, 53, 3095–3109, https://doi.org/10.1109/TGRS.2014.2368831, 2015.

495

Jones, S. B., Blonquist Jr., J. M., Robinson, D. A., Rasmussen, V. P., and Or, D.: Standardizing Characterization of Electromagnetic Water Content Sensors: Part 1. Methodology, Vadose Zone Journal, 4, 1048–1058, https://doi.org/10.2136/vzj2004.0140, 2005.

Klotzsche, A., Jonard, F., Looms, M., van der Kruk, J., and Huisman, J.: Measuring Soil Water Content with Ground Penetrating Radar: A Decade of Progress, Vadose Zone Journal, 17, 180–202, https://doi.org/10.2136/vzj2018.03.0052, 2018.

500

Laloy, E. and Vrugt, J. A.: High-dimensional posterior exploration of hydrologic models using multiple-try DREAM(ZS) and high-performance computing, Water Resources Research, 48, https://doi.org/10.1029/2011WR010608, 2012.

Léger, E., Saintenoy, A., and Coquet, Y.: Hydrodynamic parameters of a sandy soil determined by ground-penetrating radar inside a single ring infiltrometer, Water Resources Research, 50, 5459–5474, https://doi.org/10.1002/2013WR014226, 2014.

Léger, E., Saintenoy, A., Tucholka, P., and Coquet, Y.: Hydrodynamic Parameters of a Sandy Soil Determined by Ground-Penetrating Radar Monitoring of Porchet Infiltrations, IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 9, 188–200, https://doi.org/10.1109/JSTARS.2015.2464231, 2016.

505

Léger, E., Saintenoy, A., Coquet, Y., Tucholka, P., and Zeyen, H.: Evaluating hydrodynamic parameters accounting for water retention hysteresis in a large sand column using surface GPR, Journal of Applied Geophysics, 182, 104–176, https://doi.org/10.1016/j.jappgeo.2020.104176, 2020.

510

Mara, T. and Tarantola, S.: Application of Global Sensitivity Analysis of Model Output to Building Thermal Simulations, Building Simulation, 1, 290–302, https://doi.org/10.1007/s12273-008-8129-5, 2008.

Moysey, S. M.: Hydrologic trajectories in transient ground-penetrating-radar reflection data, GEOPHYSICS, 75, WA211–WA219, https://doi.org/10.1190/1.3463416, 2010.

Mualem, Y.: A new model for predicting the hydraulic conductivity of unsaturated porous media, Water Resources Research, 12, 513–522, https://doi.org/10.1029/WR012i003p00513, 1976.

515

Richards, L. A.: CAPILLARY CONDUCTION OF LIQUIDS THROUGH POROUS MEDIUMS, Physics, 1, 318–333, https://doi.org/10.1063/1.1745010, 1931.

Robinson, D. A., Campbell, C. S., Hopmans, J. W., Hornbuckle, B. K., Jones, S. B., Knight, R., Ogden, F., Selker, J., and Wendroth, O.: Soil Moisture Measurement for Ecological and Hydrological Watershed-Scale Observatories: A Review, Vadose Zone Journal, 7, 358–389, https://doi.org/10.2136/vzj2007.0143, 2008.

520

Saintenoy, A. and Hopmans, J. W.: Ground Penetrating Radar: Water Table Detection Sensitivity to Soil Water Retention Properties, IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 4, 748–753, https://doi.org/10.1109/JSTARS.2011.2171920, 2011.

Saintenoy, A., Schneider, S., and Tucholka, P.: Evaluating Ground-Penetrating Radar use for water infiltration monitoring, Vadose Zone Journal, 7, 208–214, https://hal.archives-ouvertes.fr/hal-00831408, 2008.

525
<figure><img src="image_47.png" /></figure>

<figure><img src="image_48.png" /></figure>

Saito, H., Kuroda, S., Iwasaki, T., Fujimaki, H., Nagai, N., and Sala, J.: Tracking Infiltration Front Depth Using Time-lapse Multi-offset Gathers Collected with Array Antenna Ground Penetrating Radar, Journal of Visualized Experiments, 2018, https://doi.org/10.3791/56847, 2018.

530 Scharnagl, B., Vrugt, J. A., Vereecken, H., and Herbst, M.: Inverse modelling of in situ soil water dynamics: investigating the effect of different prior distributions of the soil hydraulic parameters, Hydrology and Earth System Sciences, 15, 3043–3059, https://doi.org/10.5194/hess-15-3043-2011, 2011.

Scholer, M., Irving, J., Binley, A., and Holliger, K.: Estimating vadose zone hydraulic properties using ground penetrating radar: The impact of prior information, Water Resources Research, 47, https://doi.org/10.1029/2011WR010409, 2011.

535 Shao, Q., Younes, A., Fahs, M., and Mara, T. A.: Bayesian sparse polynomial chaos expansion for global sensitivity analysis, Computer Methods in Applied Mechanics and Engineering, 318, 474–496, https://doi.org/10.1016/j.cma.2017.01.033, 2017.

Sobol', I.: Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates, Mathematics and Computers in Simulation, 55, 271–280, https://doi.org/10.1016/S0378-4754(00)00270-6, the Second IMACS Seminar on Monte Carlo Methods, 2001.

540 van Genuchten, M. T.: A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils, Soil Science Society of America Journal, 44, 892–898, https://doi.org/10.2136/sssaj1980.03615995004400050002x, 1980.

Vereecken, H., Huisman, J., Bogena, H., Vanderborght, J., Vrugt, J., Hopmans, J., and Vereecken, C.: On the value of soil moisture measurements in vadose zone hydrology: A review, Water Resources Research - WATER RESOUR RES, 44, https://doi.org/10.1029/2008WR006829, 2008.

545 Vrugt, J. A. and Bouten, W.: Validity of First-Order Approximations to Describe Parameter Uncertainty in Soil Hydrologic Models, Soil Science Society of America Journal, 66, 1740–1751, https://doi.org/10.2136/sssaj2002.1740, 2002.

Vrugt, J. A., Stauffer, P. H., Wöhling, T., Robinson, B. A., and Vesselinov, V. V.: Inverse Modeling of Subsurface Flow and Transport Properties: A Review with New Developments, Vadose Zone Journal, 7, 843–864, https://doi.org/10.2136/vzj2007.0078, 2008.

550 Wiener, N.: The Homogeneous Chaos, American Journal of Mathematics, 60, 897–936, http://www.jstor.org/stable/2371268, 1938.

Younes, A., Mara, T. A., Fajraoui, N., Lehmann, F., Belfort, B., and Beydoun, H.: Use of Global Sensitivity Analysis to Help Assess Unsaturated Soil Hydraulic Parameters, Vadose Zone Journal, 12, vzj2011.0150, https://doi.org/10.2136/vzj2011.0150, 2013.

555 Younes, A., Delay, F., Fajraoui, N., Fahs, M., and Mara, T.: Global sensitivity analysis and Bayesian parameter in- ference for solute transport in porous media colonized by biofilms, Journal of Contaminant Hydrology, 191, 1–18, https://doi.org/10.1016/j.jconhyd.2016.04.007, 2016.

Younes, A., Mara, T., Fahs, M., Olivier, G., and Ackerer, P.: Hydraulic and transport parameter assessment using column infiltration experi- ments, Hydrology and Earth System Sciences, 21, 2263–2275, https://doi.org/10.5194/hess-21-2263-2017, 2017.

560 Younes, A., Zaouali, J., Lehmann, F., and Fahs, M.: Sensitivity and identifiability of hydraulic and geophysical parameters from streaming potential signals in unsaturated porous media, Hydrology and Earth System Sciences, 22, 3561–3574, https://doi.org/10.5194/hess-22-3561-2018, 2018.

Zhang, L., Niu, Y., Zhang, H., Han, W., Guang, L., and Xingshuo, P.: Maize Canopy Temperature Extracted From UAV Thermal and RGB Imagery and Its Application in Water Stress Monitoring, Frontiers in Plant Science, 10, https://doi.org/10.3389/fpls.2019.01270, 2019.