Predicting Soil Infiltration and Horizon Thickness for a Large-Scale Water Balance Model in an Arid Environment

Prediction of soil characteristics over large areas is desirable for environmental modeling. In arid environments, soil characteristics often show strong ecological connectivity with natural vegetation, specifically biomass and/or canopy cover, suggesting that the soil characteristics may be predicted from vegetation data. The objective of this study was to predict soil infiltration characteristics and horizon (soil layer) thickness using vegetation data for a large-scale water balance model in an arid region. Double-ring infiltrometer tests (at 23 sites), horizon thickness measurements (58 sites) and vegetation surveys (35 sites) were conducted in a 30 km ˆ 50 km area in Western Australia during 1999 to 2003. The relationships between soil parameters and vegetation data were evaluated quantitatively by simple linear regression. The parameters for initial-term infiltration had strong and positive correlations with biomass and canopy coverage (R2 = 0.64  ́ 0.81). The horizon thickness also had strong positive correlations with vegetation properties (R2 = 0.53  ́ 0.67). These results suggest that the soil infiltration parameters and horizon thickness can be spatially predicted by properties of vegetation using their linear regression based equations and vegetation maps. The background and reasons of the strong ecological connectivity between soil and vegetation in this region were also considered.


Introduction
Prediction of soil characteristics over large areas is desirable for environmental modeling, precision agriculture and natural resources management.In arid environments, soil permeability is one of the most important characteristics because it governs infiltration and runoff characteristics of the landscape.Horizon (soil layer) thickness above underground hardpans or basement rocks is also significant because hardpans and basement rocks are impermeable and constrain water movement.These soil characteristics are essential as input parameters or boundary conditions for water balance analysis by rainfall-runoff modeling.Rainfall-runoff models are important to predict the movement of precious water resources in arid landscapes.The accuracy of the input parameters used in such models has a significant effect on model performance [1].A lack of data for the models is a major limitation for simulation of regional hydrologic processes [2].However, it is well known that soil characteristics often show high spatial variability.Direct measurements of soil characteristics are time-consuming, labor-intensive and costly to characterize for areas >100 km 2 .
Stemming from land system surveys of the 1950s to 1970s [3][4][5], two main approaches are widely used to predict the spatial variability of soil characteristics on large scales: geostatistics [6] and environmental correlation [7].In this paper, we focus on the environmental correlation, specifically using correlations between soil and vegetation properties.An environmental correlation approach evaluates the spatial distribution of soil characteristics by quantifying the relationship between soil and environmental attributes.Terrain attributes derived from digital elevation models (DEM) are the most major predictive factor in this approach [8,9].Explicit, quantitative, and usually simple empirical models are developed to describe the relationships between georeferenced soil samples and predictive factors (e.g., [10,11]).Multiple linear regression has been most widely used in the modeling processes (e.g., [12][13][14][15]).In addition to the topographical factors, climatic, organic, lithologic and temporal factors can be used as predictors [16].
Land cover is one of the useful indicators of soil properties and class because the natural vegetation and its dynamics should represent some kind of equilibrium relation with soil type [16][17][18].Recently, vegetation data such as vegetation structure and canopy coverage can be obtained relatively easily even over large areas by remote sensing techniques.These data have been used as predictive factors in combination with terrain attributes and other factors in several studies [19][20][21].For example, vegetation indices from Landsat Thematic Mapper and Enhanced Thematic Mapper Plus (Landsat TM and ETM+) were used to predict soil classes [22], soil drainage classes [23], soil hydraulic properties [2] and soil profile depth, total phosphorus and total carbon [24].Zhu [25] used canopy coverage from Landsat TM to predict soil classes.Park and Vlek [26] used a vegetation map to predict soil properties.Vegetation indices from SPOT images were also used to predict soil organic matter [27] and soil organic carbon [28].In these studies, vegetation data are just one category of the predictive factors, however, successful deployment of this data suggests that the vegetation attributes may play an important role in the environmental correlation approach.
To predict soil infiltration properties and horizon thickness, vegetation factors may be particularly useful because these soil characteristics are highly correlated to vegetation due to the co-evolution of landscapes [29].A considerable number of studies have been reported on the strong relationship between infiltration parameters and vegetation data (e.g., [30][31][32][33][34][35][36]).Regarding the horizon thickness, Lamotte et al. [37] showed that the depth and distribution of hardpans were associated with differences in vegetation.Pracilio et al. [38] also reported that tree height was limited by the depth to hardpans in southwestern Australia.Therefore, these soil characteristics may be simply predicted from vegetation data only if quantitative and strong correlations can be found between the soil characteristics and vegetation attributes for large areas.
The objective of this study is to investigate the prediction of soil infiltration characteristics and soil horizon thicknesses using vegetation data for areas of the order of hundreds to thousands of km 2 in an arid environment.We conducted infiltration tests (at 23 sites), horizon thickness measurements (58 sites) and vegetation surveys (35 sites) in an established research area (30 km ˆ50 km) in Western Australia.The relationships between obtained soil parameters and vegetation data were evaluated quantitatively by simple linear regression analyses.The infiltration properties and horizon thickness in the whole research area were predicted from regression equations and vegetation maps derived by remote sensing.

