Soil hydraulic material properties and subsurface architecture from time-lapse GPR

Quantitative knowledge of effective soil hydraulic material properties is essential to predict soil water movement. ground-penetrating radar (GPR) is a non-invasive and non-destructive geophysical measurement method to monitor the hydraulic processes precisely. Previous studies showed that the GPR signal from a fluctuating groundwater table is sensitive to the soil water characteristic and the hydraulic conductivity function. In this work, we show that this signal is suitable to accurately estimate the subsurface architecture and the associated effective soil hydraulic material properties with inversion 5 methods. Therefore, we parameterize the subsurface architecture, solve the Richards equation, convert the resulting water content to relative permittivity with the complex reflective index model (CRIM), and solve Maxwell’s equations numerically. In order to analyze the GPR signal, we implemented a new heuristic event detection and association algorithm. Using events instead of the full wave regularizes the inversion as it allows to focus on the relevant measurement signal. Starting from an ensemble of Latin hypercube drawn initial parameter sets, we sequentially couple the simulated annealing algorithm with the 10 Levenberg–Marquardt algorithm. We apply the method to synthetic as well as measured data from the ASSESS test site and show that the method yields accurate estimates for the soil hydraulic material properties as well as for the subsurface architecture by comparing the results to references derived from time domain reflectometry (TDR) and subsurface architecture ground truth data.


Introduction
Quantitative understanding of soil water movement is in particular based on accurate knowledge of the subsurface architecture and the hydraulic material properties.As direct measurements are time-consuming and near to impossible at larger scales, soil hydraulic material properties are typically determined with indirect identification methods, such as inversion (Hopmans et al., 2002;Vrugt et al., 2008).Time domain reflectometry (TDR, e.g., Robinson et al., 2003) is a standard method to acquire the required measurement data because it monitors the hydraulic processes accurately.Yet, being an invasive method, the TDR sensors disturb the soil texture of interest and typically require the maintenance of a local measurement station.Hence, it is difficult to apply the method at larger scales or to transfer the sensors to another field site.Ground-penetrating radar (GPR, e.g., Daniels, 2004;Neal, 2004) is an established non-invasive method for subsurface characterization and has the potential to Carlo (MCMC) methods (e.g., Scholer et al., 2011;Thoma et al., 2014;Jonard et al., 2015) and data assimilation approaches (e.g., Tran et al., 2014;Manoli et al., 2015;Rossi et al., 2015) have been successfully applied so far.
The GPR signal has to be processed automatically for parameter estimation.Many full waveform inversion approaches directly use the resulting Green's function (e.g., Lambot et al., 2009;Busch et al., 2012;Jadoon et al., 2012) in the cost function.Using the full antenna signal may lead to many local minima prohibiting a reliable identification of the global minimum (e.g., Bradford et al., 2014).In contrast, filtering the radargram with convolution approaches to determine travel time and amplitude of a limited number of events may even allow the application of efficient locally-convergent algorithms (e.g., Buchner et al., 2012).
In homogeneous material, the transition zone above the groundwater table exhibits a smooth variation of the relative permittivity.As the resulting GPR reflection is a superposition of a series of infinitesimal contributions along the transition zone, the detailed form of this reflection is sensitive to the variation of the relative permittivity.For simplicity, we refer to this reflection as transition zone reflection.Dagenbach et al. (2013) showed that this reflection is sensitive to the hydraulic material parameterization model.Bradford et al. (2014) measured the transition zone reflection of a drainage pumping test in a fluvial area with a antenna center frequency of 200 MHz and estimated hydraulic material properties.Klenk et al. (2015) employed numerical forward simulations and experiments using GPR antennas with higher antenna center frequency (400 and 600 MHz) for a more detailed explanation of the transition zone reflection during imbibition, relaxation, and drainage.They also concluded that the transition zone reflection is sensitive on hydraulic material properties.
In this work, we use the transition zone reflection together with reflections at material interfaces to determine the subsurface architecture and the corresponding hydraulic material properties.Therefore, the ASSESS test site was forced with a fluctuating groundwater table ensuring large hydraulic dynamics.The time-lapse measurement data was acquired with a single channel on-ground bistatic antenna operating at a center frequency of 400 MHz.Similar to Buchner et al. (2012), we solve Maxwell's equations in 2D and employ a new semi-automatic heuristic approach to extract travel time and amplitude of relevant reflections.This allows the optimization procedure to focus on the relevant information in the radargram and decreases the number of local minima.We draw an ensemble of initial parameter sets with the Latin hypercube algorithm.These parameter sets serve as initial parameters for the simulated annealing algorithm which is sequentially coupled with the Levenberg-Marquardt algorithm.We show that this procedure allows to accurately estimate the subsurface architecture and the associated effective hydraulic material properties for synthetic and measurement data.and C).The hydraulic state can be manipulated with a groundwater well (white square at 18.3 m) and is monitored with three GPR antennas (1, 2, and 3), a tensiometer (black square, at 4.0 m), and 32 TDR sensors (black dots).The gravel layer at the bottom ensures a rapid water pressure distribution over the site.An L element (left wall, at 0.4 m) and compaction interfaces (white lines) were introduced during the construction.Additionally to those visualized, GPR evidence indicates additional compaction interfaces (Fig. 13).Roman numbers (I)-(VII) indicate material interfaces referred to in the text.
Note the different scales on the horizontal and vertical axes.
temperature and TDR sensors.A geotextile separates the sand from an approximately 0.1 m thick gravel layer below, which ensures a rapid water pressure distribution and is the only connection of the groundwater well to the rest of the test site.Below this gravel layer, a basement layer partially consisting of reinforced concrete confines the site.As the test site is built into a former fodder silo, a concrete L element serves as additional wall.In order to stabilize the material during the construction, it was compacted.In addition to those shown in Fig. 1, GPR measurements indicate even more compaction interfaces (Fig. 13).

Representation
We follow Bauser et al. (2016) and define the representation of a system as a set consisting of: dynamics (mathematical description), subscale physics (material properties), forcing (superscale physics), and states.

