The two-layered (0 - 50 and 50 - 250 mm) surface horizon hydraulic parameters of three dryland floodplain soil-types under aquafer water management in Postmasburg, Northern Cape Province of South Africa were estimated with HYDRUS-1D model. Time dependent water infiltration measurements at 30 and 230 mm depths from simulated rainfalls on undisturbed 1 m2 small plots with intensities of 1.61 (high), 0.52 (medium) and 0.27 (low) mm·min-1, were minimised using a two-step inversion. Firstly, separate optimisation of the van Genuchten-Mualem model parameters for the two surface-horizon layers and secondly, simultaneous optimisation for the joint two-layered horizon with first step optimal parameters entered as initial values. The model reproduced transient water-infiltration data very well with the Nash-Sutcliffe model efficiency coefficient (NSE) of 0.99 and overestimated runoff (NSE; 0.27 to 0.98). The upper surface horizon had highly optimised and variable parameters especially θs and Ks. Optimal Ks values from higher soil surface bulk-density (≥1.69 g·cm-3) were lower by at least one order of magnitude to double ring infiltrometers and water infiltration properties were different (P < 0.05) for the high rainstorm due to raindrop impact and surface crusting. Optimal α and n parameter values corresponded well with texture of the Addo (Greysols), Augrabies (Ferralsols) and Brandvlei (Cambisols) soil types. However, θs and Ksshowed greater sensitivity to model output and exerted greater influence on dryland floodplain water-infiltration and runoff characteristics. Increasing rainfall simulation period to attain near-surface saturated conditions and inclusion of surface ponding data in the inverse problem could considerable improve model prediction of hydro-physical parameters controlling surface-subsurface water distribution in fluvial environments.
Soil hydraulic properties controlling infiltration and runoff play an important part in capturing and distributing water resources in dry riverbed and floodplains. These fluvial environments are strategic sites for groundwater recharge and water-resource development. Modelling surface and subsurface water-flow requires knowledge of soilhydraulic parameters. However, sedimentation and “amphibious” conditions characterising fluvial depositional environments make soil surface hydraulic properties to be highly variable   .
Soil hydraulic properties, which describe water-flow in variably saturated media, include saturated hydraulic conductivity (Ks), unsaturated hydraulic conductivity (K), and the soil water retention curve (SWRC). The SWRC represents the relationship between water content (θ) and metric suction (h), and mathematically represented by various pore-size distribution models (     ). These analytic models use parameter based closed-form equations to describe θ-h relationship and are used to predict other more difficult properties, K(h) or K(θ) in particular. Application of analytic models has been actively studied for several decades, but relatively few studies exist in which hydraulic parameters of fluvial soil types have been estimated using data from in situ experiments.
Soil types developed from colluviation and sedimentation vary from highly permeable to impermeable  . Fluvial processes, geometry and topography influence the nature and type of alluvium deposits in floodplains  . Gravel and sands alluvium affect floodplain-segments associated with fast-running water while alluvium of silt and clays affects segments with slow-moving water. Higher permeability of course textured alluvium deposits is linked to rapid and shallow groundwater recharge  . Low relief terrain of most arid floodplains reduces flow velocities that favour water infiltration and widespread recharge. However, settling and accumulation of translocated fine particles from other positions by runoff and erosion can remarkably reduce permeability of the soil surface even of course-textured soil types  . Hot and dry conditions of arid climates make floodplain inundation events rare and far apart resulting to bare soil and sparse vegetation. Under these conditions, compaction and sedimentary crust formation characterise soil surface conditions  .
A soil crust often consists of two parts. The first is an upper skin seal, 0.1 mm thick and forms under the influence of raindrop impact, splash, slaking, swelling and sedimentation. The second is a 2-mm thick deeper region of washed-in dispersed fine particles  . Permeability of the soil beneath the crust is 200 times higher than the washed in layer and about 2000 times higher than the skin seal  . A 0.1 mm thick surface crust may reduce infiltration rate (IR) from 800 cm day to 70 cm day   . Subsequently, the matric flux potential and hydraulic gradient could decrease K(h) by an order of several magnitudes. Direct measurement of a soil crust hydraulic properties is not practical due to its dual thin layers and stringent procedures requiring highly specialised equipment. A standard procedure uses data from tension disk infiltrometers and pre-installed mini-tensiometers beneath a soil crust to estimate K(θ), pressure head (h) and sorptivity.
Alternatively, due to the delicate nature of a crust, steady-state infiltration measurements is commonly used to indirectly estimate hydraulic properties of a well-established soil crust with a constant thickness. Touma et al.  combined a rainfall simulation experiment and a single-ring infiltrometer test on the same site before and after removing the crusted layer to calculate the crust steady state IR. Alagna et al.,  adopted the same approach, but substituted the minidisk tension infiltrometer with a rain simulator and used pedo-transfer functions to estimate steady state IRof the soil beneath the crusted layer. These studies relied on saturated crust hydraulic properties and limited findings to steady state IR or Ks. In this study, saturated and variably saturated flow parameters are estimated using a numerical model with an inverse parameter optimisation algorithm  and is coupled with transient soil-surface water infiltration measurements from small-plot rainfall simulation experiments. Usefulness and advantages of the inverse method and hydraulic parameter estimation from transient soil-water measurements has been well demonstrated over several decades     .
The aim of this study was to determine the near-surface soil hydro-physical properties of the dominant soil types found at a dryland floodplain used to monitor groundwater recharge in the Northern Cape of South Africa. Referred as a drainage or dryland floodplain, the area has several boreholes (Figure 1(c)) drilled into shallow aquifers that are artificially recharged by pumping water from open cast mine pits of a nearby iron ore mine. Since the launching of the aquifer recharge program in 2012, data from groundwater measurements showed remarkable improvements. However, lack of soil surface hydro-physical data made it difficult to account for the contribution of direct water infiltration from rainstorm events to the subsurface water balance.Physical models such as the Richards  equation and Green Ampt,  model require accurate estimation of soil hydraulic properties to compute water-flow in variable saturated soils. The study only sought to near surface hydraulic parameters that are better predictors of water infiltration and runoff from the studied floodplain segment. Parameter sensitivity was assumed to depict a measure of influence on soil-water infiltration and runoff characteristics of the studied floodplain soil types. The specific objectives were therefore to: 1) estimate near-surface soil hydraulic parameters of three floodplain soil-types from transient water infiltration data using a two-step inverse parameter optimisation approach  , and 2) evaluate sensitivities between parameters of the van Genuchten  analytic model for rainwater infiltration and runoff simulation results.
Figure 1. Geographical location (o) of the Postmasburg town in the Northern Cape Province, South Africa (a), Aerial photograph of the floodplain with location soil types (b), and boreholes infrastructure of the aquifer recharge program (c).
2. Material and Methods
2.1. Site Description
The study area was located at the dryland floodplains under wild life and aquafer water management of the Anglo American Kolomela Iron Ore mine, situated 30 km south of the Postmasburg town, Northern Cape Province of South Africa (Figure 1(a)). The area had a common slope of 1% with three dominant soil types (Figure 1(b)). These included the Addo, Augrabies and Brandvlei soil forms  , which also referred respectively as Greysols, Ferralsols and Cambisols by the International Union of Soil Science (IUSS) Working Group World Resource Base (WRB)  . The Food Agriculture Organisation (FAO)WRB,  for soil resources generally referred to these floodplain soils types as fluvisols. A summary of soils physical and chemical properties is shown in Table 1. The Addo and Brandvlei soil forms found at altitude of 1252 m while the Augrabies at altitude of 2060 m.
A long-term rainfall data constituting monthly minimum, maximum and averages over a 98- year period of the Postmasburg area was summarised in Table 2. The wettest months are January to March and the driest months are June to August. Rainfalls with amounts of 130 mm and 104.5 mm were the highest recorded rainstorms and occurred in the months of February and January, respectively. Needle like vegetation and sedges were common in the Addo while clusters of tussock grasses and bushy plants were observed in the Brandvlei and Augrabies soil forms. Vegetation was withered and dormant with the ground virtually bare.
2.2. Experimental Design and Treatments
2.2.1. Rainfall Simulation Experiment
A mobile field rainfall simulator was used to generate soil water infiltration and runoff data. The simulator constituted of oscillating sprinkler nozzlewith adjustable height in a closed compartment to protect operations on windy days (Figure 2(a)). It is also fitted with water pump; pressure gauges and intensity regulators. At the floor of the compartment is a metal runoff frame of 1 m × 1 m area that when inserted on the ground at 10 cm soil depth formulated the experimental plots (Figure 2(b)). Plots were prepared close to representative soil profiles from each of the three soil types. A gutter fitted on the sloping side of the frame connects with an outlet pipe for runoff collection. The simulator created three rainstorms of 0.27 mm・min−1 (low), 0.52 mm・min−1 (medium) and 1.61 (high) mm・min−1 intensities. When calibrating the rainfall simulator, many discharge irregularities and inconsistencies confounded intensities higher than 1.62 mm・min−1 or 80 mm・hr−1 and lower than 0.27 mm・min−1 or 16 mm・hr−1. Selected rainfall intensities intheir increasing order had simulated times of 56, 50 and 40 minutes to obtain corresponding accumulative amounts of 15, 26 and 65 mm (Table 3). To ensure that the correct rainfall amount for a particular intensity a time-based calibrated automated rain gauge was placed inside the simulation plot (Figure 2(b)).
Table 1. Summary of soil types physical and chemical properties.
Table 2. Long term (98 years) monthly rainfall (mm) for the Postmasburg area.
Table 3. Characteristics of simulated rainstorm treatments.
Figure 2. A Hofrey rainfall simulator (a) mounted on a small plot with a 1 m × 1 m metal frame connected to a gutter with an outlet runoff collecting pipe (b).
2.2.2. Soil Water Infiltration and Runoff Measurement
Time dependent soil water content and runoff measurement obtained from simulated three rainstorm treatments replicated three times from the three soil types. Installation of the mobile simulator metal frame and driven to the ground up to 10 cm depth on undisturbed land representative of soil surface conditions in situ. Centrally of each plot a 1 m long capacitance soil water measuring (DFM) probe was vertically installed such that water-measuring sensors are at 3 cm, 23 cm, 43 cm, 63 cm and 83 cm depths (Figure 2(b)). Programmed capacitance continuous water meter probes measured soil water content at minute intervals during rainstorm simulation. Soil water content measurement taken prior and post rainfall simulation used as initial and final soil water contents, respectively. Gravimetric soil water content and bulk density values employed to calibrate capacitance water sensors. Runoff measurement involved time to runoff and runoff amounts for every after 5 min intervals until simulation cut-off time. Manual measurement of runoff yield obtained from the gutter with a measuring cylinder.
2.3. Theoretical and Experimental Considerations
2.3.1. Numerical Flow Model
The Richard equation  is a numerical model that describes flow in variably saturated soils and for one-dimensional vertical infiltration written as follows:
where ∂θ/∂t is the surface water flux, z is the down wide direction; H is the soil-water pressure head relative to atmospheric pressure (H = h + z), h being matric suction potential and S is the sink. HYDRUS-1D software  modifies the one dimensional vertical flow equation by using cos α instead −1 representing downward gravity gradient. The α is the angle between the flow direction and the vertical axis with α = 0˚ for vertical flow and 90˚ for horizontal flow. HYDRUS-1D code numerically solves the water flow equation using the Galerkian-type linear finite element scheme.
2.3.2. Soil Hydraulic Functions
The single porosity model of van Genuchten (1980)  described soil hydraulic functions in terms of soil-water retention θ(h), parameters. It also predict unsaturated hydraulic conductivity function from the statistical pore-size distribution model  shown in Equations (2) and (3).
where θr and θs are residual and saturated water contents (mm・mm−1), respectively; α is the air entry value also referred as bubble pressure [mm], n is the pore size distribution parameter [-], Ks is the saturated hydraulic conductivity, and l is a pore-connectivity parameter assumed to 0.5 [-]. The condition m = 1 − 1/n should be satisfied with the air entry value of −2 cm.
Initial estimate of the van Genuchten  model water retention and hydraulic conductivity parameters of Maulem  were predicted using the Rosetta pedo-transfer model of Schaap et al.  constituted in the HYDRUS-1D software  . Rosetta estimate retentive curve with good statistical comparability to known retention curves of other media with similar properties to the medium’s particle-size distribution and other soil properties  .
2.4. Inverse Modelling
The HYDRUS-1D model have an inverse modelling capability and was required to optimise unknown parameters for van Genuchten  model characterising rainstorm infiltration in the water flow equation (Equation (1)). The objective function constituted residuals between observed and predicted soil water contents at different depths and times and was minimised using the Levenberg-Marquardt nonlinear minimization method (   ). Mathematically, the objective function may be defined as
where ɸ(β) is the objective function of the parameter vector, β; θi* and θi are measured and predicted soil water contents, respectively; z is the depth; t is the time and N is the number of observations available; W is the weight of specific measurement by standard deviation.
The surface horizon constituted the flow domain discretised into upper (0 - 50) and lower (50 - 250 mm) layers with observation points at 30 mm and 230 mm depths, respectively. HYDRUS-1D model predicts well infiltration properties near the surface but because of neglecting the effect of entrapped air on water infiltration, it often underestimates the infiltration depth  . The two-layered horizon was discretised into 101 nodes with specified nodal density of 0.01 and 1 for the upper and lower layers. The high density of the near surface superficial layer was aimed to capture rapid changes in water contents by the advancing wetting front. Infiltration was initiated by atmospheric boundary condition with surface runoff and was maintained for period of the simulated rainstorm. Rainstorms intensity and zero values of evapotranspiration were specified. Free drainage was prescribed for the lower boundary condition and initial conditions given in terms of water contents measurements. Time discretization was in minutes with 20 maximum number of iterations and initial time step of 0.01.
To improve parameter identification and uniqueness for upper and lower layers the surface horizon, a two-step inverse parameter optimisation  was conducted. The first step estimated soil hydraulic optimal parameter values of the upper and lower surface horizon layers as separate independent layers. The second step estimated hydraulic parameters of a two-layered surface horizon simultaneously with optimal values from the first step as initial estimates. The Rosetta pedotransfer program contained in HYDRUS-1D estimated initial parameter values using soil type sand, silt and clay fractions and bulk density (Table 1). Time dependent soil-water contents measurements at 30 mm and 230 mm were used in the objective function for the first and second step inversion. HYDRUS calculates infiltration-excess runoff in the absence of ponding by the system-independent atmospheric boundary conditions. The surface water-infiltration flux was limited by the following two conditions  :
where, q maximum potential infiltration flux under the current atmospheric conditions, h is the matric suction potential at soil surface, z is depth, K(h) is unsaturated hydraulic conductivity, and hA and hS are, respectively, minimum and maximum matric suction potential allowed under the prevailing soil conditions.
2.5. Sensitivity Analysis
The objective function minimised in inverse parameter optimisation should provide sufficient information about the unknown parameters to be identified  . This is only possible if the objective function is sensitive enough to track changes in the unknown parameters. The ability of the objective function to capture these variations presumably within the vicinity of true parameter values can be evaluated using sensitivity analysis. Various methods exist including the coefficient method or finite difference, a sensitivity equation method and a variation method  . Sensitivities between different parameters for selected model outputs was compared using the unit less coefficient method as normalised by Simunek and van Genucthten  with 1% change effected on parameters written as:
where eij, is the change in the auxiliary variable ai corresponding to 1% change in parameter βj. Thereby β is the parameter vector, while ej is the jth unit vector. The parameter vector included the optimised θr, θs, α, n and Ks van Genuchten-Mualem parameters. Sensitivities conducted involved soil water content, infiltration rates, accumulative infiltration, runoff rate and accumulative runoff. A high sensitivity suggested a well-defined minimum and the parameter can be optimised with greater certainty once the global minimum is verified.
2.6. Statistical Analysis
The HYDRUS-1D model performance assessed with the Nash-Sutcliffe model efficiency coefficient (NSE)  and root mean square error (RMSE). The NSEwas calculated by the expression:
where θ1 is the predicted soil water content; θ*1 is the observed soil water content; θave is the average soil water content of all the observed events, and N is the number of observations i.e., the number of measured events.
The RMSE is widely used to measure agreement between the observed data and model prediction and represented by the expression:
The Duncan’s multiple range test (DMRT) was used to compare infiltration rates and accumulation infiltration of the three soil types and simulated rainstorms. The means were ranked from the highest to lowest values and ranks were compared to shortest significant ranges (Rp). The Rp was computed as follows:
where rp, tabular values of the significant studentized ranges at 0.05, standard error of the mean difference and s is the variance  .
3. Results and Discussion
3.1. Soil Water Content
Time dependent soil water contents at 30 mm and 230 mm depth measured from simulated (high, medium and low) rainstorms were minimised in the objective function and results for the model fit at 30 mm depth were presented in Figure 3. The model was able to reproduce the water contents for the three soil types very well with NSE not falling below 0.99. The lower layer had constant water content values suggesting the infiltration-wetting front was not detected or unable to reach this depth (results not shown). In all soil types, soil water content decreased with rainstorm size with the Addo and Brandvlei having earliest and latest response times regardless of rainstorm size. The good fit by the HYDRUS-1D code was not surprising because the model’s ability to reproduce measured infiltration data either on bare or vegetated surfaces is well documented (     ). Although a good fit is desirable for the minimization process, it is not sufficient for parameter identification because data used in the objective function directly or indirectly should relate to the unknown parameters (   ). Soil water content, pressure head and infiltration rate as well as accumulative infiltration data are acceptable for inverse parameter estimation (   ).
Figure 3. Observed (obs) and predicted (pred) soil water content at 3 cm depth from the Addo (a), Augrabies (b) and Brandvlei (c) soil types. Level of agreement represented with Nash-Sutcliffe model efficiency (NSE = 0.99) for the high, medium and low simulated rainstorms.
3.2. Soil-Surface Water Infiltration
Figure 4 and Figure 5 showed respective water infiltration rates and accumulative infiltration of studied soils and rainstorms estimated by HYDRUS-1D inverse parameter optimisation procedure. Final infiltration rates and total infiltration rates were summarised alongside runoff parameters in Table 5. Final infiltration rates and accumulative infiltration decreased with rainstorm intensity and size in all soil types. This finding was consistent with the observation that higher rainstorm amounts and intensity resulted to higher final infiltration rate (    ).
Initial infiltration rates declined sharply under the high rainstorm suggesting that the 1.62 mm・min−1 intensity was higher than soils infiltration capacity. The Addo and Augrabies approached steady infiltration rates after 20 minutes recording final infiltration rates and total infiltration of 0.71 mm・min−1 and 34.8 mm, and 0.91, mm・min−1 and 36.5 mm respectively. The lower infiltration values of the Addo was explained by the total fine silt plus clay content of 47.7% of the surface horizon compared to the 20.6% from the Augrabies. The Brandvlei had a constant initial infiltration rate for the first 5 min before declining to a final infiltration rate of 1.14 mm・min−1 and total infiltration of 50.51 mm. The Brandvlei higher infiltration capacity was consistent to the lower total fine silt plus clay content of 17.7%. Silt plus clay fraction were easily mobilised by high-energy raindrops into formation of surface crusts and seals that can reduce infiltration capacity by several orders of magnitudes (    ). Although the high rainstorm had significantly (P < 0.05) higher final infiltration rates, these values were lower when compared to Ks from double ring infiltrometers (Table 4). Despite double rings overestimations, the lower than expected steady infiltration rates confirmed that rainfall amount and intensity influenced steady state infiltrability especially when intensity excess steady infiltration rate  .
The medium rainstorm had constant initial infiltration rates in the first 10 min from the Addo and Augrabies, and 15 min from the Brandvlei, before declining to significantly similar (P < 0.05) final infiltration rates of 0.38, 0.14 and 0.38 mm・min−1, respectively. Corresponding accumulative infiltration were 22, 15.2 and 22.5 mm, respectively. The Augrabies produced lower infiltration rates under the medium rainstorm, which were also indifferent when compared to the low rainstorm. This result was attributed to the higher surface bulk density (1.71 g・cm−3) that supported low porosity and permeability. Compacted or crusted surfaces of 0.1 mm thick were observed to reduce final infiltration rates by more than 10 times  . Soils with surface crust of 1 mm thickness had final infiltration rate of 0.34 mm・mm−1 and crust thicker than 2 mm did not exceed 0.16 mm・min−1  .
Figure 4. Soil surface infiltration rates of the Addo (a), Augrabies (b) and Brandvlei (c) soil types for high, medium and low simulated rainstorms. Statistical mean differences for the Duncan Multiple Range Test (P < 0.05) represented by different letter symbols.
Figure 5. Accumulation infiltration of the Addo (a), Augrabies (b) and Brandvlei (c) soil types for the high, medium and low simulated rainstorms. Statistical mean differences for the Duncan Multiple Range Test (P < 0.05) represented by different letter symbols.
Table 4. Soils saturated hydraulic conductivity (Ks) derived from double ring infiltrometers.
Note: superscript letters, statistical different means depicted by different letters based on least significant Difference (LSD) mean separation test P ≤ 0.05; SD, standard deviation; CV coefficient of variation.
Table 5. Selected infiltration and runoff properties as predicted by inverse solution.
Note: superscript letters, statistical different means depicted by different letters based on the Duncan Multiple Range Test (DMRT) mean separation test at 95% confident intervals (P < 0.05). Superscripts letters also applicable for accumulative runoff (not shown).
Under the low rainstorm, initial infiltration rates were constant in the first 20, 25 and 35 min from the Addo, Augrabies and Brandvlei, respectively. Correspondingly, final or steady infiltration rates and total infiltration were 0.21, 0.15 and 0.18 mm・min−1 and 13.4, 12.7 and 14 mm, respectively. The longer constant initial infiltration rates showed that the low rainstorm intensity of 0.27 mm・min−1 was close to soil-types infiltration capacity, hence steady infiltration rates were significantly similar (P < 0.05). However, this result was not surprising because as rainstorm decreases the wetting rate and infiltration curve decreases slowly resulting to lower final infiltration rates (   ). Interestingly, the Addo despite a finer texture had a higher final infiltration rate: a phenomenon that can be attributed to mineralogy and distribution of exchangeable cations. The surface horizon had exchangeable cations of calcium (7650 mg・kg−1) and magnesium (1340 mg・kg−1) recognised as good flocculators versus exchangeable sodium (60 mg・kg−1) with high dispersivity. Such a distribution favoured a stable than dispersive surface aggregates, which supported infiltration especially under less disruptive rainstorm intensity. Stern et al.  observed that final infiltration rates greater than 0.13 mm・min−1 were associated with stable aggregates while less than 0.08 mm・min−1 characterised dispersive soils. In this study, final infiltration rate was greater than 0.13 mm・min−1 from all soil types suggesting that distribution of exchangeable cation did not favour dispersive infiltration reduction. However, in less dispersive soils, surface crusting and rainfall characteristics controlled final infiltration rates  . Prevalence of finer aeolian sand fractions under arid drylands climates favoured surface crust formations in South Africa soils (   ).
3.3. Water Infiltration-Unsaturated Hydraulic Conductivity
Figure 6 shows predicted van Genuchten-Mualem model K(θ) of the Addo, Augrabies and Brandvlei upper surface horizons. Corresponding Ks values determined with double ring infiltrometersis presented in Table 4. Saturated hydraulic conductivity between soil types were significantly different with the least (0.05 mm・min−1) and highest (2.02 mm・min−1) from the Addo and Brandvlei, respectively. Unsaturated hydraulic conductivity of soil types showed greater variation among rainstorms at near saturation. Corresponding K(θ) values at near saturation of the Addo, Augrabies and Brandvleifor high, medium and low rainstorms are 0.66, 0.23 and 0.11 mm・min−1, 0.09, 0.03 and 0.02 mm・min−1, and 0.53, 0.16 and 0.07 mm・min−1, respectively. The result showed that K(θ) function at near saturation increased with rainstorm amount and intensity suggesting volume of hydraulic active pores increases with rainstorm size until rainstorm intensity excess steady infiltration rate (    ). The K(θ) values at near saturation of the Augrabies and Brandvlei under high rainstorm were respective lower by 1.41 and 1.49 mm・min−1 when compared to corresponding Ks values. Such a reduction in near saturation conductivity is explained by the breaking down and slaking of soil surface aggregates by raindrop energy impact, which can reduce Ks exponentially  . Nciizah and Wakindiki  and Fohree, et al.  also observed that physical properties of surface crusts affecting K(θ) during and between rainstorms were always changing. This phenomenon explained the inter-changes of K(θ) functions between high and medium rainstorms from the Addoand, between medium and low rainstorms from the Augrabies and Brandvlei soil forms. Apart from these interchanges, K(θ) under all rainstorms merged at lower water contents suggesting rainstorms affected soil infiltration properties at near saturation water contents.
3.4. Soil Surface Runoff
In addition to soil-types infiltration, Table 5 presented observed and predicted time to runoff, final runoff rate and accumulative runoff. Even though runoff data was not part of the objective function, HYDRUS-1D model generated infiltration-excess runoff using Neumans type system-dependent conditions (Equations (5) and (6)). Figures 7-9 presented a 1:1 line comparing observed and predicted runoff rates and accumulative runoff for the Addo, Augrabies and Brandvlei, respectively. The model predicted poorly (NSE ≤ 0.27) rainstorm runoff rates for all soil types. However, prediction of accumulative runoff was satisfactory with NSE ranging from 0.53 to 0.98. Model prediction was highest and lowest for medium and lower rainstorms, respectively. In all cases, the model overestimated runoff rates and total runoff production, a phenomenon attributed to the exclusion of ponding conditions.
Due to shorter simulation periods (≤56 min), ponding was assumed to be negligible at small plots scale. Presence of dry and dormant vegetation on plots that can influence surface water storage and delay time to runoff  were not specified as part of the upper boundary conditions, hence the model runoff overestimation. Nevertheless, spatiality of vegetal cover and surface roughness in-situ can make ponding parameter(s) difficult to determine  . Agreement between observed and predicted accumulative runoff was sufficient though to describe soil-type surface runoff properties using optimised parameters. However, this should be done with caution only applicable for the prevailing soil surface condition.
(a) (b) (c)
Figure 6. Unsaturated hydraulic conductivity estimated from inverse solution of the Addo (A), Augrabies (B) and Brandvlei (B) soil types for High, Medium and Low rainstorm treatments.
Figure 7. Measured and predicted runoff rate (a) and accumulative runoff (b) of the Addo for High (H), Medium (M) and Low (L) rainstorms.
Figure 8. Measured and predicted runoff rate (a) and accumulative runoff (b) of the Augrabies for High (H), Medium (M) and Low (L) rainstorms.
Figure 9. Measured and predicted runoff rate (a) and accumulative runoff (b) of the Brandvlei for High (H), Medium (M) and Low (L) rainstorms.
The Brandvlei had longest time to runoff (4.3 to 38 min) for all rainstorms while the Addo and Augrabies had shortest times for the high and lower rainstorms, respectively. Total runoff increased with opportunity time making the Addo and Augrabies to have highest runoff for the respective high (27.5 mm) and lower (2.41 mm) rainstorms. This result collaborated with the higher total fine silt plus clay (47.7%) fraction of the Addo and bulk density (1.76 g・cm−3) of the Augrabies; observed earlier in this study as factors that reduced infiltration. Soils with high quartz fraction like in the Augrabies have poor bonding properties to protect aggregates against raindrop impact and hence, are highly dispersive and susceptible to surface crusting  . Total clay plus fine silt was lowest (17.7%) from the Brandvlei suggesting clogging of hydraulic pores by fine sediments was too low to affect infiltration a sentiment that was shared by Lado et al.  . In addition, excluding course and medium fractions, the Brandvlei had total sand fraction of 76% suggesting it had more uniform particle size distribution a property that discouraged translocation of dispersed fine articles and deposition into voids. Such a property is common among the wind-blown sands of arid and semi-arid regions of South Africa and in the absence of impervious underlying horizon, these soils supported deep infiltration and drainage (   ).
3.5. Soil Surface Hydraulic Parameters
Table 6 presented initial and optimised parameters of the Genucthen-Mualem model for the respective Addo, Augrabies and Brandvlei soil forms at 30 mm and 230mm depths. Inverse parameter optimisation was able to minimise objective function from the Addo, Augrabies and Brandvlei with RMSE not higher than 0.01. Initial and optimised parameters were considerable different among soil types and rainstorms with the exception of θr in the Addo and Brandvlei. The Addo and Bransdvlei were initially described as the least (Ks = 0.05 min−1) and most (Ks = 0.23 mm・min−1) permeable, respectively. However, optimised Ks values described the Addo as the most permeable (Ks; 0.11 to 0.78 mm・min−1) and the Augrabies as least permeable (Ks; 0.2 to 0.52 mm・min−1). This result showed that parameters derived from theoretical pedo-transfer functions were ill posed for describing soil-water dynamic processes in-situ. Soil heterogeneity and occurrence of superficial crusted surface layers and variability in atmospheric-surface boundary conditions are among the common reasons (    ). Nevertheless, parameters estimated from laboratory soil-water retention curves and pedo-transfer functions provided reliable initial parameter estimates foroptimisation of unknown parameters (   ).
Table 6. Initial soil hydraulic parameters obtained from Rosetta-code and final optimised parameters obtained by HYDRUS-1D inverse solution for the upper and lower surface horizon under high, medium and low rainstorm treatments.
Optimised parameters varied considerably for upper and lower surface horizons. The Ks fitted for upper layer was always lower when compared to the lower surface horizon. This result showed that the soil surface was susceptible to compaction and raindrop impact especially under higher rainstorm intensities  . Lower Ks values of the upper surface horizon corresponded with lower θs for the Augrabies and Brandvlei soil forms under the high rainstorm in particular. Higher surface bulk density of not less than 1.67 g・cm−3 and disruptive power of high intensity raindrop were the reasons. Despite the initially higher α parameters, which represented the coarseness of the Augrabies (0.004) and Brandvlei (0.003), the optimised α values for the former was remarkably lower (0.002). The lowest α parameter depicted the Augrabies as a fine textured soil because of higher bulk density of 1.76 g・cm−3, which reduced permeability and porosity for the upper surface horizon in particular. Larger n parameter represented greater pore-size distribution uniformity. The Addo had lowest n parameter (1.674 to 1.890) and were always lower in the upper compared to lower surface horizon. This result corresponded to the lower bulk density (1.52 g・cm−3) and higher total fine silt and clay (47.7%). Bonding and cementing properties of clays improved aggregate stability and inter-aggregate porosity  ; a phenomenon that supported higher θs parameter for the Addo’s upper horizon especially under the higher intensity rainstorms. Aggregate stability against dispersion in the Addo was also supported by large exchangeable magnesium and calcium compared to exchangeable sodium  .
Larger n parameter (2.801 to 3.857) characterised the Brandvlei followed by the Augrabies (2.680 to 2.919). Greater uniformity depicted by larger n values for the Brandvlei collaborated with the 76% total sand fraction excluding course and medium sand fraction. In addition to larger n parameter, larger α parameter (0.007 to 0.14) described the Brandvlei surface horizon suggesting it was a uniformly course textured soil. This description supported the higher infiltration and lower runoff of the Brandvlei from all rainstorm treatment.
3.6. Parameter Sensitivity
Sensitivities of soil water content (θ), infiltration rate (I), accumulative infiltration (Z), runoff rate (Ror) and accumulative runoff (Acc. Ro) to 1% change in optimised parameters for the medium rainstorm were presented in Figures 10-12 for the Addo, Augrabies and Brandvlei soil forms, respectively. Similar results were observed under the high and low rainstorms (data not shown).Sensitivities of model output to 1% change in parameters were noticeable at two infiltration stages; firstly, passing of wetting front at 30 mm depth resulting to a sharp rise in θ and secondly, when θ approached near saturation values. Table 7 also showed model output sensitivity to 1% change in hydraulic parameters as inferred from the RMSE. Deviation of RMSE from the optimised value due to 1% change in parameters was limited to θs and Ks parameters had greater influence on model output. Considerable deviation was from the Addo and Brandvlei with the latter higher RMSE suggesting model output showed greater parameter sensitivity.
Table 7. Effect of 1% change in parameters on models output for the medium rainstorm as inferred from the Root Mean Square Error (RMSE).
Figure 10. Addo soil type sensitivity coefficients (SCs) of soil water content (θ), infiltration rate, I(t), accumulative infiltration (Z), runoff rate (Ror) and accumulative runoff (Acc Roac) at 3 cm depth to 1% change in parameters θr, θs, alpha, n and Ks.
Figure 11. Augrabies sensitivity coefficients (SCs) of soil water content (θ), infiltration rate, I(t), accumulative infiltration (Z), runoff rate (Ro) and accumulative runoff (Acc Ro) at 3 cm depth to a 1% change in parameters θr, θs, alpha, n and Ks.
Sensitivities of soil water content to 1% change in parameters were highly variable among soil types. In the Addo, soil water content sensitivity to all parameters was limited to the first 15 minutes with well-defined peaks at 10 minutes of 0.02 and 0.01 for respective θr and θs after which became constant. Soil water content from the Augrabies showed greater sensitivity to θs and Ks after 10 minutes reaching respective peaks of 0.005 and 0.007 at 15 and 30 minutes, consecutively. Brandvlei soil water contents after 15 minutes showed sensitivity to all with the least being n parameter. Only θs had a well-defined peak at 18 min of 0.006 while α, θr and Ks parameters increased with time reaching respective peaks of 0.009, 0.008 and 0.007, before assuming nearly constant trends. Sensitivity of soil water content to Ks appeared to increase with infiltration time especially for the Augrabies and Brandvlei, soils forms, and both clay content of less than 20% in all their horizons.
Figure 12. Blandvlei’s sensitivity coefficients (SCs) of soil water content (θ), infiltration rate, I(t), accumulative infiltration (Z), runoff rate (Ro) and accumulative runoff (Acc Ro) at 3 cm depth to a 1% change in parameters θr, θs, alpha, n and Ks.
Sensitivities of infiltration rates to 1% change of parameters was characterised by defined peaks after detection of the wetting front with some parameters increasing sensitivity with time. In the Addo, except n parameter, parameters had distinct peaks within 12 to 20 minutes after which all parameter sensitivity to infiltration rate were diffused to nearly constant rates of less than 0.004. Soil water content displayed a different parameter sensitivity for the Augrabies and Brandvlei. Infiltration rate showed greater sensitivity to θs and Ks that increased after the 25 min time mark to reach sensitivities of 0.18 and 0.23, respectively. From the Brandvlei, all parameters had definite peaks at 18 minutes with the α, θr and Ks assuming sensitivities that increased with time after the 20 min time mark suggesting that increasing infiltration time measurement would have improved these parameters information. However, research has showed that sensitivity of infiltration rate to hydraulic parameters was limited to the period when of detecting wetting front to the time when gravity began to influence infiltration  . This conditions appeared to apply in the Addo which was least affected by surface compaction and raindrop impact. However, in the Augrabies and Brandvlei sensitivity of infiltration rate to parameters extended over longer period especially for θs and Ks. Clogging of pores and densification due to raindrop impact on soils with higher (≥72%) sand fraction was the reason. Parameter sensitivity showed by runoff rates was similar to that of infiltration rates for the respective soil types. This observation suggested that estimated soil hydraulic parameters from prescribed near-surface time-variable soil-moisture and atmospheric boundary conditions with surface-runoff represented the partitioning of rainwater infiltration and runoff fluxes.
Sensitivities of accumulative infiltration from soil types increased with time following rapid increase in soil water content. For the Addo sensitivities increased after 12 min from all parameters reaching 0.189 for α parameter by the end of the infiltration. Despite Ks being the least sensitive from the Addo, it was the most sensitive along with θs in the Augrabies increasing after the 30 min time mark to reach 0.447 the highest realised from this study. This trend would have continued further with time suggesting that increasing duration of the infiltration experiment would have improved parameter identification especially θs and Ks parameters. Simunek and van Genuchten  shared the same sentiment and indicated that benefit of information gained by extending duration of the infiltration experiment was limited to the occurrence of steady state conditions. In the Brandvlei sensitivity of accumulative infiltration to θs and Ks had well-defined peaks of 0.26 at the 8 min time mark. However, after the 25 min mark sensitivity to α and θr parameters increased with time a phenomenon that collaborated with the notion that measurements taken further with time would improve information about unknown parameters  . Parameter sensitivities of accumulative runoff showed similar trends to that of accumulative infiltration. This was not surprising given soil hydraulic parameters that encouraged infiltration had opposite effect on runoff.
Surface horizon hydraulic parameters controlling water infiltration and runoff of the van Genuchten-Mualem analytic model were estimated from three floodplain soil types using HYDRUS-1D model. A two-step inversion approach was used to estimate optimal parameter values for a two-layered surface horizon discretised at 0 - 50 mm and 50 - 250 mm depths. Inversion of time dependent soil water contentinfiltration measurements from the high, medium and low rainfall simulated experiments were satisfactory from all soil types with the Nash-Sutcliffe model efficiency coefficient (NSE) of not less than 0.99. Overestimation of runoff rates (NSE; 0.27) and accumulative runoff (0.53 ≤ NSE ≤ 0.98) was suggestive that inclusion of transient water infiltration data was insufficient for the water infiltration-runoff inverse problem. Inclusion of additional data to define the inverse problem such as ponding and surface water-storage parameters would require more computational time but considerably improve inverse solution and agreement between measured and predicted runoff data.
The upper surface horizon had highly optimised and variable parameters especially θs and Ks. Higher soil-surface bulk density (≥1.69 g・cm−3) had optimal Ks values lower by at least one order of magnitude to the double ring infiltrometers. This finding showed soil surface compaction and sedimentary crust to be important factors of influence for saturated hydraulic parameters, which also determined floodplain soil-type infiltration and runoff characteristics. Considering longer rainfall simulation, periods to attain saturation and steady state conditions would therefore improve model prediction of floodplain soil-types water infiltration and runoff. Significant differences (P < 0.05) in soil-types infiltration rates and accumulative infiltration under the simulated high rainstorm confirmed that surface crusting and sealing was an important factor in dryland floodplain water infiltration-runoff characteristics. Optimal parameter values were typical of fine textured, compacted sands and uniformly course textured for the respective Addo (Greysols), Augrabies (Ferralsols), Brandvlei (Cambisols) soil types.
The authors wish to acknowledge the management and staff of Kolomela Iron Ore mine (Anglo-American) for the technical, logistic and financial support of this research project.
List of Notation and Symbols
The following symbols are used in this paper:
ai = auxiliary variable;
α = air entry value also referred as bubble pressure;
β = parameter vector;
βj = specific parameter;
P < 0.05 = probability at 5% level of significance
≥ = greater than or equal to;
≤ = less than or equal to;
∂ = partial differential;
ɸ = objective function;
cosα = angle between the flow direction and the vertical axis;
ej = jth unit vector;
h = metric suction potential;
hA = minimum pressure head;
hS = maximum pressure head;
H = is the soil-water pressure head relative to atmospheric pressure;
K(h) = matric suction based unsaturated hydraulic conductivity;
K(θ) = water content based unsaturated hydraulic conductivity;
Ks = saturated or steady state hydraulic conductivity;
l = pore-connectivity parameter assumed;
m = empirical parameter
N = number of observations;
n = pore size distribution parameter;
θ = volumetricwater content;
θave = average soil water content;
θr = residual soil water content;
θs = saturated soil water content;
θi* = observed soil water contents;
θi = predicted soil water contents;
q = water flux through a specified flow domain also depicted;
Rp = shortest significant ranges;
rp = tabular values of significant studentized ranges;
S = sink;
s = variance;
= standard error of the mean difference;
t = time;
W = weight by standard deviation;
z = depth; and,
Z = vertical down wide direction or gravitational potential.