Research Area and Measurement Sites
The Leonora research area (lat 28 ˝38 1 13" S, long 120 ˝59 1 15" E) was located about 600 km northeast from Perth, the capital city of Western Australia, covering approximately 1500 km 2 (30 km in an east-west direction and 50 km in a north-south direction) (Figure 1).This research area was originally established to assess arid zone afforestation and to estimate the effect of carbon fixation in arid lands [39].A large-scale water balance model and a tree growth simulator have also been developed to estimate water and carbon flow in the research area [40,41].Thus, information on soil characteristics for the whole research area was needed for these models and also to establish boundary conditions for the models [42].The Leonora research area (lat 28°38′13″ S, long 120°59′15″ E) was located about 600 km northeast from Perth, the capital city of Western Australia, covering approximately 1500 km 2 (30 km in an east-west direction and 50 km in a north-south direction) (Figure 1).This research area was originally established to assess arid zone afforestation and to estimate the effect of carbon fixation in arid lands [39].A large-scale water balance model and a tree growth simulator have also been developed to estimate water and carbon flow in the research area [40,41].Thus, information on soil characteristics for the whole research area was needed for these models and also to establish boundary conditions for the models [42].The area is arid: average annual rainfall measured in the area (1917-1998) is 211 mm with a standard deviation of 109 mm [43].Although the topographical gradient is less than 1% for almost the entire research area, there are many wadis (dry riverbeds containing water only during times of heavy rain).When intensive rainfall occurs, the wadis function as drainage paths for runoff to the large and shallow salt lake at the lowest elevation.The surface soil contains abundant quantities of ferric oxide and is colored red due to laterization [44].The feature of the soil in this area is the presence of the Wiluna hardpan [45], which is extremely hard, impermeable and consists of laminated layers of silica with ferruginous and calcareous cement, and which is often up to 10 m thick [46].According to a fluorescent X-ray analysis, these hardpans in the research area consist of 45% Al2O3, 39% SiO2 and 12% Fe2O3.
The research area has been used for extensive grazing by one farm family, but all of the forest and woodland in the research area is considered natural forest and woodland [47].The vegetation of the area is a mixture of Acacia woodland, Acacia shrub land (Acacia aneura with small shrubs, e.g., Eremophilla sp.) and bare ground [47,48].About half of the research area is bare ground or low open woodland, where canopy coverage is less than 10%, and the other half is woodland and open forest, where canopy coverage varies between 10% and 50% [49].The dominant species is Acacia aneura; Eucalyptus camaldulensis open forest and woodland are adjacent to some wadis; and salt-tolerant shrubs and bushes are observed around the salt lake.Abe et al. [50] reported that the biomass distribution is related with the slope, position of the wadi and the distance from the wadi.
Measurements on soil infiltration, horizon thickness and vegetation were conducted at a maximum of 58 sites of size 50 m × 50 m.The soil surveys were intensively conducted in November The area is arid: average annual rainfall measured in the area (1917-1998) is 211 mm with a standard deviation of 109 mm [43].Although the topographical gradient is less than 1% for almost the entire research area, there are many wadis (dry riverbeds containing water only during times of heavy rain).When intensive rainfall occurs, the wadis function as drainage paths for runoff to the large and shallow salt lake at the lowest elevation.The surface soil contains abundant quantities of ferric oxide and is colored red due to laterization [44].The feature of the soil in this area is the presence of the Wiluna hardpan [45], which is extremely hard, impermeable and consists of laminated layers of silica with ferruginous and calcareous cement, and which is often up to 10 m thick [46].According to a fluorescent X-ray analysis, these hardpans in the research area consist of 45% Al 2 O 3 , 39% SiO 2 and 12% Fe 2 O 3 .
The research area has been used for extensive grazing by one farm family, but all of the forest and woodland in the research area is considered natural forest and woodland [47].The vegetation of the area is a mixture of Acacia woodland, Acacia shrub land (Acacia aneura with small shrubs, e.g., Eremophilla sp.) and bare ground [47,48].About half of the research area is bare ground or low open woodland, where canopy coverage is less than 10%, and the other half is woodland and open forest, where canopy coverage varies between 10% and 50% [49].The dominant species is Acacia aneura; Eucalyptus camaldulensis open forest and woodland are adjacent to some wadis; and salt-tolerant shrubs and bushes are observed around the salt lake.Abe et al. [50] reported that the biomass distribution is related with the slope, position of the wadi and the distance from the wadi.Measurements on soil infiltration, horizon thickness and vegetation were conducted at a maximum of 58 sites of size 50 m ˆ50 m.The soil surveys were intensively conducted in November and December in 1999 and 2002.The vegetation surveys were conducted from 1999 to 2003.All sites were classified into three vegetation categories by the dominant species at the sites: Eucalyptus (EP), Acacia aneura (AA), or bush or bare ground (BB) sites.The locations of sites on which measurements were taken were chosen by preliminary ground inspections and interpretation of aerial photographs in order to satisfy the following conditions: (i) the ecological and topographical conditions within each site were homogeneous: and (ii) there was a wide range in the biomass (Figure 1).In addition, the distances of the sites from the nearest wadis were measured by visual interpretations of aerial photographs in order to categorize the sites into "around wadis" (within 100 m of the nearest wadis) and "others".The reason for the categorization was that the areas around wadis had a characteristic ecosystem.The distribution of wadis was able to be determined also from terrain analysis of a DEM [50].

