Efficient estimation of effective hydraulic properties of stratal undulating surface layer using time-lapse multi-channel GPR

Multi-scale soil architectures in shallow subsurface are widespread in natural and anthropogenic depositional environments, and acquisition of the surface stratal structure and hydrological properties are essential in quantifying water cycling. Geophysical methods like ground-penetrating radar (GPR) can provide quantitative information like soil architecture and spatiotemporal soil water content distribution for the shallow layer. Concerning the informative multi-dimensional water flow in the surface layer with an undulating bottom at the plot scale, this study assesses the feasibility of efficiently estimating soil hydraulic properties using a few time-lapse multi-channel GPR observations, namely soil water storage and layer thickness of the surface layer, at reclamation land near an old river channel. We show that effective hydraulic properties of the surface layer can be obtained with a small number of time-lapse GPR measurements during a rainfall event. Additionally, we analyze the effect of some key factors controlling the informative lateral water redistribution on the results of the proposed approach using synthetic simulations.

1 Introduction

The exchanges of water and energy fluxes between the land surface and the atmosphere are strongly influenced by geomorphic features like landforms and soils (e.g., Corenblit et al., 2011). Surface morphometric features can be comparatively conveniently obtained from remote sensing at the sub-meter scale via satellites or drones. In contrast, the acquisition of subsurface soil features like textural and structural information is not so efficient, although they are essential to quantify soil functions in various ecosystems (e.g., Rihani et al., 2010). Particularly, soils in fluvial depositional environments exhibit a hierarchical stratal architecture at the scale of meters to hundreds of meters, and modeling subsurface water and mass transport is difficult due to large uncertainties in representation of the hydraulic properties (Huggenberger and Aigner, 1999). Apart from these natural depositional systems, land reclamation also forms a surface layer overlying undulating landforms. The surface layer is usually characterized by a flat surface and an undulating bottom, and the internal lateral water flow can be non-negligible in wetting conditions. Model predictions of water flow and nutrient transport in the surface layer are essential to agriculture. However, it is challenging for the acquisition of such soil architecture information and soil hydraulic properties with the commonly used approaches using point measurements, as they are costly and time-consuming.

For in situ estimation of hydraulic properties, inverse methods are commonly used at multiple scales, depending on observed state variables (Vrugt et al., 2008). For the one-dimensional (1-D) in situ monitoring profile, soil hydraulic parameters for each layer can be reasonably estimated based on a 1-D soil hydrological model and a time series of point observations. For larger-scale observations, e.g., from the field scale to catchment scale, a stronger assumption of predominant vertical water flow is implicitly set from the 3-D soil hydrological model. Alternatively, an assumption of the homogeneous soil column for the shallow measurement depth (< 5 cm) is used for the satellite remote-sensing observations (e.g., Mohanty, 2013). Overall, these methods ignore lateral water redistribution and usually suffer from insufficient data information concerning the soil architecture and the soil water dynamics. Yet this information is essential to all inverse parameter estimation methods (Bandara et al., 2013).

Quantitative observations of soil architecture and soil water content have become viable during the past decade using geophysical methods. In particular, ground-penetrating radar (GPR) has been used for efficiently imaging the shallow soil water content distribution (e.g., Huisman et al., 2003; Weihermüller, et al., 2007) together with soil architecture (Gerhards et al., 2008; Bradford, 2008; Pan et al., 2012a; Klenk et al., 2016). In view of informative multi-dimensional water flow contained in the time-lapse spatial observations of GPR measurements, the common problem of insufficient data information from 1-D observations is mitigated for inverse modeling of soil water dynamics in the vadose zone soils. This comes at the cost of increasing the complexity of inverse modeling and the need for spatial observations of the soil architecture as well as the soil water content. In addition, since the information of small-scale heterogeneities within the layers cannot be quantified for the GPR measurements using the reflection approach, the inverse modeling approaches only focus on the stratal soil properties and effective soil water flow at the plot scale.

