Effect of hysteresis on water flow in the vadose zone under natural boundary conditions, Siloam Village case study, South Africa

Flow in subsurface formations is usually described based on quantity of water in the intergranular medium. This could be either operated as saturated or unsaturated porous media. Typical of a natural state occurrence in a vertical profile is a zone of aeration whereby the pores contain both gases and water overlying a zone of saturation completely filled with water. The foregoing, also called the vadose zone or the unsaturated zone acts as a filter to unwanted substances that might originate from the ground surface such as pesticides, fertilisers, and hazardous waste. The dynamics of water in the unsaturated zone is directly linked to the hydrologic cycle by partitioning of water at the land surface and regulating movement of water to and from groundwater. As such, it effectively controls the interrelationships between precipitation, runoff, infiltration, evapotranspiration, and groundwater recharge (Heinse & Link 2013). This explains why studies that model water flow and solute transport in the unsaturated zone are increasingly becoming an issue of major concern in terms of water resources planning and management and groundwater contamination (Rumynin 2011). These models, usually represented by mathematical expressions, provide a rational and scientific basis for catchment management decisions related to water flow and solute transport in the vadose zone.

Over the past few decades, a large number of numerical models have been developed to compute water flow and solute transport in the vadose zone. These models usually require as forcing functions, meteorological input data, while some have included complex physical processes such as hysteresis, which have been shown to significantly influence water flow and solute transport in variably saturated flow (Kaluarachchi & Parker 1987; Jaynes 1992). Drying and wetting processes driven by water potential differences often occur simultaneously at different depths based on ambient field meteorological conditions. Hysteresis is the relationship that exists between the matric potential and water content which depends on the wetting and drying history of the soil (Hillel 1980). This phenomenon was long brought forward by Haines (1930) who recognised hysteresis as an important aspect of soil water dynamics. He was able to demonstrate that hysteresis largely depends on the geometric non-uniformity of the individual pores of soils. Since then a variety of models have been developed to describe the hysteretic behaviour of soils and have been broadly classified into empirical (e.g. Poulovassilis 1962; Gillham et al. 1976; Hoa et al. 1977) and theoretical (e.g. Topp 1971; Mualem 1974; Poulovassilis & Kargas 2000). The theoretical approach assumes the soil to be a complex continuum of dependent or independent pore domains while the empirical models, on the contrary, are a much simplified version of the domain models (Dohnal et al. 2006).

A number of studies have brought forward the argument that hysteresis is negligible for certain soil types or flow conditions (Kool & Parker 1987; Novác & Gallová 1998). Royer & Vachaud (1975) rated the effect of hysteresis on soil water movement as inferior relative to the influence of spatial variability of soil hydraulic properties, while Mantoglou & Gelhar (1987) considered changes in soil hydraulic conductivity as a result of spatial variability to be attributed to hysteresis. Beese & van der Ploeg (1976), in their study, were able to show significant improvements in soil moisture prediction dynamics when the effect of hysteresis was included. However, their study was done under undisturbed homogeneous soil with natural boundary conditions. Kerkides et al. (2006) were able to show that hysteresis largely affects vadose zone flow and depends on the choice of method used to determine the hysteretic properties of the porous media.

This study, aims to quantify the impact of hysteresis on model predictions of soil moisture under natural boundary conditions when simulating temporal changes in meteorological input data for undisturbed soil profiles. The outcome of the study will contribute to the discussion on the importance of hysteresis in soil water flow models.

With the vertical coordinate z (⁠⁠) being positive downwards, the matric potential head ψ (⁠⁠), and the hydraulic conductivity k (⁠⁠) which depends on the volumetric water content θ (⁠⁠). At saturation and otherwise the nonlinear relation between θ and k must be defined. This one-dimensional partial differential equation describes the movement of soil moisture through unsaturated porous media under appropriate boundary and initial conditions.

The relationship between the matric potential (ψ) and volumetric water content (θ) and that of the matric potential (ψ) and hydraulic conductivity (k) constitute the soil water retention curve (SWRC) and hydraulic conductivity function respectively. The two complete characteristic curves from saturation to dryness and vice versa form the main branches of the hysteretic soil moisture characteristics (Izady et al. 2008). Partially wetting and drying conditions, as it moves from one main branch to another, are described as scanning curves, as shown in Figure 1. There are two main groups of hysteresis models, the first one describes the scanning curves with expressions similar to the main wetting and drying curves while the second explains the shape of the scanning curves by means of the physical properties of the soil (Lehmann et al. 1998). A comprehensive review on the hysteresis process can be found in Mualem (1974), Kool & Parker (1987), Luckner et al. (1989), Jaynes (1992) and Šimŭnek et al. (1999).