Infiltration Data
Double-ring infiltrometer tests were performed at 23 sites (4 EP, 10 AA and 9 BB sites) (Figure 1).All other variables shown in Sections 2.2.2 and 2.2.3 were also recorded at these sites.At each site, one point was chosen that gave a satisfactory representation of the surface soil condition of the site.An inner ring, 0.30 m in diameter and 0.30 m deep, was carefully driven 0.10 m into the soil at the chosen point.Each infiltration run continued for 90 min by which time a steady state rate had been attained.Soil profile surveys to measure depth to hardpans were also conducted at all sites at which measurements were made.The horizons were generally dry during the measurement period and the variation of the volumetric water contents was small (0.02-0.08 m 3 ¨m´3 ).The effect of the difference in the soil moisture on infiltration rate was therefore small.At one site with little vegetation, infiltration tests were performed at three locations to evaluate the spatial variability in the soil permeability within the site.
Two commonly used infiltration equations were fitted to the infiltration data collected during tests: the Kostiakov [51] and Philip [52] equations.Cumulative infiltration (I, mm) is expressed by the Kostiakov equation as: where t (min) is the time that water has been in contact with the soil, and k and m are constants.The Kostiakov equation is empirical; its parameters do not have clear physical meanings but it is simple and easy to define.Note that if t = 1 in Equation ( 1), then I = k.This means k is 1 minute of infiltration.
In contrast, the Philip equation is based on physical parameters and given by: where S (mm¨min ´1{2 ) accounts for the sorption of water at the beginning of the infiltration process and A (mm¨min ´1) accounts for the long-term (steady-state) infiltration rate.We introduced a parameter I t (mm), which is the measured cumulative infiltration for any period of time (t).As above, parameters k, m, S, A and I t were obtained as soil infiltration characteristics at each of the sites.
According to visual interpretation of the aerial photographs, six of the 23 sites were within 100 m of the nearest wadis.Analyses excluding these six sites were also performed.A subscript _ew was added to the infiltration parameters when the six sites were excluded in the analyses (e.g., k _ew ).

Horizon Thickness Data
Horizon thickness was measured at 58 sites (10 EP, 25 AA and 23 BB sites).The surface horizon above the underground hardpan was dug out by engine augers and the horizon thickness (T h , cm) was determined by measuring the depths of the augered holes.The horizon thicknesses sometimes Water 2016, 8, 96 5 of 15 varied within a short distance thus they were measured at 10-15 locations within each site area, and the average horizon thickness was adopted as the T h data at the site.The maximum measurement range of T h was limited to 115 cm by the length of the drill.Therefore, at 10 sites (seven EP and three AA sites) where the horizon thickness was over 115 cm, T h was tentatively assigned a value of 115 cm.Analyses excluding these 10 sites were also performed.

Vegetation Data
The stand biomass (B m , kg¨m ´2) and canopy coverage (C c , m 2 ¨m´2 ) were measured or estimated at every site.The stand biomass was defined as the sum of the biomass of the individual trees above ground that comprised stand per the unit ground area.The canopy coverage was defined as the ratio of the vertical projection area of the forest and woodland canopy per unit area [47].C c was obtained relatively easily from both vegetation surveys and remote sensing analyses.In contrast, B m was more difficult to estimate than C c but we considered that it was a more efficacious measure of the amount of vegetation.
At 35 of the 58 sites, vegetation surveys were conducted originally for the remote sensing analyses.Various vegetation parameters were measured for all trees such as tree height, diameter at breast height and canopy silhouette area.Each tree biomass was calculated by the allometric equation for each tree species using a destructive sampling method [53].From these results, the values of B m and C c were calculated for each site.To extrapolate the site based B m and C c values to the entire study area, satellite images (Landsat TM and ETM+; 16 scenes from June 1998-December 1999; 30 m resolution) and aerial photographs (on 10-11 November 1999; digitized ortho color photograph consisted of Red, Green and Blue bands, 0.25 m quadrate) were used.As results of analyses for both images, distribution maps of B m and C c were made for the entire area.The details of the estimation methods adopted are described by Taniguchi et al. [54] and Suganuma et al. [47,49].An example of the B m distribution map developed by satellite image analyses is shown in Figure 2 [54].At 23 sites at which vegetation surveys were not conducted (Figure 1), B m and C c were estimated from the aerial photographs that had higher resolution and accuracy than the satellite images.The C c values were estimated by visual interpretation of the aerial photographs (e.g., R 2 = 0.99 [49]).The B m values were estimated using high correlation between B m and C c for each vegetation type (e.g., R 2 = 0.92 [47]).
Water 2016, 8, 96 5 of 15 range of Th was limited to 115 cm by the length of the drill.Therefore, at 10 sites (seven EP and three AA sites) where the horizon thickness was over 115 cm, Th was tentatively assigned a value of 115 cm.Analyses excluding these 10 sites were also performed.

Vegetation Data
The stand biomass (Bm, kg•m −2 ) and canopy coverage (Cc, m 2 •m −2 ) were measured or estimated at every site.The stand biomass was defined as the sum of the biomass of the individual trees above ground that comprised stand per the unit ground area.The canopy coverage was defined as the ratio of the vertical projection area of the forest and woodland canopy per unit area [47].Cc was obtained relatively easily from both vegetation surveys and remote sensing analyses.In contrast, Bm was more difficult to estimate than Cc but we considered that it was a more efficacious measure of the amount of vegetation.
At 35 of the 58 sites, vegetation surveys were conducted originally for the remote sensing analyses.Various vegetation parameters were measured for all trees such as tree height, diameter at breast height and canopy silhouette area.Each tree biomass was calculated by the allometric equation for each tree species using a destructive sampling method [53].From these results, the values of Bm and Cc were calculated for each site.To extrapolate the site based Bm and Cc values to the entire study area, satellite images (Landsat TM and ETM+; 16 scenes from June 1998-December 1999; 30 m resolution) and aerial photographs (on 10-11 November 1999; digitized ortho color photograph consisted of Red, Green and Blue bands, 0.25 m quadrate) were used.As results of analyses for both images, distribution maps of Bm and Cc were made for the entire area.The details of the estimation methods adopted are described by Taniguchi et al. [54] and Suganuma et al. [47,49].An example of the Bm distribution map developed by satellite image analyses is shown in Figure 2 [54].At 23 sites at which vegetation surveys were not conducted (Figure 1), Bm and Cc were estimated from the aerial photographs that had higher resolution and accuracy than the satellite images.The Cc values were estimated by visual interpretation of the aerial photographs (e.g., R 2 = 0.99 [49]).The Bm values were estimated using high correlation between Bm and Cc for each vegetation type (e.g., R 2 = 0.92 [47]).4 and 7a).The stand biomass was estimated using SAVI (Soil Adjusted Vegetation Index) from satellite images (Landsat TM and ETM+; 16 scenes from June 1998-December 1999; 30 m quadrate) [54].4 and 7a).The stand biomass was estimated using SAVI (Soil Adjusted Vegetation Index) from satellite images (Landsat TM and ETM+; 16 scenes from June 1998-December 1999; 30 m quadrate) [54].