There are at least two popular approaches for using GPR measurements in the inverse estimation of soil hydraulic properties. One approach is the coupled inverse modeling of hydrological processes and GPR measurements. The model for the simulation of the propagation of the electromagnetic wave in the soils is computationally expensive but necessary for full waveform inversion (e.g., Lambot et al., 2009; Busch et al., 2012; Jadoon et al., 2012) or other evaluation approaches (e.g., Buchner et al., 2012; Jaumann and Roth, 2018). The other approach directly uses the evaluated soil water content and depth from GPR measurements in the inverse hydrological modeling, similar to the above-mentioned 1-D inversion using point observations, e.g., time-domain reflectometer (TDR) measurements. In this study, we propose using pre-evaluated GPR data together with 2-D inverse hydrological modeling to estimate soil hydraulic properties of the surface soil layer with underlying undulating structures. We use the optimization procedure described by Jaumann and Roth (2017) to estimate the effective hydraulic properties and demonstrate that this approach allows us to efficiently estimate effective hydraulic properties of the surface layer at reclamation land. Controlling factors of this approach for practical applications are discussed subsequently.

Table 1 Summary of the equations used in the muPhi model (Ippisch et al., 2006) for simulating two-dimensional Darcian water flow in a variably saturated isotropic medium.

2 Scheme of the proposed hydraulic parameter estimation

2.1 Simulation model and setup

The muPhi model (Ippisch et al., 2006) is used for simulating two-dimensional Darcian water flow in a variably saturated isotropic medium. It employs the Richards equation (Richards, 1931) and the soil hydraulic functions of van Genuchten (1980) and Mualem (1976) for the soil water characteristic θ (h), usually in terms of the water saturation Θ (–) and the hydraulic conductivity function K_w (h), as shown in Table 1.

In this work, the water dynamics in a two-layer vadose zone with an undulating interface during a rainfall event is simulated. The model domain is discretized with rectangular grids, where the material for the top layer is characterized by five parameters, $p_{\mathrm{1}}=\left\{\mathit{\alpha },n,{\mathit{\theta }}_{\mathrm{r}},{\mathit{\theta }}_{\mathrm{s}},K_{\mathrm{s}}\right\}$, as shown in Table 1. Since soil porosity of the top layer is relatively easy to obtain by soil coring, θ_s is assumed to be known beforehand in the evaluation of the GPR observations of the mean soil water content and reflector depth of the interface (e.g., Pan et al., 2012a). With this assumption, there are four unknown parameters, $p_{\mathrm{1}}=\left\{\mathit{\alpha },n,{\mathit{\theta }}_{\mathrm{r}},K_{\mathrm{s}}\right\}$, for the top layer. Material for the bottom layer is sand and the hydraulic parameters, $p_{\mathrm{2}}=\left\{{\mathit{\alpha }}^{\prime },n^{\prime },{{\mathit{\theta }}^{\prime }}_{\mathrm{r}},{{\mathit{\theta }}^{\prime }}_{\mathrm{s}},{K^{\prime }}_{\mathrm{s}}\right\}$, are set to the values listed in Table 2 for all studied cases. A Neumann no-flow boundary condition is implemented at both sides, and two types of boundary conditions are applied at the upper boundary in consideration of numerical convenience. For the period with strong precipitation events, a Neumann condition is applied for the upper boundary, while a Dirichlet boundary condition is applied for the left period with a time series of outflux. For the bottom boundary, a Dirichlet boundary condition is applied with a fixed water pressure of −0.4 m. The equilibrium state is assumed as an initial condition at the lower boundary.

2.2 Parameter estimation

Given a time series of multi-channel GPR measurements, the surface soil structure and corresponding soil water dynamics were obtained using the evaluation algorithm proposed by Gerhards et al. (2008). Then, the effective hydraulic parameters were estimated accordingly. The framework of the 2-D inversion procedure is shown in Fig. 1.

Figure 1 The framework of efficient estimation of effective hydraulic properties of the surface layer with an undulating interface in a two-layer transect. Provided that the water dynamic range is captured by a small number of time-lapse multi-channel GPR observations, final parameter sets are inversely estimated based on the Levenberg–Marquardt algorithm and proper convergence criteria (see Sect. 2). Latin hypercube-sampled initial parameter sets for each ensemble member are used in the 2-D inversion. The final median parameter set is derived from the 68 % estimated ensemble parameter sets with minimal χ².