Dynamics
The Richards equation (Richards, 1931), is the standard model to describe the propagation of the volumetric water content θ (−) and the matric head h m (m) in space and time t (s).The solution of this partial differential equation requires the specification of material properties, namely the soil water characteristic θ(h m ) and the hydraulic conductivity function K w (θ), which are (i) highly non-linear, (ii) varying over many orders of magnitude, (iii) showing hysteretic behavior, (iv) impossible to determine a priori, and (v) very expensive to measure directly.The unit vector in z-direction e z indicates the direction of gravity, typically pointing downwards.vertical lines), the position of the groundwater table was measured manually in the groundwater well and automatically with the tensiometer (Fig. 1).Notice that the difference between them is proportional to the driving force of water flow in the gravel layer.

Subscale physics
We choose the Brooks-Corey parameterization (Brooks and Corey, 1966) for the soil water characteristic θ(h m ), because it describes the materials in ASSESS appropriately (Dagenbach et al., 2013).Neglecting hysteresis, this parameterization may be inverted for θ r ≤ θ ≤ θ s , leading to including parameters representing a saturated water content θ s (−), a residual water content θ r (−), a scaling parameter h 0 (m) related to the air entry pressure (h 0 < 0 m) and a shape parameter λ (−) related to the pore size distribution (λ > 0).
Inserting the Brooks-Corey parameterization into the hydraulic conductivity model of Mualem (1976) yields the parameterization for the hydraulic conductivity function where K 0 (m s −1 ) is the saturated hydraulic conductivity and τ (−) a fudge factor.

Forcing
The ASSESS site was forced with a fluctuating groundwater table leading to three characteristic phases (Fig. 2): (i) initial drainage phase, (ii) multistep imbibition phase, and (iii) multistep drainage phase.We neglect evaporation in the following, because the experiment took place at the end of November and the weather was cloudy with 2-7 • C air temperature.The last precipitation was measured approximately 10 days before the experiment.More details about the experiment are given in Jaumann and Roth (2017).In this work, we only focus on the initial drainage and multistep imbibition phase.

State
The hydraulic state was monitored with GPR as well as with measurements of the position of the groundwater table in the groundwater well and at the position of the tensiometer.We used three shielded bistatic single channel 400 MHz GPR antenna pairs (Ingegneria dei Sistemi S.p.A., Italy).These antennas are referred to as antenna 1, 2, and 3, respectively.The measurement resolution was set to 2048 samples for 60 ns.In order to analyze the initial state of the test site, a multi-channel common offset measurement was acquired with antennas 1 and 2. The internal separation of the transmitter and receiver of these antennas is 0.14 m.During the experiment, the antennas were used to measure three time-lapse radargrams.In this work, we focus on the quantitative evaluation of the time-lapse data from GPR antenna 3 (Fig. 1).These data are analyzed in detail in Sect.3.3.
Additionally, a mean soil temperature (T s = 8.5 • C) and a mean direct current conductivity (σ dc = 0.003 S m −1 ) was estimated from TDR related measurements available in ASSESS.
The observation operator required to compare the hydraulic state with the GPR measurement data involves the solution of the time-dependent Maxwell's equations in linear macroscopic isotropic media.These equations quantify the propagation of the electromagnetic field consisting of the electric field E and the magnetic field B (Jackson, 1999): The dielectric permittivity ε = ε 0 ε r , magnetic permeability µ = µ 0 µ r , and direct current conductivity σ dc are generally spatially variable and represent the electromagnetic properties of the subsurface.Here, we neglect dispersive effects (∂ε r /∂ω = 0) as well as the imaginary part of the dielectric permittivity (ε r ∈ R).The relative magnetic permeability is assumed to be that of vacuum (µ r = 1).The source current density J is applied at the position of the transmitter antenna.
The relative permittivity of the subsurface ε r = ε r,b is calculated from the water content distribution θ resulting from the Richards equation using the complex refractive index model (CRIM) (Birchak et al., 1974): with the geometry parameter α = 0.5 (Roth et al., 1990).In order to apply the CRIM, the porosity φ, the relative permittivity of water ε r,w , the relative permittivity of air ε r,a , and the relative permittivity of the soil matrix ε s have to be known.The relative permittivity of air ε r,a was set to 1. Assuming that the sand matrix consists mainly of Quartz (SiO 2 ) grains, the relative permittivity of the soil matrix ε r,s was set to 5 (Carmichael, 1989).The porosity φ is assumed to be equal to the saturated water content θ s (Eq.2) which is estimated from the data.Following Kaatze (1989), we parameterize the dependency of the relative permittivity of water ε r,w on the soil temperature T s ( • C) with ε r,w (T s ) = 10 1.94404−Ts•1.991•10−3 .
In this work, the required direct current conductivity σ dc of the subsurface is assumed to be constant in the whole architecture.

GPR Analysis
Similar to Buchner et al. (2012), we extract the signal travel time t and the according amplitude A at M samples of the GPR signal (events) with a heuristic approach.This allows us to focus on the phenomena that are represented in the model and to exclude events of, e.g., reflections originating from compaction interfaces or confining walls.However, this procedure demands an automatic event association algorithm which associates events extracted from the measured signal with events extracted from the simulated signal.Thus, the evaluation method presented in this section consists of four steps: (i) signal processing, (ii) event detection, (iii) event selection, and (iv) event association.