Data Analysis and Development of Prediction Equation
Water 2016, 8, 96 6 of 15

Data Analysis and Development of Prediction Equation
Liner regression analysis was used to quantify relationships between the infiltration parameters or horizon thickness and vegetation data.In this paper, soil characteristics were predicted from vegetation data, thus k, m, S, A, I t and T h were used as the dependent variables (y) and B m and C c were used as the independent variables (x).In addition to regression parameters, some statistical measures were provided to evaluate the model performance: mean (M) and standard deviation (SD) for each observed and predicted variable, mean bias error (MBE) and root mean square error (RMSE) [55].
The dependent variables highly correlated with the vegetation variables were predicted using their linear regression based equations.Prediction equations of time series of infiltration were also developed by substituting the regression equations for each infiltration constant into the original infiltration equations (Equations ( 1) and ( 2)).The detailed development is shown in the results section (Section 3.2).

Relationships between Infiltration Parameters and Vegetation Data
The cumulative infiltrations varied greatly from site to site, suggesting a high spatial variability in the soil permeability in the research area (Figure 3).At a final cumulative infiltration (90 min), the maximum value was 685 mm, minimum 37 mm, mean 226 mm and standard deviation 155 mm.Both Equations ( 1) and ( 2) were in close agreement with the measured cumulative infiltration at every site; The coefficients of determination of the fittings were within a range of 0.99 to 1.00, and the RMSE between the measured and fitted cumulative infiltrations were within a range of 0.45 mm to 10.0 mm.At one site where the infiltration tests were performed at three locations, the final cumulative infiltrations (I 90 ) of three tests were 180, 193 and 214 mm and the standard deviation was 17 mm, suggesting the spatial variability in the soil permeability was small within the site.
Water 2016, 8, 96 6 of 15 Liner regression analysis was used to quantify relationships between the infiltration parameters or horizon thickness and vegetation data.In this paper, soil characteristics were predicted from vegetation data, thus k, m, S, A, It and Th were used as the dependent variables (y) and Bm and Cc were used as the independent variables (x).In addition to regression parameters, some statistical measures were provided to evaluate the model performance: mean (M) and standard deviation (SD) for each observed and predicted variable, mean bias error (MBE) and root mean square error (RMSE) [55].
The dependent variables highly correlated with the vegetation variables were predicted using their linear regression based equations.Prediction equations of time series of infiltration were also developed by substituting the regression equations for each infiltration constant into the original infiltration equations (Equations ( 1) and ( 2)).The detailed development is shown in the results section (Section 3.2).

