the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Coupled hydrogeophysical inversion of an artificial infiltration experiment monitored with groundpenetrating radar: synthetic demonstration
Rohianuu Moua
Nolwenn Lesparre
JeanFrançois Girard
Benjamin Belfort
François Lehmann
Anis Younes
In this study, we investigate the use of groundpenetrating radar (GPR) timelapse monitoring of artificial soil infiltration experiments. The aim is to evaluate this protocol in the context of estimating the hydrodynamic unsaturated soil parameter values and their associated uncertainties. The originality of this work is to suggest a statistical parameter estimation approach using Markov chain Monte Carlo (MCMC) methods to have direct estimates of the parameter uncertainties. Using the GPR time data from the moving wetting front only does not provide reliable results. Thus, we propose to use additional information from other types of reflectors to optimize the quality of the parameter estimation. Water movement and electromagnetic wave propagation in the unsaturated zone are modeled using a onedimensional 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 K_{s}, the saturated and residual water contents θ_{s} and θ_{r}, and the Mualem–van Genuchten shape parameters α 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 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 GPR data to the hydrodynamic soil parameters are different between moving and fixed reflectors, whereas fixed reflectors at various depths have similar sensitivities. K_{s} has a similar and strong influence on the data of both types of reflectors. Concerning the other parameters, for the wetting front, only θ_{s} and α have an influence, and only at long time steps since the total variance is zero at the very beginning of the experiment. On the other hand, for the fixed reflectors, the total variance is not zero at the very start and the parameters θ_{s}, θ_{r}, α and n can have an influence from the very beginning of the infiltration. Results of parameter estimation show that the use of calibration data from the moving or fixed reflectors alone does not enable a good identification of all soil parameters. With the moving reflector, the error between the estimated mean value and the exact target value for θ_{r} and α is 9 % and 45 %, respectively, and less than 3 % for the other parameters. The best reduction of the size of the parameter distribution is obtained for n, with a posterior distribution 9 times smaller than the prior one. For the others, this reduction ratio varies between 1 and 5. For the fixed reflectors, the estimated mean values are far from the target values for α, θ_{r} and n, representing for a reflector located at 120 cm 15 %, 27 %, and 121 %, respectively. On the other hand, when both data are combined, all soil parameters can be well estimated with narrow confidence intervals. For instance, when using both data from the moving wetting front and a fixed reflector located at 120 cm for calibration, the estimated mean values of the errors of all parameters are less than 5 %. Moreover, all parameter distributions are well reduced, with a maximum reduction for K_{s}, leading to a posterior distribution being 46 times smaller than the prior one, and the worst but still satisfactory being for θ_{r} for which the posterior distribution is 8 times smaller than the prior one. The methodology was applied to fine, medium, and coarse sands with very good results, particularly for the finest soil. The thickness of the unsaturated zone was also tested (0.5, 1, and 2 m) and a better estimation of the hydrodynamic parameters is obtained when the water table is deeper. In addition, the height of water applied in the infiltrometry test influences the speed of the test without affecting the performance of the proposed method.
 Article