The study site, Siloam Village, falls under the delineated portion of the quaternary catchment A80A of Nzhelele River Catchment. It is located on the leeward side of the Soutpansberg Mountains in the northern region of Limpopo Province, South Africa. It has a drainage area of 26 km²; an elevation gradient of 800–860 m above mean sea level; temperatures fluctuate between summer and winter with records of up to 40°C during summer and 16°C in winter and the mean annual evaporation varies between 1,300 and 1,400 mm. It has average rainfall of 350–400 mm/annum which is predominantly received during the summer season (October to March). The area is extensively covered with areas of human settlement and patches of land used for subsistence farming, which is primarily the maize crop.

The soil moisture dynamics was monitored at five locations by means of continuous logging Neutron Probes (DFM, South Africa). However, analysis was based only on three probe sites (20922, 12699, and 20918) chosen based on their spatial distribution shown in Figure 2. At each location, every installed probe is equipped with a nest of six sensors measuring soil moisture at depth increments of 30 cm up to a depth of 180 cm at hourly intervals, though these measuring depths did not correspond to the soil horizons in the area. Meteorological variables such as rainfall, temperature, wind speed and net radiation were observed at the site with a Vantage Pro2 weather console and these observations were used to compute hourly potential evapotranspiration fluxes using the Penman-Monteith method (Monteith 1981).

The area consists of volcanic assemblages at the base which is about 400 m thick and overlain by argillaceous and arenaceous sediments (Brandl 1999). A vertical soil profile for each of the three sites was characterised after a granulometric analysis which revealed three layered profiles, used as input for HYDRUS-1D (Table 1). HYDRUS-1D code is coupled with a Rosetta Dynamically Linked Library (Schaap et al. 2001), which implements pedotransfer functions (PTFs) to predict van Genuchten (1980) water retention parameters and the saturated hydraulic conductivity (⁠⁠) in a hierarchical manner from soil textural class to produce the initial estimates in Table 2. The bulk density was calculated from a volume-mass relationship from each core sample recovered at respective depths (Grossman & Reinsch 2002).

The total depth of each profile was 180 cm, covering the length of the soil moisture probes. The model domain was discretised into 30 nodes with unequal distances at the layer boundaries. The model boundary conditions were set as atmospheric with surface layer and deep drainage for upper and lower boundaries respectively. Vertical drainage across the lower boundary of each soil profile was defined according to Hopmans & Stricker (1989). Parameters for this equation were obtained from Leterme et al. (2012) except for the position of groundwater level which was set constant throughout the simulation period at 21 m, based on the findings of Arrey (2015). The soil moisture values observed at the beginning of the simulation periods at each site were used to specify the initial conditions for the simulation.

Allen et al. (1998) provided the length of growth stages and values of k_c and LAI for maize and grass crops. It is important to note that land cover during simulation period is typically made of patches of fallow land with settlement and gardens of maize plant and grass. The atmospheric boundary conditions for the entire simulation period summarised in Figure 3 show the hourly values of precipitation.

In order to obtain the optimal soil hydraulic parameters for each simulation scenario described in the next section, an inverse modelling was performed. The objective function used in our study for the parameter optimisation process consisted of the water content at multiple depths (0–30, 30–60, 60–90, 90–120, 120–150, 150–180 cm). Minimisation of the objective function was accomplished using the Levenberg-Marquardt nonlinear minimisation method for each depth increment for the three sites (Marquardt 1963).

Two major scenarios were implemented in order to assess the effect of hysteresis for unsaturated zone flow. The first scenario assumed no hysteresis in the soil hydraulic properties (moisture retention curve and hydraulic conductivity function) implying the main drying curve is identical to the main wetting curve while the second scenario involved hysteresis in both Equations (3) and (4) while assuming ⁠. Hence for the hysteretic scenario, the air entry parameter in the main drying curve was twice that in the main wetting curve, a common assumption used (Kool & Parker 1987; Nielsen & Luckner 1992; Šimŭnek et al. 1999). Flow was first simulated with the estimated hydrodynamic soil hydraulic properties given in Table 2. Thereafter, it was simulated with the optimised soil hydraulic parameters after inversion, given in Table 4.

Running HYDRUS-1D with the estimated hydraulic parameter estimates from Rosetta (Table 2) resulted in a significant amount of residuals between observed and simulated data. The model was therefore calibrated for the soil hydraulic properties using measured soil moisture data for 1,700 hrs (14:00 6th June 2013 to 17:00 16th August 2013) as shown in Figure 4. Initially, all the van Genuchten parameters were subjected to calibration except for the pore connectivity parameter ⁠. Several possible parameterisations were considered, varying the number of soil layers and type of hydraulic parameters. The best overall parameterisation of the calibration process shown in Table 3 was obtained based on diagnostic information from the output of HYDRUS-1D routines on the model fit and convergence behaviour of the algorithm as well as visual inspection of the model fit to observed data.