Relationships between Infiltration Parameters and Vegetation Data
The cumulative infiltrations varied greatly from site to site, suggesting a high spatial variability in the soil permeability in the research area (Figure 3).At a final cumulative infiltration (90 min), the maximum value was 685 mm, minimum 37 mm, mean 226 mm and standard deviation 155 mm.Both Equations ( 1) and ( 2) were in close agreement with the measured cumulative infiltration at every site; The coefficients of determination of the fittings were within a range of 0.99 to 1.00, and the RMSE between the measured and fitted cumulative infiltrations were within a range of 0.45 mm to 10.0 mm.At one site where the infiltration tests were performed at three locations, the final cumulative infiltrations (I90) of three tests were 180, 193 and 214 mm and the standard deviation was 17 mm, suggesting the spatial variability in the soil permeability was small within the site.Table 1 shows the output of the resulting linear regression analysis relating vegetation data to infiltration parameters and statistical measures.Several infiltration parameters correlated highly with vegetation data.The constants k in the Kostiakov equation and S in the Philip equation showed particularly high and positive correlations with Bm and Cc (Figure 4).Both k and S are measures of the infiltrability at the beginning of the infiltration process and results indicate that surface soil permeability was strongly correlated with the vegetation.With an increase in the t value of It, the correlation coefficients between It and Bm or Cc decreased (Table 1).Thus, the rate of infiltration into shallower soil has a stronger correlation with vegetation.No significant correlation was obtained between the constant A of Philip's equation and Bm or Cc.This result suggests that permeability of horizons low in the profile is not correlated with the vegetation because the constant A accounts for the long-term infiltration rate.Table 1 shows the output of the resulting linear regression analysis relating vegetation data to infiltration parameters and statistical measures.Several infiltration parameters correlated highly with vegetation data.The constants k in the Kostiakov equation and S in the Philip equation showed particularly high and positive correlations with B m and C c (Figure 4).Both k and S are measures of the infiltrability at the beginning of the infiltration process and results indicate that surface soil permeability was strongly correlated with the vegetation.With an increase in the t value of I t , the correlation coefficients between I t and B m or C c decreased (Table 1).Thus, the rate of infiltration into shallower soil has a stronger correlation with vegetation.No significant correlation was obtained Water 2016, 8, 96 7 of 15 between the constant A of Philip's equation and B m or C c .This result suggests that permeability of horizons low in the profile is not correlated with the vegetation because the constant A accounts for the long-term infiltration rate.
When the six sites within 100 m of the nearest wadis were excluded from the regression analyses, the R 2 values between I 60_ew or A _ew and the vegetation data substantially increased compared with those of I 60 or A (Table 1).The initial infiltrability, k _ew and S _ew , also showed a strong correlation with vegetation data, suggesting that the soil permeability from the surface to depth is strongly correlated with the vegetation except close to wadis.

