Estimating hydraulic conductivity of internal drainage for layered soils in situ

The soil hydraulic conductivity ( K function) of three layered soils cultivated at Paradys Experimental Farm, near Bloemfontein (South Africa), was determined from in situ drainage experiments and analytical models. Pre-ponded monoliths, isolated from weather and lateral drainage, were prepared in triplicate on representative sites of the Tukulu, Sepane and Swartland soil forms. The first two soils are also referred to as Cutanic Luvisols and the third as Cutanic Cambisol. Soil water content (SWC) was measured during a 1200 h drainage experiment. In addition soil physical and textural data as well as saturated hydraulic conductivity (Ks) were derived. Undisturbed soil core samples of 105 mm with a height of 77 mm from soil horizons were used to measure soil water retention curves (SWRCs). Parameterization of SWRC was through the Brooks and Corey model. Kosugi and van Genuchten models were used to determine SWRC parameters and fitted with a RMSE of less 2 %. The SWRC was also used to estimate matric suctions for in situ drainage SWC following observations that laboratory and in situ SWRCs were similar at near saturation. In situ K function for horizons and the equivalent homogeneous profiles were determined. Model predictions based on SWRC overestimated horizons K function by more than three orders of magnitude. The van Genuchten–Mualem model was an exception for certain soil horizons. Overestimates were reduced by one or more orders of magnitude when inverse parameter estimation was applied directly to drainage SWC with HYDRUS-1D code. Best fits ( R2 ≥ 0.90) were from Brooks and Corey, and van Genuchten–Mualem models. The latter also predicted the profiles’ effective K function for the three soils, and the in situ based function was fitted with R2 ≥ 0.98 irrespective of soil type. It was concluded that the inverse parameter estimation with HYDRUS-1D improved models’ K function estimates for the studied layered soils.