The initial hydraulic parameters for each ensemble member were generated based on the Latin hypercube algorithm. The Richards equation solver (muPhi; Ippisch et al., 2006) was used to simulate the spatiotemporal soil water dynamics, and an optimization procedure is implemented. Generally, the optimal soil hydraulic parameters, $p=\mathit{\{}\mathit{\alpha },n,{\mathit{\theta }}_{\text{r}},K_{\text{s}},({\mathit{\alpha }}^{\prime },n^{\prime },{{\mathit{\theta }}^{\prime }}_{\text{r}},{{\mathit{\theta }}^{\prime }}_{\text{s}},{K^{\prime }}_{\text{s}}\left)\mathit{\right\}}$, were determined with an objective function by minimizing the differences between observed and simulated water storages, $l_{\text{obs}}(x,t),l_{\text{mod}}(x,p,t)$, at location x:

$$\begin{array}{}\text{(1)}& {\mathit{\chi }}^{\mathrm{2}}\left(p\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _t^N\sum _x^M{\displaystyle \frac{\left[l_{\mathrm{obs}}\right(x,t)-l_{\mathrm{mod}}(p,x,t){]}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{obs}}}},\end{array}$$

where M is the number of grid cells in the x direction and N is the number of time series of observations. The Levenberg–Marquardt algorithm as implemented in Jaumann and Roth (2017) was used to minimize χ²(p). Convergence requires 5 to 40 iterations, depending on the existence of structural representation errors in the GPR observations. Depending on the data, two different criteria are used in this study. For the synthetic data without considering structural errors, the “optimal stopping” criterion is applied with a number of iterations of up to 40. In contrast, for the measured data with structural errors, the “early stopping” criterion is applied with only five iterations to avoid overfitting to the structural errors. All the inversions were conducted in a cluster via parallel computation. The final estimated parameter set was extracted from the ensemble members at the end.

3 Applications

In this section, we apply the approach introduced in Sect. 2 to a field dataset (Sect. 3.1). Afterwards, two controlling factors which can be used to improve the performance of the approach are investigated with synthetic studies (Sect. 3.2).

3.1 Field study: Daheigang dataset

The field test site is located near the village of Daheigang at 35^∘2.1^′ N, 114^∘33.8^′ E, in Fengqiu County, China. As an aeolian–fluvial depositional environment in the Yellow River floodplain, complex soil architecture is widespread in this area. After land reclamation in the middle of 1980s, aeolian–fluvial landforms are rarely seen except in nearby stripped dunes in the woods. The surface soils are dominated by the Ochric Aquic Cambisol and Ustic Sandic Entisol (Research Group of Chinese Soil Taxonomy System, 1995). This provides a suitable test site for the proposed approach.

Figure 2 The underlying interface of soil architecture derived from GPR measurements (adapted from Pan et al., 2012) and the selected transect (magenta curve along the dashed cutting line) perpendicular to the longitudinal troughs. The bottom lines are a projection of the 25 GPR survey lines.

The GPR data used in this study were introduced in Pan et al. (2012a) but are shortly described in the following for the convenience of the reader. The GPR survey was conducted using an IDS (Ingegneria dei Sistemi S.p.A., Italy) multi-channel GPR system, where two antennas operating at a central frequency of 400 MHz were connected in tandem. The setup was employed using three different antenna separations, S1 = S2 = 0.14 m, S3 = 1.94 m and S4 = 1.66 m. Given the investigated reflector depth around 1.0 m, this setup has an accuracy of 0.05 m in depth (Pan et al., 2012b). Five two-dimensional GPR surveys with two antennas were repeated along prefixed parallel lines with a 1.5 m interval spacing on 22, 23, 25, 27 and 29 May 2011 after a heavy rainfall event which followed 10 d without any rain. These GPR measurements were recorded with a resolution of 0.05 m along the acquisition line. The measuring wavelet data are the result of stacking 12 scans. Each measurement point was recorded for a time window of 80 ns, discretized in 1024 samples.

