Introduction

The unsaturated soil zone refers to the area between the land surface and the groundwater table, within which the water content is lower than at full saturation, and the pressure is lower than atmospheric pressure (Hillel et al. 2005; Stephens et al. 2018). This soil zone is of great importance both for the vertical transport of water in the course of groundwater recharge or the saturation of dikes and for the retention of substances in the course of adsorption and degradation processes (Alam et al. 2021; Dillon et al. 2009).

The hydraulic properties of unsaturated soil, specifically the relationships between volumetric water content and matric potential, can be described by soil hydraulic characteristic curves such as the soil water retention curve (SWRC) (Hillel et al. 2005; Plate and Zehe 2008; Schanz 2007; Vereecken et al. 2016). The relationship between soil water tension and soil moisture is characteristic of the pore size distribution and is a crucial basis for simulating the water flow in the unsaturated soil zone as it contains important information related to the water that can be retained in the pores (Hopmans et al. 1999; Gupta et al. 1979; Le Bourgeois et al. 2016).

The determination of the SWRC can be done by different direct methods in laboratory and field scale. The most commonly used tests are the evaporation method (saturated soil samples are drained by evaporation, sample weight and matrix potential are recorded (Peters and Durner 2008; Schindler and Müller 2006; Wind 1966)), the overpressure method (a positive pressure is applied to the air phase compared to the atmospheric water pressure, the soil samples are dewatered (Delage & Cui 2008)), dew point method (determination of water potential by equilibrating water phase in the soil sample and vapor pressure in the air (Maček et al. 2013)) and suction plate method (saturated soil samples placed on porous plates are dewatered to a certain negative pressure by means of hanging water column (Osinski et al. 2016)).

A parameterization of the determined SWRC is performed using empirical models. For example, a well-known model for parameterization of sandy soils was developed by Brooks and Corey (Brooks and Corey 1964). However, the most widely used empirical-statistical approach for determining retention parameters is that of van Genuchten (Van Genuchten 1980; Abkenar et al. 2019) according to the following equation:

$$\theta \left( \psi \right)\; = \;\theta r + \frac{{\left( {\theta s + \theta r} \right)}}{{[1 + (\alpha *\psi )^{n} ]^{m} }}$$
(1)

\(\theta \left( \psi \right)\) = volumetric water content depending on tension (m3 m−3).

\(\psi\) = tension (m).

\(\theta_{r}\) = residual water content (m3 m−3).

\(\theta_{s}\) = saturated water content (m3 m−3).

n = van Genuchten shape parameter, measure of pore size distribution (-).

\(\alpha\) = van Genuchten shape parameter, air entry parameter (m−1).

m = transformation parameter (-).

The parameters determined on the basis of the SWRC are used for any calculations of water movements in the unsaturated soil zone. For example, these parameters are the starting point for numerical simulations of the time variable water content and matric potential in the unsaturated soil zone (Ket et al. 2018; Qiao et al. 2019). The performance of the numerical simulations relies on the accuracy of the information from hydrological parameters. The better these are characterized, the more efficient the calculation and the more accurate the simulation result.

However, the determined SWRC often provides only local information about the hydraulic properties of the unsaturated soil and is not representing the whole soil system. Especially in laboratory and field investigations on a larger scale, heterogeneity’s over the complete unsaturated soil body are often existing due to the variability of soil pore structure. (Iiyama 2016). Therefore, the local determination of the soil water retention curves often leads to an over- or underestimation of the parameters describing the matric potential—water content relationship. In case of using for numerical simulations, this parameter uncertainty leads initially to worse fit between measured and simulated data and following to greater effort for calibration (Ghanbarian et al. 2010).

Therefore, it is more optimal to obtain these input parameters directly from the system to be simulated and to overcome this issue by using more global information about the soil parameters. For this reason, an experiment was started to show how retention parameters determined on the basis of local and global measurements differ and how they influence the simulation. The following questions should be answered by the results of this investigation:

  • How big is the variance between the hydraulic soil parameters estimated on the basis of local and global determined soil water retention function in an unsaturated soil system?

  • Is there a difference depending on the soil type?

  • What is the influence of these parameters on the calibration procedure?