(3293 KB)  Fulltext XML
 BibTeX
 EndNote
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, which 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 nonlinear behavior and is characterized by two fundamental hydraulic relationships, namely, the water retention and the hydraulic 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 α and n.
Different approaches can be applied to estimate the unsaturated soil parameters. In soil physics, the reference method relies on laboratory measurements conducted on soil core samples. Such experiments can use various techniques such as thermogravimetry or tensiometry, but common practices rely on measurements of hydraulic fluxes (Vereecken et al., 2008). Laboratory measurements can provide direct measurements of the soil hydraulic properties or state variables and can therefore be of great accuracy at the column scale. On the other hand, they are 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 laborintensive for deep or largescale investigations (Binley et al., 2015). Furthermore, the conservation of collected samples can be challenging because of issues of compaction and changes in porosity.
At the field scale, the soil hydraulic properties and state variables can be estimated using numerous approaches. Measurements of the soil water content, water pressure, and hydraulic conductivity can show significant variations because of their sensitivity to different hydrological processes. As a consequence, such measurements are convenient for the estimation of soil parameters of the subsurface at the field scale by inverse modeling. Measuring techniques for soil hydraulic properties and state variables can be classified into two categories based on whether the measuring devices have to be in direct contact with the soil or not. In the former, when the measuring devices must be in direct contact with the soil, measurements can present a spatial support around the micro (mm–cm) and local scale (cm–m) with water content sensing techniques (using thermal or electromagnetic sensors, e.g., capacitance or time domain reflectometry, Jones et al., 2005; Belfort et al., 2019), water pressure measurements with tensiometers (Cassel and Klute, 1986) or psychrometers (Rawlins and Campbell, 1986), and hydraulic conductivity measurements with permeameters (Kodešová et al., 1998) and infiltrometers (Muntz et al., 1905). 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 for a good characterization of the state of the soil. For measurements with sensors, however, their installation is often laborious, timeconsuming, and destructive (Huisman et al., 2003; Dal Bo et al., 2019). Furthermore, their reliability requires an accurate calibration (Robinson et al., 2008).
Other techniques use noninvasive devices that do not have to be in direct contact with the soil, such as remote sensing and hydrogeophysical methods. Remote sensing techniques use devices that operate remotely and relatively far from the ground, such as unmanned aerial vehiclebased 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 contactbased 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).
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 geophysical to hydraulic properties or states requires the use of robust petrophysical relationships to provide reliable estimates of hydraulic parameters.
Today, GPR is widely used in the field of hydrogeophysics. Different techniques have been reviewed and discussed by Huisman et al. (2003) and Klotzsche et al. (2018). Indeed, 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 (Klotzsche et al., 2018). 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) enable deeper penetrations (until a maximum depth between 1 m and 3 m in most organic media). The temporal variability of the soil water content can be characterized from timelapse GPR measurements. In this case, the GPR method is applied in a static approach, where, instead of classically imaging the spatial variation of the properties of the subsurface, the device is set immobile and captures how the properties of the soil change over time. GPR data can be collected during artificial hydraulic processes (e.g., infiltration, runoff, drainage, imbibition) that can provide valuable information on the flow characteristics. Compared to other hydraulic processes, artificially forced infiltration is particularly fast. It also induces a rapidly evolving transient hydraulic perturbation. Timelapse GPR is characterized by a high spatial and temporal resolution and is therefore well adapted for monitoring such a fast hydraulic process. Artificial infiltration process is also easy to establish since it only requires the application of a positive water pressure head on the soil surface. Hence, timelapse GPR monitoring of artificial infiltration experiments is usually effortless and timesaving. Furthermore, except in the case of borehole investigations, the GPR device can be laid on or raised above the surface. For these reasons, timelapse GPR monitoring of artificial infiltration is fast as well as easy to apply and repeat at multiple locations, and, when used on or above the soil surface, it is nondestructive. Therefore, it is one of the cheapest approaches that fits well into the context of mapping the heterogeneity of unsaturated soil parameters at a small catchment scale.
Various studies have investigated the monitoring of different types of flow processes with timelapse GPR in the context of evaluating the soil hydraulic states, hydraulic properties, or unsaturated soil parameters (e.g., Saintenoy et al., 2008; Moysey, 2010; Scholer et al., 2011; Busch et al., 2013; Tran et al., 2014; Jonard et al., 2015; Jaumann and Roth, 2018; Léger et al., 2014; Léger et al., 2016; Léger et al., 2020).
At the laboratory scale, Léger et al. (2020) have monitored imbibitiondrainage experiments using a singleoffset surface GPR. Jaumann and Roth (2018) conducted similar experiments but at the test site scale, where they showed reasonable results when estimating the soil unsaturated parameters and the subsurface architecture. As already pointed out, however, this hydraulic process can take longer than an infiltration to reach a steady state and is also practically harder to conduct at the field scale. Busch et al. (2013) calibrated the Mualem–van Genuchten parameters of their model by monitoring natural precipitation and evapotranspiration events at the field scale. Such slow hydraulic processes can last several days or months and are therefore not suitable for easy and fast characterization. Infiltration experiments have been conducted at the laboratory scale by Moysey (2010), where they considered the GPR twoway travel time (TWT) from various sources of reflection. They showed that the Mualem–van Genuchten shape parameter n is the most poorly constrained among all unsaturated soil parameters. On the field, infiltration processes have been monitored with borehole GPR (Scholer et al., 2011), singleoffset (Léger et al., 2014; Léger et al., 2016) and multioffset (Saito et al., 2018) surface GPR, or offground GPR (Jadoon et al., 2008, Jadoon et al., 2012; Jonard et al., 2015). For practicality, surface GPR is preferred over offground and borehole GPR, the latter also being destructive by nature. Saito et al. (2018) used a more complex multioffset and multichannel surface GPR to directly monitor the wetting front progression. The monochannel multioffset technique is usually not suited for monitoring experiments with high temporal variability, as the offset must be adjusted between each measurement. The multichannel technique has the advantage of being multioffset and is, therefore, able to simultaneously determine the propagation speed and the depths of reflectors.
In the present study, we are interested in using a quick, easytoapply, and cheap fieldscale method to characterize the unsaturated soil parameters. To this end, timelapse GPR monitoring of artificial infiltration is a wellsuited protocol. It is similar to ring infiltrometry methods but with additional information from GPR measurements. In the literature, the work of Léger et al. (2014) is the closest one to consider this protocol for parameter estimation. The authors demonstrated the relevance of such a methodology in evaluating the hydraulic parameters of sandy soil. They 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) used an optimizationbased inversion algorithm which did not make it possible to assess the reliability of the estimated values since the uncertainty associated with the calibrated parameters was not 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, which was satisfactory enough for them to obtain such remarkable results. The original work presented here aims to extend the actual state of the art 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;

