Identifying a parameterisation of the soil water retention curve from on-ground GPR measurements

We show the potential of on-ground GroundPenetrating Radar (GPR) to identify the parameterisation of the soil water retention curve, i.e. its functional form, with a semi-quantitative analysis based on numerical simulations of the radar signal. An imbibition and drainage experiment has been conducted at the ASSESS-GPR site to establish a fluctuating water table, while an on-ground GPR antenna recorded traces over time at a fixed location. These measurements allow to identify and track the capillary fringe in the soil. The typical dynamics of soil water content with a transient water table can be deduced from the recorded radargrams. The characteristic reflections from the capillary fringes in model soils that are described by commonly used hydraulic parameterisations are investigated by numerical simulations. The parameterisations used are (i) full van Genuchten, (ii) simplified van Genuchten with m = 1− 1 n and (iii) Brooks–Corey. All three yield characteristically different reflections, which allows the identification of an appropriate parameterisation by comparing to the measured signals. We show that for the sand used here, these signals are not consistent with the commonly used simplified van Genuchten parameterisation withm = 1− 1 n .


Introduction
Parameterising the soil hydraulic properties is essential for modelling and thus predicting water movement in soils.The strategies for estimating the corresponding hydraulic parameters can be divided into two main lines: (i) laboratory methods using soil samples and (ii) inversion of field measurements.In general both strategies concentrate on using inverse methods, where a given set of hydraulic parameters is adjusted such that modelled and corresponding measured data are in optimal agreement.The laboratory methods allow for a rather precise estimation of hydraulic properties of soil samples with methods ranging from Multi-Step Outflow experiments (e.g.Van Dam et al., 1994;Vereecken et al., 1997) to evaporation experiments (e.g.Simunek et al., 1998b;Schindler et al., 2010) and to combinations of the two methods (e.g.Iden et al., 2011).The main issue with lab methods is the applicability of the thus gained descriptions to the field for which the soil samples are to be representative.This gave rise to strategy (ii) of estimating hydraulic properties directly at the field scale.It has been demonstrated by a variety of methods, which involved artificial forcing (three of them compared by Simunek et al., 1998a) or the analysis of long time series with natural forcing (e.g.Wollschläger et al., 2009).A main aspect in both strategies is the requirement of a suitable parameterisation, which is capable of representing the soil of interest.Mostly used in this context are the parameterisations of van Genuchten (1980) and Brooks (1966) for the soil water characteristic, typically combined with the parameterisation of Mualem (1976) for the soil hydraulic conductivity.Russo (1988) compared the suitability of different parameterisations from an experimental perspective.Ippisch et al. (2006) showed from a theoretical point of view that the parameter space has to be limited, if the mathematical formulations are to represent physical feasible systems.Hence, the identification of an appropriate parameterisation is essential for estimating hydraulic properties.
GPR is a powerful non-invasive measurement instrument.It is already widely used as a method to investigate architecture and average water content of soils.In the last years, the precision of estimating these quantities could be increased by optimised inversion methods with quite different approaches.
For instance, Buchner et al. (2012) developed a more precise method for the inversion of time and amplitude information of on-ground GPR and Lambot et al. (2004b) established a method using amplitude and phase information from monostatic off-ground GPR.Choosing a different way, van der Kruk et al. (2010) demonstrated that low-velocity layers induced by precipitation events act as a wave-guide and allow for very precise estimations of thickness and water content of these layers.With this GPR emerges as a valuable tool in soil hydrology to monitor the dynamics of water content.This was demonstrated with a 1.5 yr time series with natural forcing by Steelman and Endres (2012).GPR is capable of monitoring the movement of infiltration fronts as shown by Saintenoy et al. (2008) for an infiltration through a tube within the sand and by Moysey (2010) for an infiltration from the surface.Also the capillary fringe can be tracked, which was demonstrated by Endres et al. (2000) during monitoring a pumping test and was investigated further by Tsoflias et al. (2001) regarding shape and amplitude of the reflections.While it could be demonstrated that GPR reflection signals can provide enough information to estimate the hydraulic parameters of a given parameterisation (most precisely by Lambot et al., 2004aLambot et al., , 2009)), we here focus on another problem in hydraulic parameter estimation already mentioned above: the identification of an appropriate parameterisation.This is feasible because of the high sensitivity of GPR reflection signals to small differences in the water content distribution and opens the door to more accurate non-invasive measurements of soil hydraulic properties of the field scale.