The GPR measurements were evaluated using the multi-channel GPR method (Gerhards et al., 2008). This approach has been used in several studies (e.g., Wollschlager et al., 2010; Westermann et al., 2010; Pan et al., 2012a) and uses the travel times from the four channels as a common midpoint (CMP) recorded during the measurement of each trace along a survey line. To obtain absolute travel times for all channels, time-zero calibration was applied (Pan et al., 2012a). Wide-angle reflection–refraction (WARR) measurements were conducted in air to get the time zero of the cross-antenna channels; offsets of the two box-internal channels were assumed to be equal to the travel time of the air-wave wavelet. Thus, reflector depth d (m) and the soil dielectric permittivity number ε_r (–) can be obtained by inversion. The soil water content can be derived from a petrophysical relationship between the depth-averaged volumetric water content θ (–) and soil dielectric permittivity number (Roth et al., 1990). To mitigate the negative correlation between d and θ, the total water storage $l=d\cdot \mathit{\theta }$ (m) is used in this study (Pan et al., 2012b).

Figure 3 Setup for the inverse modeling. (a) Observed water flux through the upper boundary. The colored dots (excluding the black one) correspond to the timing of soil water storage observations in (c). (b) The soil architecture of the transect. (c) Mean soil water content of the upper layer from time-lapse GPR observations on 22, 23, 25, 27 and 29 May.

3.1.1 Setup of the parameter estimation

Data analysis from Pan et al. (2012a) yields an interpolated 3-D undulating architecture and water content observations. A qualitative interpretation of the water content difference between the soil above the trough and hill indicates that lateral flow could happen at the weather condition shown in Fig. 2a. To minimize computational effort, a 2-D model of the water dynamics in a transect (magenta section along the dashed line in Fig. 2) perpendicular to the longitudinal troughs is used so that the 2-D water flow could approximate the actual 3-D water dynamics. Observations of reflector depth and water storage over the transect are calculated using the original survey lines (parallel lines at the bottom in Fig. 2) via 3-D linear interpolation. Considering the resolution of 0.05 m along the GPR survey lines and a 1.5 m in between the survey lines, the spacing interval of the transect is set to 0.1 m. Given a 1 m top layer with a mean soil water content of 15 %, the uncertainties of the depth measurement and water storage are set to 0.05 and 0.0075 m, respectively.

For the hydraulic model of the transect, the geometry (16.82 m × 1.3 m) is discretized with a resolution of 0.10 m × 0.05 m (Fig. 3b). The upper Neumann boundary flux (Fig. 3a) is calculated by subtracting evaporation from the natural infiltration flux. The rainfall observations are obtained from rain gauge measurements. The evaporation observations are estimated according to the empirical relationship between reference evaporation and pan evaporation proposed by Allen et al. (1998). Two different empirical coefficients, 0.7 and 0.63, are used for the two consecutive periods in view of the removal of land cover of wheat on the morning of 22 May 2011. The lower Dirichlet boundary is constantly set to a pressure of −0.4 m with respect to a water table of 1.7 m below the ground surface, inferred from drilling. The initial equilibrated condition is used due to a long preceding dry period. The soil water dynamics in the domain is simulated over a period of 18 d. Five time-lapse 2-D water storage observations (Fig. 3c) are used for the inverse modeling.

Compared to 1-D inverse estimation of hydraulic properties, one major advantage of the proposed 2-D approach is using a small number of time-lapse 2-D GPR observations, which increases the available information about the soil water dynamics. However, apart from the measurement precision, the proposed approach is sensitive to some structural errors presented in the field study. First, the conceptual error can also originate from the simplification of a 3-D flow to a 2-D flow. Second, the GPR-derived water storage observations might contain some structural errors due to the varying accuracy of multi-channel GPR evaluation with different ratios of antenna separation to reflector depth, small-scale soil heterogeneity and the interpolation of unevenly spaced data. Our approach, not representing these structural errors, may lead to overfitting when the stopping criterion is too small, i.e., Δχ²≤1, or χ²(p)<χ²(true). Hence, a two-step procedure is applied in this study. Initially, we used the optimal stopping criterion, only considering measurement precision, and it resulted in a large number of iterations as well as overfitting. After an analysis of the convergence behavior of the parameter distribution, we identified that overfitting could be avoided by using the early stopping criterion and evaluating the output of the fifth iteration. Note that the optimal stopping criterion should be applied once the structural errors are alleviated in the field study.

