Articles | Volume 27, issue 23
Research article
07 Dec 2023
Research article |  | 07 Dec 2023

Coupled hydrogeophysical inversion of an artificial infiltration experiment monitored with ground-penetrating radar: synthetic demonstration

Rohianuu Moua, Nolwenn Lesparre, Jean-François Girard, Benjamin Belfort, François Lehmann, and Anis Younes

In this study, we investigate the use of ground-penetrating radar (GPR) time-lapse 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 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 Ks, 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. Ks 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 Ks, 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.

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, 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 non-linear 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 (Mualem1976, van Genuchten1980) 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 labor-intensive for deep- or large-scale 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 Klute1986) or psychrometers (Rawlins and Campbell1986), 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, time-consuming, 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 non-invasive 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 vehicle-based 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 contact-based 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 Brevik2014), direct current resistivity (de Jong et al.2020), nuclear magnetic resonance (Costabel and Günther2014), 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 time-lapse 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. Time-lapse 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, time-lapse GPR monitoring of artificial infiltration experiments is usually effortless and time-saving. Furthermore, except in the case of borehole investigations, the GPR device can be laid on or raised above the surface. For these reasons, time-lapse 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 non-destructive. 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 time-lapse GPR in the context of evaluating the soil hydraulic states, hydraulic properties, or unsaturated soil parameters (e.g., Saintenoy et al.2008; Moysey2010; Scholer et al.2011; Busch et al.2013; Tran et al.2014; Jonard et al.2015; Jaumann and Roth2018; 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 imbibition-drainage experiments using a single-offset 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 two-way 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), single-offset (Léger et al.2014; Léger et al.2016) and multi-offset (Saito et al.2018) surface GPR, or off-ground GPR (Jadoon et al.2008, Jadoon et al.2012; Jonard et al.2015). For practicality, surface GPR is preferred over off-ground and borehole GPR, the latter also being destructive by nature. Saito et al. (2018) used a more complex multi-offset and multi-channel surface GPR to directly monitor the wetting front progression. The 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 of being multi-offset and is, therefore, able to simultaneously determine the propagation speed and the depths of reflectors.

Figure 1Test 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, respectively, of the GPR system; effective saturation Se (b) and reflection coefficient r (c) profiles with depth.


In the present study, we are interested in using a quick, easy-to-apply, and cheap field-scale method to characterize the unsaturated soil parameters. To this end, time-lapse GPR monitoring of artificial infiltration is a well-suited 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 optimization-based 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 time-lapse 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 Test case description and numerical solution

2.1 Test case description

In this work, we conduct a synthetic study on the time-lapse 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 one-dimensional 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 double-ring 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 Dirichlet-type 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 time-lapse TWT signal for reflection caused by the infiltration wetting front is noted TWTf and that from the two immovable diffracting points R50 and R120 are, respectively, noted TWT50 and TWT120.

2.2 The mathematical model

2.2.1 Unsaturated flow model

Water infiltration in unsaturated or saturated soils is governed by the one-dimensional Richards equation (Richards1931):

(1) θ t = z K ( θ ) h z - 1 ,

where h (cm) is the pressure head; z is the depth (cm), taken as positive in the downward direction; t is the time (s), θ (cm3 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 Mualem1976 and van Genuchten (1980):

(2)Se(h)=θ(h)-θrθs-θr=[1+(α|h|)n]-m if h<01 if h0(3)K(h)=Ks×Se(h)L[1-(1-Se(h)1/m)m]2 if h<0Ks if h0,

where Se(h) (–) is the effective saturation, θs and θr (cm3 cm−3) are the saturated and residual water contents, respectively, Ks (cm s−1) is the saturated conductivity, m=1-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 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 ϵ (–) of a three-phase (water–solid–air) medium to its water content by:

(4) ϵ ( z , t ) = θ ( z , t ) ϵ w + ( 1 - ϕ ) ϵ s + [ ϕ - θ ( z , t ) ] ϵ a ,

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 non-magnetic 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) (Annan2003):

(5) V = c ϵ ,

where c≈30cm 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 one-dimensional scenario (the offset between the antennas is null) and discretize the domain into N cells i, centered at a depth zi, with element boundaries at zi-1/2 and zi+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:

(6) TWT ( z i - 1 / 2 ) = 2 j = 1 i - 1 | l j | V j ,

in which |lj| (cm) is the length of the element j above i, and Vj (cm ns−1) is the GPR propagation speed in the element j.