Prediction of Infiltration Parameters from Vegetation Data
The initial infiltrability (k, S and I 5 ) can be predicted with high accuracy from B m and C c over the whole area using the linear regression coefficients shown in Table 1.For example, S (sorptivity: mm¨min ´1{2 ) can be expressed as a function of B m using S (Bm) = 2.98 + 0.96B m from Figure 4, where S (Bm) is the estimated sorptivity from B m .Using this equation, Figure 2 is regarded as the distribution map of S for the whole research area.The long-term infiltrability (I 60 and A), which was poorly correlated, was not predicted for the whole area, however, it can be predicted as well as the initial infiltrability if the area of within 100 m wadis is eliminated.I 60 is an important long-term infiltrability parameter because it is used to model hourly runoff generation when hourly rainfall data are available.I 60_ew showed a strong correlation with B m and C c .Therefore, the prediction of hourly runoff becomes possible by the concurrent use of the hourly rainfall data, the regression equation for I 60_ew , and distribution maps of wadis and B m or C c .
An example of the development of the prediction equation of time series infiltration, using the relationship between k, m and B m is given below.Both constants k and m of the Kostiakov equation are closely correlated with B m .The regression equations between k and B m , and m and B m from Table 1 are k = 2.53 + 1.16B m and m = 0.86 ´0.020B m .Substituting these into Equation (1) gives the estimated cumulative infiltration, I K(Bm) (mm), as: Simultaneous use of Equation ( 3) and the B m distribution map (Figure 2) enables the prediction of cumulative infiltration at any elapsed time and at any place in the area.The value of I K(Bm) increases with increases in t and B m (Figure 5).However, at high values of t and B m , I K(Bm) decreases with an increase in B m .The decrease in I K(Bm) with the increase in B m is incongruous with the observed positive correlation between infiltrability and vegetation.
Water 2016, 8, 96 9 of 15 Simultaneous use of Equation ( 3) and the Bm distribution map (Figure 2) enables the prediction of cumulative infiltration at any elapsed time and at any place in the area.The value of IK(Bm) increases with increases in t and Bm (Figure 5).However, at high values of t and Bm, IK(Bm) decreases with an increase in Bm.The decrease in IK(Bm) with the increase in Bm is incongruous with the observed positive correlation between infiltrability and vegetation.A similar procedure was used for the other infiltration parameters.The relationships between k_ew and m_ew for Equation (1), and S_ew and A_ew for Equation ( 2) are shown below.The prediction equations are presented as, respectively: where IK(Bm)_ew and IP(Bm)_ew are the estimated cumulative infiltrations in mm.The subscript K and P denote that the bases of the equations are Kostiakov and Philip equations, respectively.These equations can be applied to the study area excluding the wadi surrounds.A similar procedure was used for the other infiltration parameters.The relationships between k _ew and m _ew for Equation (1), and S _ew and A _ew for Equation ( 2) are shown below.The prediction equations are presented as, respectively: I P(Bm)_ew " p2.69 `0.91B m qt 1/2 `p1.05 `0.17B m qt, where I K(Bm)_ew and I PpBm)_ew are the estimated cumulative infiltrations in mm.The subscript K and P denote that the bases of the equations are Kostiakov and Philip equations, respectively.These equations can be applied to the study area excluding the wadi surrounds.
The prediction accuracy of Equations ( 3)-( 5) was evaluated by the coefficient of determination (R 2 ) between the estimated cumulative infiltration and I t (measured cumulative infiltration) for all sites.Figure 6 presents the variations in the R 2 values with time.The R 2 value between I K(Bm) and I t was initially high (R 2 > 0.6, t < 10 min), but decreased rapidly with time because of: (i) the low correlation between long-term infiltrability and vegetation; and (ii) the decrease in I K(Bm) with the increase in B m in the high ranges of t and B m (Figure as previously mentioned.We propose that Equation (3) should be used only for the prediction of initial-term infiltrations, such as for the prediction of infiltrations and runoff generations for brief rainfall events (t < 20 min).In contrast, the R 2 value between I K(Bm)_ew or I PpBm)_ew and I t was more than 0.7 for all time periods, suggesting that even long-term infiltration can be predicted by Equations ( 4) or ( 5) with an acceptable degree of accuracy in the area outside the wadi surrounds.There was no significant difference between Equations ( 4) and ( 5) in prediction accuracy.The RMSE between the measured cumulative infiltration and I K(Bm) , I K(Bm)_ew and I PpBm)_ew at 60 min were 83.5, 35.6 and 36.0 mm, respectively.These values were slightly smaller than the RMSE values of B m -I 60 and B m -I 60_ew (Table 1), suggesting the developed time series equations predicted the cumulative infiltration more accurately than the original regression equations.

Prediction of Horizon Thickness from Vegetation Data
The horizon thickness had strong positive correlations with Bm (R 2 = 0.67; n = 58; p < 0.001) and Cc (R 2 = 0. 56; n = 58; p < 0.001) (Figure 7).Th can be expressed as a function of Bm using Th(Bm) = 25.40 + 6.00Bm from Figure 7a, where Th(Bm) is the estimated horizon thickness from Bm in cm.Using this equation, Figure 2

Prediction of Horizon Thickness from Vegetation Data
The horizon thickness had strong positive correlations with B m (R 2 = 0.67; n = 58; p < 0.001) and C c (R 2 = 0. 56; n = 58; p < 0.001) (Figure 7).T h can be expressed as a function of B m using T h(Bm) = 25.40 + 6.00B m from Figure 7a, where T h(Bm) is the estimated horizon thickness from B m in cm.Using this equation, Figure 2 is regarded as the distribution map of T h for the whole research area.When 10 sites, the horizon thicknesses of which were greater than 115 cm were excluded in regression analyses, T h showed strong correlation with B m (R 2 = 0.58; n = 48; p < 0.001) and C c (R 2 = 0.53; n = 48; p < 0.001), meaning the excluding these sites does not affect the correlation relationship.Eight of these 10 sites and another two sites were located within 100 m of the nearest wadi.If the total of 12 irregular sites, T h > 115 cm or around wadis, were excluded in regression analyses, T h still showed strong correlation with B m (R 2 = 0.55; n = 46; p < 0.001) and C c (R 2 = 0.53; n = 46; p < 0.001).
equation, Figure 2 is regarded as the distribution map of Th for the whole research area.When 10 sites, the horizon thicknesses of which were greater than 115 cm were excluded in regression analyses, Th showed strong correlation with Bm (R 2 = 0.58; n = 48; p < 0.001) and Cc (R 2 = 0.53; n = 48; p < 0.001), meaning the excluding these sites does not affect the correlation relationship.Eight of these 10 sites and another two sites were located within 100 m of the nearest wadi.If the total of 12 irregular sites, Th > 115 cm or around wadis, were excluded in regression analyses, Th still showed strong correlation with Bm (R 2 = 0.55; n = 46; p < 0.001) and Cc (R 2 = 0.53; n = 46; p < 0.001).