To increase our understanding of the performance of the proposed approach in the field study, two additional synthetic studies were set up. First, since the model is forced at the surface, the inversion is more sensitive to the parameters of the first layer than to the parameters of the underlying layer. Thus, an initial synthetic study was conducted to investigate the effect of the parameters of the bottom layer on that of the first layer. For this inversion we used the same settings as the field study, but adding unknown hydraulic parameters for the bottom layer. The required synthetic observations of soil water storage were generated using the forward simulations at the same time as field measurements and adding Gaussian distributed random errors of GPR observations. Hence, the above-mentioned structural errors were not included in the inversion and the optimal stopping criterion was used. Then, another synthetic study was conducted to investigate the effect of the above-mentioned structural errors on the parameter estimation of the first layer. All the settings are the same as the previous synthetic study, but using fixed parameters for the bottom layer.

Table 2 Hydraulic parameters of the two-layer soil and their allowed ranges for parameter estimation. The true values for the upper layer are obtained from Rosetta (Schaap et al., 2001) based on measured soil texture information, and the values for the lower layer are cited from Carsel and Parrish (1988).

Figure 4 Effective hydraulic properties estimated with the proposed approach for the field study (a, b) and the synthetic study (c, d). The grey curves in all subplots represent estimated hydraulic functions from ensemble members using the early stopping criterion. The black curves in (a) and (b) are derived from the median parameters of the best 68 members, while the black curves in (c) and (d) are the same members but using the optimal stopping criterion. The red curves in (c) and (d) are derived from the true parameter given in Table 2. The histograms in (a) and (c) (blue histogram; the same in Figs. 6 and 7) mark the covered range of the mean water content evaluated from the GPR observations.

3.1.2 Results and discussion

In the previously described initial synthetic study (Sect. 3.1.1), the parameters for both layers are estimated based on the synthetic data of the first layer. The resulting parameters for the first layer match the true parameter well, but the resulting parameter set for the bottom layer is biased dramatically. Here, we used a measure $S_j=\frac{\partial {\mathit{\chi }}^{\mathrm{2}}}{\partial p_j}$ to analyze the sensitivity of the inversion to the specific parameters and found that the sensitivity of the parameters of the first layer typically surpasses that of the bottom layer by 2 orders of magnitude. Therefore, we decided to reduce the computation cost by fixing the parameter of the bottom layer, although this introduces a small structural error. Hereafter, all the inversions only estimate the parameters for the first layer and use the fixed parameters listed in Table 2.

Using the setup described in Sect. 3.1.1, the estimated hydraulic properties for the upper layer from the field study and the synthetic study are shown in Fig. 4. Figure 4a and b present the resulting water retention curves and hydraulic conductivity curves, respectively, using the estimated parameters from the field study based on the early stopping criterion. We show those 68 best ensemble members with a minimum χ², accounting for 68 % of the 100 ensemble members. Given the observed water content range (18 %–27 %; histogram in Fig. 4a), the water retention curve is mainly constraint over a small water content range, and large uncertainties appear at both ends of the soil water retention curve, in particular at the dry end. Since we set the initial parameters with the Latin hypercube approach, the parameter space is equidistantly sampled. Due to the low sensitivity of the hydraulic conductivity and the limited iterations of the early stopping criterion, the uncertainty of the resulting hydraulic conductivity is high. The overall uncertainty of the resulting hydraulic conductivity function is relatively high due to the small dynamic water range and measurement errors.

