Abstract

A one-dimensional vadose zone model was used to simulate flow under natural boundary conditions. The effects of hysteresis and temporal variability of meteorological conditions were evaluated. Simulations were performed in HYDRUS-1D code for the period April 2013–January 2014 (6601 hours) at three different locations in a delineated portion of the sub-quaternary catchment A80A of Nzhelele with different soil textures. Soil hydraulic characteristics were estimated in a Rosetta library dynamically linked to the HYDRUS-1D model which is based on the numerical solution of a one-dimensional Richard's equation. Analysis of the simulation results suggests that ignoring hysteresis for soils of similar textural class does not lead to any significant deviation of the model predicted soil moisture, unlike for soils with different textural classes.

INTRODUCTION

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.

THEORY OF FLOW THROUGH THE VADOSE ZONE

In this paper, water flow in porous media is described by use of the Richard's equation in Equation (1).  
formula
(1)

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).

Figure 1

Schematic representation of the SWRC (Lehmann et al. 1998).

Figure 1

Schematic representation of the SWRC (Lehmann et al. 1998).

MATERIALS AND METHODS

Study area

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 km2; 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).

Figure 2

Geology and location of monitoring site in delineated catchment.

Figure 2

Geology and location of monitoring site in delineated catchment.

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).

Table 1

Soil textural properties for sampled sites

Depth (cm) Sand, % Silt, % Clay, % Bulk density g/cm3 
Site 1 
 0–30 44.22 3.03 52.76 1.1 
 30–80 53.85 11.1 35.05 0.8 
  >80 78.3 18.8 2.9 1.8 
Site 2 
 0–30 40.72 4.06 55.22 1.1 
 30–120 30.76 3.85 65.39 1.7 
  >120 93.1 5.1 1.8 1.56 
Site 3 
 0–30 82.2 11.6 6.2 1.55 
 30–120 73.4 9.1 17.5 1.53 
  >120 93.1 5.1 1.8 1.56 
Depth (cm) Sand, % Silt, % Clay, % Bulk density g/cm3 
Site 1 
 0–30 44.22 3.03 52.76 1.1 
 30–80 53.85 11.1 35.05 0.8 
  >80 78.3 18.8 2.9 1.8 
Site 2 
 0–30 40.72 4.06 55.22 1.1 
 30–120 30.76 3.85 65.39 1.7 
  >120 93.1 5.1 1.8 1.56 
Site 3 
 0–30 82.2 11.6 6.2 1.55 
 30–120 73.4 9.1 17.5 1.53 
  >120 93.1 5.1 1.8 1.56 
Table 2

