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

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 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.

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 (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 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 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, 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 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 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; 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 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 1 Test case and experimental device illustration at an advanced time step (a). R50 and R120 are fixed reflectors considered in this experiment. Tx and Rx refer to the transmitter and receiver antennas, respectively, of the GPR system; effective saturation S_e (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 TWT_f and that from the two immovable diffracting points R50 and R120 are, respectively, noted TWT₅₀ and TWT₁₂₀.

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 (Richards, 1931):

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial \mathit{\theta }}{\partial t}}={\displaystyle \frac{\partial }{\partial z}}\left[K\left(\mathit{\theta }\right)\left({\displaystyle \frac{\partial h}{\partial z}}-\mathrm{1}\right)\right],\end{array}$$

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³ cm⁻³) is the actual water content, and K(θ) (cm s⁻¹) is the hydraulic conductivity which is a function of water content. The initial condition is a hydrostatic pressure distribution corresponding to a water table at 100 cm depth. The boundary condition at the top of the domain is a fixed Dirichlet condition of 10 cm maintained during the experiment. The boundary condition at the bottom is a piezometric head fixed at −100 cm which corresponds to the water table position (Fig. 1).

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

$$\begin{array}{}\text{(2)}& {\displaystyle }{S}_{\mathrm{e}}\left(h\right)={\displaystyle \frac{\mathit{\theta }\left(h\right)-{\mathit{\theta }}_{\mathrm{r}}}{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}}}=\left\{\begin{array}{ll}[\mathrm{1}+(\mathit{\alpha }\left|h\right|{)}^n{]}^{-m}& \text{ if }h<\mathrm{0}\\ \mathrm{1}& \text{ if }h\ge \mathrm{0}\end{array}\right.\text{(3)}& {\displaystyle }K\left(h\right)=\left\{\begin{array}{ll}{K}_{\mathrm{s}}\times S_{\mathrm{e}}\left(h{)}^L\right[\mathrm{1}-(\mathrm{1}-S_{\mathrm{e}}(h{)}^{\mathrm{1}/m}{)}^m{]}^{\mathrm{2}}& \text{ if }h<\mathrm{0}\\ K_{\mathrm{s}}& \text{ if }h\ge \mathrm{0},\end{array}\right.\end{array}$$

where S_e(h) (–) is the effective saturation, θ_s and θ_r (cm³ cm⁻³) are the saturated and residual water contents, respectively, K_s (cm s⁻¹) is the saturated conductivity, $m=\mathrm{1}-\mathrm{1}/n$, α (1 cm⁻¹), 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:

$$\begin{array}{}\text{(4)}& \sqrt{\mathit{ϵ}(z,t)}=\mathit{\theta }(z,t)\sqrt{{\mathit{ϵ}}_w}+(\mathrm{1}-\mathit{\varphi })\sqrt{{\mathit{ϵ}}_s}+[\mathit{\varphi }-\mathit{\theta }(z,t\left)\right]\sqrt{{\mathit{ϵ}}_a},\end{array}$$

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⁻¹) (Annan, 2003):

$$\begin{array}{}\text{(5)}& V={\displaystyle \frac{c}{\sqrt{\mathit{ϵ}}}},\end{array}$$

where c≈30 cm ns⁻¹ 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 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:

$$\begin{array}{}\text{(6)}& \text{TWT}\left(z_{i-\mathrm{1}/\mathrm{2}}\right)=\mathrm{2}\sum _{j=\mathrm{1}}^{i-\mathrm{1}}{\displaystyle \frac{\left|l_j\right|}{V_j}},\end{array}$$

in which $\left|l_j\right|$ (cm) is the length of the element j above i, and V_j (cm ns⁻¹) is the GPR propagation speed in the element j.

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

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:

$$\begin{array}{}\text{(7)}& r\left(z_{i-\mathrm{1}/\mathrm{2}}\right)={\displaystyle \frac{\mathit{ϵ}\left(z_i\right)-\mathit{ϵ}\left(z_{i-\mathrm{1}}\right)}{\mathit{ϵ}\left(z_i\right)+\mathit{ϵ}\left(z_{i-\mathrm{1}}\right),}}\end{array}$$

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 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, TWT₅₀ and TWT₁₂₀, are calculated at each time step using Eq. (6) (dashed and dotted curves in Fig. 2).

Figure 2 Hydrogeophysical model responses for K_s=0.08 cm s⁻¹, θ_s=0.4 (–), θ_r=0.07 (–), α=0.145 cm⁻¹, n=2.68 (–). TWT_f corresponds to the TWT signal for the wetting front, while TWT₅₀ and TWT₁₂₀ are the TWT signals for fixed objects located at 50 and 120 cm below the surface, respectively.

The time-lapse signal TWT_f, induced by wave reflection on the wetting front because of the sharp water content variation at the front position, is calculated in two steps. First, we search 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).

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

Note that TWT₅₀ and TWT₁₂₀ 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₅₀ and TWT₁₂₀, 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⁻². This reflects a fully saturated domain with an almost uniform water content distribution (solid curve Fig. 2). An explanation of the computation of all TWT signals is summarized in Fig. 3.

3 Global sensitivity analysis of TWT signals

3.1 GSA method