Figure 5 Impacts of errors on the inverse estimation in the field study and in the synthetic study. (a) Structural pattern of standardized residuals in the field study indicates notable structural errors, apart from Gaussian distributed error at each point. (b) Gaussian distributed standardized residuals in the synthetic study.

Apart from the influence of the water dynamic range, the ignored structure errors also play an important role in the inversion for the field study. As shown in Fig. 5, the different structures of the standardized residuals in the field study and in the synthetic study indicate that notable structural errors exist in the field study. Once excluding structural errors in the synthetic study, the resulting hydraulic properties in Fig. 4c and d are much better than the field study. The bands of the resulting hydraulic curves based on the optimal stopping criterion are much more narrow than those based on the early stopping criterion. The structural errors lead to different histograms in Fig. 4a and c and result in uncertainties of the resulting hydraulic parameters that are a little bit larger in comparison to the synthetic study when using the same criteria. This is also the reason why overfitting occurs when using the optimal stopping criterion in the field study.

In summary, the barely passable performance of the proposed approach in the field study is not only attributed to the narrow observed water dynamic range but also the structural errors. However, better performance can be achieved by improving the experiment design. One solution is to deploy the GPR transect measurements directly over the undulating structure with steeper slopes, which results in a wider water content range by internal water redistribution. Another solution is to conduct the GPR measurements to capture a wider water content range, i.e., record data before a rainfall event. The two controlling factors are further verified in the following synthetic study.

Figure 6 Undulating soil architectures used in synthetic studies for investigating the performance of the proposed approach with different slope gradients. Three interfaces with different amplitudes are (a) S1, 0.25 m; (b) S2, 0.5 m; and (c) S3, 0.75 m.

Figure 7 Effects of the slope gradient on the performance of the proposed approach. The top, middle and bottom panels compare the results for the soil architectures in Fig. 6. Similar to the lower panel of Fig. 4 for the synthetic study, the true parameter set is shown together with the 20 best ensemble members and the median of these parameter sets. The uncertainty of the resulting hydraulic properties decreases along with the sequentially increase in the slope gradient and water content range of the histograms from S1 to S3.

3.2 Synthetic study: investigation of controlling factors

3.2.1 Setup of the parameter estimation

In this section we used synthetic studies to verify the proposed approach and assess the effects of the slope gradient on the accuracy of the approach. Three two-layer architectures (S1, S2 and S3 in Fig. 6) are employed with a domain of 6.28 m × 2 m and are discretized with a resolution of 0.04 m × 0.02 m. The internal layer interface is given as a sine-wave shape with different slope gradients (amplitude: 0.25, 0.5 and 0.75 m), and soil hydraulic properties of both layers are the same as the synthetic case discussed in Sect. 3.1. Given the same boundary conditions and initial station, parameter estimations are conducted using the optimal stopping criterion described in Sect. 2.2. In addition, the impact of the water dynamic range on the parameter estimation is also investigated by adding one more measurement of GPR observations before the heavy rainfall event (Fig. 3a).

Figure 8 The same as Fig. 7, but using wider soil water content range by adding one more GPR measurement before the rainfall (black dot in Fig. 3a). Excellent performances of the approach are achieved in S2 and S3, while only a little improvement happens in S1 compared to Fig. 7. Apart from the soil water content range, the slope-induced lateral water redistribution is also essential to the proposed approach.

3.2.2 Results and discussion

The hydraulic functions estimated from the three undulating architectures are shown in Fig. 7. The left panels compare the estimated water retention curves of S1, S2 and S3, and the right panels compare the estimated hydraulic conductivity curves. In each plot, the curves represent the best 20 estimates with minimal χ², accounting for 67 % of the 30 ensemble members. Results from the panels show that the uncertainty band of the estimates decreases sequentially from S1 to S3. In other words, the higher the slope gradient, the better the estimation. The histograms in Fig. 7a, c and e show an increase in water content range from top to bottom, where the ranges are 19 % to 25 %, 19 % to 26 % and 18 % to 30 %, respectively. This is mainly ascribed to the increasing intensity of lateral water redistribution, which results in a wider water dynamic range.