Experimental setup
The ASSESS-GPR test site is a cuboid (19 m × 4 m × 1.9 m) consisting of different well-known sand layers which are uniform in the short horizontal direction (Buchner et al., 2012).A gravel layer is placed at the bottom 0.1 m to allow for uniform imbibition and drainage via a pumping tube where the water table height is measured as well.The tube has a diameter of 0.47 m and is only perforated over the extend of the gravel.Since the hydraulic conductivity of the gravel can be assumed to be large compared to the one of the sands, the water is expected to imbibe uniformly over the complete area.
In the following we focus on a 4 m long section of the site that consists of two layers plus another compaction horizon formed during the construction of the testbed (Fig. 1).
During the experiment, a stationary GPR antenna was installed at the surface and was set to record a trace every 3 s.We employed a shielded bi-static antenna (Ingegneria dei Sistemi S.p.A., Italy) with a nominal centre frequency of 400 MHz.The transmitter-receiver distance within the antenna box is 0.14 m.
Starting at a water table height of 0.47 m, the water table was altered by imbibing a total amount of 3800 L within 2 h, corresponding to an imbibition rate of 0.0238 m h −1 .Afterwards the system equilibrated for 2 h prior to draining approximately the same amount of water, i.e. 3868 L within 1.5 h (0.0242 m h −1 ).With the given cross-sectional area of the tank, 80 m 2 , and by assuming a porosity between 0.3 and 0.4, this led to a maximum amplitude of the water table oscillation of 0.12-0.16m.GPR traces were recorded during imbibition and drainage and additionally 10-15 min before and after that (Fig. 2).

Empirical results
We consider time-series of GPR traces obtained from stationary antennas and also refer to them as "radargrams".The individual traces went through minimal postprocessing which consisted of just a dewow filter, which is given as the remainder of a runmean hat-filter with a width of 5 ns.

Imbibition
Figure 2a shows an excerpt of the radargram recorded during the imbibition period, focussing on travel times between 16 to 26 ns to monitor the dynamics of interest.
We expect reflections from boundaries between different materials and structural differences.We can observe both in the experiment, the latter in the very special case of a compaction.At the beginning (in equilibrium state) there are two sharp reflections, one from the upper layer boundary at 0.90 m depth (A), the other one from the compaction horizon at 1.20 m depth (B1).The reflection wavelets consist of three significant extrema plus a forth weaker one only visible for high intensities.
Some minutes after starting the imbibition the reflection at (B1) is separating into two reflections, one moving upwards (C1) and one downwards (B2), corresponding to shorter and longer travel times, respectively.The upcoming reflection (C1) has a significantly different shape compared to a reflection from a layer boundary.The wavelet only consists of two significant extrema.This reflection can only originate from the uprising capillary fringe since a reflection from a static layer boundary would move downwards during the imbibition, since an increase in water content above the layer boundary is resulting in a lower propagation velocity.Such a behaviour can be observed at (B2) for the reflection from the compaction horizon at 1.2 m depth.In this case one can also observe a phase flip of the wavelet which can be explained by the change of hydraulic properties due to the compaction (Sect.3.3).At (D) the top of the capillary fringe is reaching the upper boundary.The capillary fringe reflection signal gets significantly stronger near the upper boundary (C2) which can be caused either by a sharpening of the capillary fringe due to a capillary barrier at the boundary or by interference with the top layer reflection.

Drainage
While lowering the water table (Fig. 2b), the shape of the capillary fringe reflection (F) shows the same behaviour as observed in the imbibition period.Due to the variation of the hydraulic conductivity function over several orders of magnitude for different water contents in a specific sand, lower water contents translate to a slower relaxation time concerning externally induced changes.This means that parts with lower water contents might not be able to follow the water movement, which results in a potential sharpening during imbibition and widening during drainage.Hence, the shape of the capillary fringe reflection consisting of two extrema cannot be caused by a sharpening of the capillary fringe and the overlying transition zone during imbibition since drainage would show the opposite behaviour and widen the capillary fringe and the transition zone.
The upper boundary reflection is denoted by (E).The capillary fringe reflection (F) is well trackable throughout the descent of the water table until the capillary fringe moves through the compaction.After that only the reflection from the lower boundary is visible which is moving upwards during drainage (G1 to G2).Again the lower boundary reflection shows a phase flip of the wavelet.
In the following we concentrate on the different shape of the capillary fringe reflection in comparison to an ordinary layer boundary reflection.To retrieve the shape of the capillary fringe reflection, we choose the time interval between 32 and 55 min during drainage, where the capillary fringe can be observed without significant interference with the reflections from the boundaries (Fig. 3a).To illustrate temporal changes and allow for an identification of the characteristic shape, we look at all traces in this time interval, divided into smaller time intervals, marked by differently coloured areas in the radargram.The deformations of the reflection wavelet can be identified as interferences with reflections from structural heterogeneities.As visible in the radargrams, these are frequently present and constantly interfere significantly with the capillary fringe reflection due to its weak amplitude.This emphasises the importance of transient measurements to retrieve the shape of the capillary fringe reflection.Nevertheless, the shape of the capillary fringe can be identified as consisting of two extrema.In contrast, the observed reflection from the compaction (Fig. 3b) has three significant extrema.However, both reflections have a comparable main wavelength.
In the following we will show, with the help of numerical simulations, that the shape of the capillary fringe reflection is closely linked to the general shape of the water retention curve.In fact, this information allows to identify a hydraulic parameterisation most likely describing the water content dynamics for the observed sand since not every parameterisation is able to reproduce the observed reflection.
Before proceeding we complete the recovery of the basic hydraulic dynamics by investigating the observed phase flip of the reflection wavelet from the compaction horizon.