Methods

Experimental setup

The overall objective of the performed dike experiments is to investigate the influence of changing water contents in the unsaturated part of dikes, mainly influenced by precipitation events, on the slope stability as well as erosion phenomena during occurring flood events. Therefore, a rectangular-shaped, Plexiglas tank with dimensions (L × W × H) of 3.4 m × 1.0 m × 1.3 m) is constructed to address the study aims (Figs. 1A and 2B). Different scenarios with different flooding and precipitation events are realized.

Fig. 1
figure 1

Sketch of the experimental setup including the positions of soil moisture sensors and pressure head sensors (A) and picture of the dike experiment exemplary with 20 cm flooding (B)

Full size image

For answering the questions regarding the soil water retention curves, the following experimental setup is used:

  • A homogeneous dike consisting of one soil material with extends (L × W × H) of 2.7 m × 1.0 m × 0.6 m is built on the bottom of the tank (Figs. 1A and 1B), and the slope gradient is 1:2.

  • Flooding in stages from one side is occurred, the water level is raised after reaching steady-state conditions, and the heights of the stages were 20, 30, 40 and 48 cm.

  • For the estimation of spatial and temporal distribution of soil moisture, the dike body is equipped with FDR probes (Hydraprobes, Co. Stevenswater, Portland, OR, USA) at seven different observation points (Fig. 1A).

  • To record qualitative changes in the seepage line, gauge pipes are installed at four positions (Fig. 1A), in which the water level is measured using water-level loggers (Levellogger, Co. Solinst, Georgetown, Canada); in addition, the water level in the inflow is measured.

Two experiments with two different soil materials (soil 1: pure sand (Ss), soil 2: weakly loamy sand (Sl2) according to the German Standards Institute (GSI) norm number 4220 (2020), for further characterization, see Table 1) are performed, to identify the influence of the fine grain content on the determination of the retention parameters.

Determination of the soil water retention curve

The retention parameters are determined using two different methods, which differ in the different size of the soil system considered for the measurement.

Local measurement

The determination of the soil water retention curve, respectively, the pF curve, is carried out by retention behavior laboratory experiments (Hyprop, Co. METER Group, Munich, Germany). The principle is based on the laboratory evaporation method according to Wind/Schindler. In this method, the unsaturated hydraulic conductivity and water retention properties are determined as a function of tension and water content of soil samples (Peters and Durner 2008; Schindler and Müller 2006; Wind 1966).

Several soil samples are taken from the dike body using a soil sample ring and subsequently saturated with water. The saturated soil samples are dewatered by evaporation, the changing sample weight of the soil sample, open to the top, is recorded by means of a balance, and the matrix potential is recorded by means of two tensiometers installed in the soil samples. The parameterization of the hydraulic characteristic curves was carried out by means of the van Genuchten model (Abkenar et al. 2019).

Global measurement

The aim of the global measurement is to determine a soil water retention curve directly from measured data of the laboratory dike experiment without disturbing the soil properties.

The water content data required for this is obtained directly through the FDR probes installed at seven locations in the dike. For the determination of the matrix potential ψ, the use of an indirect measurement is necessary. The binding of water to the soil matrix, the matric potential ψ, represents the pressure head difference between the soil water and the free water table and thus corresponds to a negative hydrostatic pressure. If the influence of evaporation and possibly other source or sinkage terms is neglected, the matric potential ψ at an observation point depends only on the vertical distance to the free water level, in this case, the seepage line. This distance corresponds in hydrostatic equilibrium to pressure in, e.g., m water column and can then be converted into a matric potential according to:

$$- \Delta \left( {OP - SPL} \right) = 98.0638*\psi .$$
(2)

Δ(OP-SPL) = vertical distance from observation point (OP) to the free water level (SPL) (m).

Ѱ = matric potential (hPa).

and

$$pF = \log_{10} \left( \psi \right).$$
(3)

pF = decadic logarithm of the matric potential (-).

Ѱ = matric potential (hPa).