investigating the influence of all soil parameters (the saturated hydraulic conductivity, the saturated and residual water contents, and the Mualem–van Genuchten shape parameters α and n) on the GPR TWT data of the three reflectors using global sensitivity analysis (GSA). The GSA allows us to estimate the soil parameter range where timelapse GPR data monitoring is sensitive to these parameters. It also provides some insights into which parameter is more sensitive at the beginning of the infiltration experiment or at the end of the infiltration;

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: Sect. 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, Sect. 4 discusses the results of the soil parameter estimation with MCMC for different scenarios including varying soil types, water table depths, and surface boundary conditions.
2.1 Test case description
In this work, we conduct a synthetic study on the timelapse GPR monitoring of artificial infiltration protocol, prior to applying it in real conditions. The idea is to perform synthetic experiments under the same conditions as real experiments so as to better understand the pertinence of the investigated protocol when used for estimating the unsaturated soil parameters. The test case considered is a hypothetical onedimensional experiment of water infiltration in a homogeneous sandy soil of 150 cm (Fig. 1a). The approach used to drive the artificial infiltration is comparable to other techniques commonly used to estimate the properties of the porous medium, such as single or doublering infiltrometry. As evidenced in other studies (e.g., Léger et al., 2014), the idea is to add information from the GPR data monitored during the infiltration to have access to more of the hydrodynamic parameters. In the present synthetic case, a constant pressure head of 10 cm is applied at the surface of the soil (i.e., a 10 cm water ponding Dirichlettype boundary 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 thickness above a saturated zone of 50 cm thickness. 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 timelapse 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}.
2.2 The mathematical model
2.2.1 Unsaturated flow model
Water infiltration in unsaturated or saturated soils is governed by the onedimensional Richards equation (Richards, 1931):
where h (cm) is the pressure head; z is the depth (cm), taken as positive in the downward direction; t is the time (s), θ (cm^{3} cm^{−3}) is the actual water content, and K(θ) (cm s^{−1}) 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).
The interdependencies of the pressure head, conductivity, and water content are described using the standard models of Mualem, 1976 and van Genuchten (1980):
where S_{e}(h) (–) is the effective saturation, θ_{s} and θ_{r} (cm^{3} cm^{−3}) are the saturated and residual water contents, respectively, K_{s} (cm s^{−1}) is the saturated conductivity, $m=\mathrm{1}\mathrm{1}/n$, α (1 cm^{−1}), and n (–) are the Mualem–van Genuchten shape parameters, and L (–) is a parameter characterizing the tortuosity of the flow paths of moving water in the interconnected pores of the soil. It is set at L=0.5 here, following the works of Mualem (1976) and van Genuchten (1980).
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 timelapse 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 ϵ (–) of a threephase (water–solid–air) medium to its water content by:
where ϕ (–) is the porosity, considered equal to the saturated water content θ_{s}, and ϵ_{w}=80, ϵ_{s}=2.5 (Léger et al., 2014), and ϵ_{a}=1 are the dielectric constants of water, silica (sand), and air, respectively.
In this work, the soil is considered as a linear and isotropic nonmagnetic medium. When working with frequencies below 1 GHz, the soil electrical conductivity can be neglected. In this case, the electromagnetic waves propagate at speed V (cm ns^{−1}) (Annan, 2003):
where c≈30 cm ns^{−1} is the speed of electromagnetic waves in air, and ϵ (–) 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 wave.
We consider a onedimensional 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\mathrm{1}/\mathrm{2}}$ and ${z}_{i+\mathrm{1}/\mathrm{2}}$. The TWT for the reflection occurring at the interface $(i\mathrm{1}/\mathrm{2})$ between the elements i−1 and i can be expressed as the sum of the vertical TWT in each element above i:
in which $\left{l}_{j}\right$ (cm) is the length of the element j above i, and V_{j} (cm ns^{−1}) 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\mathrm{1}/\mathrm{2})$ is defined by:
where ϵ(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 at the wetting front makes it easily detectable. This statement is true in the presented test case and for any parameter value taken from the prior distributions tested (Table 1). On the contrary, the water table may be hidden due to the softer change in the reflection coefficient at the capillary fringe (Bano, 2006; Saintenoy and Hopmans, 2011).
Note that one could easily consider a nonperpendicular incidence of the GPR wave at the interface, introducing the incidence angle in Eq. (7). Nevertheless, the offset between the TX and RX antenna for a 800 MHz GPR system is typically around 10 cm. By simple trigonometry, the incidence angle is 5^{∘} at 50 cm depth, 2^{∘} at the deeper reflector, the reflection coefficient is then very close to the normal incidence, and Eq. (7) is considered in the following. Considering more closely the physics of the radar wave emission and propagation in porous media, if one needs to consider precisely the wave amplitude (such as in full waveform inversion), one should consider the radiation pattern of the antenna. The latter is linked to the dielectric contrast at the surface and the antenna characteristics. Whether one needs to calculate it precisely or one should consider specific acquisition configurations to handle this effect (generally normalization of the signal by a reference signal), we choose to keep our approach as simple as possible, to be applicable to any system easily, without the need for calibration. Hence, we have chosen to consider only the travel time, and not to use the amplitude. Note that the latter is more sensitive to attenuation (electrical conductivity) than to the hydraulic parameters.
2.3 The numerical model
The variation of the water content in the soil during the infiltration is computed using the WAMOS1D code (Belfort et al., 2018). The model describes the water movement in the porous medium using the Richards equation (Eq. 1), and the constitutive relationships between the pressure, the hydraulic conductivity, and the volumetric water content given by Eqs. (2) and (3). The domain of 150 cm depth is discretized with uniform elements 1 cm thick with homogeneous properties. Such discretization enables an appropriate model precision and a low enough computation time. The WAMOS1D 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 ϵ using the petrophysical relationship in Eq. (4) and into a GPR wave propagation speed profile V using Eq. (5). Then, the timelapse 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 timelapse 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 for the wetting front position ${z}_{i\mathrm{1}/\mathrm{2}}^{\ast}$, 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 ${\mathrm{TWT}}_{\mathrm{f}}=\text{TWT}\left({z}_{i\mathrm{1}/\mathrm{2}}^{\ast}\right)$ 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 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 the simulation time. The water table is assumed to be reached when the maximum reflection coefficient of Eq. (7) is below 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.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 variancebased sensitivity analysis, allowing for 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}, θ_{s}, θ_{r}, α, 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 $\mathit{X}=\mathit{\{}{x}_{\mathrm{1}},{x}_{\mathrm{2}},\mathrm{\dots},{x}_{p}\mathit{\}}$ that yields a random response y(X), the two variancebased sensitivity measures, also called Sobol indices (Sobol', 2001), are

the firstorder sensitivity index:
$$\begin{array}{}\text{(8)}& {S}_{\mathrm{i}}={\displaystyle \frac{\text{Var}\left[E\left[y\right(\mathit{X}\left)\right{x}_{i}]\right]}{\text{Var}\left[y\right(\mathit{X}\left)\right]}}\in [\mathrm{0},\mathrm{1}]\end{array}$$ 
the total sensitivity index:
$$\begin{array}{}\text{(9)}& {\mathrm{ST}}_{\mathrm{i}}={\displaystyle \frac{E\left[\text{Var}\left[y\right(\mathit{X}\left)\right{x}_{i}]\right]}{\text{Var}\left[y\right(\mathit{X}\left)\right]}}\in [\mathrm{0},\mathrm{1}],\end{array}$$
where ${x}_{i}=\mathit{X}\setminus {x}_{i}$ is the set of all parameters except x_{i}, E() and $E(..)$ are the expectation and the conditional expectation operators, respectively, and Var() and Var$(..)$ are the variance and the conditional variance, respectively. The firstorder index S_{i} quantifies the contribution of the parameter x_{i} alone to the total variance of y(X), while ST_{i} also includes all interactions of x_{i} with the other parameters x_{−i}.
To perform a variancebased 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(X) as a set of orthonormal multivariate polynomials of a degree not exceeding D:
where $\mathit{\beta}={\mathit{\beta}}_{\mathrm{1}},{\mathit{\beta}}_{\mathrm{2}},\mathrm{\dots},{\mathit{\beta}}_{p}\in {\mathbb{R}}^{p}$ is a pthdimensional index, s_{β} are the PC coefficients, and Ψ_{β} are the generalized polynomial chaos of degree $\left\mathit{\beta}\right={\sum}_{i=\mathrm{1}}^{p}{\mathit{\beta}}_{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). This combination of parameters investigated in the GSA is exhaustive and allows one to consider a large panel of soil types. Note that while simulations with values of n between 1 and 1.5 would allow us to investigate a wider range of porous media, they also take much longer to end.
The number of coefficients for a full PCE representation is $P=(p+D)\mathrm{!}/p\mathrm{!}D\mathrm{!}$. 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 for the best fit (in the least square sense) of the PCE surrogate model to the hydrogeophysical model for the M realizations. To work with lowdiscrepancy sets, the M realizations correspond to sets of input parameters sampled from their prior probability distributions, using quasirandom Sobol sequences (Shao et al., 2017). Because each parameter varies in its own range and has a proper unit, the parameter prior intervals are normalized to $[\mathrm{1},\mathrm{1}]$ during PCE computation. We illustrate the principle of the construction of the PCE with our hydrogeophysical model in Fig. 3.
A PCE is constructed at each time step for all model responses (TWT_{f}, TWT_{50}, and TWT_{120}) 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_{50} 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^{2}). 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.
3.2 GSA results
The temporal distribution of the output variance of the three TWT signals (TWT_{f}, TWT_{50}, and TWT_{120}) is represented 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.
TWT_{f} has a different behavior from the TWT signals of fixed reflectors TWT_{50} and TWT_{120} (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 5 times more significant for TWT_{120} than for TWT_{50}. This is in agreement with the physics since the zone of the porous medium affecting the GPR wave is more important for the TWT_{120} signal than for the TWT_{50}. In addition, the period of influence of the unsaturated parameters (θ_{r}, α, n) is also more important for TWT_{120} than for TWT_{50} 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_{f} and TWT_{120} signals.
3.2.1 GSA of the TWT_{f} signal
TWT_{f} variance is zero at the beginning of the infiltration (Fig. 4a), which means that the TWT_{f} signal is not affected by the initial conditions. Indeed, the infiltration wetting front and the TWT_{f} signal start at zero for all parameter sets. Then, the variance of the signal increases until a maximum of 60 ns^{2}, reached at around 3 min. After that, the variance decreases, but maintains a significant value of around 25 ns^{2} (Fig. 4a). Concerning parameter sensitivities, at the beginning, the TWT_{f} signal is mainly affected by K_{s}. The influence of this parameter decreases over time and reaches zero for long times when steadystate conditions (corresponding to a fully saturated soil) are reached. The parameter θ_{s} has a moderate influence on the TWT_{f} signal. Its influence is not observable at short times since unsaturated conditions occur. Overall, the most influential parameter on the TWT_{f} 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_{f} is artificially maintained when the infiltration wetting front reaches the water table, while physically the signal disappears. The effects of the parameters θ_{r} 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. To estimate this interaction, we plot the difference between the total (ST_{i}) and the firstorder (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}≈S_{i}), except for K_{s} and α (Fig. 5a). Hence, the interaction between these two parameters affects the variance of the TWT_{f} signal as represented by the blank region between the variance curve and the shaded area (Fig. 4a).
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. Figure 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 us to determine 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, and therefore we represent them at the three time steps (t_{1} = 1 min, t_{2} = 5 min, and t_{3} = 200 min) highlighted with dashed 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, the following can be noted:

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} does not 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}.