Introduction
Soil profile physical and permeability properties determine the rate and extent of water flow especially in soils with contrasting horizons.At near saturation, water flow is more rapid because the majority of soil pores conduct water.The amount of water that drains away is depicted as deep drainage losses while the water content at which internal drainage allegedly ceases is referred as field capacity or drainage upper limit (DUL) (Hillel, 2004).These are very important components of the soil water balance function, which is applied in many agricultural and environmental issues.Reliable knowledge about water flow and storage at near saturation relies on accurate estimates of soil hydraulic properties.
Soil hydraulic properties serve as inputs in Richard's flow equation and include saturated hydraulic conductivity (K s ), unsaturated hydraulic conductivity (K), the soil water retention curve (SWRC) also depicted as water content (θ ), and the matric suction (h) relationship.The K can be expressed as function of θ or h.Direct measurements of these properties are difficult, given the equilibrium and maintenance of stringent initial and boundary conditions required over several stages during transient experiments (Dirksen, 1999).Despite these difficulties there are direct methods that are still commonly used for laboratory and in situ experiments.These include the hanging water column (Haines, 1930) and pressure cell plate methods (Richards, 1941) for the SWRC, and various forms of infiltration evaporation and outflow laboratory techniques (Gardener, 1956;Wind, 1968;Klute, 1986;Klute and Dirksen, 1986).For in situ conditions, possible methods are the plane of zero flux (Rose et al., 1965), constant flux vertical time domain reflectometry (Parkin et al., 1995) and the instantaneous profile method (Rose et al., 1965).The instantaneous profile method of Rose et al. (1965) is also considered the reference method (Vichuad and Dane, 2002).
Consequently, the instantaneous profile method has been modified and validated to improve the precision of measurements in situ.Soil water content (SWC) and h measurements are noisy because of increased spatial variability.Differentiating the measurements often produces mediocre K functions in soils with restrictive layers (Fluhler et al., 1976;Tseng and Jury, 1993).Several authors recommend subjective smoothing of the data prior to differentiation (Ahuja et al., 1980;Libardi et al., 1980;Luxmoore et al., 1981) or the adoption of the fixed or gravity gradient (Sisson et al., 1980(Sisson et al., , 1988;;Sisson, 1987).The fixed gradient depicts K in the expression dK/dθ to be equal to depth over time (z/t).The expression is modified by directly substituting analytical functions for K and reformulating data analysis in terms of the gravity-drainage optimisation code UNIGRA (Sisson and van Genuchten, 1992).This expression was extended to layered soils by scaling water content to assume equivalent homogeneous soil profiles (Shouse et al., 1991;Durner et al., 2008).By so doing, the K functions for the different layers are linearized into a single effective property, and the effect of spatial variability is minimized.Other forms of linearization include simple arithmetic, weighted or geometric statistical average schemes, as well as stochastic means (Wildenschild, 1996;Baker, 1998;Belfort and Lehman, 2005).
Alternatively, in situ K functions can be determined indirectly by applying transient experimental data to the inverse technique (Hopmans et al., 2002;Kosugi et al., 2002).Given the increased availability of computer models for solving the Richard equation and analytical expressions that describe K and θ-h functions by using a few parameters, hydraulic properties can be estimated simultaneously with a single transient experiment.The hydraulic parameters are usually based on the SWRC, because it is easily measured and can be estimated using a parameter optimisation technique (Kool et al., 1987;Hopmans and Simunek, 1999).Several studies applied inverse modelling and parameter estimation of hydraulic functions directly from in situ drainage transient experiments (Dane and Hruska, 1983;Romano, 1993;Zijilstra and Dane, 1996;Musters and Bouten, 1999;Dikinya, 2005).Agreement between in situ and predicted hydraulic functions was generally satisfactory, even though the K function was highly variable.The K function variability was attributed not only to spatial variability but also to the computational and convergence efficiency of the model in which many parameters are simultaneously optimised (Zijlstra and Dane, 1996).Consequently, lack of parameter uniqueness and lower boundary conditions limitations are common deficiencies in inverse modelling of soils with contrasting horizons (Romano, 1993).Narrow SWC range depicted by internal drainage experiments especially from poorly drained soils can result in ill-posed inverse solution and parameter estimation of soil hydraulic properties (Zijlstra and Dane, 1996;Simunek et al., 1998;Hopmans and Simunek, 1999).
The HYDRUS-1D software simulates water, heat and solute movement in one-dimensional variably saturated media and has the inverse method and parameter estimations, which can be applied to a wide range of in situ conditions (Simunek et al., 1998b;Simunek and Hopmans, 2009;Jiang et al., 2010).The code numerically solves the Richard equation by Galerkin-type linear finite element schemes and is equipped with the Marquardt-Levenberg type parameter optimisation technique (Simunek et al., 1998b).This is a local optimisation algorithm that requires initial estimate of the unknown parameters to be optimised.Local optimisers have been shown not to be powerful enough to handle topographic complexities of the objective function such as those emanating from lack of a well-defined global minimum or having several local minima in the parameter space (Simunek and Hopmans, 2002;Vrugt and Bouten, 2002).More computationally intensive and robust external global techniques have been developed and can be interlinked with HYDRUS (Vrugt et al., 2003;Wohling et al., 2008;Zhu et al., 2007).When using the Marquardt-Levenberg technique, it is recommended to identify a limited number of parameters and to solve the inverse problem repeatedly using different initial estimates of the optimised parameters, and then select those parameter values that minimised the objective function (Hopmans et al., 2002;Simunek et al., 2012).Nevertheless, robustness of HYDRUS-1D code to apply the inverse method and parameter estimation directly from transient data for soils with contrasting horizons has been well demonstrated (Sonnleitner et al., 2003;Sumunek et al., 2008;Montzka et al., 2011;Rubio and Poyatos, 2012).
Sub-standard estimates of layered soil hydraulic properties cultivated at Paradys Experimental Farm of the University of the Free State in South Africa have been a concern for the past decade.These are marginal soils with saprolite or high swelling clay content, either in the B or C horizons.Several studies have been confronted with challenges when installing, calibrating and reading tensiometer from these soils (Fraenkel, 2008;Chimungu, 2009;Bothma, 2009).In this regard a case study was undertaken to characterise hydraulic properties of undisturbed soils.The objective of this study therefore was to improve the prediction of in situ K function for individual soil horizons and equivalent homogeneous soil profiles.

Experimental site and classification
Three soil types cultivated at the Paradys Experimental Farm (29 • 13 25 S, 26 • 12 08 E; altitude 1417 m) of the University of the Free State located south of Bloemfontein, South Africa, were selected.These included Tukulu, Sepane and Swartland (Soil Classification Working Group, 1991).Tukulu and Sepane are also referred as Cutanic Luvisols and Swartland as Cutanic Cambisols (World Reference Base for Soil Resources, 1998).Three excavated soil profile pits at a depth of 1 m on each representative soil form (Fig. 1) were used for pedological and physical classification of diagnostic horizons.

Soil sampling
Two forms of soil samples were taken from the sites.Disturbed soil samples were taken from each excavated soil profile diagnostic horizon for textural and chemical analysis as proposed by the Non-Affiliated Soil Analysis Work Committee (1990).Undisturbed samples from the soil profile horizons were taken for bulk density and soil water retention curve measurement.A core sampler with an inner diameter of 105 mm and a height of 77 mm, mounted on a manually operated hydraulic jack, was used to take samples from the horizons at the end of the internal drainage experiment.For the determination of gravimetric soil water content, soil samples were oven dried at 105 • C for 24 h.

Saturated hydraulic conductivity
Saturated hydraulic conductivity (K s ) for the individual profile layers of the three soils was measured using double-ring infiltrometers as described by Scotter et al. (1982).Soil profile pits were excavated in a stepwise manner to allow the fitting of both rings with diameters of 400 and 600 mm at a depth of 20 mm.The falling head over a distance of 10 mm depth was used to determine K s with every fall recorded by means of a timer and a calibrated floater.After steady state was recorded for three consecutive times, the K s constant value (mm h −1 ) was computed using the Jury et al. (1991) formula given as where L is the depth of the soil layer in question (mm), b o the initial depth of total head above the soil column, b 1 the depth that the falling head is not allowed to fall below (mm), and t 1 the time taken for b o to fall to b 1 (in hours). 1

Soil profile drainage curves
Tukulu and Sepane monoliths with a 4 m × 4 m surface area and 1 m depth were prepared in triplicate.There were difficulties in excavating the Swartland as a result of the dolerite and saprolite rock, and the monoliths were reduced to 1.2 m × 1.2 m area and 0.5 m depth.To capture soil water content measurement within and just below the monolith, the neutron access tubes were installed at a depth of 1.1 m on the central area of each monolith in a V-shaped arrangement.Given the shallow depth of the Swartland, soil water sensors (DFM capacitance probes) were installed at a depth of 0.6 m with the water sensors positioned at 100, 300 and 550 mm.Polythene plastic was used to isolate side walls with slurry used to seal the sides from the surface.A ridge around the monoliths was also used to keep away surface runoff.In the absence of tensiometer measurements, monoliths were pre-ponded for three consecutive days to ensure wetting of the soil profile to near saturation.On the third day, each monolith surface was covered with a polythene plastic sheet to protect the trial from weather elements.Neutron access tubes and probes were inserted through openings in the plastic sheet and sealed with tape.Immediately after sealing, soil water content measurements were taken at the centre of each profile horizon and then daily for 50 days.Corresponding measurements for the A, B and C horizons for the Tukulu were at 150, 450 and 725 mm, Sepane at 150, 500 and 800 mm, respectively.

Laboratory-based experiments
The SWRC for the three soil profile horizons was determined with a laboratory desorption experiment.At the end of the drainage experiment, undisturbed soil samples were obtained from monoliths.These were first de-aired with a vacuum chamber pump set at −70 kPa for 48 h at room temperature.De-aired water was then introduced to saturate samples by capillarity for 24 h.Samples were then desorbed through the following series of pressure heads, 0 to −10 kPa, −10 to −100 kPa, and −100 to −1500 kPa.The first phase of desorption involved the hanging water column method, described by Dirksen (1999).At every step, interval samples were weighted before and after equilibration.The desorption chamber for −100 to −1500 kPa was designed to take samples of smaller volume, and thus the samples were disturbed and packed in 2000 mm 3 PVC tubes at the measured bulk densities.Reducing the sample volume also improved experimental time measurements and had little effect on the quality of the desorption data because, at high matric suction, range desorption mainly occurred in the soil matrix.Measured θ at h level was plotted to produce the SWRC.

Mathematical description of the soil water characteristic curve
The Brooks andCorey (1964), van Genutchen (1980) and Kosugi (1996) parametric models were used to describe the SWRC for the selected three soil profiles diagnostic horizons.These models describe the θ-h relationship from the expression representing pore size distribution (Kosugi, 1996) of many soils written as where S e is effective saturation, θ s and θ r are the respective saturated and residual values of the volumetric water content, θ (mm mm −1 ), h is the matric suction (mm), m equals 1, while α and n are the shape and pore size distribution parameters, respectively.The Brooks and Corey (1964) reduced Eq.( 2) into the following general equation: where α is the inverse of air-entry value, and the rest is as defined previously.This expression allows a zero slope to be imposed on SWRC as h equals air-entry value.S e equals unity when h ≥ −1/α.The van Genutchen (1980) and Kosugi (1996) model assumed the following respective expressions: (5) where, for the van Genutchen (1980) model, the condition m = 1 − 1/n should be satisfied with the air-entry value of −2 cm.For the Kosugi model, symbol α instead of h o and n instead of σ are adopted for uniformity reasons by some computer optimisation programs such as RECT (van Genuchten et al., 1991) and HYDRUS-1D (Simunek et al., 1998b(Simunek et al., , 2008)).
Model description of the experimental SWRC was carried out using the RECT program that constituted all three parametric models.Saturated soil water content (θ s ) was initially equated to total porosity (f ) as defined by the expression in Eq. 7 for pre-saturated undisturbed soil samples, and θ r was assumed to be equal to SWC for desorbed samples at −1500 kPa.
where ρ b is dry bulk and ρ s is particle density.The initial estimates of saturated and residual soil water contents were then optimised together with the α and n values determined using the Rosetta Lite pedotransfer software (Schaap et al., 2001).The K s initial estimate was determined from the double-ring experiment and therefore was unchecked for optimisation.

Estimation of internal drainage tensiometry data
The instantaneous profile method determines matric suction from tensiometers installed at various depths corresponding to soil water measurement points.This standard procedure was slightly modified following preliminary failure of tensiometry instruments to provide reliable calibration.Tensiometry data for the internal drainage experiment were then inferred from the θ-h relationship of the SWRC under the assumption that SWRCs from laboratory and in situ experiments are similar at near saturation irrespective of the structural effect and air entrapment (Bouma, 1982;Wessolek et al., 1994;Morgan et al., 2001).The van Genuchten (1980)

Calculation of unsaturated hydraulic conductivity in situ
The instantaneous profile method (Hillel et al., 1972;Marion et al., 1994) was used to determine the K(θ) function for the three soils' internal drainage boundary conditions.Changes in water storage between time intervals at different depths (z) corresponding to soil profile horizons were computed into drainage flux q(z, t) (mm h −1 ), which was then fitted to the following mass balance expression such that where K (θ) is unsaturated hydraulic conductivity (mm h −1 ), and dh is the change in the estimated matric suction (mm) between the neighbouring horizons, z (mm), which is the thickness of the horizon layer in question.The positive unity value on the hydraulic gradient component represents the effect of gravity with change with profile depth (dz/dz) and positive for downward flow.The calculated K functions together with the corresponding K s from the soil profile horizons were plotted on a semi-log scale.

Predicting unsaturated hydraulic conductivity
Firstly, unsaturated hydraulic conductivity was predicted from models based on the knowledge of SWRC and K s parameters.The conductivity functions of Brooks and Corey (1964), van Genuchten (1980) integrated with the Mualem (1976) expression and Kosugi (1996) were used to predict the K functions of the three diagnostic horizons respectively given as (11) where symbols are as previously defined.The RECT program was used to predict the K function simultaneously with the optimisation of the SWRC parameters.
Secondly, the inverse option of the HYDRUS-1D code was used to predict unsaturated hydraulic conductivity for individual horizons and for the average soil profile.Transient internal drainage SWC-time data were used in the objective function with soil hydraulic parameters optimised from the SWRC and in situ based K s entered as initial estimate in the inverse problem.Separate inverse solutions were run for the single porosity Brooks and Corey (1964), van Genuchten- Mualem (1980), and Kosugi (1996) models.
For the layered profile inverse solution, the graphical profile was discretized into three layers and observation points located at centre blocks corresponding to in situ profile horizons and SWC measurements.A constant flux and a free drainage were selected for the upper and lower boundary conditions, respectively.Initial conditions were set in water content measured at the onset of the drainage experiment.Given that the Marquardt-Levenberg-type parameter optimisation technique is only applicable to identify a limited number of unique parameters, no more than three parameters were optimised for each horizon.The θ r and θ s were among the first set of parameters to be checked alongside the l exponent parameter.Hydraulic parameters of soil profile horizons were optimised simultaneously during the application of the inverse solution.
The HYDRUS-1D was also used to estimate unsaturated soil hydraulic properties for equivalent homogeneous soil profiles of the Tukulu, Sepane and Swartland.Geometric mean average scheme as defined by Barker (1995Barker ( , 1998) ) was used to determine the representative effective profile drainage curves and pertinent hydraulic parameters.This average scheme provided estimates that were as accurate as the more sophisticated stochastic mean (Wildenschild, 1996;Abbasi et al., 2004).
Soil water contents measured during drainage for each horizon were averaged to give an effective profile drainage curve that was in turn used to compute effective water fluxes.Estimated matric suctions from horizons were not averaged, but effective matric suction gradient was calculated using the values of the surface and underlying horizons that bordered the flow domain.The effective flux and hydraulic gradient are then fitted in Eq. 9 to approximate the in situ effective K function.The K s values of individual horizons were also linearized using the same average scheme to estimate effective K s .The effective K function was presented in a semi-log scale with the effective K s being the first on the plot.
Effective SWC-time data were also used in the objective function during the optimisation process of inverse effective parameter estimation.Other effective parameters estimated by averaging were the K s and θ s .To improve model prediction, θ r for the most restricting layer was used in the initial estimates.The same was also done for the α and n parameters because their high non-linearity discouraged the use  of simple averages (Wildenschild, 1996).In the simulation of equivalent homogenous profile, the flow domain had one observation point in the central position.Upper and lower boundary conditions were similar to those applied for a layered profile.

Sensitivity analysis
Sensitivity analysis of the optimised parameters was also carried out to identify the parameters whose variation has a large effect on the model output.Sensitivity coefficients (SCs) for SWC were calculated using the influence method as described by Simunek and van Genuchten (1996): where SC(z, t, b j ) is the soil water content change at time t and depth z due to a variation of the parameter b j .In this study each parameter was varied by 10 % of its optimised value.Thereby b is the parameter vector, while e j is the j th  unit vector.This function depicts sensitivity coefficients that depict the behaviour of the objective function at a particular location in a parameter space.In this regard a high sensitivity means that the minimum is well defined, and that one can estimate the parameters with greater certainty once the global minimum is identified.
The sensitivity analysis of soil water content to parameters of the Brooks and Corey (1964), van Genuchten- Mualem (1980), and Kosugi (1996) models was carried out only for the Tukulu soil assuming a 1200day drainage experiment.Sensitivity analysis of the effective soil water content to equivalent homogeneous soil profile parameters was only performed on the three soil types using the van Genuchten-Mualem model.

Statistical analysis
Measured and optimised drainage, unsaturated hydraulic conductivity, as well as the pertinent hydraulic parameters constituted the major findings.The coefficient of determination (R 2 ), root mean square error (RMSE) and the index of agreement or D-index as proposed by Willmot et al. (1985) were the statistical tools used to quantify the quality of fit and variability between measured and fitted data.3 Results and discussions

Soil water retention curve
Figure 2 shows the experimental and fitted SWRC for the Tukulu, Sepane and Swartland diagnostic horizons, whose soil physical properties are summarised Table 1.The corresponding hydraulic parameters and statistical indicators are presented in Table 2. From these results, it is clear that the shape of the SWRC varied with the horizons' textural and structural properties and that the model's fit was satisfactory (R 2 ≥ 0.93) for the three soil forms.Variations in SWRC with soil physical properties illustrated the importance of texture and structure on soil water release and storage.The "S" shape of the SWRC for the sandy textured orthic and neocutanic horizons was well defined.In the clay-rich (≥ 35 %) prismacutanic and pedocutanic horizons, the SWRC diffused to almost a straight line.The orthic A horizon from three soil forms had the highest sand fraction (≥ 80 %), and average SWC for the SWRC ranged from 0.34 to 0.12 mm mm −1 .Despite small differences in θ s between sandy and clay textured horizons, θ r showed remarkable variability with change in clay content with a range of 0.19 to 0.26 mm mm −1 observed for a clay content range of 35 to 48 %.These findings are similar to those made from sandy and clayey horizons by various authors (Wilson et al., 1997;Wildenchild et al., 2001;Fraenkel, 2008;Chimungu, 2009).Sandy soils are well known for their large volume of macro-pores that drain readily at near saturation as a result of the small air-entry value that was approximated at −1 kPa.However, the narrow pore size distribution increased matric suction to as high as −100 kPa and steepened hydraulic gradients.The SWRC for the Tukulu C and Sepane B and C horizons was consistent with high surface area and micro-pore volume of clay soils, which induced slow water release due to strong ionic adsorption and capillarity at an air-entry value as high as −1.5 kPa.
Although the models' fit had a good coefficient of determination (R 2 ) ranging from 0.92 to 0.99, discrepancies between measured and fitted data were observed.Most models showed a poor fit at near saturation (0 to −1 kPa) and at very high matric suction (−100 to −1500 kPa).Accurate measurement of matric suction at near saturation is a common challenge in desorption experiments.At near saturation flow through macro-pores is difficult to control and is less sensitive to changes in matric suctions.Entrapped air also reduces θ s in the range of 0.85 to 0.9 of porosity (Kosugi et al., 2002).Dependence of bulk density on SWC for swelling and shrinking clays was the primary source of discrepancy, especially for θ-h relationships at higher matric suctions.This phenomenon explained the poor fit of θ r at −1500 kPa for all models.Nevertheless, the fitted curve was able to agree with measured data points and shape, with the most consistent fit provided by the van Genuchten (1980) model.This confirmed previous studies that found the van Genuchten model to fit the SWRC of a wide range of soils accurately.
The Brooks and Corey model had the poorest fit especially at near air-entry value.This model also produced a poor fit in fine textured soils and undisturbed core soil samples (Kosugi et al., 2002).That the model imposes a zero slope on the SWRC near the air-entry point could explain the poor fit.Additionally, the measurement of θ-h relationships at saturation above 85 % was impractical because of the general disconnection of the gas phase at this SWC range (Brooks and Corey, 1999).

Estimating matric suction from parameterised retention curves
Figure 3 shows the estimated matric suction where parameterised SWRC of van Genuchten (1980) was fitted directly to transient SWC measured during the internal drainage experiment.The estimated matric suction showed consistency with soil profile physical properties and water gradient depicted by the Tukulu, Sepane and Swartland horizons.Decreasing matric suction head with depth for the Tukulu and Sepane soil profiles supported the presence of the prismacutanic C horizons that restricted drainage to near-saturated conditions.Consequently, the Tukulu and Sepane profiles had a matric suction range not higher than −1000 mm (−10 kPa) for the 1200 h drainage experiment suggesting that the restrictive C horizon impaired the overall drainage of the soil profile.Similar observations were made by Freankel ( 2008) and Chimungu (2009) for the same soil types.Greater spatiality was observed in the Swartland soil profile with the highest matric suction of −1200 mm (−12 kPa) for the saprolite C horizon.
Even though the validity of the estimates cannot be detected, the results and procedures that were followed provided a reasonable account of the internal drainage process.Estimated matric suction ranged from 0 to −1200 mm (−12 kPa) and was within the 0 to −33 kPa range proposed by Ratliff et al. (1983) for a number soils that drain to field capacity.In variably structured soils, −10 kPa is often used as a hypothetical boundary for separating drained structural pores and water-filled micro-pores (Marshall, 1959;Kutilek, 2004).Various work from local and international drainage experiments recorded suctions around −10 kPa from variably structured soils (Hensley et al., 2000;Sonnleitner et al., 2003;Nhlabatsi, 2011;Adhanom et al., 2012).The use of undisturbed core samples three times larger than the area sensitive to tensiometer ceramic cup qualifies this procedure, even though estimates were made from parameters based on SWRC.The fit of the estimated matric suction was also supported by the narrow SWC near saturation, depicted in drainage experiments, and required no extrapolation outside the experimental data.In addition, at this wet range it was difficult to measure the θ-h relationship accurately, especially under in situ conditions.

Comparison of in situ and predicted K function
Estimated matric suction from SWRC fitted with van Genuchten (1980) model parameters was used to determine the matric suction gradients (dh/ z).These together with drainage fluxes were fitted to Eq. 9 to calculate in situ K function for Tukulu, Sepane and Swartland diagnostic horizons.The K function was also predicted from Brooks and Corey (1964), Kosugi (1996) and van Genuchten- Mualem (1980) 3. Optimised parameters from inverse modelling and the corresponding statistics are presented in Table 4.The results showed that the fit from inverse modelling produced a better fit compared to that of the SWRC parameters irrespective of model and soil type.For In situ K function, the curves were characterised by steep gradient over narrow SWC ranges, especially from soil profile horizons with high clay content (> 35 %).For a change in SWC of 0.02 to 0.03 mm mm −1 , the K (θ ) values declined from saturation by three and four orders of magnitude from the Tukulu and Sepane, respectively.For the Swartland, a change in SWC of 0.1 to 0.2 mm mm −1 initiated a decline in K (θ ) of about four orders of magnitude from saturation.The gentle slope of the K functions from the Swartland was consistent with the low clay content (< 22 %) and the presence of saprolite rock in the C horizon.Similar observations of clay soils were made by Freankel (2008) and Nhlabatsi (2011).Given the poor drainage properties of Tukulu and Sepane, it was proposed that a zero drainage flux be assigned at the bottom of these soil profiles for soil water balance studies (Hensley et al., 2000).For retention-based models, predictions based on SWRC overestimated the K function of soils horizons for the drainage SWC range.Overestimates were pronounced at lower SWC with three to four orders of magnitude observed from the Tukulu and Sepane horizons, respectively.The Tukulu C horizon was the exception, where the best fit among these models was observed from the van Genuchten-Mualem (R 2 = 0.94) and Kosugi (R 2 = 0.96) models that over-and underestimated K function by one order of magnitude.The Swartland profile was better fitted by all models (R 2 > 0.66) with the best fit produced by the van Genuchten-Mualem model for the B horizon (R 2 > 0.92).The Brooks and Corey (1964) model overestimated the K function irrespective of soil profile textural and structural formation.
Although the models fitted the experiment SWRC data very well, there was a strong disagreement between the in situ and predicted K function, especially at lower SWC.Similar observations were made in various studies (Dane and Hruska, 1983;Zavattaro and Grignani, 2001;Abbasi et al., 2003;Dikinya, 2005;Adhanom et al., 2012).Poor representation of field conditions by laboratory measurements is acknowledged to be the primary reason especially for layered soils (Sonnleitner et al., 2003).Discrete soil columns used in desorption experiments are devoid of layers, and optimised parameters will tend to agree more with homogeneous and well-drained soils compared to structured soils (van Genuchten, 1980;Knopman and Voss, 1987;Dikinya, 2005).This analogy is supported by the overall better fit observed in the Swartland profile horizons.The better fit from the Tukulu C horizon could be explained by the limited discrepancy between SWRC-based θ r (0.26 mm mm −1 ) and the lowest SWC or drainage upper limit (DUL) (0.31 mm mm −1 ) from the drainage experiment.Therefore the optimisation of SWRC parameters for field conditions is essential for better predictions.
For inverse modelling, optimisation of hydraulic parameters for in situ conditions was carried out by fitting transient drainage data into the HYDRUS-1D inverse solution for the Brooks and Corey, van Genuchten-Mualem and Kosugi models.Figure 8 shows sensitivity analysis for soil water content to the models parameters (θ r , θ s , α, n and K s ).The most sensitive parameter in the van Genuchten-Mualem model was K s and θ s in the Brooks and Corey and Kosugi models, irrespective of horizons suggesting that these parameters were of critical importance to the minimisation of the objective function.Thus the K s parameter was treated as known from the double-ring experiments.During the optimisation process, models were able to reproduce the drainage curves very well with a coefficient of determination of no less than 0.90 (Fig. 4).Brooks and Corey, and van Genuchten-Mualem models had an overall better convergence (R 2 = 0.98) compared to Kosugi (R 2 = 0.93) irrespective of soil type.The most fitted parameters for the Tukulu and Sepane were the θ r , θ s and α, and for the Swartland the α and n.The similarities between the Tukulu and Sepane could be attributed to the poorly drained prismacutanic C horizon shared by these soil profiles compared to the Swartland with an underlying saprolite layer.Constant flux and free drainage of the respective upper and lower boundary conditions were applied in the Tukulu and Swartland numerical solutions.The Sepane models converged readily when the constant water content was selected for lower boundary conditions, suggesting that this soil had the most restrictive properties -hence, the K s value of 1.9 mm h −1 for the underlying horizon.
Quality of fit between in situ and predicted K functions from inversely optimised parameters was improved by more than one order of magnitude irrespective of model and soil type.The van Genuchten-Mualem model better fitted the Tukulu (R 2 ≥ 0.97) while the Brooks and Corey model fitted the Sepane (R 2 ≥ 96).The Swartland K function was fairly predicted by all models although the van Genuchten-Mualem produced the best estimates for the A and B horizons.This tendency for model performance to be soil-specific was not unique to this study.Many studies have shown that various models are likely to fit experimental data and that the van Genuchten- Mualem (1980) model was among the most robust models (Russo, 1988;Chen et al., 1999;Mallants et al., 1996;Simunek et al., 2008).However, because of the exponential mathematical background of the van Genuchten-Mualem model, it often shows better fit on weakly structured soils.This sentiment is confirmed when the models produced better fit for all the soils' A horizons.The good fit that the Brooks and Corey model obtained from the Tukulu and Sepane structured horizons could be attributed to the high air-entry point associated with clay soils (van Genuchten, 1980;Brooks and Corey, 1999;Kosugi et al., 2002).It is therefore not surprising that optimised parameters for the same soil profile varied among horizons.Interestingly, the optimised α and n values of the different horizons were nearly the same for the Tukulu and Sepane, especially for the Brooks and Corey, and van Genuchten-Mualem models.

Unsaturated hydraulic conductivity for equivalent homogeneous soil profiles
Sensitivity coefficients for effective SWC on optimised parameters of van Genuchten-Mualem model for the Tukulu, Sepane and Swartland soils are presented in Fig. 9. Figure 10 shows the fitting of the geometric mean drainage curve with HYDRUS-1D in the objective function using the van Genuchten-Mualem model during the estimation of the overall profile hydraulic parameters.In all the soils the model was able to reproduce effective drainage curves very well (R 2 = 0.98).Calculated and predicted K functions are compared in Fig. 11 with corresponding optimised parameter and statistical indicators shown in Table 5.
Sensitivity coefficients show significant loops for almost all parameters particularly in the Tukulu and Sepane.Maximum sensitivity coefficients were associated with the Swartland (SC ≤ 10) with n being the most sensitive parameter.High sensitivity was distributed near the start, middle and end of the 1200-day internal drainage for the Swartland, Tukulu and Sepane, respectively.This observation can be an indication that the sensitivity analysis was able to provide useful information for all the soil types' parameter optimisation process.Results show a strong agreement (R 2 ≥ 0.98) between the estimated and predicted effective K function for the three soil types.During the linearization of profile hydraulic properties, the θ r was the most optimised parameter followed by the α and n parameters.This was expected given that the profile drainage curve assumed an effective SWC and θ-h relationship.Interestingly, the optimised θ r was almost equal to the lowest SWC of the effective drainage curve suggesting that model predictions can be improved if there were minimum discrepancy between initially estimated θ r and the lowest SWC of the experimental data.Similar observations were also made by van Genuchten (1980).Over and above the θ r , the α and n parameters were observed to be very sensitive to changes in experimental data that are used in the objective function (Sonnleitner et al., 2003;Saito et al., 2009).The Sepane-optimised θ r value of 0.30 mm mm −1 compared to the 0.27 mm mm −1 of the Swartland confirmed earlier observations that the former had the highest clay content.This result shows that effective parameters were consistent with profile physical properties and thus can be used with reasonable confidence to predict K function for an equivalent homogenous soil profile.Some researchers qualified the use of effective parameters on the basis not only of reducing the enormous data required but also of improving convergence of the inverse solution (Santini and Romano, 1992;Abbasi et al., 2003Abbasi et al., , 2004)).

Conclusions
This study estimated unsaturated hydraulic conductivity of Tukulu, Sepane and Swartland soils from parametric models using information from saturated hydraulic conductivity, laboratory soil water retention and in situ internal drainage curves.The Tukulu and Sepane shared a prismatic C horizon rich in clay content (≥ 45 %) compared to the Swartland horizons that had less than 22 % clay content.The in situ based unsaturated hydraulic conductivity was determined with the standard instantaneous profile method slightly modified to allow estimation of the θ-h relationship from parameterised soil water retention curve.The soil water retention curves were parameterised with the Brooks and Corey (1964), Kosugi (1996) and van Genuchten (1980) models using the RECT code.These models fitted the measured retention curves well with RMSE of less than 2 % and R 2 of no less than 0.98, and the most consistent was the van Genuchten model.
Direct predictions of K from retention parameters produced overestimates of more than three orders of magnitude, especially at lower soil water content.The only exception was the van Genuchten-Mualem model, which produced estimates around one order of magnitude for the Tukulu C and Swartland B horizons.This result confirmed that hydraulic parameters from laboratory-measured soil water retention curves were generally ill posed for predicting in situ K conditions.Estimation of soil horizons K functions was improved by one or more orders of magnitude with inverse parameter estimation applied directly to drainage transient soil water content measurement using HYDRUS-1D.The Brooks and Corey, and the van Genuchten-Mualem models produced the best K estimates (R 2 ≥ 0.90) irrespective of soil type and horizon material.Further improvement was observed when in situ K function was predicted from effective soil hydraulic properties with R 2 of no less than 0.98 in all soil types.Based on this result it can be concluded that the prediction of in situ K function can be remarkably improved by inverse parameter estimation for individual soil horizons and equivalent homogenous soil profiles.

Figure 5 Fig. 5 .
Figure 5 Models predictions of in-situ hydraulic conductivity based on the soil water retention curve (a) and inverse modelling (b) for the Tukulu diagnostic horizons.

Figure 6 Fig. 6 .
Figure 6 Models predictions of in-situ hydraulic conductivity based on the soil water retention curve (a) and inverse modelling (b) for the Sepane diagnostic horizons.

Figure 7 Fig. 7 .
Figure 7 Models predictions of in-situ hydraulic conductivity based on the soil water retention curve (a) and inverse modelling (b) for the Swartland diagnostic horizons.

Fig. 8 .
Fig. 8. Sensitivity coefficients (SCs) of soil water content (θ) to parameters of the van Genuchten-Mualem (i), Brooks and Corey (ii) and Kosugi (iii) models for the Tukulu A, B and C horizons.
models using SWRC and parameters and inverse modelling.Fitted drainage curves used in the objective function during hydraulic parameter optimisation with HYDRUS-1D are shown in Fig. 4. The resulting in situ and predicted K functions are plotted in Figs. 5 to 7. Statistical indicators from SWRC-based K function are summarised in Table

11Figure 11 Fig. 11 .
Figure 11 Comparison of in situ and fitted hydraulic conductivity of equivalent homogeneous soil profiles for the Tukulu (a), Sepane (b) and Swartland (c) soil forms.

Table 1 .
Summary of the physical characteristics of the three soil types.
K s = Saturated hydraulic Conductivity; ( ) optimised values considered in this paper.

Table 2 .
Fitting models' hydraulic parameters of the SWCC for the Tukulu, Sepane and Swartland soil.

Table 3 .
Statistical measure of fit for conductivity-based parametric models on in situ K coefficient for the Tukulu, Sepane and Swartland soil horizons.

Table 4 .
Optimised models parameters with statistical indicators for the prediction of in situ hydraulic conductivity functions of the Tukulu, Sepane and Swartland soil horizons using HYDRUS-1D.
Italic indicates checked parameters during the optimisation process.

Table 5 .
Optimised parameters of van Genuchten-Mualem model for equivalent homogenous soil profiles.