The performance of the model for all of the three sites was satisfactory as shown in their computed correlation coefficient (R²) and RMSE values given in Table 5. Values for R² were greater than 0.5 which are considered as acceptable (Santhi et al. 2001). The error indices (RMSE) for all three sites were close to 0 which also indicates a good model fit (Moriasi et al. 2007). When the results of the optimisation looked promising, the simulations were repeated using different parameter estimates to check if the same results would be achieved. If two or more parameter sets obtained a similar fit to the observed data, the one with the fewer fitted parameters was chosen. The best parameterisation was found to have three layers for site 1, and two layers each for site 2 and site 3 as shown in Table 4. The inverse algorithm failed to converge for sites 2 and 3 using three layers.

The observed and simulated water content for the top soil at site 1, shown in Figure 4, displays a good match at the 0–30 cm depth. Similar behaviour is observed at 30–60 cm depth except for the periods between 1,500 and 2,583 hours where the water content is underestimated. The observed fluctuations could be possibly attributed to numerical instability. At the 60–90 cm depth, the model fluctuated between 1,100 to 1,200 and 1,400 to 1,700 hours but continued with a good fit for the rest of the simulation period. The subsoil had an even better fit of the simulated water content with the observed water content, with the exception of the 120–150 and 150–180 cm depths where there was slight underestimation of the soil water content at the beginning of simulation. These discrepancies are possibly due to the distinct hydraulic characteristics of each layer causing instabilities before equilibrium is achieved. A study conducted by Leão & Perfect (2010) to model water movement in horizontal columns using the fractal theory, found that another soil hydraulic parameter, soil sorptivity, varied as the content of coarse-grained material increased in porous media. This increases the pressure head gradient in a layered soil profile leading to nonlinearity in infiltration rate. Lehmann et al. (1998) in their study also attributed these differences to local heterogeneities which could lead to inaccurate determination of hydraulic properties at low matric potential values or uncertainties in measurements from sensing devices.

With the fitted parameters, HYDRUS-1D was next used to predict the root zone dynamics of soil moisture fluxes during the entire simulation period at the three sites for both hysteretic and non-hysteretic scenarios, as shown in Figure 5. At site 1 and site 2, both hysteretic and non-hysteretic flow had similar flux rate patterns throughout the simulation period. These two sites have similar soil class with clay dominating the top soil and sand underlying it. These results indicate that not incorporating hysteresis in soils of the same taxonomy lead to a slightly lower flux rate compared to using hysteretic flow.

This similarity could be attributed to the presence of a similar soil stratification for both profiles with the presence of fine grain dominating the top soil. This is in line with the results of the studies by Dohnal et al. (2006), and Novác & Gallová (1998) who considered spatial variability of soil hydraulic properties to be more influential than hysteresis on soil water movement. Unlike sites 1 and 2, site 3 had a completely different soil stratification (coarse topsoil) make up and an expected deviation from the hysteretic and non-hysteretic flux rate pattern was observed. A significant increase in flux rate was observed during hysteretic flow, unlike during non-hysteretic flow. In other words, the draining process occurs faster when hysteresis is considered than when neglected. These results are similar to that of Elmaloglou & Soulis (2013), who showed that the distribution process in soil water movement progresses more quickly when hysteresis is considered than when not considered for stratified soils with a coarser topsoil. However, their study was simulated under surface drip irrigation from equidistant line sources. Kaluarachchi & Parker (1987), on the other hand, found hysteresis to be more sensitive to surface and boundary conditions than to the heterogeneous nature of soils while concluding that the effect of hysteresis on vadose zone flow is problem specific. Similarly, Bashir et al. (2015) found initial wetting and drying conditions to be a critical factor in predicting vadose zone flow rate accurately but also cautioned about the importance of an appropriate temporal resolution for the climate data set which in their case was yearly. This aspect of initial wetting and drying conditions and their effect on infiltration rate was earlier reported by Kerkides et al. (2006) and agrees with the use of wetting branches for a better representation of soil moisture profiles. Although we are in agreement that initial wetting conditions provide the most appropriate means of replicating field conditions when hysteresis is considered during infiltration, there seems to be no clear indication of how the layering structure in stratified soils and their taxonomic group alters vadose zone flow.

The results of this study comprised a series of forward and inverse simulations (and only a small fragment is presented in this paper). The results reaffirm the significance of hysteresis in subsurface flow as well as the spatial variability of soil hydraulic properties. The results of this study show that neglecting hysteresis in the soil hydraulic properties for soils of the same class is less significant than doing so for soils of different classes.

This study investigated the impact of hysteresis on vadose zone flow in natural boundary conditions. Analysis of the simulation results suggests that hysteresis in the retention curve and hydraulic conductivity is not a dominant factor for soils of the same class. In contrast, soils of different classes (texture) might have a substantial deviation in the model prediction of the soil water dynamics of the vadose zone between hysteretic and non-hysteretic flows.

According to the results obtained, it can be concluded that hysteresis significantly influences soil water movement in the vadose zone for a semi-arid surface boundary condition in stratified soil profiles of different textural class. The study therefore suggests hysteresis should be considered in similar studies with different soil classes in order to achieve adequate results.

This research was funded by the Research and Publications Committee (RPC) which is under the Directorate of Research and Innovation of the University of Venda, South Africa.