Table 1Prior intervals of the unsaturated soil parameters for both GSA and Bayesian estimation.

Download Print Version | Download XLSX

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:

(7) r ( z i - 1 / 2 ) = ϵ ( z i ) - ϵ ( z i - 1 ) ϵ ( z i ) + ϵ ( z i - 1 ) ,

where ϵ(zi) 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 (Bano2006; Saintenoy and Hopmans2011).

Note that one could easily consider a non-perpendicular 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 WAMOS-1D 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 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 ϵ using the petrophysical relationship in Eq. (4) and into a GPR wave propagation speed profile V using Eq. (5). Then, the time-lapse TWT signals for the fixed objects, TWT50 and TWT120, are calculated at each time step using Eq. (6) (dashed and dotted curves in Fig. 2).

Figure 2Hydrogeophysical model responses for Ks=0.08cm s−1, θs=0.4 (–), θr=0.07 (–), α=0.145cm−1, n=2.68 (–). TWTf corresponds to the TWT signal for the wetting front, while TWT50 and TWT120 are the TWT signals for fixed objects located at 50 and 120 cm below the surface, respectively.


The time-lapse signal TWTf, 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 zi-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 TWTf=TWT(zi-1/2) from Eq. (6) (solid curve Fig. 2).

Figure 3Summary of the working process of the forward hydrogeophysical model and how it is used to build the PCE surrogate model.


Note that TWT50 and TWT120 signals are induced by fixed objects; thus, these signals exist regardless of the position of the infiltration front. On the other hand, TWTf is induced by the infiltration wetting front whose position varies over time. Besides, contrarily to TWT50 and TWT120, the TWTf 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 TWTf 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 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 Tarantola2008). Among the various forms of GSA, a variance-based 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 (Ks, θs, θr, α, n) and the output variables are the TWT signals (TWTf, TWT50, TWT120).

Given a model with a set of p independent random parameters X={x1,x2,,xp} that yields a random response y(X), the two variance-based sensitivity measures, also called Sobol indices (Sobol'2001), are

  • the first-order sensitivity index:

    (8) S i = Var E [ y ( X ) | x i ] Var [ y ( X ) ] [ 0 , 1 ]
  • the total sensitivity index:

    (9) ST i = E Var [ y ( X ) | x - i ] Var [ y ( X ) ] [ 0 , 1 ] ,

where x-i=Xxi is the set of all parameters except xi, 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 first-order index Si quantifies the contribution of the parameter xi alone to the total variance of y(X), while STi also includes all interactions of xi with the other parameters xi.

To perform a variance-based GSA, a practical approach (to save computational time) is to use polynomial chaos expansion (PCE; Wiener1938). The PCE approach consists in developing any signal y(X) as a set of orthonormal multivariate polynomials of a degree not exceeding D:

(10) y ( X ) = | β | D s β Ψ β ( X ) ,

where β=β1,β2,,βpRp is a pth-dimensional index, sβ are the PC coefficients, and Ψβ are the generalized polynomial chaos of degree |β|=i=1pβ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)!/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 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 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). Because each parameter varies in its own range and has a proper unit, the parameter prior intervals are normalized to [-1,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 (TWTf, TWT50, and TWT120) 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 TWT50 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 ns2). 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.

Table 2Variance of TWTf, TWT50, and TWT120 signals at t= 50, 150, and 2000 s calculated with the PCE surrogate model and with the hydrogeophysical model.

Download Print Version | Download XLSX