The GSA method evaluates how the outputs of a model are influenced by the variation of the input parameters (Mara and Tarantola, 2008). Among the various forms of GSA, a variance-based sensitivity analysis, allowing 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₅₀, TWT₁₂₀).

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 variance-based sensitivity measures, also called Sobol indices (Sobol', 2001), are

• the first-order 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 first-order 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 variance-based GSA, a practical approach (to save computational time) is to use polynomial chaos expansion (PCE; Wiener, 1938). The PCE approach consists in developing any signal y(X) as a set of orthonormal multivariate polynomials of a degree not exceeding D:

$$\begin{array}{}\text{(10)}& y\left(\mathit{X}\right)=\sum _{\left|\mathit{\beta }\right|\le D}{s}_{\mathit{\beta }}{\mathrm{\Psi }}_{\mathit{\beta }}\left(\mathit{X}\right),\end{array}$$

where $\mathit{\beta }={\mathit{\beta }}_{\mathrm{1}},{\mathit{\beta }}_{\mathrm{2}},\mathrm{\dots },{\mathit{\beta }}_p\in {\mathbb{R}}^p$ is a pth-dimensional 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 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 $[-\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₅₀, and TWT₁₂₀) 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₅₀ 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²). 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 2 Variance of TWT_f, TWT₅₀, and TWT₁₂₀ signals at t = 50, 150, and 2000 s calculated with the PCE surrogate model and with the hydrogeophysical model.

Figure 4 Time distribution of the variance of TWT_f (a), TWT₅₀ (b) and TWT₁₂₀ (c). The shaded area under the variance curve represents the partial marginal contributions of the uncertain parameters; the blank region between the shaded area and the variance curves represents the contribution of interactions between the parameters. The marginal effects shown in Fig. 6 are represented at three time steps, t₁ = 1 min, t₂ = 5 min, and t₃ = 200 min, highlighted here (dotted black lines).

3.2 GSA results

The temporal distribution of the output variance of the three TWT signals (TWT_f, TWT₅₀, and TWT₁₂₀) 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₅₀ and TWT₁₂₀ (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₁₂₀ than for TWT₅₀. This is in agreement with the physics since the zone of the porous medium affecting the GPR wave is more important for the TWT₁₂₀ signal than for the TWT₅₀. In addition, the period of influence of the unsaturated parameters (θ_r, α, n) is also more important for TWT₁₂₀ than for TWT₅₀ 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₁₂₀ signals.

Figure 5 Difference between the total (ST_i) and the first-order (S_i) sensitivity indices for all parameters for the TWT_f (a) and the TWT₁₂₀ (b) 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², reached at around 3 min. After that, the variance decreases, but maintains a significant value of around 25 ns² (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 steady-state 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 first-order (S_i) sensitivity indices for all parameters (Fig. 5a). This difference reflects the interaction between the parameters over time. Interactions between parameters are negligible for all parameters (ST_i≈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).

Figure 6 Marginal effects of the unsaturated soil parameters K_s, θ_s, θ_r, α, and n on the TWT_f and TWT₁₂₀ signals at three different times, t₁ = 1 min, t₂ = 5 min, and t₃ = 200 min, highlighted in Fig. 4.

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 min, t₂ = 5 min, and t₃ = 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 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₂ = 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₃ = 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 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₃ = 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⁻¹.

• 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⁻¹ and n < 2. This means a wide range of soil types.

3.2.2 GSA of the TWT₁₂₀ signal

The variance of the TWT₁₂₀ signal is nonzero at the beginning of the experiment which means that the TWT₁₂₀ 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 first-order 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₁₂₀ 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₁₂₀ varies almost linearly with θ_s even at the beginning (t₁ = 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₁₂₀ increases.

• The van Genuchten parameter α is highly sensitive. However, contrarily to the TWT_f signal where α is highly sensitive at long times (t₃ = 200 min), the sensitivity of α for the TWT₁₂₀ signal is high at short times (t₁ = 1 min). For long times, the influence of α disappears since the soil becomes fully saturated. The negative slope of the curve of the TWT₁₂₀ 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₁₂₀ 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₁₂₀ signal is sensitive to n at the beginning of the experiment (t₁ = 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₁₂₀ 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₁₂₀). TWT₁₂₀ 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 Bouten, 2002; Vrugt et al., 2008). The statistical calibration is performed for a GPR-monitored infiltration experiment in order to address the following questions:

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

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

3. What is the optimal set of TWT measurements that yields good reliability of all unsaturated soil parameters?

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 3 Reference 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)

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⁻¹, ${\mathit{\theta }}_{\mathrm{s}}^{\ast }=\mathrm{0.4}$, ${\mathit{\theta }}_{\mathrm{r}}^{\ast }=\mathrm{0.07}$, ${\mathit{\alpha }}^{\ast }=\mathrm{0.145}$ cm⁻¹, $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₅₀, and TWT₁₂₀ 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₅₀, and TWT₁₂₀ calibration signals, illustrated before noise corruption in Fig. 2, increase almost linearly until reaching a plateau. For the TWT₅₀ 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₁₂₀ 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 wetting-front TWT_f signal are used for the calibration. The second and third scenarios use only the TWT₅₀ and TWT₁₂₀ signal, respectively, obtained from reflection on the fixed reflector R50 and R120. The fourth scenario uses both data of TWT_f and TWT₁₂₀ as fitting data. The last scenario investigates the benefit of adding a fixed reflector by using data of the TWT_f, TWT₅₀, and TWT₁₂₀ 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⁻¹ 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⁻¹ 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₅₀ or TWT₁₂₀ 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₁₂₀ signal allows for a much better estimate of both K_s and θ_s as their CIs are smaller than the ones estimated with TWT₅₀. 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₅₀ and TWT₁₂₀ 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₁₂₀ 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 MCMC 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 TWT_f, TWT₅₀ and TWT₁₂₀ signals, show performances very similar to scenario 4. Additional information from TWT₅₀ 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.

Figure 8 MCMC solutions using TWT_f, TWT₅₀, and TWT₁₂₀ 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 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₁₂₀ 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.

Figure 9 MCMC 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 TWT_f and TWT₁₂₀ 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 4 Hydraulic 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).

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 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₁₂₀ 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₅₀ or TWT₁₂₀ 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.