Figure 8 shows the results from the estimates, using the same setup as Fig. 7 but extending the observed water dynamic range by adding one more GPR observation. Compared to Fig. 7, performances of the approach in the three architectures are all improved. This is attributed to the extended water content range, which controls the performance of the inverse estimation. The histograms in Fig. 8a, c and e present two peaks and a much wider water dynamic range in comparison to Fig. 7, and the ranges all increased from 13 % to 26 %, 13 % to 26 % and 11 % to 31 %, respectively. However, the improvements of the three architectures are distinct. The most significant improvement happened in the architecture S2, with the intermediate undulating amplitude. The improvement in S1 is not as much as in S2, although the soil water content range in S1 is similar to S2. This is attributed to the non-negligible slope-induced lateral water redistribution in S2, since the inversion is not only influenced by the water content range but also the structural information change resulting from the slope-induced lateral water redistribution. The bell-shaped peaks in Figs. 7a and 8a indicate little structural information in S1. As a result, only the extended water content range has a little effect in S1. In contrast, both factors affect the inverse estimation for S2 and S3, although the added effect is not so significant for S3.

Results of both synthetic studies in Figs. 7 and 8 show that excellent performance of the approach can be achieved when the information of the water content range and slope-induced lateral water redistribution is enough to constrain the inversion. Particularly, the structural information resulting from the slope-induced lateral water redistribution is essential to the proposed approach. To improve the performance of the approach, other factors contributing to lateral flow, such as the soil hydraulic properties and the intensity and duration of the precipitation, should be considered in designing field experiments.

4 Discussion

For practical application, we describe some necessary conditions for this approach to work. First of all, a surface layer with an undulating bottom is required, which leads to notable lateral water redistribution. Following this concept, applying this approach to a three-dimensional architecture could be possible but needs further computational and experimental efforts. Secondly, for rain-based cases, prerequisite conditions like large rainfall intensity and high soil permeability are necessary to ensure a proper water dynamic range for GPR data acquisition. But the wetting range of soils under natural conditions typically regulated by precipitation controls the applicability of the inverse modeling approach (e.g., Steenpass et al., 2011; Scharnagl et al., 2011). Finally, to capture the lateral water redistribution, timing for the time-lapse GPR observations is essential.

Except for the measurement uncertainties in reflector depth and water content, other errors which were ignored in this study influence the application of this method. Three types of errors are summarized as follows. (1) The uncertainty in upper-boundary fluxes should be considered when relevant observations are not well conducted. (2) The structural errors in GPR observations from the multi-channel GPR evaluation and the selection of the investigated transect can be reduced with better GPR setup and survey design. (3) Serious conceptual errors might arise when ignoring small-scale heterogeneous soil properties within the layer. As pointed out in the Introduction, the proposed approach focuses on the stratal soil properties and effective soil water flow at the scale of meters to tens of meters. The layer with strong small-scale heterogeneities is unfavorable in this approach.

5 Conclusions

The surface layer with a flat surface and an undulating bottom at the scale of meters to hundreds of meters is widespread in naturally and anthropogenically affected fluvial depositional environments. This study demonstrates an approach for efficiently estimating effective stratal hydraulic properties of the surface layer. It employs multi-channel GPR to capture soil architecture and time-lapse soil water content, a 2-D simulation of soil hydrology, and an optimization procedure. The approach is applied to reclamation land with a surface layer over buried undulating landforms. A small number of time-lapse GPR observations were conducted after a heavy rain event, and effective hydraulic parameters were obtained using 2-D inverse modeling. The performance of the proposed approach for the given field data is mainly limited by the observed narrow water content range and structural errors. Further synthetic studies show that better performance can be achieved, e.g., by adding one more GPR observation to extend the water dynamic range.

Application of the demonstrated approach mainly relies on the slope gradient of the undulating structure and the lateral water redistribution. Other factors contributing to lateral flow, such as soil hydraulic properties and the intensity and duration of the precipitation, also influence the performance of this approach. Overall, the major advantages of this approach include (i) non-destructive observations, (ii) a bigger scale of the effective soil hydraulic properties and (iii) efficiency in field applications.