Figure 4Time distribution of the variance of TWTf (a), TWT50 (b) and TWT120 (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, t1= 1 min, t2= 5 min, and t3= 200 min, highlighted here (dotted black lines).


3.2 GSA results

The temporal distribution of the output variance of the three TWT signals (TWTf, TWT50, and TWT120) 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.

TWTf has a different behavior from the TWT signals of fixed reflectors TWT50 and TWT120 (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 TWT120 than for TWT50. This is in agreement with the physics since the zone of the porous medium affecting the GPR wave is more important for the TWT120 signal than for the TWT50. In addition, the period of influence of the unsaturated parameters (θr, α, n) is also more important for TWT120 than for TWT50 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 TWTf and TWT120 signals.

Figure 5Difference between the total (STi) and the first-order (Si) sensitivity indices for all parameters for the TWTf (a) and the TWT120 (b) signals.


3.2.1 GSA of the TWTf signal

TWTf variance is zero at the beginning of the infiltration (Fig. 4a), which means that the TWTf signal is not affected by the initial conditions. Indeed, the infiltration wetting front and the TWTf signal start at zero for all parameter sets. Then, the variance of the signal increases until a maximum of 60 ns2, reached at around 3 min. After that, the variance decreases, but maintains a significant value of around 25 ns2 (Fig. 4a). Concerning parameter sensitivities, at the beginning, the TWTf signal is mainly affected by Ks. 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 θs has a moderate influence on the TWTf signal. Its influence is not observable at short times since unsaturated conditions occur. Overall, the most influential parameter on the TWTf 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 TWTf 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 (STi) and the first-order (Si) 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 (STiSi), except for Ks and α (Fig. 5a). Hence, the interaction between these two parameters affects the variance of the TWTf signal as represented by the blank region between the variance curve and the shaded area (Fig. 4a).

Figure 6Marginal effects of the unsaturated soil parameters Ks, θs, θr, α, and n on the TWTf and TWT120 signals at three different times, t1= 1 min, t2= 5 min, and t3= 200 min, highlighted in Fig. 4.


To evaluate further the effect of the unsaturated soil parameters on the TWTf, 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 (t1= 1 min, t2= 5 min, and t3= 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:

  • Ks is highly influential at the beginning of the experiment. At t1= 1 min, the TWTf signal varies almost linearly with Ks. Indeed, at the beginning of the experiment, when Ks increases, the wetting front is more advanced, thus, the GPR wave propagates at a lower speed and the TWTf signal increases. At t2= 5 min, the TWTf signal is sensitive only for small Ks values. Indeed, for high Ks values, the soil is fully saturated and the perturbation of the high value of Ks does not change the TWTf signal. At t3= 200 min, the soil is fully saturated for almost all Ks values, thus, the TWTf signal becomes insensitive to Ks.

  • θs has no influence at the first times (t1= 1 min) since unsaturated conditions occur. For long times, the TWTf signal is very sensitive to θs with an almost linear behavior. Indeed, when the soil is fully saturated, the dielectric permittivity and thus the TWTf signal is almost proportional to θs.

  • The sensitivity of TWTf 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 TWTf signal increases).

  • The van Genuchten parameter α is highly influential notably for long times (t3= 200 min). A small variation of the parameter α can induce a strong variation of the TWTf signal. Notably, the sensitivity of α is very high for α≤0.05cm−1.

  • The sensitivity of TWTf 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 TWTf signal and, as a consequence, it is expected to be poorly identifiable from the TWTf data. TWTf is sensitive to all Ks, θ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 TWT120 signal

The variance of the TWT120 signal is nonzero at the beginning of the experiment which means that the TWT120 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 Ks. 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 Ks 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 TWTf 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 Ks, α, 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 TWT120 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 TWTf signal, Ks is highly sensitive, especially for t= 1 min and 5 min.

  • The saturated water content θs is very influential for all times. The TWT120 varies almost linearly with θs even at the beginning (t1= 1 min), since the fixed reflector is located in the lower saturated region.

  • As for the TWTf 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 TWT120 increases.

  • The van Genuchten parameter α is highly sensitive. However, contrarily to the TWTf signal where α is highly sensitive at long times (t3= 200 min), the sensitivity of α for the TWT120 signal is high at short times (t1= 1 min). For long times, the influence of α disappears since the soil becomes fully saturated. The negative slope of the curve of the TWT120 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 TWT120 signal decreases.

  • Surprisingly, and contrarily to the TWTf signal which showed a flat curve for the marginal effect of the parameter n for all parameter values and at all investigated times, the TWT120 signal is sensitive to n at the beginning of the experiment (t1= 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, TWT120 shows similar sensitivity for Ks, θs and θr but slightly higher than TWTf. For α, they show complementarity, which makes the procedure very efficient for α < 0.05 using all the infiltration experiments (early time for TWTf and late time for TWT120). TWT120 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.

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 Bouten2002; 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?

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.

Table 3Reference values used to build the synthetic calibration data for the different scenarios (estimated mean values in bold, size of the posterior confidence intervals between parentheses, and ratio of prior to posterior intervals in italics)

Download Print Version | Download XLSX

In the sequel, the MCMC method is performed with the DREAM(ZS) (DiffeRential Evolution Adaptive Metropolis) software (Laloy and Vrugt2012; Vrugt2016). 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 (Ks, θ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, Ks=0.08cm s−1, θs=0.4, θr=0.07, α=0.145cm−1, n=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 TWTf, TWT50, and TWT120 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.5ns. This error corresponds to an uncertainty of 1 ns, which is realistic in the instance of an 800 MHz GPR antenna.

The TWTf, TWT50, and TWT120 calibration signals, illustrated before noise corruption in Fig. 2, increase almost linearly until reaching a plateau. For the TWT50 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 TWT120 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 TWTf, 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 wetting-front TWTf signal are used for the calibration. The second and third scenarios use only the TWT50 and TWT120 signal, respectively, obtained from reflection on the fixed reflector R50 and R120. The fourth scenario uses both data of TWTf and TWT120 as fitting data. The last scenario investigates the benefit of adding a fixed reflector by using data of the TWTf, TWT50, and TWT120 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 Rubin1992)) criterion is verified and the chains are stable and not autocorrelated.

The results in Table 3 for scenario 1 using only data of the TWTf signal for the estimation of the unsaturated soil parameters show the following:

  • An accurate estimation of Ks, 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 TWTf 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 TWTf (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 TWTf signal yields a mean estimated value n=2.75±0.47 which is close to the reference value n=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 TWTf 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 Ks, θ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 TWT50 or TWT120 signal for the calibration, shows the following:

  • The parameters Ks 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 Ks, and 19 for θs, of their intervals of variation. We note that the TWT120 signal allows for a much better estimate of both Ks and θs as their CIs are smaller than the ones estimated with TWT50. It is especially true for Ks 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 TWT50 and TWT120 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 TWTf and TWT120 signals, show the following:

  • Parameters Ks, θ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 7MCMC solutions using scenarios 1 (a1–a5), 3 (b1–b5), and 4 (c1–c5) for calibrating the hydrogeophysical model. The histograms are built from the posterior distributions. The estimated mean values are represented as a dashed black line and compared to the exact target value (solid red line). The displayed parameter intervals correspond to the prior upper and lower limits of Table 1.


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 TWTf, TWT50 and TWT120 signals, show performances very similar to scenario 4. Additional information from TWT50 helped to reduce slightly the posterior intervals of Ks, θs, and α, which in that case show a reduction of 49, 44, and 10 times their prior intervals, respectively.

Figure 8MCMC solutions using TWTf, TWT50, and TWT120 signals for the calibration of the hydrogeophysical model. The diagonal plots represent the inferred posterior parameter distributions, showing the estimated mean value (dashed black line) and the target value (solid red line). The off-diagonal represents the pairwise correlations between parameters.


The results of MCMC for this last scenario are shown in Fig. 8 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 Ks and θr (r=-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 TWTf and TWT120 signals yields a much poorer identification of the parameter n. For instance, 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.

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.

Figure 9MCMC solutions considering a water table at 50 cm, 1 m, and 2 m depth (from top to bottom) for calibrating the hydrogeophysical model. The histograms are built from the posterior distributions. The estimated mean values are represented as a dashed black line and compared with the exact target value (solid red line). The displayed parameter intervals are equal from each parameter (each column).


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 TWTf and TWT120 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.

Table 4Hydraulic parameters for three types of sand. X shows the target values while the estimation results are represented below for each sand (estimated mean values in bold, size of the posterior CIs between parentheses, and ratio of prior to posterior intervals in italics).

Download Print Version | Download XLSX

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.

5 Conclusions

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: time-lapse GPR monitoring of artificial infiltration experiments. Water infiltration into an initially unsaturated sandy soil was simulated using a one-dimensional 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 TWTf 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 TWTf signal is highly sensitive to Ks and α and moderately sensitive to θs. A low sensitivity was observed for θr, whereas the parameter n was insensitive. The TWT120 signal of the fixed reflector located at 120 cm depth is highly sensitive to Ks, θ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 TWTf 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 Ks, θs, and n were correctly identified (with narrow posterior intervals). Further, using only data of the TWT50 or TWT120 signals for the calibration enabled also good identification of Ks 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 TWTf with at least one fixed reflector. In this case, the four parameters Ks, θ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.

Code availability

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

Data availability

All data presented in this paper are available from the authors upon request.

Author contributions

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 co-authors. RM prepared the manuscript with reviewing and editing contributions from all co-authors.

Competing interests

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 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 Alsacalcul/Big Data. Computing time was provided by the HPC-UdS.

Financial support

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 (grant no. ANR-17-CE06-0012-01).

Review statement

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, 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,, 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 1D-Model Calibration, Geosciences, 8, 72,, 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,, 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,, 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,, 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 time-lapse surface GPR data to estimate hydraulic properties of a layered subsurface, Water Resour. Res., 49, 8480–8494,, 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,, 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,, 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,, 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,, 2020. a

Doolittle, J. A. and Brevik, E. C.: The use of electromagnetic induction techniques in soils studies, Geoderma, 223–225, 33–45,, 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,, 2015. a

Edemsky, D., Popov, A., Prokopovich, I., and Garbatsevich, V.: Airborne Ground Penetrating Radar, Field Test, Remote Sens.-Basel, 13, 667,, 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 laboratory-scale porous media, Water Resour. Res., 47, W02521,, 2011. a, b

Gelman, A. and Rubin, D. B.: Inference from Iterative Simulation Using Multiple Sequences, Stat. Sci., 7, 457–472, (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,, 2003. a, b, c

Jadoon, K. Z., Slob, E., Vanclooster, M., Vereecken, H., and Lambot, S.: Uniqueness and stability analysis of hydrogeophysical inversion for time-lapse ground-penetrating radar estimates of shallow soil hydraulic properties, Water Resour. Res., 44, W09421,, 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 Time-Lapse Ground-Penetrating Radar Data, Vadose Zone J., 11, vzj2011.0177,, 2012. a

Jaumann, S. and Roth, K.: Soil hydraulic material properties and layered architecture from time-lapse GPR, Hydrol. Earth Syst. Sci., 22, 2551–2573,, 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 Ground-Based Active and Passive Microwave Remote Sensing, IEEE T. Geosci. Remote Sensing, 53, 3095–3109,, 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,, 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,, 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.: High-dimensional posterior exploration of hydrologic models using multiple-try DREAM(ZS) and high-performance computing, Water Resour. Res., 48, W01526,, 2012. a

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 Resour. Res., 50, 5459–5474,, 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 Ground-Penetrating Radar Monitoring of Porchet Infiltrations, IEEE J. Sel. Top. Appl., 9, 188–200,, 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,, 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,, 2008. a

Moysey, S. M.: Hydrologic trajectories in transient ground-penetrating-radar reflection data, GEOPHYSICS, 75, WA211–WA219,, 2010. a, b

Mualem, Y.: A new model for predicting the hydraulic conductivity of unsaturated porous media, Water Resour. Res., 12, 513–522,, 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,, 1986. 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 Watershed-Scale Observatories: A Review, Vadose Zone J., 7, 358–389,, 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,, 2011. a

Saintenoy, A., Schneider, S., and Tucholka, P.: Evaluating Ground-Penetrating Radar use for water infiltration monitoring, Vadose Zone J., 7, 208–214, (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 Time-lapse Multi-offset Gathers Collected with Array Antenna Ground Penetrating Radar, JoVE-J. Vis. Exp., 2018,, 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,, 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,, 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,, 2017. a, b

Sobol', I.: Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates, Math. Comput. Simulat., 55, 271–280,, 2001. a, b

Tran, A. P., Vanclooster, M., Zupanski, M., and Lambot, S.: Joint estimation of soil moisture profile and hydraulic parameters by ground-penetrating radar data assimilation with maximum likelihood ensemble filter, Water Resour. Res., 50, 3131–3146,, 2014. a

van Genuchten, M. T.: A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils, Soil Sci. Soc. Am. J., 44, 892–898,, 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,, 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,, 2016. a

Vrugt, J. A. and Bouten, W.: Validity of First-Order Approximations to Describe Parameter Uncertainty in Soil Hydrologic Models, Soil Sci. Soc. Am. J., 66, 1740–1751,, 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,, 2008. a

Wiener, N.: The Homogeneous Chaos, Am. J. Math., 60, 897–936, (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,, 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,, 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,, 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,, 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,, 2019. a

Short summary
Hydraulic properties of soil include the ability of water to move through the soil and the amount of water that is held in the soil in dry or wet conditions. In this work, we further investigate a protocol used to evaluate such hydraulic properties. We propose a modified version of the protocol, with which we show (i) how the data obtained with this protocol are influenced by the soil hydraulic properties and (ii) how one can use it to estimate these properties.