Impact of the compaction
The observed phase flip of the reflection wavelet can be explained by looking at the change of hydraulic properties due to the compaction.
A compaction of the sand leads to a more dense packing of the grains.This obviously results in a lower porosity leading to a lower saturated water content.Additionally, the generally smaller pore spaces lead to a stronger capillary rise in the compacted sand.The impact on the water retention curve is shown in Fig. 4a conceptionally by using the Brooks-Corey parameterisation (Brooks, 1966).This has a well-defined impact on the water content distribution around the compaction for different heights of the water table.Figure 4b illustrates this for a sand profile showing the discussed properties for two different characteristic water tables with a compaction at 1.2 m depth.The reflection then originates from the jump in water content at the compaction horizon leading to a jump in the permittivity with the same sign.Following Fresnel's equation for an incident signal perpendicular to the layer boundary and neglecting conductivity, this translates to a phase flip between different signs in the jump.For a lower water table, the stronger capillary rise in the compacted sand leads to a higher water content at the compaction horizon compared to the upper sand.For reasonably higher water tables only the saturated water content determines the jump in water content and thus permittivity at the compaction horizon.Since the permittivity is higher in the upper sand due to the higher water content, one gets a jump in permittivity with an opposite sign, resulting in the observed phase flip.

Hydraulic parameterisations
To parameterise the hydraulic properties of porous media, several parameterisations are available.In our case the water content distribution is of interest.The following parameterisation equations state the volumetric water content θ in dependence of the pressure head h, which corresponds to the height over the water table in equilibrium.θ s and θ r denote the saturated and residual water content, respectively.To investigate the possibility of reproducing the observed shape of the capillary fringe reflection, we study the three commonly used hydraulic parameterisations, namely the van Genuchten parameterisation, its simplified version and the Brooks-Corey parameterisation.
The full van Genuchten parameterisation (van Genuchten, 1980) is given by with scale parameter α and shape parameters n, m.
By fixing m = 1 − 1 n this is simplified to which is the most commonly used parameterisation in soil hydrology.
Furthermore the Brooks-Corey parameterisation (Brooks, 1966) is given by with scale parameter h 0 and shape parameter λ.

Numerical simulation
Since the observed behaviour of the capillary fringe reflection is similar for both imbibition and drainage, an investigation of the stationary water content is sufficient.All stationary water content profiles are calculated using θ s = 0.35 and θ r = 0.05.This decision represents typical values rather than an optimal estimate, which is also not the intention of this study.The given values just lie in a realistic range for a sand.Furthermore a 2 m 1-D sand profile is assumed with a water table at 0.6 m in all cases.
As an input for the GPR simulation, a water content profile is converted to a permittivity profile via the CRIM formula using ε air = 1, ε water = 80 and ε matrix = 5.
Since the focus of this study lies on the shape of the reflection signal and dispersive effects can be neglected in the frequency range used, electric conductivity is set to zero.Introducing it would only change the total amplitude and not affect the analysis.
We used the FDTD solver Meep (described in Oskooi et al., 2010).All simulations are carried out in a 2 m × 2 m × 2.7 m domain including a 0.5 m PML (perfectly matched layer) at each boundary.For z ≤ 2 m the permittivity is given by the calculated profiles, for z > 2 m, above the sand, the permittivity is set to ε = 1 for air.The permittivity is assumed to be constant in each horizontal direction.
The transmitter antenna is represented by a point source transmitting a Ricker wavelet with a centre frequency of f = 400 MHz, polarised in E y -direction.It should be noted that this is not exactly the same source wavelet as observed in the experiment but sufficient to reproduce the observed behaviour.The field in E y -direction is observed over time in 0.2 m distance to represent the receiver antenna.Both points are located 0.02 m above the ground to avoid non physical coupling due to the finite transition width between sand and air coming from the averaging procedure of the permittivity by the algorithm.The spatial discretisation is set to 5 mm to guarantee a good resolution of the capillary fringe and minimise numerical dispersion.
We now look at the characteristic reflections which can be observed in the modelled data when using the different parameterisations.A special focus is placed on the ability to reproduce the characteristics of the experimental data.
To allow a deeper understanding of the presented results, the general formation of a reflection from a continuous permittivity profile is described as conceptualised by the authors.