Rosetta estimated soil hydraulic parameters for the three sites

  Layer θr (cmcm−3θs (cmcm−3α (cm−1n ( − ) Ks (cm/hr) 
SITE 1 0.1096 0.5612 0.0298 1.2732 1.95 
0.0997 0.633 0.0275 1.328 5.85 
0.034 0.2983 0.055 1.6275 1.68 
SITE 2 0.1106 0.5636 0.0294 1.2674 1.82 
0.0948 0.3874 0.0257 1.158 0.23 
0.0488 0.3677 0.0343 3.2081 23.2 
SITE 3 0.0447 0.3742 0.0386 1.8802 5.06 
0.0591 0.3936 0.0274 1.4679 1.69 
0.0488 0.3677 0.0343 3.2081 23.2 
  Layer θr (cmcm−3θs (cmcm−3α (cm−1n ( − ) Ks (cm/hr) 
SITE 1 0.1096 0.5612 0.0298 1.2732 1.95 
0.0997 0.633 0.0275 1.328 5.85 
0.034 0.2983 0.055 1.6275 1.68 
SITE 2 0.1106 0.5636 0.0294 1.2674 1.82 
0.0948 0.3874 0.0257 1.158 0.23 
0.0488 0.3677 0.0343 3.2081 23.2 
SITE 3 0.0447 0.3742 0.0386 1.8802 5.06 
0.0591 0.3936 0.0274 1.4679 1.69 
0.0488 0.3677 0.0343 3.2081 23.2 

Numerical modelling

Water flow and root water uptake were simulated using the finite-element numerical model HYDRUS-1D (Šimŭnek et al. 2005). A comprehensive model description including recent modifications can be found at Šimŭnek et al. (2008). Assuming thermal gradient is negligible enough to not influence water flow, that the air phase does not interact with water flow, and that flow is predominantly vertical, the governing equation for water flow is the Richards equation (Equation (1)) including the sink term, S (), described as the volume of water removed from a unit volume of soil per unit time due to plant water uptake, as given in Equation (2).  
formula
(2)
The sink term was specified in terms of a potential uptake rate and a stress factor according to Feddes et al. (1978).  
formula
(3)
where is a dimensionless water stress response function that prescribes the reduction in uptake that occurs due to drought stress. For this study, root water uptake due to water stress was described using the model introduced by Feddes et al. (1978).  
formula
(4)
where are threshold parameters such that uptake is at the potential rate when the pressure head is between and , and when pressure head values are between (or and ), root water decreases or increases linearly with h. Specific values for maize and grass were obtained from the database available in HYDRUS-1D (Šimŭnek et al. 2005). In HYDRUS-1D, is a function of the potential evapotranspiration and users can specify two different and , respectively as shown in Equation (5).  
formula
(5)
The soil hydraulic properties for the unsaturated soil were described using the set of closed form analytical solutions of van Genuchten-Mualem single porosity constitutive relationships (Mualem 1976; van Genuchten 1980).  
formula
(6)
 
formula
(7)
where saturated water content (); = residual water content (); = saturated hydraulic conductivity (); h is the pressure head ; and the empirical coefficients are air entry parameter (), pore size distribution, pore connectivity. In this study, we applied to reduce the number of free parameters based on the work of Mualem (1976). is the effective saturation given by:  
formula
(8)
In this study, a simplified scaling procedure of the soil hydraulic properties is available in HYDRUS1-D to describe the spatial variability of the soil hydraulic properties by means of a set of three linear scaling transformations which relate each individual soil hydraulic characteristics and to the reference characteristics and . These three independent scaling factors are embodied in HYDRUS through the following general equations (Vogel et al. 1991).  
formula
(9)
 
formula
(10)
 
formula
(11)
where , , and are mutually independent scaling factors for the water content, pressure head, and hydraulic conductivity, respectively. However, for this study, scaling was applied only to the hydraulic conductivity equation (Equation (9)) at each soil profile using a scaling factor computed from Equation (12) since it was the only parameter considered to differ from the main wetting and drying curves.  
formula
(12)
where is the arithmetic mean of saturated hydraulic conductivity obtained from Carsel & Parrish (1988) which is the default parameter of HYDRUS for van Genuchten functions and is the initial saturated conductivity obtained from Table 2. The computed scaling factor for the hydraulic conductivity for each layer in all profiles is shown in Table 3.
Table 3

Computed scaling factor for hydraulic conductivity for the three sites

Site Layer  
0.10 
0.02 
17.68 
0.11 
0.87 
1.28 
2.88 
8.63 
1.28 
Site Layer  
0.10 
0.02 
17.68 
0.11 
0.87 
1.28 
2.88 
8.63 
1.28 

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.

Meteorological data

Hourly precipitation, evaporation and transpiration during the simulation period which lasted for 6,601 hours (23:00 25th April 2013 to 23:00 25th January 2014) were used as time variable boundary conditions. Meanwhile the upper boundary condition involved hourly rainfall intensities and potential evapotranspiration rates. Potential evapotranspiration rates were calculated by the Penman-Monteith method. Using the calculated reference evapotranspiration , the potential evapotranspiration was derived from Allen et al. (1998).  
formula
(13)
is reference evapotranspiration in hourly time steps, is crop specific coefficient that characterises plant water uptake and evaporation relative to the reference crop (maize) and was set to be uniform spatially across the entire modelling site and only varied temporally with growing season.
With obtained, potential evaporation was calculated using Equation (14) as described in Kroes & Van Damm (2003) and Pachepsky et al. (2004).  
formula
(14)
where , is the radiation extinction coefficient and leaf area index in hourly time steps.
Knowing and from Equations (13) and (14), the potential transpiration was obtained using Equation (15).  
formula
(15)

Allen et al. (1998) provided the length of growth stages and values of kc 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.

Figure 3

Summary of modelled soil surface boundary conditions for all three sites.

Figure 3

Summary of modelled soil surface boundary conditions for all three sites.

Modelling hysteresis

This version of HYDRUS-1D (4.0) incorporates hysteresis in the flow model by employing the extended model of Kool & Parker (1987) which considers hysteresis not only on the retention curve but also in the hydraulic conductivity function (Vogel et al. 1996). The adopted procedure for modelling hysteresis in the retention curve required knowledge of both the main drying (d) and wetting (w) and curves using the parameters () and () in the modified form of the unsaturated hydraulic conductivity function with a pressure head range (Vogel et al. 2001). The initial condition for the hysteretic model was specified with the main wetting branch for the single porosity van Genuchten-Mualem model. Furthermore, it was assumed that, the only parameters differing from each other from the main wetting and drying curves at saturation were given the following restrictions:  
formula
(16)
 
formula
(17)
 
formula
(18)

Inverse modelling

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).

The performance of the optimisation procedure was expressed by means of the root mean square error (RMSE) and correlation coefficient (R2) which were chosen over the other popular criterion based on the fact that they can be used to compare data series with different numbers of measurements (Dohnal et al. 2006).  
formula
(19)
 
formula
(20)
where n is the number of measurements, y and are the observed and predicted responses respectively, and is the average value of observations.

Simulation scenarios

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.

Table 4

Optimised soil hydraulic parameters for the three sites

Site Layer θr (cmcm−3θs (cmcm−3α (cm−1n 
0.312 0.5845 0.051 1.749 
0.1318 0.7378 0.0677 1.265 
0.4282 0.615 0.0252 1.15 
0.2264 0.5164 0.011 1.5545 
0.511 0.715 0.0906 3.34 
0.0135 0.5642 0.0441 1.4725 
0.252 0.615 0.0376 3.14 
Site Layer θr (cmcm−3θs (cmcm−3α (cm−1n 
0.312 0.5845 0.051 1.749 
0.1318 0.7378 0.0677 1.265 
0.4282 0.615 0.0252 1.15 
0.2264 0.5164 0.011 1.5545 
0.511 0.715 0.0906 3.34 
0.0135 0.5642 0.0441 1.4725 
0.252 0.615 0.0376 3.14 

RESULTS AND DISCUSSION

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.

Figure 4

Observed water content data measured at various depths for site 1 along with simulated final fitted HYDRUS simulations (1, 2, 3, 4, 5, 6, represent 30, 60, 90, 120, 150, and 180 cm depths respectively).

Figure 4

Observed water content data measured at various depths for site 1 along with simulated final fitted HYDRUS simulations (1, 2, 3, 4, 5, 6, represent 30, 60, 90, 120, 150, and 180 cm depths respectively).

The performance of the model for all of the three sites was satisfactory as shown in their computed correlation coefficient (R2) and RMSE values given in Table 5. Values for R2 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.

Table 5

Goodness of fit measures for simulated and observed data

Site Depth (cm) R2 RMSE 
Profile 1 30 0.87 0.0040 
60 0.89 0.0015 
90 0.82 0.0035 
120 0.78 0.0038 
150 0.90 0.0046 
180 0.92 0.0041 
Profile 2 30 0.97 0.0110 
60 0.63 0.0057 
90 0.91 0.0060 
120 0.95 0.0065 
150 0.66 0.0003 
180 0.83 0.0021 
Profile 3 30 0.93 0.0043 
60 0.64 0.0051 
90 0.84 0.0034 
120 0.91 0.0040 
150 0.85 0.0043 
180 0.98 0.0020 
Site Depth (cm) R2 RMSE 
Profile 1 30 0.87 0.0040 
60 0.89 0.0015 
90 0.82 0.0035 
120 0.78 0.0038 
150 0.90 0.0046 
180 0.92 0.0041 
Profile 2 30 0.97 0.0110 
60 0.63 0.0057 
90 0.91 0.0060 
120 0.95 0.0065 
150 0.66 0.0003 
180 0.83 0.0021 
Profile 3 30 0.93 0.0043 
60 0.64 0.0051 
90 0.84 0.0034 
120 0.91 0.0040 
150 0.85 0.0043 
180 0.98 0.0020 

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.

Figure 5

Comparison of bottom flux for hysteretic and non-hysteretic flows at monitoring sites.

Figure 5

Comparison of bottom flux for hysteretic and non-hysteretic flows at monitoring sites.

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.

CONCLUSIONS

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.

ACKNOWLEDGEMENT

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.

REFERENCES

REFERENCES
Allen
,
R. G.
,
Pereira
,
L. S.
,
Raes
,
D.
&
Smith
,
M.
1998
Crop Evapotranspiration. Guidelines for Computing Crop Water Requirements. Irrigation and Drainage. Paper No. 56
,
FAO
,
Rome
,
Italy
.
Arrey
,
I. A.
2015
A Modelling Approach to Estimate Groundwater Recharge From Infiltration in the Unsaturated Zone: Siloam Village Case Study, South Africa
.
Unpublished Master Dissertation
,
Department of Hydrology and Water Resources, University of Venda
,
South Africa
.
Bashir
,
R.
,
Sharma
,
J.
&
Stefaniak
,
H.
2015
Effect of hysteresis of soil-water characteristic curves on infiltration under different climatic conditions
.
Canadian Geotechnical Journal
53
(
2
),
273
284
.
Beese
,
I.
&
Van Der Ploeg
,
R. R.
1976
Influence of hysteresis on moisture flow in an undisturbed soil monolith
.
Soil Science Society of America
40
(
4
),
480
484
.
Brandl
,
G.
1999
Soutpansberg Group. Catalogue of South African Lithostratigraphic Units, SA Committee for Stratigraphy, Council for Geosciences
, pp.
639
641
.
Carsel
,
R. F.
&
Parrish
,
R. S.
1988
Developing joint probability distributions of soil water characteristics
.
Water Resources Research
24
,
755
769
.
Dohnal
,
M.
,
Dušek
,
J.
&
Vogel
,
T.
2006
The impact of the retention curve hysteresis on prediction of soil water dynamics
.
Journal of Hydrology and Hydromechanics
54
,
258
268
.
Elmaloglou
,
S.
&
Soulis
,
K. X.
2013
The effect of hysteresis on soil water dynamics during surface trickle irrigation in layered soils
.
Global Nest Journal
15
(
3
),
351
365
.
Feddes
,
R. A.
,
Kowalik
,
P. J.
&
Zaradny
,
H. X.
1978
Simulation of Field Water Use and Crop Yield
.
Centre for Agricultural Publishing and Documentation, Wageningen
,
The Netherlands
.
Gillham
,
R. W.
,
Klute
,
A.
&
Heermann
,
D. F.
1976
Hydraulic properties of porous medium: measurements and empirical representation
.
Soil Science Society of America
40
,
203
207
.
Grossman
,
R. B.
,
Reinsch
,
T. G.
2002
Bulk density and linear extensibility
. In:
Methods of Soil Analysis, Part 4, SSSA Book Series
,
Vol. 5
(
Dane
,
J.
&
Topp
,
C.
, eds).
American Society of Agronomy
,
Madison, WI
, pp.
201
228
.
Heinse
,
R.
&
Link
,
T. E.
2013
Vadose zone processes: a compendium for teaching interdisciplinary modelling
.
Journal of Contemporary Water Research and Education
152
(
1
),
22
31
.
Hillel
,
D.
1980
Fundamentals of Soil Physics
.
Academic Press
,
New York, NY
.
Hoa
,
N. T.
,
Gaudu
,
R.
&
Thirriot
,
C.
1977
Influence of the hysteresis effect on transient flows in saturated–unsaturated porous media
.
Water Resources Research
13
(
6
),
992
996
.
Hopmans
,
J. W.
&
Stricker
,
J. N. M.
1989
Stochastic analysis of soil water regime in a watershed
.
Journal of Hydrology
105
,
57
84
.
Izady
,
A.
,
Ghahraman
,
B.
&
Davari
,
K.
2008
Hysteresis: phenomenon and modelling in soil-water relationship
.
Iran Agricultural Research
27
,
1
2
.
Jaynes
,
D. B.
1992
Estimating hysteresis in the soil water retention function
. In:
Proceedings of the International Workshop on Indirect Methods for Estimating the Hydraulic Properties of Unsaturated Soils
(
van Genuchten
,
M. Th.
,
Leij
,
F. J.
&
Lund
,
L. J.
, eds).
University of California
,
Riverside, CA
, pp.
219
232
.
Kaluarachchi
,
J. J.
&
Parker
,
J. C.
1987
Effects of hysteresis with air entrapment on water flow in the unsaturated zone
.
Water Resources Research
23
(
10
),
1967
1976
.
DOI: 10.1002/ird.266
.
Kerkides
,
P.
,
Kargas
,
G.
&
Argyrokastritis
,
I.
2006
The effect of different methods used for hysteretic K(H) determination on the infiltration simulations
.
Irrigation and Drainage
55
,
403
418
.
Kroes
,
J. G.
&
Van Damm
,
J. C.
2003
Reference Manual SWAP: Version 3.0.3. Rep. 773
.
Alterra Green World Research
,
Wageningen
,
The Netherlands
.
Leão
,
T. P.
&
Perfect
,
E.
2010
Modelling water movement in horizontal columns using fractal theory
.
Revista Brasileira de Ciȇncia do Solo
34
,
1463
1468
.
Lehmann
,
P.
,
Stauffer
,
F.
,
Hinz
,
C.
,
Dury
,
O.
&
Flühler
,
H.
1998
Effects of hysteresis on water flow in a sand column with a fluctuating capillary fringe
.
Journal of Contaminant Hydrology
33
,
81
100
.
Leterme
,
B.
,
Jacques
,
D.
,
Mallants
,
D.
,
Hooker
,
P.
,
De Craen
,
M.
&
Van den Hoof
,
C.
2012
Long-Term Climate Change and Effects on Disposal Facility, Geosphere, and Biosphere, Report NIROND-TR 2009-07E, Technical report, SCK•CEN, Mol
,
Belgium
.
Luckner
,
L.
,
van Genuchten
,
M. T.
&
Nielsen
,
D. R.
1989
A consistent set of parametric models for the two-phase flow of immiscible fluids in the subsurface
.
Water Resources Research
25
(
10
),
2187
2193
.
Mantoglou
,
A.
&
Gelhar
,
L. W.
1987
Stochastic modelling of large-scale transient unsaturated flow systems
.
Water Resources Research
29
,
3769
3774
.
Marquardt
,
D. W.
1963
An algorithm for least squares estimation of nonlinear parameters
.
Journal of the Society for Industrial and Applied Mathematics
11
,
431
444
.
Monteith
,
J. L.
1981
Evaporation and surface temperature
.
Quarterly Journal of the Royal Meteorological Society
107
,
1
27
.
Moriasi
,
D. N.
,
Arnold
,
J. G.
,
Van Liew
,
M. W.
,
Bingner
,
R. L.
,
Harmel
,
R. D.
&
Veith
,
T. L.
2007
Model evaluation guidelines for systematic quantification of accuracy in water shed simulations
.
American Society of Agricultural and Biological Engineers
50
(
3
),
885
900
.
Mualem
,
Y.
1974
A conceptual model of hysteresis
.
Water Resources Research
10
(
3
),
514
520
.
Nielsen
,
D. R.
,
Luckner
,
L. M.
1992
Theoretical aspects to estimate reasonable initial parameters and range limits in identification procedures for soil hydraulic properties
. In:
Proceedings of the International Workshop on Indirect Methods for Estimating the Hydraulic Properties of Unsaturated Soils
(
van Genuchten
,
M. T.
,
Leij
,
F. J.
&
Lund
,
L. J.
, eds).
US Salinity Laboratory, Agricultural Research Service, US Dept. of Agriculture
,
Riverside, CA
, pp.
147
160
.
Novác
,
V.
&
Gallová
,
C.
1998
Mathematical modelling of annual soil water content changes and sensitivity of the model to soil hydraulics characteristics
.
Journal of Hydrology and Hydromechanics
46
,
432
446
(in Slovak)
.
Pachepsky
,
Y. A.
,
Smettem
,
K. R. J.
,
Vanderborght
,
J.
,
Herbst
,
M.
,
Vereecken
,
H.
,
Wosten
,
J. H. M.
2004
Reality and fiction of models and data in soil hydrology
. In:
Unsaturated Zone Modelling
(
Feddes
,
R. A.
,
de Rooij
,
G. H.
&
van Dam
,
J. C.
, eds).
Kluwer Academic Publishers
,
Dordrecht
,
The Netherlands
.
Poulovassilis
,
A.
1962
Hysteresis of pore water, an application of the concept of independent domains
.
Soil Science Society of America Journal
93
(
6
),
405
412
.
Poulovassilis
,
A.
&
Kargas
,
G.
2000
A note on calculating hysteretic behaviour
.
Soil Science Society of America Journal
64
(
6
),
1947
1950
.
Royer
,
J. M.
&
Vachaud
,
G.
1975
Field determination of hysteresis in soil-water characteristics
.
Soil Science Society of America
39
,
221
223
.
Rumynin
,
V. G.
2011
Subsurface Solute Transport Models and Case Histories
.
Springer
,
New York
.
Santhi
,
C. J.
,
Arnold
,
J. R.
,
Williams
,
W. A.
,
Dugas
,
R.
,
Srinivasan
,
L.
&
Hauck
,
M.
2001
Validation of the SWAT model on a large river basin with point and non-point sources
.
Journal of the American Water Resources Association
37
(
5
),
1169
1188
.
Schaap
,
M. G.
,
Leij
,
F. J.
&
van Genuchten
,
M. T.
2001
ROSETTA: a computer program for estimating soil hydraulic parameters with hierarchical pedotransfer functions
.
Journal of Hydrology
251
,
163
176
.
Šimŭnek
,
J.
,
Šejna
,
M.
&
van Genuchten
,
M. T.
1999
The HYDRUS-2D software package for simulating two-dimensional movement of water, heat, and multiple solutes in variably saturated media. Version 2.0, IGWMC – TPS - 53. Golden, International Ground Water Modeling Center, Colorado School of Mines
, p.
251
.
Šimŭnek
,
J.
,
van Genuchten
,
M. T.
&
Šejna
,
M.
2005
The HYDRUS-1D software package for simulating the one-dimensional movement of water, heat, and multiple solutes in variably-saturated media. Version 3.0. HYDRUS Software Ser. 1. Dep. of Environmental Sciences
,
Univ. of California
,
Riverside
.
Šimŭnek
,
J.
,
van Genuchten
,
M. T.
&
Šejna
,
M.
2008
Development and applications of the HYDRUS and STANMOD software packages and related codes
.
Vadose Zone Journal
7
,
587
600
.
Topp
,
G. C.
1971
Soil–water hysteresis: the domain theory extended to pore interaction conditions
.
Soil Science Society of America Proceedings
35
,
219
225
.
Vogel
,
T.
,
Císlerová
,
M.
&
Hopmans
,
J. W.
1991
Porous media with linearly variable hydraulic properties
.
Water Resources Research
27
,
2735
2741
.
Vogel
,
T.
,
Huang
,
K.
,
Zhang
,
R.
&
Van Genuchten
,
M. T.
1996
The HYDRUS Code for Simulating one-Dimensional Water Flow, Solute Transport, and Heat Movement in Variably-Saturated Media, Version 5.0, Research Report No. 140
,
US Salinity Lab., ARS, USDA
,
Riverside
,
CA
.
Vogel
,
T.
,
Van Genuchten
,
M. T.
&
Císlerová
,
M.
2001
Effect of the shape of soil hydraulic functions near saturation on variably-saturated flow predictions
.
Advances in Water Resources
24
,
133
144
.