These conditions are fulfilled in the dike experiment carried out with different, constant flooding levels so that steady-state conditions are formed after sufficient time according to the prevailing flooding. The level of the seepage line is obtained directly through the gauge pipes installed at four positions in the dike.

Numerical investigation-set up

A two-dimensional model has been realized while using PcSiWaPro (Gräber et al. 2006), applying Richards equation (Richards 1931) for unsaturated flow and Mualem–van Genuchten retention and relative hydraulic conductivity model (Mualem 1976; Schaap & Genuchten 2006).

The numerical model is conceptualized based on the given geometry of the tank experiment (Fig. 1A). For the two-dimensional numerical model (3488 calculation nodes with 6649 elements), a one-side flooding at different water levels (20, 30, 40 cm) was considered for the dike body. For this purpose, a time-dependent potential head as boundary condition (defined according to the water level) on the “water side” and a seepage surface (discharge possible from a pressure head of 0 m) on the “air side” is implemented. The impermeable bottom is represented by a no-flow boundary condition. Existing observation points (4 × pressure head, 7 × water content) were also integrated into the model so that every 2 min measured data for pressure head and water content can be compared with the simulated data. Together with the balanced inflow and outflow, these serve as the basis for model calibration and validation.

Within the calibration process, an inverse analysis based on observed pressure heads and water contents during the flooding of the dike is performed by PcSiWaPro internally to optimize the initial soil parameters determined by local and global measurement (Table 1). This serves for the minimization of a suitable objective function, which expresses the discrepancy between the observed values and the predicted system responses. The system response is represented here by the numerical solution of the Richards equation, augmented with single porosity hydraulic model (Mualem–van Genuchten) and used boundary conditions. Initial estimates are then iteratively improved during the minimization process until the desired degree of precision is obtained. For the optimization of the objective function, a weighted least-squares approach based on the Levenberg–Marquardt nonlinear method is used in PcSiWaPro (Marquardt 1963).

Results and discussion

Determination of the soil water retention curve

For the two soils used for the experiments, the determination of the soil water retention curve occurs both by means of local and global measurement. The results obtained show a clear difference between the local and global determined soil water retention curve and the related van Genuchten parameters, respectively (Fig. 2A, B).

Fig. 2
figure 2

Local and global determined soil water retention curves—measured (dots) and by means of the van Genuchten model (lines) for soil1 (A) and soil 2 (B)

Full size image

In particular, the different values of n and ɑ in the van Genuchten function lead to a shift of the curves in the diagram. The residual water contents also show bigger differences, with the global values being lower than the local values for both soils (Table 1). In case of soil 2, volumetric water contents smaller than 0.22 m3 m-3 cannot be observed in the flooding experiment, so that data points for the calculation of the global soil water retention curve are not available for this range. The reason for this is the relatively high water saturation up to the area of the upper observation points due to the strongly pronounced capillary rise of water already for the lowest flooding stage. For soil 1, volumetric water contents and data points, respectively, up to the range of 10 percent could be measured due to the coarser soil texture accompanied by a much lower capillary rise.

Table 1 Hydraulic parameters of the used soil materials—determined by local and global measurement
Full size table

Simulation results

Hydraulic soil parameters originating from local and global measurement—soil 1

The visual comparison of measured and simulated water contents in the dike body (Fig. 3) is clearly showing that the use of the hydraulic soil parameters determined by the global measurement compared to those determined by the local measurement leads to a better agreement between measured (green line) and simulated (global—red line, local—blue) water contents. Furthermore, it can be found that both for the temporal change of the variable water contents and the level of the water content itself the fit of the simulated to the measured water contents is better in most observation points by using the global determined soil parameters.

Fig. 3
figure 3

Observed (green) and simulated water content for soil 1 at the observation points H1 (A), H3 (B), H4 (C) and H7 (D) by using basic hydraulic soil parameters originating from local (blue) and global (red) measurement as well as after calibration process (black)

Full size image