Formation of a reflection from a continuous permittivity profile
A given permittivity profile can be described as a composition of infinitesimal thin layers.The reflection signal then results from an infinite number of reflections and transitions originating from the layer boundaries.
Using the linearity of Maxwell's equations, it is possible to express the incoming wavelet as a superposition of circular monochromatic waves.Depending on the scale of the wave given by the wavelength, the wave is reflected differently from a given feature of the profile.
A wave with a small wavelength compared to the extent of a continuous feature will not be reflected since partial reflections interfere destructively.A wave with a much longer wavelength on the other hand will experience this feature as a highly localised one and will be reflected accordingly.With this understanding we investigate and explain the characteristic reflections from permittivity profiles calculated with the different hydraulic parameterisations.

Simplified van Genuchten parameterisation
In Fig. 5 we show the different permittivity profiles (right panel) and their resulting modelled reflection signals (left panel).For comparison also a sharp transition is shown (blue).
The common characteristics of a profile parameterised by the simplified van Genuchten parameterisation are illustrated with α = 4 m −1 and n = 6 (green).Also here the decision for the parameters should not be seen as an estimation attempt but a realistic parameter guess for this sand.The permittivity profile shows a continuous transition throughout the profile not distinguishably divided in capillary fringe and transition zone.The consequences for the corresponding characteristic  reflection compared to a reflection from a sharp transition are obvious: the main wavelength significantly increases and the signal has a much weaker amplitude.This is caused by the fact that only the low frequency components of the wavelet get reflected due to the width of the transition zone.By increasing n (red), the transition zone gets sharpened, corresponding to a more localised feature and therefore resulting in an increased amplitude and a smaller main wavelength since also higher frequency components get reflected.Nevertheless, the reflected wavelet always consists of three significant extrema, whereas the experiment only shows two significant extrema.
Only for the very special case of α ≥ 15 m −1 and n of around 2 (cyan), the observed signal shape consisting of only two extrema can be reproduced.However, these parameters do not represent a realistic sand since α indicates very large pores, whereas n is more representative for a rather finetextured material.
Summarizing the results, the simplified van Genuchten parameterisation is not suitable for reproducing the given experimental data and would lead to an erroneous result in a parameter estimation using this parameterisation and capillary fringe reflection data.

Brooks-Corey parameterisation
Figure 6 shows the results of employing the Brooks-Corey parameterisation with h 0 = 0.25 m and different λ compared to a sharp transition.The value of h 0 is chosen such that it allows an easy comparison with the full van Genuchten parameterisation (Sect.5.4).The permittivity profile always shows a kink at height h 0 above the water table, followed by a continuous transition above.This translates to a clear separation between capillary fringe and transition zone with the air entry point at h 0 .The corresponding modelled reflections show two significant extrema with a period comparable to the reflection from a sharp transition (blue).Since the kink at the air entry point is a spatially high localised feature, every frequency component gets reflected in the same way while only smaller frequency components get reflected from the transition zone.As a result of interference, the reflection wavelet only consists of a pronounced later part in terms of travel time compared to a reflection from a sharp transition.By lowering λ the general characteristic remains but the signal undergoes a stronger damping due to the larger transition zone.Comparing the results with Fig. 3, the modelled reflections match the experimental data showing two extrema.Going back to the radargrams (Fig. 2), these extrema correspond to the later part of a reflection from a boundary as seen in (C1), for example, when the capillary fringe reflection is separating from the layer boundary reflection.
This shows that reflections from sands parameterised by the Brooks-Corey approach reproduce the characteristic signal changes observed in the experimental data.To emphasise the fact that the different characteristic reflection compared to the simplified van Genuchten parameterisation is a direct consequence of the sharp air entry point, we investigate the full van Genuchten parameterisation.It enables us to control the sharpness around the air entry point by an additional parameter.