Discussion
The six sites within 100 m of the nearest wadis had relatively high Bm and Cc values, and showed high initial infiltration rates (Figure 4).From the result of soil profile surveys, horizons with low permeability were observed at three of these six sites.The depths of these horizons were about 20-50 cm from the soil surface.Although these horizons were not hardpans, they were extremely compacted and hard; the values for the long-term infiltration parameters were strongly affected by them.It seemed that these horizons of low permeability served as transient river beds around wadis

Discussion
The six sites within 100 m of the nearest wadis had relatively high B m and C c values, and showed high initial infiltration rates (Figure 4).From the result of soil profile surveys, horizons with low permeability were observed at three of these six sites.The depths of these horizons were about 20-50 cm from the soil surface.Although these horizons were not hardpans, they were extremely compacted and hard; the values for the long-term infiltration parameters were strongly affected by them.It seemed that these horizons of low permeability served as transient river beds around wadis during large runoff events.In contrast, there were no horizons of low permeability in the other three sites within 100 m of wadis and the soil permeability was high even at depth.These results indicate that the spatial variability of permeability was particularly high in lower horizons around wadis.
The constant m in the Kostiakov equation had a negative and strong correlation with B m and C c (Figure 4).This constant does not have clear physical meaning, but when the m value is small, the infiltration rate rapidly decreases with time (Equation ( 1)).Thus, m values were particularly small at the three sites around wadis where there are horizons of low permeability.These m values enhanced the negative slope of the regression line (b) between m and B m or C c (Table 1 and Figure 4).Therefore, the slope values of m _ew for both B m and C c were smaller than the slope values of m for B m and C c .The prediction equation of time-series infiltration was also affected by the steep negative gradient due to the small m values; this was the reason why the predicted I K(Bm) by Equation (3) decreased with an increase in B m at high values of t and B m (Figure 5).BB (bush or bare ground) sites were concentrated in the lower range of k and S, and the lower range of B m and C c (Figure 4).In contrast, EP (Eucalyptus) sites concentrated in the upper range of k and S, and B m and C c , respectively.AA (Acacia aneura) sites were distributed along regression lines, which means the stand biomass and canopy coverage of Acacia aneura spatially varied and correlated with the surface soil permeability.
It appeared that the main reason for the strong correlation between horizon thickness and vegetation was the restriction of the rooting depth for plants in the vicinity of hardpans [38,56].Furthermore, vegetation may also affect the horizon thickness as a natural selection process and for its advantage [57][58][59].For example, rich vegetation supplies surface soils through litter falls and micro-organisms, in contrast, poor vegetation cover often causes runoff and erosion which reduces the horizon thickness.BB sites were concentrated in the lower ranges of T h and B m or C c , whereas EP sites were concentrated in higher ranges (Figure 7).These tendencies were similar to the distributions of BB and EP shown in the relationships between initial infiltrability and vegetation data (Figure 4).AA sites were dispersed widely in thin to deep soil areas; Acacia aneura may be able to grow regardless of the limitation of depth to the hardpan.This would be one of the reasons why Acacia aneura was the dominant species in this area where hardpans distribute in shallow underground.
It seemed that the strong correlation between soil infiltration, horizon thickness and vegetation was caused by their ecological connectivity through water movement.Such associations may be particularly pronounced in arid lands because water is the main driving force of vegetation distribution and hence ecosystems in arid environments (e.g., [60][61][62]).The regression equations presented in this paper are region-specific.We contend that the same concept of prediction of these soil characteristics would be applicable for other arid lands that have similar ecosystems.
In most cases, above ground stand biomass showed slightly higher correlations with soil characteristics than canopy coverage (Table 1).We propose that this was because the stand biomass obtained in this study included the allometric equation for each tree species in the calculation, and hence it was a more robust measure of the current condition and amount of vegetation than canopy coverage.However, canopy coverage can be obtained more easily than biomass.The choice of either stand biomass or canopy coverage as the predictor will depend on the required accuracy and the resources required to collect the field data.
Strong statistical relationships are proven between vegetation indices calculated from satellite imagery such as NDVI (Normalized Difference Vegetation Index) or SAVI (Soil Adjusted Vegetation Index) and biomass or canopy coverage.Vegetation indices can be obtained more easily for even larger areas than our current study site; they may be able to be effective soil characteristic predictors.Further studies are required to explore other predictors based on remote sensing.Improvement of field measurements is another task.To determine depths to deep hardpans, electric resistivity surveys may be useful [37].We used the double-ring infiltrometer for the convenience of use and transportation; however, infiltrability is often overestimated under ponding conditions compared with dripping conditions (e.g., [63,64]).Quantification of spatial variability of infiltrability within sites may also be important, specifically at high vegetation sites.Conducting multiple infiltration tests within a site using rainfall simulators will contribute to increase in accuracy of runoff prediction.