To evaluate the model performance besides visual comparison, two error indicators are calculated to assess the quality of the conformity between observed and simulated water contents (Golmohammadi et al. 2014) (Table 2). Root-mean-square error (RMSE), one of the commonly used error indicators (Gupta et al. 2009; Wöhling et al. 2013), shows how well the simulated data represent the pattern of the observed data. In principal applications, the lower the RMSE, the better the model performance. One of the drawbacks of using RMSE is that it does not differentiate between over- and underestimated values. If there are some values which possess large error due to any faulty instrument or any wrong technical procedure, it will trigger an overall large RMSE value (Chai and Draxler 2013; Willmott and Matsuura 2005).

Table 2 Model fit to observation regarding water content by using hydraulic soil parameters originating from local, global and calibration measurements—soil 1
Full size table

In addition, NSE (Nash–Sutcliffe) error indicator is used to assess the quality of the conformity between observed and simulated water contents. It is a proven method to depict model performance especially for hydrological related operations (Knoben et al. 2019; Bai et al. 2009; Krause et al. 2005; Pachepsky et al. 2016). The optimal value of NSE is 1, indicating that both the simulated and observed values are 100% in accordance. In the case of values > 0 it shows a coherence between the simulated and observed values and can be considered as good-to-moderate fit. If the NSE < 0, it means the simulated values are no way near the observe values (Knoben et al. 2019; Houska et al. 2014; Schaefli and Gupta 2007).