Full van Genuchten parameterisation
In addition to α and n, the full van Genuchten parameterisation includes a second shape parameter m.By setting n × m constant and increasing n, the sharpness at the air entry point can be increased without changing the behaviour for large heights.This is done for α = 4 and n × m = 5, starting with n = 6 (simplified van Genuchten, blue) (Fig. 7).
As the curvature at the air entry point gets sharpened in the permittivity profile (green and red), the shape of the reflection changes accordingly.The later part of the reflection wavelet gets stronger while the main wavelength decreases and becomes comparable to the source wavelet.The full van Genuchten parameterisation approaches the Brooks-Corey parameterisation for n → ∞, λ = n × m and h 0 = 1 α .Still, there remains a difference between the two at a value of n = 60.
This shows that for a large value of n compared to n × m, the full van Genuchten parameterisation also reproduces the observed reflection behaviour of the experimental data with two significant extrema.Furthermore, it can be stated that the key feature in reproducing the observed reflection is the sharpness around the air entry point.Additionally, the amplitude is very sensitive to even slight changes in this sharpness.Although we are not able to make a quantitative statement this last point gives a hint to the limitations of the method.The reflection amplitudes in Fig. 7 suggest a high probability that a sand following the van Genuchten parameterisation with a continuous transition around the air entry might not have a satisfactory signal-to-noise ratio due to the small amplitudes of the reflection.This fact also gives the limitations concerning a quantitative estimation of hydraulic parameters.In all cases the reflection originating from the transition above the air entry has a very small amplitude, presumably in the range of the deviations over time shown in Fig. 3b.

Conclusions
It was shown that an on-ground 400 MHz GPR system can provide valuable information about the basic shape of the capillary fringe without any further complex postprocessing.Simulations show that this can be used to select an appropriate hydraulic parameterisation for the observed sand which is an essential but often neglected aspect of estimating hydraulic properties.This is possible since the different commonly used parameterisations discussed here show significantly different reflections.
While the simplified van Genuchten parameterisation with m = 1 − 1 n and the Brooks-Corey parameterisation show completely different characteristic reflections, the full van Genuchten parameterisation can be understood as a continuous transition between them.This comes at the cost of an additional parameter, however.We can identify the shape at the air entry point as a key feature for the characteristic reflections.For the provided data it is shown that a sharp air entry with a transition zone above is required to reproduce the reflections.Therefore, the commonly used simplified van A. Dagenbach et al.: Identifying a soil water retention curve parameterisation with GPR Genuchten parameterisation is not suitable for reproducing the observed reflections since its profiles show no visible air entry point but a continuous transition throughout the profile.The full van Genuchten parameterisation and the Brooks-Corey parameterisation are able to reproduce the data well and therefore qualify as appropriate parameterisations.They include a sharp air entry point (Brooks-Corey) or they are able to model a sharp air entry (full van Genuchten).Nevertheless, the amplitude can vary significantly with a slight change of this sharpness by the full van Genuchten parameterisation.

Fig. 1 .
Fig. 1.Sketch of geometry of interest shown in vertical crosssection: two different sand layers plus another compaction horizon (dashed line).

Fig. 2 .
Fig. 2. Radargrams recorded during (a) imbibition and (b) drainage.The radargrams show the reflection from the layer boundary at 0.9 m depth (A and E), the reflection from the compaction horizon at 1.2 m depth (B and G) and the capillary fringe reflection (C and F).

Fig. 3 .
Fig. 3. Comparison of reflections from capillary fringe and compaction: (a) capillary fringe reflection taken from the traces marked in the radargram and the resulting mean value.The reflection is sensitive to interference with reflections from heterogeneities, nevertheless a shape with two significant extrema is observable.(b) Reflection from compaction taken out from the trace marked in the radargram.A shape with three significant extrema is observable.

Fig. 5 .
Fig. 5. Simplified van Genuchten parameterisation (Eq.3): permittivity profiles (left panel) and corresponding modelled reflections (right panel) with amplification factor if used.Starting from a realistic set of parameters (green), the impact of n is shown compared to a sharp transition (blue), while the cyan curve shows an unrealistic parameter set reproducing the measured signal.

Fig. 7 .
Fig. 7. Full van Genuchten parameterisation (Eq.2): permittivity profiles (left panel) and corresponding modelled reflections (right panel).Starting from parameters equal to the simplified van Genuchten parameterisation with α = 4 m −1 and n = 6 (blue), n is increased while n × m = const to sharpen the curve around the air entry.Furthermore, a profile parameterised by the Brooks-Corey parameterisation with h 0 = 1 α and λ = n × m (cyan) is shown for comparison.