θ_{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 θ_{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 θ_{s}.

The sensitivity of TWT_{f} to θ_{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 with the physics of the process (when θ_{r} increases, the speed of the electromagnetic wave decreases, and the TWT_{f} signal increases).

The van Genuchten parameter α is highly influential notably for long times (t_{3} = 200 min). A small variation of the parameter α can induce a strong variation of the TWT_{f} signal. Notably, the sensitivity of α is very high for α≤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. TWT_{f} is sensitive to all K_{s}, θ_{s} and θ_{r} values tested (see Table 1), but is also sensitive to α < 0.05 cm^{−1} and n < 2. This means a wide range of soil types.
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 θ_{s}, θ_{r}, α, and n. The most influential parameter at the beginning of the experiment is the parameter α. Over time, the effect of this parameter decreases, whereas the effect of θ_{s} increases. For long times, θ_{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, after which 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 firstorder Sobol indices is negligible for all parameters except after 1 min for the parameters K_{s}, α, and θ_{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 the following:

As for the TWT_{f} signal, K_{s} is highly sensitive, especially for t = 1 min and 5 min.

The saturated water content θ_{s} is very influential for all times. The TWT_{120} varies almost linearly with θ_{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, θ_{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 θ_{r} increases, the water content increases and, hence, the TWT_{120} increases.

The van Genuchten parameter α is highly sensitive. However, contrarily to the TWT_{f} signal where α is highly sensitive at long times (t_{3} = 200 min), the sensitivity of α for the TWT_{120} signal is high at short times (t_{1} = 1 min). For long times, the influence of α disappears since the soil becomes fully saturated. The negative slope of the curve of the TWT_{120} signal as a function of α observed at the beginning of the experiment is consistent with the physics of the process. Indeed, when α 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≥3.5. Finally, TWT_{120} shows similar sensitivity for K_{s}, θ_{s} and θ_{r} but slightly higher than TWT_{f}. For α, they show complementarity, which makes the procedure very efficient for α < 0.05 using all the infiltration experiments (early time for TWT_{f} and late time for TWT_{120}). TWT_{120} is more sensitive to n, but only until values of 3.5. The GSA study shows that the monitoring of the infiltration using both the TWT from the infiltration front and (at least) a fixed reflector displays a significant sensitivity for a wide range of soil types (see in Table 1 the hydraulic parameters range tested). The next section demonstrates the use of this sensitivity to calibrate these parameters.
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 GPRmonitored infiltration experiment in order to address the following questions:

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

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?

What is the optimal set of TWT measurements that yields good reliability of all unsaturated soil parameters?
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 for 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 make it possible to evaluate the reliability of the parameters via uncertainty quantification.
In the sequel, the MCMC method is performed with the DREAM_{(ZS)} (DiffeRential Evolution Adaptive Metropolis) software (Laloy and Vrugt, 2012; Vrugt, 2016). This software samples the posterior probability density function (pdf) by running multiple Markov chains simultaneously for global exploration of the parameter space. The prior distributions of the parameters are the same as in the GSA (Table 1). The DREAM_{(ZS)} then automatically tunes the scale and orientation of the proposal distribution until we get the posterior target pdf. A MATLAB toolbox of the DREAM_{(ZS)} algorithm is available for Bayesian inference in fields ranging from physics, chemistry, and engineering to ecology, hydrology, and geophysics. The vector of unknowns is formed by the five unsaturated soil parameters (K_{s}, θ_{s}, θ_{r}, α, n). Compared to the GSA, which enables investigation of a large panel of soil types, the parameter estimation is demonstrated on a single synthetic case. A reference solution is generated by simulating the hydrogeophysical problem formed by the system of Eqs. (1)–(6) using the following reference parameter values, ${K}_{\mathrm{s}}^{\ast}=\mathrm{0.08}$ cm s^{−1}, ${\mathit{\theta}}_{\mathrm{s}}^{\ast}=\mathrm{0.4}$, ${\mathit{\theta}}_{\mathrm{r}}^{\ast}=\mathrm{0.07}$, ${\mathit{\alpha}}^{\ast}=\mathrm{0.145}$ cm^{−1}, ${n}^{\ast}=\mathrm{2.68}$, as shown in Table 3. These parameter values correspond to those of a sandy porous medium present in an experimental platform where we test the protocol under real conditions. The modeled TWT_{f}, TWT_{50}, and TWT_{120} signals used as synthetic 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 σ=0.5 ns. This error corresponds to an uncertainty of 1 ns, which is realistic in the instance of an 800 MHz GPR antenna.
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 reaches 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 depth of a fully saturated porous medium. For TWT_{f}, the plateau value is also reached when the infiltration front reaches 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 depth.
The reliability of the unsaturated soil parameters is assessed for five different scenarios of measurement sets. In the first scenario, only data of the wettingfront 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.
In the five scenarios, the MCMC sampler uses three parallel chains and a total number of 50 000 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). The posterior histograms and the derived statistics are obtained from the last 12 500 simulations, as mentioned earlier, for which the Gelman–Rubin (Gelman and Rubin, 1992)) criterion is verified and the chains are stable and not autocorrelated.
The results in Table 3 for scenario 1 using only data of the TWT_{f} signal for the estimation of the unsaturated soil parameters show the following:

An accurate estimation of K_{s}, the most sensitive parameter (Fig. 4a), is obtained with a CI of 0.037 cm s^{−1} and a variation interval reduced by 4.

A fair estimate of the parameters θ_{s} with a standard deviation of 0.031 (–) and a reduction of the interval of variation by 5 is obtained. This result is relatively surprising as this parameter did not show a strong influence on TWT_{f} sensitivity (Fig. 4a).

The parameter θ_{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.

There is a poor estimation of α, while the sensitivity analysis showed it has a strong influence on TWT_{f} (Fig. 4a). Its CI is large, with a value of 0.167 cm^{−1} and its posterior interval size is half the prior one.

The TWT_{f} signal yields a mean estimated value $n=\mathrm{2.75}\pm \mathrm{0.47}$ which is close to the reference value ${n}^{\ast}=\mathrm{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 (Figs. 4a and 6a5). n values between 1.5 and 3 could represent silty loam, sandy loam, or sand. Therefore, it is not obvious to consider this result as a good identification. We have performed another estimation (the results are not presented here) with a target value of the n parameter equal to 6, which is located in the low sensitivity region of the parameter. In this case, the inversion led to a relatively good estimation (estimated mean value of 6.97). However, the parameter was poorly identified since its posterior interval was still large, although it was a bit smaller than the prior parameter interval.
In summary, using only data of the TWT_{f} signal as conditioning information for the hydrogeophysical model calibration yielded well estimated mean values of the parameters, close to the reference values for all unsaturated soil parameters. However, the examination of the associated uncertainties showed that only K_{s}, θ_{s}, and n are correctly identified (with narrow posterior intervals with respect to the prior ones). This highlights the importance of statistical calibration methods for highly nonlinear problems for investigating 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_{50} or TWT_{120} signal for the calibration, shows the following:

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

The soil parameters θ_{r}, α, and n, although sensitive (Fig. 4b and c), cannot be identified from the TWT_{50} and TWT_{120} signals since their posterior intervals are as large as, or at best 2 times smaller than, their prior intervals.
The results of scenario 4, which combines data of TWT_{f} and TWT_{120} signals, show the following:

Parameters K_{s}, θ_{s}, and n are very well identified, with very narrow posterior intervals showing a strong reduction of 46, 37, and 17 times of their prior intervals, respectively.

Parameters θ_{r} and α are reasonably well estimated with mean values very close to their reference and intervals of variation reduced by 8 and 9, respectively.
Figure 7 shows the posterior histograms obtained from scenarios 1, 3, and 4. For all parameters, the displayed intervals correspond to the prior upper and lower limits of Table 1.
Finally, the results of the last scenario, which combines data of TWT_{f}, TWT_{50} and TWT_{120} signals, show performances very similar to scenario 4. Additional information from TWT_{50} helped to reduce slightly the posterior intervals of K_{s}, θ_{s}, and α, which 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. 8 where diagonal plots depict the inferred posterior parameter distributions and the offdiagonal scatterplots represent the pairwise correlations in the MCMC draws. Almost bellshaped posterior distributions are obtained for all unsaturated soil parameters. Negligible correlations are observed between the parameters, except moderate correlations observed between K_{s} and θ_{r} ($r=\mathrm{0.78}$) and between n and θ_{r} (r=0.64) .
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_{f} and TWT_{120} signals yields a much poorer identification of the parameter n. For instance, using scenario 5 with a reference value ${n}^{\ast}=\mathrm{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.
To complete the numerical study, the protocol was tested varying the boundary conditions. One may wonder how much the thickness of the vadose zone would impact the calibration of the hydraulic parameters. For that purpose, three scenarios were considered, varying the water table depth from 50 cm to 1 and 2 m, and assuming a hydrostatic initial profile. Results of the MCMC calibration depicted in Fig. 9 show that the five parameters are even better estimated when the water table is deeper. We explain this result as follows: when the vadose zone is thicker, then the initial water content profile highlights a larger variation with depth, which is perturbed when the infiltration propagates. In the shallowest case, with a 50 cm deep water table, α and n could not be recovered because the water content (which directly affects the radar propagation) in the vadose zone is already close to the saturation conditions. One should note that in this case, we maintained the fixed reflector at 120 cm depth, which means it is above the water table in the 2 m case. In this latter case, a deeper fix reflector would enhance the result (as seen in previous scenarios), but in field conditions, a deeper fixed reflector would be harder to be reliably detected. Nevertheless, we show here that it is not necessary to have a reflector below the water table to obtain an accurate calibration. We also vary the surface boundary conditions by using different heights of water ponding at the surface, which would practically mean varying the height of the infiltration ring, from 5 to 10 and 20 cm. As one would expect, it only affects the duration of the infiltration experiment without impacting on the accuracy of the hydraulic parameter calibration. The infiltration duration could be divided by 2 if the pressure head is doubled. This is worth noting, especially for a medium with low permeability, in order to speed up the experiment.
Last, the efficiency of the protocol was numerically tested on three types of soils used in the experimental platform SCERES in Strasbourg (Bohy et al., 2006). A coarse, a medium, and a fine sand are considered (see Table 4). The permeability varies over 2 orders of magnitude, and θ_{r} and α are, respectively, multiplied by 4 and 10, and n varies from 2.0 to 2.7. The results of the calibration using the TWT_{f} and TWT_{120} summarized in Table 4 show that the parameters for the three materials are well estimated, even better when the sand is finer. All five hydraulic parameters are recovered here, considering a water table at 1 m depth.
These results evidence that the GPR signal data of both the wetting front and a fixed reflector can provide very different but complementary information for the identification of the unsaturated soil parameters. They also point to the significant benefit 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 leads to 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.
The aim of the present study was to optimize a cheap method used at the field scale to characterize the hydraulic parameters of a porous medium. To this end, we investigated a particular protocol: timelapse GPR monitoring of artificial infiltration experiments. Water infiltration into an initially unsaturated sandy soil was simulated using a onedimensional hydrogeophysical model. GPR time signals were 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, was used to assess the effect of the unsaturated soil parameters (saturated hydraulic conductivity, saturated and residual water contents, and Mualem–van Genuchten shape parameters α and n) on the different TWT signals. Statistical calibration of the unsaturated soil parameters was 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 α and moderately sensitive to θ_{s}. A low sensitivity was observed for θ_{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}, θ_{s}, and α, and moderately sensitive to θ_{r}. The van Genuchten parameter n has a high sensitivity for n<3.5 and a poor sensitivity for n≥3.5. The GSA study shows that the monitoring of the infiltration using both the TWT from the infiltration front and (at least) a fixed reflector has a significant sensitivity for a wide range of soil types.
The reliability of the unsaturated soil parameters was assessed for five different scenarios of TWT measurement sets. When only data of the TWT_{f} signal were 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}, θ_{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 enabled also good identification of K_{s} and θ_{s} with a significant decrease in their intervals of variation. The best results, in terms of parameter reliability, were obtained with the combination of TWT_{f} with at least one fixed reflector. In this case, the four parameters K_{s}, θ_{s}, θ_{r}, and α were very well identified with very narrow posterior intervals. The van Genuchten parameter n was estimated with a low uncertainty but its estimation degraded in the low sensitivity region n≥3.5. We note that the deeper reflectors provide more information as the inversion of the signal furnishes parameters with lower uncertainty. Then using two or three reflectors in addition to the wetting front signal does not consequently reduce the uncertainty of the parameters. The procedure was applied for three types of soil ranging from coarse to fine sand and the results of MCMC simulations show its efficiency. The best estimate was obtained for the finest material. In field conditions, one could expect a higher clay content, which would decrease the electrical resistivity and then would attenuate the GPR signal, limiting the penetration depth of the radar wave. Our numerical study shows that using a thicker water level in the infiltration ring to apply a greater pressure head could speed up the protocol without having any impact on the MCMC results. Furthermore it appears that a deeper water table makes the calibration protocol more efficient and accurate. A limitation is observed for very shallow water tables (e.g., 50 cm) where the van Genuchten parameters α and n could not be estimated because the vadose zone is already almost saturated.
The results of this study highlight the significant 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 in assessing the influence of the parameters on the output signals and the necessity of performing statistical calibration to assess the reliability of model parameters by evaluating not only estimated parameter values but also their associated uncertainties.
The hydrological code WAMOS1D, the coupled hydrogeophysical model, and the global sensitivity analysis tools can be provided by the authors upon request.
All data presented in this paper are available from the authors upon request.
JFG and NL formulated the aims of the research. BB, FL, and AY provided the hydrological code. AY provided the resources for sensitivity analysis and Bayesian estimation studies. RM established the hydrogeophysical coupled model and ran all experiments under the supervision of JFG and NL, and validation of all coauthors. RM prepared the manuscript with reviewing and editing contributions from all coauthors.
The contact author has declared that none of the authors has any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
The authors would like to acknowledge the HighPerformance 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 Alsacalcul/Big Data. Computing time was provided by the HPCUdS.
The doctoral position of the first author is cofunded by the GrandEst Region and the University of Strasbourg. This research work is partly funded by the French National Research Agency through the Exciting project (grant no. ANR17CE06001201).
This paper was edited by Marnik Vanclooster and reviewed by two anonymous referees.
Annan, P.: Ground penetrating radar: Principles, procedures and applications, Sensors and Software Inc, https://geolportal.sdsu.edu/jiracek/sage/documents/Sensors and Software GPR Manual.pdf (last access: 5 December 2023), 2003. a
Bano, M.: Effects of the transition zone above a water table on the reflection of GPR waves, Geophys. Res. Lett., 33, L13309, https://doi.org/10.1029/2006GL026158, 2006. a
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 1DModel Calibration, Geosciences, 8, 72, https://doi.org/10.3390/geosciences8020072, 2018. a
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, 2232, https://doi.org/10.3390/w11112232, 2019. a
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 Resour. Res., 51, 3837–3866, https://doi.org/10.1002/2015WR017016, 2015. a, b
Birchak, J., Gardner, C., Hipp, J., and Victor, J.: High dielectric constant microwave probes for sensing soil moisture, P. IEEE, 62, 93–98, https://doi.org/10.1109/PROC.1974.9388, 1974. a
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 timelapse surface GPR data to estimate hydraulic properties of a layered subsurface, Water Resour. Res., 49, 8480–8494, https://doi.org/10.1002/2013WR013992, 2013. a, b
Cassel, D. K. and Klute, A.: Water Potential: Tensiometry, Chap. 23, in: Methods of Soil Analysis, edited by: Klute, A., John Wiley & Sons, Ltd, https://doi.org/10.2136/sssabookser5.1.2ed.c23, 563–596, 1986. a
Costabel, S. and Günther, T.: Noninvasive Estimation of Water Retention Parameters by Observing the Capillary Fringe with Magnetic Resonance Sounding, Vadose Zone J., 13, vzj2013.09.0163, https://doi.org/10.2136/vzj2013.09.0163, 2014. a
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/10.1016/j.catena.2019.04.023, 2019. a
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, 5313, https://doi.org/10.3390/s20185313, 2020. a
Doolittle, J. A. and Brevik, E. C.: The use of electromagnetic induction techniques in soils studies, Geoderma, 223–225, 33–45, https://doi.org/10.1016/j.geoderma.2014.01.027, 2014. a
Dusek, J., Dohnal, M., Snehota, M., Sobotkova, M., Ray, C., and Vogel, T.: Transport of bromide and pesticides through an undisturbed soil column: A modeling study with global optimization analysis, J. Contam. Hydrol., 175–176, 1–16, https://doi.org/10.1016/j.jconhyd.2015.02.002, 2015. a
Edemsky, D., Popov, A., Prokopovich, I., and Garbatsevich, V.: Airborne Ground Penetrating Radar, Field Test, Remote Sens.Basel, 13, 667, https://doi.org/10.3390/rs13040667, 2021. a
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 laboratoryscale porous media, Water Resour. Res., 47, W02521, https://doi.org/10.1029/2010WR009639, 2011. a, b
Gelman, A. and Rubin, D. B.: Inference from Iterative Simulation Using Multiple Sequences, Stat. Sci., 7, 457–472, http://www.jstor.org/stable/2246093 (last access: 7 March 2023), 1992. a
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 J., 2, 476–491, https://doi.org/10.2136/vzj2003.4760, 2003. a, b, c
Jadoon, K. Z., Slob, E., Vanclooster, M., Vereecken, H., and Lambot, S.: Uniqueness and stability analysis of hydrogeophysical inversion for timelapse groundpenetrating radar estimates of shallow soil hydraulic properties, Water Resour. Res., 44, W09421, https://doi.org/10.1029/2007WR006639, 2008. a
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 TimeLapse GroundPenetrating Radar Data, Vadose Zone J., 11, vzj2011.0177, https://doi.org/10.2136/vzj2011.0177, 2012. a
Jaumann, S. and Roth, K.: Soil hydraulic material properties and layered architecture from timelapse GPR, Hydrol. Earth Syst. Sci., 22, 2551–2573, https://doi.org/10.5194/hess2225512018, 2018. a, b
Jonard, F., Weihermüller, L., Schwank, M., Jadoon, K. Z., Vereecken, H., and Lambot, S.: Estimation of Hydraulic Properties of a Sandy Soil Using GroundBased Active and Passive Microwave Remote Sensing, IEEE T. Geosci. Remote Sensing, 53, 3095–3109, https://doi.org/10.1109/TGRS.2014.2368831, 2015. a, b
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 J., 4, 1048–1058, https://doi.org/10.2136/vzj2004.0140, 2005. a
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 J., 17, 180052, https://doi.org/10.2136/vzj2018.03.0052, 2018. a, b, c
Kodešová, R., Gribb, M. M., and Šimůnek, J.: Estimating soil hydraulic properties from transient cone permeameter data, Soil Sci., 163, 436–453, 1998. a
Laloy, E. and Vrugt, J. A.: Highdimensional posterior exploration of hydrologic models using multipletry DREAM(ZS) and highperformance computing, Water Resour. Res., 48, W01526, https://doi.org/10.1029/2011WR010608, 2012. a
Léger, E., Saintenoy, A., and Coquet, Y.: Hydrodynamic parameters of a sandy soil determined by groundpenetrating radar inside a single ring infiltrometer, Water Resour. Res., 50, 5459–5474, https://doi.org/10.1002/2013WR014226, 2014. a, b, c, d, e, f, g
Léger, E., Saintenoy, A., Tucholka, P., and Coquet, Y.: Hydrodynamic Parameters of a Sandy Soil Determined by GroundPenetrating Radar Monitoring of Porchet Infiltrations, IEEE J. Sel. Top. Appl., 9, 188–200, https://doi.org/10.1109/JSTARS.2015.2464231, 2016. a, b
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, J. Appl. Geophy., 182, 104176, https://doi.org/10.1016/j.jappgeo.2020.104176, 2020. a, b
Mara, T. and Tarantola, S.: Application of Global Sensitivity Analysis of Model Output to Building Thermal Simulations, Build. Simul.China, 1, 290–302, https://doi.org/10.1007/s1227300881295, 2008. a
Moysey, S. M.: Hydrologic trajectories in transient groundpenetratingradar reflection data, GEOPHYSICS, 75, WA211–WA219, https://doi.org/10.1190/1.3463416, 2010. a, b
Mualem, Y.: A new model for predicting the hydraulic conductivity of unsaturated porous media, Water Resour. Res., 12, 513–522, https://doi.org/10.1029/WR012i003p00513, 1976. a, b, c
Muntz, A., Faure, L., and Laine, E.: Études sur la perméabilité des terres, faites en vue de l'arrosage, Ann. de la Direction de l'Hydraulique, f33, 45–53, 1905. a
Rawlins, S. L. and Campbell, G. S.: Water Potential: Thermocouple Psychrometry, Chap. 24, in: Methods of Soil Analysis, edited by: Klute, A., John Wiley & Sons, Ltd, 597–618, https://doi.org/10.2136/sssabookser5.1.2ed.c24, 1986. a
Richards, L. A.: CAPILLARY CONDUCTION OF LIQUIDS THROUGH POROUS MEDIUMS, Physics, 1, 318–333, https://doi.org/10.1063/1.1745010, 1931. a
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 WatershedScale Observatories: A Review, Vadose Zone J., 7, 358–389, https://doi.org/10.2136/vzj2007.0143, 2008. a, b
Saintenoy, A. and Hopmans, J. W.: Ground Penetrating Radar: Water Table Detection Sensitivity to Soil Water Retention Properties, IEEE J. Sel. Top. Appl., 4, 748–753, https://doi.org/10.1109/JSTARS.2011.2171920, 2011. a
Saintenoy, A., Schneider, S., and Tucholka, P.: Evaluating GroundPenetrating Radar use for water infiltration monitoring, Vadose Zone J., 7, 208–214, https://hal.archivesouvertes.fr/hal00831408 (last access: 7 March 2023), 2008. a
Saito, H., Kuroda, S., Iwasaki, T., Fujimaki, H., Nagai, N., and Sala, J.: Tracking Infiltration Front Depth Using Timelapse Multioffset Gathers Collected with Array Antenna Ground Penetrating Radar, JoVEJ. Vis. Exp., 2018, https://doi.org/10.3791/56847, 2018. a, b
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, Hydrol. Earth Syst. Sci., 15, 3043–3059, https://doi.org/10.5194/hess1530432011, 2011. a
Scholer, M., Irving, J., Binley, A., and Holliger, K.: Estimating vadose zone hydraulic properties using ground penetrating radar: The impact of prior information, Water Resour. Res., 47, W10512, https://doi.org/10.1029/2011WR010409, 2011. a, b
Shao, Q., Younes, A., Fahs, M., and Mara, T. A.: Bayesian sparse polynomial chaos expansion for global sensitivity analysis, Comput. Method. Appl. M., 318, 474–496, https://doi.org/10.1016/j.cma.2017.01.033, 2017. a, b
Sobol', I.: Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates, Math. Comput. Simulat., 55, 271–280, https://doi.org/10.1016/S03784754(00)002706, 2001. a, b
Tran, A. P., Vanclooster, M., Zupanski, M., and Lambot, S.: Joint estimation of soil moisture profile and hydraulic parameters by groundpenetrating radar data assimilation with maximum likelihood ensemble filter, Water Resour. Res., 50, 3131–3146, https://doi.org/10.1002/2013WR014583, 2014. a
van Genuchten, M. T.: A Closedform Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils, Soil Sci. Soc. Am. J., 44, 892–898, https://doi.org/10.2136/sssaj1980.03615995004400050002x, 1980. a, b, c
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 Resour. Res., 44, W00D06, https://doi.org/10.1029/2008WR006829, 2008. a, b, c
Vrugt, J. A.: Markov chain Monte Carlo simulation using the DREAM software package: Theory, concepts, and MATLAB implementation, Environ. Model. Softw., 75, 273–316, https://doi.org/10.1016/j.envsoft.2015.08.013, 2016. a
Vrugt, J. A. and Bouten, W.: Validity of FirstOrder Approximations to Describe Parameter Uncertainty in Soil Hydrologic Models, Soil Sci. Soc. Am. J., 66, 1740–1751, https://doi.org/10.2136/sssaj2002.1740, 2002. a
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 J., 7, 843–864, https://doi.org/10.2136/vzj2007.0078, 2008. a
Wiener, N.: The Homogeneous Chaos, Am. J. Math., 60, 897–936, http://www.jstor.org/stable/2371268 (last access: 29 August 2022), 1938. a
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 J., 12, vzj2011.0150, https://doi.org/10.2136/vzj2011.0150, 2013. a
Younes, A., Delay, F., Fajraoui, N., Fahs, M., and Mara, T.: Global sensitivity analysis and Bayesian parameter inference for solute transport in porous media colonized by biofilms, J. Contam. Hydrol., 191, 1–18, https://doi.org/10.1016/j.jconhyd.2016.04.007, 2016. a
Younes, A., Mara, T., Fahs, M., Grunberger, O., and Ackerer, P.: Hydraulic and transport parameter assessment using column infiltration experiments, Hydrol. Earth Syst. Sci., 21, 2263–2275, https://doi.org/10.5194/hess2122632017, 2017. a
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, Hydrol. Earth Syst. Sci., 22, 3561–3574, https://doi.org/10.5194/hess2235612018, 2018. a
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, Front. Plant Sci., 10, 1270, https://doi.org/10.3389/fpls.2019.01270, 2019. a
 Abstract
 Introduction
 Test case description and numerical solution
 Global sensitivity analysis of TWT signals
 Bayesian soil parameter estimation from the TWT signals
 Conclusions
 Code availability
 Data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 Test case description and numerical solution
 Global sensitivity analysis of TWT signals
 Bayesian soil parameter estimation from the TWT signals
 Conclusions
 Code availability
 Data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References