The results can confirm the findings obtained from the visual observation (Table 2). In all observation points, except for point H7, the determined error indicators show a better agreement between measured and simulated water contents when using the global determined soil parameters. Especially for the observation points H1 to H5, a significant difference between the simulation results based on the different hydraulic soil parameters can be demonstrated. Looking at the RMSE values, the fit can be rated as good, when compared to the general accuracy of the FDR probes (± 0.025 cm3 cm-3 as given by the manual (Stevens Water 2021). The simulated water contents based on the initial hydraulic soil parameters determined by the local method are much further away from agreement with the measured water contents.

The calibration performed based on the hydraulic soil parameters determined by the local and global method results in an improved fit of the simulated to the measured water contents in almost all observation points (Table 2). However, the improvement for almost all observation points (H1 to H5) is much higher for the model based on the local determined hydraulic soil parameters. The bigger the deviation of the simulated from the measured water values, the greater the change in the error indicators due to the performed calibration.

Thus, reduction in RMSE up to 84 percent can be recorded compared to the initial simulation. In contrast, the reduction in RMSE by the performed calibration for the model based on global determined hydraulic soil parameters is only up to 57 percent (H7), whereby for most of the observation points, the reduction is much smaller (H1–H6).

The final hydraulic soil parameter set is presented in Table 3; noticeable differences can be identified compared to the local and global laboratory data. The results show quite clearly that the initial global determined hydraulic soil parameters were quite close to the best fit hydraulic soil parameters than the initial local determined hydraulic soil parameters (Table 3).

Table 3 Final hydraulic soil parameter set for soil 1 based on calibration performance
Full size table

Hydraulic soil parameters originating from local and global measurement—soil 2

When comparing the simulation results for the water content based on the hydraulic soil parameters determined by the global (Fig. 4, red line) and local (Fig. 4, blue line) measurement, the differences to the measured values (Fig. 4, green line) are very small for both methods compared to soil 1.

Fig. 4
figure 4

Observed (green) and simulated water content for soil 2 at the observation points H1 (A), H3 (B), H4 (C) and H6 (D) by using basic hydraulic soil parameters originating from local (blue) and global (red) measurement as well as after calibration process (black)

Full size image

This also applies for both the temporal change of the variable water contents and the level of the water content itself, where the goodness of the fit of the simulated to the measured water contents is similar for both methods.

The results of the statistical analysis can confirm these findings (Table 4). In all observation points, except for point H3, the determined error indicators are showing similar agreement between measured and simulated water contents when using the local and global determined hydraulic soil parameters.

Table 4 Model fit to observation regarding water content by using hydraulic soil parameters originating from local, global and calibration measurements—soil 2
Full size table

Even the opposite effect can be observed for point H4 and H5, the fit is here a little bit better when using the determined local hydraulic soil parameters.

The calibration performed based on the hydraulic soil parameters determined by the local and global method also results in an improved fit of the simulated to the measured water contents in all observation points (Table 4). There is a similar adaptation of the simulated to the measured water contents both for the water contents determined on the basis of the local and the global method. Thus, reduction in RMSE up to 68, respectively, 66 percent (H7) can be recorded compared to the initial simulation, whereby for most of the observation points, the reduction is smaller (H2, H3—global, H4—local, H5, H6).

In case of soil 2, it cannot be demonstrated that the initial global determined hydraulic soil parameters are quite closer to the best fit soil hydraulic parameters than the initial local determined hydraulic soil parameters (Table 5). Noticeable differences to the final hydraulic soil parameter set are similar for both methods of determination.

Table 5 Final hydraulic soil parameter set for soil 2 based on calibration performance
Full size table

Discussion

The results are indicating that a determination of the soil water retention function based on a local sampling cannot always sufficient represent the soil hydraulic parameters of a larger soil system. Especially in soils with lower content of fine particles, such as for soil 1, the difference between the locally determined hydraulic soil parameters and the hydraulic soil parameters characterizing the system under consideration best is relatively high.

Reason for this observation is a varying heterogeneity of the soil pore structure including the continuity of the pores and their arrangement, which have a crucial influence on the water movement in the unsaturated soil zone. The local soil heterogeneities influencing the compaction as well as the sampling of the soil are decreasing with an increasing proportion of fine particles in the soil material what is leading to a reduction in the heterogeneity of hydraulic soil parameters (Kuckelkorn 2005). For this reason, determined soil water retention function based on a local sampling can represent better the hydraulic soil parameters of a larger soil system in case of soils with higher content of fine particles due to small deviation of the hydraulic soil parameters from the parameters characterizing the system under consideration best. In case of soils with low content of fine particles, the greater heterogeneity of the soil pore structure leads to an over- or underestimation of the parameters describing the hydraulic properties of the soil when determining soil water retention function based on local sampling.

Since the hydraulic soil parameters are often the starting point for numerical simulations, it means a higher effort for adapting the simulated to the measured values in terms of calibration when using soil hydraulic parameters determined by the local method, especially for soils with lower content of fine particles. More simulation time and a higher demand for computing power is necessary due to the bigger deviation of the hydraulic soil parameters from the best fit parameters. Using the parameters determined by the global method cuts down the need of computation power and simulation time, especially for large-scale models, because the hydraulic soil parameters are closer to the best-fit parameters.

Conclusions

The findings of the study presented are showing that in case of increased local heterogeneity of the soil pore structure, occurring more often in soils with a low content of fine particles, it leads to uncertainties or incorrect assumption of hydraulic soil parameters when determining them based on local sampling. This is proven by the simulation results and the performed statistical analysis obtained on the basis of the different soil hydraulic parameters. The use of the hydraulic soil parameters determined by the local measurement leads to much worse fit between simulated and measured values compared to those determined by the global measurement. These differences between local and global determined hydraulic soil parameters cannot be observed for soils with higher content of fine particles due to the reduced local soil heterogeneities in such soils.

Based on the study results, it is preferable to use hydraulic soil parameters based on the global method for soils with increased local heterogeneity of the pore structure, because they are describing better the soil characteristics in a larger soil system. In particular, when using the hydraulic soil parameters for following numerical simulations inclusive calibration process, their use can reduce computation power and simulation time, in particular for large-scale models. In case of soils with higher proportion of fine particles, the use of the hydraulic soil parameters determined by the local method is sufficient due to the relatively good agreement of the hydraulic soil parameters determined by both methods.

Repeated measurements with the same boundary conditions have shown that the measured water contents and pressure heads are subject to fluctuations. Also the settling of the soil, what is changing the compaction and the hydraulic properties of the soil over the time is influencing the results. Therefore, the parameterization of the soils is underlying to some uncertainties. For this reason, further experiments should be performed to validate the results and to discuss the observed deviations between individual experiments in the course of an uncertainty or sensitivity analysis.

In addition, the influence of the fine particles on the determined hydraulic soil parameters obtained on the basis of local measurement should be checked by conducting further experiments with other soil types. Furthermore, it is certainly useful to apply the gained knowledge when obtaining data from other large soil systems.