Conclusions
Soil infiltration characteristics and horizon thickness were predicted using vegetation data for a 1500 km 2 site in an arid environment.The parameters for the initial-term infiltration and horizon thickness had strong and positive correlations with the biomass and canopy coverage over the entire area studied.The parameters for the long-term infiltration also had strong, positive correlations with vegetation when the data obtained around wadis was excluded.We conclude that the infiltration parameters and horizon thickness above hardpans can be predicted by use of the linear regression equations with vegetation data.The prediction equations of time series infiltration were developed by substituting the regression equations for each infiltration parameter into the original infiltration equations.Concurrent use of these equations and the distribution map of vegetation (and wadis) enabled the prediction of the cumulative infiltration at any time for place in the research area.We contend that the same concept of prediction of soil characteristics would be applicable for other arid lands that have similar ecosystems.

Figure 1 .
Figure 1.Map of the research area and the distribution of the measurement sites.The darker the line is, the larger the catchment area of the wadi is.

Figure 1 .
Figure 1.Map of the research area and the distribution of the measurement sites.The darker the line is, the larger the catchment area of the wadi is.

Figure 2 .
Figure 2. Distribution map of stand biomass and the distributions of sorptivity and horizon thickness predicted by the regression equations (Figures4 and 7a).The stand biomass was estimated using SAVI (Soil Adjusted Vegetation Index) from satellite images (Landsat TM and ETM+; 16 scenes from June 1998-December 1999; 30 m quadrate)[54].

Figure 2 .
Figure 2. Distribution map of stand biomass and the distributions of sorptivity and horizon thickness predicted by the regression equations (Figures4 and 7a).The stand biomass was estimated using SAVI (Soil Adjusted Vegetation Index) from satellite images (Landsat TM and ETM+; 16 scenes from June 1998-December 1999; 30 m quadrate)[54].

Figure 3 .
Figure 3. Cumulative infiltrations as a function of time at 23 sites.

Figure 3 .
Figure 3. Cumulative infiltrations as a function of time at 23 sites.

Figure 4 .
Figure 4. Relationships between infiltration parameters (k, m, S and A) and vegetation data (Bm and Cc).BB = bush or bare ground site; AA = Acacia aneura site; EP = Eucalyptus site.

Figure 4 .
Figure 4. Relationships between infiltration parameters (k, m, S and A) and vegetation data (B m and C c ). BB = bush or bare ground site; AA = Acacia aneura site; EP = Eucalyptus site.

R 2 =
Correlation coefficient.a and b = Regression coefficients (intercept and slope, respectively).n = Number of sites.M = mean.SD = standard deviation.MBE = mean bias error.RMSE = root mean square error.I 5 , I 10 , I 30 and I 60 = Measured cumulative infiltration at 5, 10, 30 and 60 min, respectively.p = P-value for a correlation coefficient test._ew = Denotes the exclusion of the sites located within a 100 m of the nearest wadis._o and _p = Denote the observed and predicted variables, respectively.

Figure 5 .
Figure 5.Estimated cumulative infiltration by Equation (3) (I K (Bm) : mm) as a function of stand biomass (B m : kg¨m ´2) and time (t: min).

Figure 6 .
Figure 6.Variations in the coefficient of determination (R 2 ) values between measured and estimated cumulative infiltrations for all sites with time.
is regarded as the distribution map of Th for the whole research area.When 10 sites, the horizon thicknesses of which were greater than 115 cm were excluded in regression analyses, Th showed strong correlation with Bm (R 2 = 0.58; n = 48; p < 0.001) and Cc (R 2 = 0.53; n = 48; p < 0.001), meaning the excluding these sites does not affect the correlation relationship.Eight of these 10 sites and another two sites were located within 100 m of the nearest wadi.If the total of 12 irregular sites, Th > 115 cm or around wadis, were excluded in regression analyses, Th still showed strong correlation with Bm (R 2 = 0.55; n = 46; p < 0.001) and Cc (R 2 = 0.53; n = 46; p < 0.001).

Figure 6 .
Figure 6.Variations in the coefficient of determination (R 2 ) values between measured and estimated cumulative infiltrations for all sites with time.

Figure 7 .
Figure 7. Relationships between (a) horizon thickness (T h ) and stand biomass (B m ) and (b) horizon thickness (T h ) and canopy coverage (C c ). BB = bush or bare ground site; AA = Acacia aneura site; EP = Eucalyptus site.

Table 1 .
Linear regression equations relating the vegetation data to the infiltration parameters and statistical measures.

Table 1 .
Linear regression equations relating the vegetation data to the infiltration parameters and statistical measures.