Signal processing
The GPR signal is processed for further evaluation according to the following steps: (i) time-zero correction, (ii) dewow filter, (iii) 2D to 3D conversion, (iv) removal of the direct and trailing signal, and (v) normalization.
As the time-zero of the GPR antennas changes over time, we pick the direct signal and subtract it from the radargram for time-zero correction.Then, a dewow filter is applied to subtract inherent low frequency wow noise of the GPR signal.Since the observation is in 3D and the simulation in 2D, we convert the simulated signal to 2.5D, meaning to 3D with translational symmetry perpendicular to the survey line and parallel to the ground surface (Bleistein, 1986).Note that ASSESS is built accordingly (Sect.2.1).In this work, the conversion is done for the frequency and the amplitude separately.First, each trace is transformed to the frequency domain with the fast Fourier transform (FFT, denoted by ).Afterwards, the amplitude is modified depending on the angular frequency ω: where i is the complex unit, . ., K }, K is number of samples per trace enlarged to the next power of two).Subsequently, all traces are transformed back to the time domain with the inverse FFT.Due to the frequency conversion and the manipulation, a high frequency noise remains on the signal which is smoothed with a fourth order Savitzky-Golay filter (e.g., Press, 2007, we employed the implementation of the 'signal' package for GNU Octave: https://octave.sourceforge.io/signal/) using a window width of 41 samples.
Accounting for energy dissipation in 3D requires additional manipulation of the amplitude.Assuming a direct ray path and horizontal reflector with the reflector distance d and mean square root of dielectric permittivity √ ε along the ray path, this is done via .The amplitude of a trace is searched for extrema with a neighborhood search algorithm.For the subsequent evaluation, the amplitude of the detected events is normalized to the maximal absolute amplitude of all events detected in the trace.The direct signal and the trailing signal of the dewow filter are set to zero in a preprocessing step (Sect.2.3.1) and events close to these signals are ignored.
Subsequently, the direct signal and the trailing signal of the dewow filter are set to zero.Finally, the each trace is normalized to its maximal absolute amplitude.Notice that the signal is renormalized later in the analysis of the GPR data (Sect.2.3.4)

Event detection
To facilitate the identification of relevant events at large signal travel times, the normalized amplitude (original amplitude) is amplified quadratically with travel time.Subsequently, the extrema are detected with a local neighborhood search.Then, we keep a predefined number of events (15) with the largest amplified absolute amplitude.If the original amplitude of an detected extremum is below a predefined amplitude threshold (0.006) it is discarded in any case.In order to correct the perturbation in travel time due to the amplification and to cope with the discrete measurement resolution, we fit a Gaussian centered at the travel time of the detected event with width of ±5 samples to the original amplitude data.The resulting amplitude and travel time of the extremum are used for the following evaluation.

Event selection
After the event detection, the measured signal is inspected manually together with the detected events.In this one-time preprocessing step events can either be deleted or added manually.This ensures that only those events enter the parameter estimation that are also represented in the model.This step is skipped for the analysis of the simulated data.

Pairwise event association
The selected events extracted from the measured data have to be associated with the detected extracted from the simulated data for the parameter estimation.Therefore, Buchner et al. (2012) tested all possible combinations of events, using the one with the minimal summed absolute travel time difference.However, this is only feasible for a small number of events.As we .Exemplary association of simulated (s) and measured (m) events with indices 1-6 and 1-7, respectively.The color of the dots indicates the sign of amplitude of the events.(a) The detected events (Fig. 3) are aggregated in clusters to minimize the number of possible event combinations.The clusters are associated such that the summed absolute travel time difference of the mean travel time of the events in the cluster is minimal.(b) The events in the clusters are associated according to consistent temporal order and amplitude sign.Hence, if (ts,1, As,1) is associated with event (tm,2, Am,2), event (ts,2, As,2) can not be associated with event (tm,1, Am,1), if tm,1 < tm,2 or sign(As,2) = sign(Am,1).Solid (dashed) arrows indicate some of the accepted (declined) association combinations.The combination with maximal number of associations and minimal summed absolute travel time difference is used for evaluation.
are not using a Gaussian convolution of the data but the data themselves, the number of events increases.Hence, testing all combinations is often prohibitively expensive.In order to exclude combinations a priori, the detected events are aggregated in clusters (Fig. 4a).Then, these clusters are associated by testing all possible combinations.We use the combination with the minimal summed absolute travel time difference.Afterwards, the events aggregated in the associated clusters are associated themselves.The applied association procedure requires the events to have an identical amplitude sign and a consistent temporal order which reflects the principle of causality (Fig. 4b).We iterate over all according combinations to find the association with the maximal number of associated events and the minimal summed absolute travel time difference.It is critical to also consider combinations where some intermediate events (e.g., (t s,2 , A s,2 ) in Fig. 4) can not be associated.After the association of the events, outliers are detected by calculating the mean and standard deviation of the travel time differences.All associations are discarded which exhibit an absolute travel time difference larger than 3 standard deviations of all absolute travel time 10 differences.Finally, the amplitude of the associated events is normalized to the maximal absolute amplitude of the associated events in each trace.

Parameter estimation
Inversion of GPR data typically requires globally convergent parameter estimation algorithms which are computationally expensive.In order to keep the parameter estimation procedure efficient, we use an iterative strategy (Fig. 5).We start the optimization procedure by drawing an ensemble of initial parameter sets with the Latin hypercube algorithm (implemented by the pyDOE package, https://github.com/tisimst/pyDOE).The most expensive operation of the forward simulation is the calculation of the observation operator, which includes the solution of Maxwell's equations (Sect.2.2.4) and the subsequent event association (Sect.2.3).Hence, as time-lapse GPR data are highly correlated in experiment time (e.g., Fig. 13), we equidistantly subsample the number of traces of the time-lapse GPR radargram and generate a data set with lower temporal resolution.We use those data to improve the the distribution of the initial parameters (preconditioning).Therefore, the drawn parameter sets are used to initialize the simulated annealing algorithm (Sect.2.4.2) which allows for a robust, fast, and easy to implemented parameter update.Subsequently, the resulting parameters serve as initial parameters for the Levenberg-Marquardt algorithm

Objective function
Assuming P parameters p π and M associations of measured events (t µ,m , A µ,m ) with simulated events (t µ,s (p), A µ,s (p)), the χ 2 objective function is given by with the constant standard deviation of the measured normalized travel times σ t,µ = σ t and of the measured normalized amplitudes σ A,µ = σ A leading to the residuals in travel time r t,µ and amplitude r A,µ .Due to the oscillating nature of the GPR signal and due to the analysis (Sect.2.3), the χ 2 function is not convex and may even be discontinuous at some points, as the number of associated events M is not necessarily constant during the minimization process.To compensate for adding and removing simulated or measured events, Buchner et al. (2012) introduced tagging.If the number of measured events is smaller than the number of the simulated events, the simulated events that are not associated are excluded.Alternatively, if there are more measured events, measured events without partner are tagged as partnerless.If a reflection event has been tagged and becomes untagged after the parameter update, the contribution of the event and its new partner to the objective function is added to the previous objective function value.If an event has not been tagged and becomes tagged after the parameter update, the contribution to the cost function is subtracted from the previous objective function value.

Simulated annealing
We choose the simulated annealing algorithm (Press, 2007) to start the minimization of the objective function (Eq.11), because this algorithm is gradient-free and updates the parameters statistically.Additionally, tagging (Sect.2.4.1) can be implemented easily and it also allows uphill steps, which can be favorable if the events are not yet associated to their appropriate reflection.
If the parameter update is drawn from the whole parameter space, the algorithm is globally convergent.However, this approach is typically inefficient.We mainly use the simulated annealing algorithm to find a parameter set that associates the events to their appropriate reflection such that the more efficient gradient-based algorithm can take over.Hence, we search the neighborhood for better parameters starting from Latin hypercube sampled initial parameters p π,0 .For each iteration i (1, . . ., I), new parameters are proposed randomly via with a mobility parameter m = 0.1, uniformly distributed random number u p ∼ U(−1, 1), and the parameter limits p π,max and p π,min .In order provide the control parameter T , which is an analogue of temperature, we choose an exponential cooling schedule choosing parameter k = 1, and accept the proposed parameter set if P i+1 > u d or else draw a new parameter set.

Levenberg-Marquardt
The Levenberg-Marquardt algorithm is implemented as described by Jaumann and Roth (2017).However, in order to successfully apply this gradient-based algorithm to GPR data, the optimization has to be regularized.Therefore, we focus with this algorithm in particular on the improvement of small residuals, because if the small residuals improve, the larger residuals are likely to also improve in subsequent iterations due to the temporal correlation of the data.Therefore, we tag events with r t,µ > 100 or r A,µ > 100.Tagged events are excluded from the optimization by setting the according entries in the Jacobi matrix (J µ,π = ∂r µ /∂p π ) to zero.The event association may also change during the perturbation of the parameters for the numerical assembly of the Jacobi matrix.This can lead to large changes in the residuals, which in turn may lead to a disturbed parameter update.Hence, corresponding entries of large changes in the residual abs (r µ (p perturbed ) − r µ (p)) > 50 are also set to zero together with entries of the Jacobi matrix that are larger than 10 4 .
We choose λ initial = 5 as initial value for λ and apply the delayed gratification method by decreasing (increasing) λ by a factor of 2 (3) if the parameter update is successful (not successful).This assures that the algorithm takes small steps such that association and the Jacobi matrix can adapt smoothly.

Application
In this section, we apply the presented methods to the acquired GPR data.Therefore, we first explain the setup of the case study and the parameter estimation procedure (Sect.3.1.1).Then, we test the method with synthetic data to understand the phenomenology of the data and capabilities of the method (Sect.3.2).Finally, we apply the method to the measured data (Sect. 3.3) and analyze the accuracy of the resulting parameters.
3.1 Setup of the case study

Implementation
The numerical solution of the Richards equation (Eq. 1) is based on µϕ (muPhi, Ippisch et al., 2006) which uses a cell centered finite volume scheme with full upwinding in space and an implicit Euler scheme in time.The nonlinear equations are linearized by an inexact Newton method with line search and the linear equations are solved with an algebraic multigrid solver.
We solve Richards equation in 1D (z dimension) on a structured grid with a resolution of ≈ 0.005 m.Generally, the boundary condition is implemented with a Neumann no-flow condition.However, during the forcing phases, we prescribe the measured groundwater table as a Dirichlet boundary condition at the position of the groundwater well.We initialize the simulation with hydraulic equilibrium based on the measured groundwater table position.The simulated water content is converted to relative permittivity via the CRIM using the mean soil temperature T s = 8.5 • C (Sect.2.2.4).
To simulate the temporal propagation of the electromagnetic signal, we solve Maxwell's equations (Sect.2.2.4) in 2D with the MIT electromagnetic equation propagation software (MEEP, Oskooi et al., 2010).The transmitter antenna is represented with an infinite dipole pointing in x dimension.Thus, we neglect any effects from the real antenna geometry (bow tie), cross coupling or antenna shielding.The antenna source current density J is given by a Ricker excitation function (first derivative of a Gaussian-shaped function) with a center frequency of 400 MHz.The receiver antenna is not represented explicitly.Instead, E x is read directly at the position of the receiver antenna.We use the antenna separation of the real GPR system (0.14 m) in the simulation.Perfectly matched layers (PML) of 0.15 m thickness serve as boundary condition.The initial electromagnetic field in the domain is zero.We use one tenth of the minimal wavelength λ w,min as upper limit for the spatial resolution ∆z: with the speed of light in vacuum c 0 , maximal frequency f max = 2•400 MHz, and ε r,max = 31.25 corresponding to θ s,max = 0.5.
Hence, we choose the numerical resolution ∆z = 0.005 m for the 2D isotropic, structured, and rectangular grid.Therefore, the simulated one dimensional relative permittivity distribution is extruded in the y dimension.The Courant number for the FDTD method is set to 0.5.
To avoid multiple reflections at the air-soil boundary, we set the relative permittivity above the soil to 3.5, which is typical for dry sand.This is justified, as no evaluation of air wave or ground wave is done and the amplitude is normalized according to the detected events.The permittivity of the basement below ASSESS is set to 23.0, based on previous simulations.The direct current conductivity of the subsurface σ dc is set to 0.003 Sm −1 (Sect.2.2.4).All electromagnetic properties are smoothed by MEEP according to Farjadpour et al. (2006).The subsurface architecture is represented with layers.The position of these layers is parameterized and can be estimated.For illustration, the setup is shown in Fig. 6.

Setup of the parameter estimation
General setup of the optimization is explained with Fig. 7.This setup is used in an iterative approach (Fig. 5), where we selected ever fewer of the traces of the time-lapse GPR data to generate a data set with decreased temporal resolution.the parameter fit range given in Table 1 determines the parameter update via p π,max and p π,max (Eq.12).After the simulated annealing algorithm, we run maximally 15 iterations of the Levenberg-Marquardt algorithm (Sect.2.4.3) which completes the precondition step.The resulting parameter sets serve as initial parameters for the Levenberg-Marquardt algorithm which is applied to the data with high temporal resolution.
In order to evaluate the performance of the ensemble members, we use the mean absolute error in travel time e MA,t , because this statistical measure is independent of the number of associated events.These are accounted for by choosing those 10 members with minimal e MA,t where at least 85% of the measured events are associated.Each of these members has locally optimal parameters.However, the exact position of these local minima typically depends on (i) the settings of the optimization algorithms, (ii) the choice of the events to be evaluated, and (iii) the random numbers drawn in the simulated annealing algorithm.There is also no guarantee that the global optimum was found by one of the ensemble members.However, the distribution of these 10 best ensemble members contains valuable information about the shape of the χ 2 surface.To account for this information, we (i) analyze the mean parameter set of the best members and (ii) use the according standard deviation to indicate the uncertainty of these parameters.Notice that the mean parameter set is not necessarily optimal.However, if uncertainty on the parameters is small, this result is typically more reliable results than the best ensemble member.
The standard deviation of the measured data, σ t ≈ 6 • 10 −4 and σ A ≈ 5 • 10 −3 for normalized travel times and amplitudes is used as the standard deviation of the residuals in the objective function (Sect.2.4.1).

Phenomenology
The phenomenology of the transition zone reflection for characteristic times during imbibition, relaxation, and drainage was discussed by Klenk et al. (2015) for typical coarse sand.Here, we focus on the temporal development of this reflection during imbibition and equilibration.Therefore, we simulated water content distribution in the 1D profile located at 17.05 m of ASSESS using parameters typical for coarse-textured sandy soils (Table 2).The results are visualized over time (Fig. 8a) and over the water content (Fig. 8b).Initialized with hydraulic equilibrium, the simulation starts with the initial drainage step (Sect.2.2.3) where the groundwater table is lowered.Hence, the material with high initial water content is desaturated.After the subsequent equilibration step, the groundwater table is raised during the subsequent imbibition step.Generally, the Brooks-Corey parameterization (Eq.2) features a sharp kink where air enters the material at the upper end of the capillary fringe.Additionally, the imbibition causes another kink in the water content distribution (at marker (2) in Fig. 8b), because the relaxation time from the hydraulic non-equilibrium is much shorter at high water contents compared to the relaxation time at low water contents.This is due to the strong non-linear dependency of the hydraulic conductivity (Eq. 3) on the water content leading the differences in hydraulic conductivity of several orders of magnitude.Hence, the transition zone is sharpened during the imbibition.During the equilibration step after the first imbibition, the transition zone smoothes.Thus, the water content increases in the material with low water content (3) and decreases in the material with high water content (4).This smoothing process strongly depends on both the soil water characteristic and the hydraulic conductivity function.Sharpening and smoothing of the transition zone are repeated consistently for the other subsequent imbibition and equilibration phases ( 5), ( 6) and ( 7), ( 8).
According to the CRIM (Sect.2.2.4), the relative permittivity distribution has the same shape as the water content distribution.
Hence, kinks in the water content distribution directly induce partial reflections of the GPR signal (Fig. 8c).Shortly after starting the imbibition, the amplitude of the reflection at the additional kink (2) increases.After passing the material interface (V), the spatial distance of the kinks increases such that the two resulting reflection wavelets (3) and ( 4) are separable.The signal in between these wavelets is a superposition of infinitesimal reflections which contain detailed information about the form of the transition zone.Notice that the reflection (3) scans the initial water content distribution, which in steady state corresponds to the soil water characteristic.With progressing equilibration, the amplitude of reflection (3) decreases as the transition zone smoothes.The GPR signal of the subsequent imbibition and equilibration phases ( 5), ( 6) and ( 7), ( 8) show similar behavior and emphasize the relatively long time scale for hydraulic equilibration of sandy materials.In summary, this numerical simulation confirms qualitatively (i) that the dynamics of the fluctuating groundwater table is sensitive to both the soil water characteristic and the hydraulic conductivity function and (ii) that the transition zone reflection leads to tractable reflections during the imbibition step.

Results and discussion
The resulting soil water characteristics for material A (Fig. 9a) exhibit a similar curvature but are shifted.Both the parameters h 0 and λ influence the shape of the desaturated transition zone.Hence, merely evaluating the shape of the desaturated part of the transition zone is not necessarily sufficient to uniquely identify both parameters leading to large correlation coefficients.
However, parameter h 0 additionally determines the extent of the capillary fringe.If the evaluation is also sensitive on the extent of the capillary fringe, h 0 can be uniquely identified which significantly decreases the correlation between h 0 and λ.Hence, we conclude that the strong correlation of the parameters h 0,C and λ C (−0.7, Fig. 10) indicates that the evaluation is more sensitive to the shape of the desaturated part of the transition zone than to the extent of the capillary fringe.
As the architecture is a layered structure where material C is located above material A (Fig. 6), the water content in material C contributes to the travel time of the other reflections.This introduces correlations of all the parameters associated with the soil water characteristic of material C to θ s,A .A high correlation of parameters indicates that the problem is not well-posed.This typically increases the number of local minima and thus the uncertainty of the parameters.
The saturated hydraulic conductivity of material A (Fig. 9b) is approximately one order of magnitude smaller than the saturated hydraulic conductivity of material C. As the 1D architecture is forced at the lower boundary, the hydraulic conductivity of material A limits the water flux into material C. Hence, the data are not sensitive on K s,C .Yet, the uncertainty of the hydraulic conductivity decreases for low water contents as the reflection at the additional kink (Sect.3.2.1) is sensitive to the hydraulic conductivity.The hydraulic conductivity function (Eq. 3) is not unique if K s is not fixed.This leads to a strong correlation of the parameters K s,C and τ C (0.6, Fig. 10).Note that the uncertainty of the saturated hydraulic conductivity of material A also influences the uncertainty of the hydraulic conductivity of material C.
The uncertainty in the soil water characteristic of material C (Fig. 9c) is largest for low water contents, as there are only few data points available.In particular, this increases the uncertainty of λ A (±0.7, Table 2).The material properties of the unsaturated -5 with the mean of these parameter sets and the true parameter set (Table 2).The plot range of the parameters is adjusted to the water content range of the data.
material A are only monitored during the first ≈ 5 h of the experiment and are independent of the largest part of the other data.
This regularizes the optimization leading to fewer local minima.Similar to material C, the parameters h 0,A and λ A are strongly correlated (−0.6).Yet, the uncertainty in h 0,A is relatively small (±0.008,Table 2) mainly because it is essentially uncorrelated to other parameters.In contrast, the parameter θ s,A is correlated to the parameters h 0,C , λ C , θ s,C , and θ r,C as wrong parameters for material C introduce changes in the total water content which can be partially balanced out by adjusting θ s,A .
The uncertainty of the saturated hydraulic conductivity of material A (Fig. 9d) is comparably small as the largest fraction of the data are influenced by this parameter.Hence, the parameters τ A and K s,A are only very weakly correlated.
The correlation coefficients (Fig. 10) also show that both the permittivity and the thickness of the gravel layer can be estimated reliably with the presented evaluation method using travel time and amplitude information of a single channel.Evaluation methods that merely exploit the signal travel time (e.g., Gerhards et al., 2008), require additional measurements to achieve this goal.
In order to further investigate the quality of the mean parameter set, we simulated the resulting water content distribution (Fig. 11a) and calculated the difference to the true water content distribution (Fig. 11b).Due to the narrow pore size distribution of the sandy material, small deviations in the parameters h 0 and λ lead to large differences in the volumetric water content above the capillary fringe (≈ ±0.04).Combined with deviations in the position of the material interface, the largest differences in volumetric water content reach up to 0.17.Still, the mean absolute deviation of the volumetric water content is 0.004.
These deviations also cause residuals in the GPR signal (Fig. 12), which are most evident for the reflection at the gravel layer (VI).The bias of its travel time shows that the total water content above the gravel layer is underestimated with the mean parameter set.This bias is essentially balanced out with the properties of the gravel layer.However, the reflection originating from the basement of ASSESS (VII) reveals residuals that decrease as soon as the groundwater table is above the initial groundwater table.This indicates (i) deviations in the initial water content distribution and (ii) that the hydraulics during the initial drainage phase is not correct.
Similar to the analysis of the deviation in water content (Fig. 11), the largest residuals in unsaturated hydraulics are found where the groundwater table is crossing the interface of materials A and C.This indicates that the interference of the according reflections still contains information which could not be exploited in the parameter estimation procedure.Apart from that, the deviations in the material properties of unsaturated material C do not lead to significant residuals in the GPR signal, although the deviations in water content are considerable.

Phenomenology
The common offset GPR measurement (Fig. 13a) reveals the initial state of ASSESS.The reflections of the material interfaces are marked with uppercase roman numbers (I, II, III, IV, V, VI, and VII) which were introduced in Fig. 1.Compaction interfaces were generated during the construction process of ASSESS.The most pronounced of them are marked with lowercase roman numbers (i, ii, iii, iv, and v).In particular, the reflection of compaction interface (iv) close to the reflection of the groundwater table (1) are difficult to separate from reflections from material interfaces.Reflections from confining walls are most visible around 1 m (W) but influence the signal for more than 2 m.The reflection of the edge of the L-element (L) is particularly prominent.As ASSESS is confined by walls at all four sides and approximately 4 m wide, the walls parallel to measurement direction also influence the measured signal.
The time-lapse GPR measurement was recorded at 17.05 m (Fig. 13b).As the groundwater table is raised, the reflection originating from the groundwater table (2) separates from the reflection of the compaction interface (iv) and its amplitude increases.
After passing the material interface, the reflection (2) splits in two separate reflections (3) and (4) due to the strong dependency of the hydraulic conductivity on the water content (Sect.3.2.1).As the transition zone is smoothing during the equilibration phase, the amplitude of reflection (3) decreases and the distance of the reflections (3) and (4) increases.During the subsequent imbibition step, the reflections are directly separated.According to the analysis of the synthetic data (Sect.3.2.1), the effects of the smoothing water content distribution are most clearly visible during the relaxation phase at the reflections ( 5) and ( 6).
However, the associated measured signals interfere with the direct wave, the ground wave, and the reflection from the compaction interface (i) which exacerbates the identification of these effects.The reflections ( 7) and ( 8 Together with the water content distribution, the time-lapse GPR data also contain information about the subsurface architecture.However, separating signal contribution from the subsurface architecture and the hydraulic dynamics is not always possible.Here, this is most prominent for the reflection of the material interface (V).Initially, the amplitude of this reflection is large because the water content in material C is near the residual water content, whereas the water content in material A is significantly higher at the material interface.As soon as both materials are water saturated, the amplitude of the material interface reflection (V) is low since the effective porosities of the two materials are similar.Thus, the amplitude of the reflected signal originating from the material interfaces may change depending on the hydraulic state.
Additional information about the subsurface architecture can be inferred from the reflection at the material interface between material A and the gravel layer (VI) and from the reflection at the material interface of the gravel layer and the concrete basement (VII).These reflections are in particular suitable to analyze the total change of water content over time.
In summary, we note that the characteristic properties of the transition zone reflection during the imbibition and equilibration steps that were identified in the simulation (Fig. 8) can also be identified in the measured data (Fig. 13).

Results and discussion
As the GPR measurements cover only a small portion of the subsurface architecture, we restricted the hydraulic representation to 1D (Sect.3.1.1).Hence, we neglect 2D effects such as lateral flow.This has to be considered during the event selection of measured data (Sect.2.3.3).Therefore, we merely focus the evaluation on the imbibition phase of the experiment as the effect of lateral flow in fluctuating groundwater table experiments is largest during drainage and close to the capillary fringe (Jaumann and Roth, 2017).
The identification of relevant events is more difficult for the measured data than for synthetic data due to additional reflections of the GPR signal which are not represented.The most important representation errors are presumed to be the (i) reduced dimensionality of the representation of the inherently three dimensional GPR antennas and ASSESS test site using a 1D hydraulic and a 2.5D electrodynamic model, (ii) neglected small-scale heterogeneity in particular associated with compaction interfaces, (iii) neglected reflections from confining walls, (iv) neglected roughness of material interfaces, (v) influences of the antenna characteristics on the GPR signal in particular including the direct wave, the ground wave, temporal fluctuation of time-zero, and the source wavelet, (vi) assumption of a constant direct current conductivity, and (vii) assumption of a constant soil temperature for the calculation of the relative permittivity of water.
The main findings concerning the mean parameters for the synthetic data (Sect.3.2.2) can also be identified for the measured data, i.e. (i) the shift in the soil water characteristic of material C, (ii) the large uncertainty of the saturated hydraulic conductivity of material C, (iii) the high uncertainty of the soil water characteristic of material A for low water contents, and (iv) the high sensitivity on K s,A .
Compared to the uncertainties based on synthetic data (Table 2), the uncertainty of the resulting mean parameters (Table 3) mostly increased due to representation errors.Yet, the parameters estimated from TDR measurements (Jaumann and Roth, 2017) are within the standard deviation of the mean parameter set except for four parameters (   together with the mean of these parameter sets (Table 3) and a reference parameter set determined from TDR data acquired during the same experiment (Jaumann and Roth, 2017).The plot range of the parameters is adjusted to the water content range of the corresponding data.
parameters are analyzed in the following.
The parameter θ r,C estimated from the GPR data is approximately a factor 3 larger than the estimated value based on the TDR data.Essentially, there are three main reasons for this.First, by evaluating the travel time of reflection (V), we use the integrated water content for inversion.This also includes the compaction interface (i) which is not represented in the model.
At the beginning of the experiment, the amplitude of this reflection is comparable to the amplitude of the reflection originating from the interface of material A and C (V).Notice that the amplitude of the reflection (i) does not vanish, but decreases when the material is saturated at the end of the experiment.This indicates that this reflection originates from changes in both the small-scale texture of the material and the stored water content at the beginning of the experiment.Hence, as this compaction interface is not represented in the model, the resulting θ r,C is increased coping for this representation error.Second, a deviation in the position of the groundwater table with reference to the antenna position at the surface can be partially adapted by changing θ r,C .As the position of the surface is subject to change over the years, the measurements of the groundwater table are referenced to a fixed point at the end of the groundwater well, leaving the exact position of the surface relative to groundwater table uncertain.According to Buchner et al. (2012), the accuracy of the ASSESS architecture when compared to GPR Compared to the evaluation of the TDR data, the resulting value for parameter τ C is a factor 2 larger for the GPR evaluation.
This parameter adjusts the hydraulic conductivity for the unsaturated material and is mainly determined with the reflections (3) and ( 5) originating from the additional kink during imbibition (Fig. 15).These reflections exhibit a small amplitude for low water contents.Yet, both reflections interfere with the rather prominent reflection originating from the compaction interface (i).Additionally, the reflection (5) also interferes with parts of the direct and the ground wave.Hence, the travel time of these reflections does hardly change leading to an underestimation of the hydraulic conductivity for low water contents.
Similar to parameter θ r,C , the parameter θ r,A yielded from the GPR evaluation is approximately a factor 3 smaller than the result from the TDR evaluation.However, this parameter can only be approximated evaluating the GPR data as they lack events that are influenced by low water content.
The resulting value for parameter K s,A is factor 2 smaller for the GPR evaluation compared to the result from the TDR evaluation.This parameter limits the flux through the lower boundary as the domain is forced with a Dirichlet boundary condition.
Hence, the parameter can be used to cope with errors in the boundary condition.Forcing ASSESS with a groundwater well instantiates a 3D flux (Jaumann and Roth, 2017).Since this is not represented, the hydraulic potential and hence the water flux is overestimated.This is compensated in the GPR evaluation by decreasing K s,A .
The estimated interface position of the material A and C corresponds well the to ground truth measurements acquired during the construction of ASSESS (Table 3).In contrast, the estimated position of the gravel layer deviates from the according ground truth measurements.However, the estimates are still within the uncertainty of the ground truth measurements when compared to GPR measurements.
The evaluation of the GPR measurement data (Fig. 15) shows that deviations from the shape of the reflected wavelet contributes to the residuum significantly.These deviations have three main origins: (i) unknown shape of the source wavelet of the GPR antenna when coupled to the subsurface, (ii) the assumption of a constant direct current conductivity for the whole subsurface, and (iii) neglected roughness of the material interfaces which influences the shape of the reflected wavelet.
Evaluating the different reflections (V) and (VII), the residuals in travel time suggest an incorrect width of the wavelet.At the beginning of the experiment, the simulated wavelet is too broad for reflection (V) whereas it is too narrow for the reflection (VII).This indicates that (i) the assumption of a constant direct electric conductivity in the whole architecture is suboptimal and that (ii) the direct electric conductivity can be assessed with GPR measurements.
The residuals for reflection (VI) are the major contribution to the cost function as different representation errors are combined.
Of all the events in the wavelet, the events with the largest travel time exhibit the highest residuals.Due to the large grain size of the gravel, the real material interface is rougher than its representation.This leads to a non-symmetric broadening of the measured compared to the simulated wavelet (Dagenbach et al., 2013).The residual of the other events in the wavelet can be attributed to a combination of representation errors associated with the boundary condition, suboptimal parameters, incorrect The resulting parameters for the measured data are mostly consistent with results from the TDR measurement data.We discussed the deviations of the parameters and basically associated them with representation errors or the lack of available measurement data.Critical representation errors comprise the neglected (i) compaction interfaces, (ii) spatial variation of the direct current conductivity, and (iii) roughness of certain material interfaces.

Data availability
The underlying measurement data are available at http://ts.iup.uni-heidelberg.de/data/jaumann-roth-2017-gpr-hess.zip Author contributions.S. Jaumann designed and conducted the experiment, developed the main ideas, implemented the algorithms, and analyzed the measurement data.K. Roth contributed with guiding discussions.S. Jaumann prepared the manuscript with contributions of both authors.
Table 3.The mean and the standard deviation are calculated using the resulting parameters from the ten best ensemble members (Sect. 3.1.2)estimated from measured data.The corresponding material functions are given in Fig. 14.The reference parameters for the materials A and C are determined from TDR data acquired during the same experiment investigated in this work (Jaumann and Roth, 2017).Note that the according standard deviations are determined formally from a single Levenberg-Marquardt run and hence are only representative for one local minimum.Also, these standard deviations are given with the understanding that they are specific to the applied algorithm and will change for different algorithm parameters.Hence, these standard deviations are in particular not suitable to compare the precision of the TDR and GPR evaluation.Notice that for the TDR evaluation the porosity of the materials is assumed to be known from core samples.
The reference parameters for the subsurface architecture are calculated from ground truth measurements acquired during the construction of ASSESS.The corresponding standard deviations are given according to Buchner et al. (2012).

Figure 1 .
Figure 1.ASSESS emulates an effective 2D geometry with three distinct kinds of sand (A, B, and C).The hydraulic state can be manipulated

Figure 2 .
Figure 2.During the experiment with three distinct phases (initial drainage, multistep imbibition, and multistep drainage -separated by Figure3.The amplitude of a trace is searched for extrema with a neighborhood search algorithm.For the subsequent evaluation, the Figure 4. Exemplary association of simulated (s) and measured (m) events with indices 1-6 and 1-7, respectively.The color of the dots

Figure 5 .
Figure 5.We choose an iterative parameter estimation procedure which (i) allows to minimize the computational cost and (ii) facilitates the implementation of tagging (Sect.2.4.1).Therefore, we precondition the Latin hypercube sampled parameter sets with a data set with reduced number of traces (low resolution data) and use optimization algorithms which compare the parameter sets sequentially.The preconditioned parameter sets for each ensemble member serve as initial parameters for the final parameter estimation based on high resolution data.The subsequent evaluation of the ensemble is based on the number of associated events M and the mean absolute error in travel time eMA,t (Sect.3.1.2).

10(
Sect. 2.4.3)concluding the preconditioning step.The preconditioned parameter sets are used as the initial parameter sets for the more expensive optimization of high resolution data set with the Levenberg-Marquardt algorithm.The details of the setup and the analysis of the parameter estimation are given in Sect.3.1.2. 10 Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2017-538Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 13 September 2017 c Author(s) 2017.CC BY 4.0 License.

Figure 6 .
Figure6.For the simulation of the GPR signal of antenna 3, we assume a layered subsurface architecture (Fig.1).The transmitter of the antenna is represented with an infinitesimal dipole (t) and the electric field is read at the position of the receiver antenna (r).An absorbing layer is used as boundary condition.

Figure 7 .
Figure 7.The available hydraulic potential hwt is measured at the position of the groundwater well x λ times t β .These measurements are used as a boundary condition for the Richards equation (Sect.2.2.1).Estimates for the soil temperature Ts and the direct current conductivity σdc are derived from TDR related measurements.The actual signal of the GPR system is proportional to x component of the electric field Ex and measured discretely at experiment time t ξ and signal travel time tτ .This signal is used for event detection and event selection (Sect.2.3).The simulated water content distribution is converted to relative permittivity distribution with the CRIM and used to solve Maxwell's equations (Sect.2.2.4 and 3.1.1).After event detection, the simulated events are assigned to measured events.The mapping of the events is used to calculate objective function value during the optimization step (Sect.2.4).Dashed arrows indicate initial preprocessing steps, whereas solid arrows indicate iterative steps required for the optimization.

Figure 8 .
Figure 8.We used hydraulic parameters representing coarse-textured sandy soils (Table 2) to generate the synthetic data for parameter estimation according to Sect.3.1.1.Subfigure (a) shows the simulated water content in color code over experiment time, whereas subfigure (b) shows the same data in a line plot emphasizing temporal development of the water content distribution.The initial water content distribution is marked with a black dashed line.Subfigure (c) shows the according simulation of the GPR signal.The imbibition leads to a characteristic transition zone reflection (2).The temporal evolution of this reflection is sensitive on the initial water content distribution, the soil water characteristic and the hydraulic conductivity function.Except for the normalization, the data are processed according to Sect.2.3, including a dewow filter and 2D to 3D conversion.In contrast to the quantitative evaluation, the radargram is normalized to the maximal absolute amplitude, facilitating the visual comparison of the traces.The markers are used consistently in this paper and are further explained in the text.

Figure 9 .
Figure 9.The resulting material parameters estimated from synthetic data are shown for the 10 best ensemble members (Sect.3.1.2)together

Figure 10 .Figure 11 .Figure 12 .
Figure 10.The correlation coefficients for the mean parameter set are analyzed in detail in the text.Notice that the porosity of the gravel (θs,G) as well as the position of the material layers (h1 and h2) can be reliably estimated from single channel GPR when evaluating both the signal travel time and amplitude.

Figure 13 .
Figure 13.We show a common offset measurement of the hydraulic state of ASSESS at the beginning of the experiment (a).The vertical line indicates the position of the time-lapse measurement shown in subfigure (b).The common offset (time-lapse) data was measured with antenna 2 (3).Hence, in particular the measured GPR signal of the direct wave and the ground wave is slightly different.Both radargrams are measured with internal channels with an antenna separation of a = 0.14 m.Except for the normalization, the data are processed according to Sect.2.3.In contrast to the quantitative evaluation, the radargram is normalized to the maximal absolute amplitude, facilitating the visual comparison of the traces.The markers are used consistently in this paper and are further explained in the text.

Figure 14 .
Figure 14.The resulting material parameters estimated from measured data are shown for the 10 best ensemble members (Sect.3.1.2)

Figure 15 .
Figure15.The evaluation the measured GPR data is separated in three parts: Subfigure (a) shows the detected (Sect.2.3.2) and selected (Sect.2.3.3)events which are used as synthetic measurement data.Except for the normalization, the data are processed according to Sect.2.3.The radargram is normalized to the maximal absolute amplitude, facilitating visual comparison of the traces.Subfigure (b) shows resulting differences in travel time and amplitude of the mean parameter set.The differences of the amplitude are given in arbitrary units which are consistent for the whole work.Subfigure (c) shows standardized residuals of the differences given in subfigure (b).Note that outliers are set onto the boundary.The markers for the reflections are used consistently in this paper.Since the hydraulic 1D model cannot represent lateral flow present in initial drainage, we neglect the measured events of the first 2 h.Additionally, outlying events which have a different sign amplitude are not considered for event association (Sect.2.3.4).
Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2017-538Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 13 September 2017 c Author(s) 2017.CC BY 4.0 License.correlations are mostly specific to the experiment type and the subsurface architecture.Using travel time and amplitude information in the evaluation allowed to estimate effective permittivity and layer depth simultaneously with a single GPR channel.