Abstract

Key Words

Deep drainage, salinity, regression kriging, SWAP

Introduction

There is continued concern about the changes in land-use in Australia and the impact on the hydrological cycle, primarily because of increased deep drainage (DD) (Walker et al. 2002). DD can be described as “potential recharge” or as a loss of water below the root zone depending on the researchers viewpoint (Asseng et al. 2001). As increased DD is recognised as one of the main drivers for dryland and irrigated salinity (Vervoort et al. 2002, Walker et al. 2002), it is important to have some estimate of the risk of such events occurring. Risk is an estimation of the expected amount of harm that will occur to an asset when a condition, such as excessive DD, occurs. For example, salinity risk is the probability that certain management practices contribute to the expression of land or water salinity in the landscape (Collins and Donaldson 2002). As such, it is an expansion of the concept of hazard, which only describes the intrinsic nature of the issue. Increasingly hazard and risk maps are being used to aid management decisions, but most of the mapping in Australia appears to be based on subjective and empirical approaches (Gilfedder and Walker, 2001) with only a few studies employing the use of models to estimate hazard or risk (Bui et al. 1996; Asseng et al. 2001 Triantafilis et al. 2003). Furthermore mapping in Australia has focussed primarily on salinity hazard, based on landscape characteristics and salt stores.

Risk of DD is of particular interest in the Northern Murray Darling Basin (NMDB), as past research and management of the water balance has mainly focussed on infiltration, runoff and soil conservation (Vervoort et al. 2004). But since irrigated cotton production is an important and growing industry in the area, DD, water-use efficiency and minimising the impact of the cotton industry on the environment are issues of growing importance. While there are a sufficient number of physically based field methods for determining DD, these methods are labour intensive and expensive (Silburn et al. 2004). An extensive range of models exists for describing soil-water movement and models are an effective way of estimating DD across an area (Walker and Zhang 2002). But, because of the one-dimensionality of most soil vegetation models and because deep drainage constitutes only a small fraction of the total water balance, little success has been made in extrapolating DD values beyond the field scale (Walker et al. 2002). Triantafilis et al. (2003) used an empirical model to derive DD risk on a field scale and Bui et al. (1996) estimated DD hazard using a pedological model and tried to validate some of the local outcomes using the mechanistic model SWIM. Asseng et al. (1997) used the model APSIM to derive DD exceedance values based on 24 rainfall stations and 5 soil types. These values were subsequently interpolated using ordinary kriging for the Southwest Australian wheatbelt. This study uses remote sensing data (such as LANDSAT) and geostatistics (Odeh et al. 1995), to interpolate the outcomes of a one-dimensional mechanistic soil-water movement model in the landscape. In addition, uncertainty due to unknown hydraulic properties is included into the model using Monte-Carlo simulations.

Methods

Study Area

For this research a 50km by 60km section of the Border Rivers Catchment was chosen. The area encompasses the main town in the region, Goondiwinidi, and the border of NSW and Queensland: the Macintyre River (150⁰19’52.4S, 28⁰21’22E). This region has summer dominated rainfall with average rainfall of 600mm/year and evapotranspiration rates which generally exceed the annual rainfall. The landscape is characterised by level plains with very low gradients.

Figure 1. Landuse in the Border Rivers catchment study area based on a classified LANDSAT-7 Image and locations of the soil sampling points (black crosses), January 2002. The Macintyre River can be seen in the middle of the image, bordered by natural vegetation (light gray), and irrigated cotton fields (dark gray), grazing and cropping (middle gray).

In order to establish the upper boundary for the one-dimensional model, it was necessary to determine the land use for the area. A LANDSAT-7 image of the area taken in January 2002 was acquired (Figure 1). LANDSAT-7 images contain data from 6 bands or channels of energy (blue, green, red, near infrared, shortwave infrared and thermal infrared). In January 2003, land use in the area was also mapped using a hand-held GPS and this data was imported into ARCGIS (ESRI Australia, Melbourne, Victoria). Through contacting landholders in the study area, and with the aid of farm maps and a copy of the LANDSAT image, land use at a field scale could be determined for January 2002. This information was combined with the LANDSAT data to perform a supervised classification using IMAGINE (Leica Geosystems, Atlanta, USA). For the land use map in this study, initial classification of obvious features and the field information meant that meaningful classes in the scene could be identified. For example cotton has a strong spectral signature in January, appearing dark green in the LANDSAT (Figure 1) and was therefore easily recognised. The main identified land uses were cotton, pasture, fallow (winter wheat) and natural vegetation (trees).

SWAP

Deep drainage (DD) was estimated using the one-dimensional finite difference model SWAP (Soil Water Atmosphere and Plant) (Van Dam et al. 1997). The SWAP upper boundary is located just above the vegetation, while free drainage at a depth of 3 m was simulated as the lower boundary. SWAP is a one-dimensional model, and does not consider lateral flow, but considering the low gradients in the landscape this was probably not a major assumption. Ten years of climate data was sourced from a weather station in the area (Goondiwindi airport). Reference evapotranspiration (ETref) was based on the Penman-Monteith equation. While rainfall data was more abundant or could have been obtained from patched and interpolated datasets such as DATADRILL (www.bom.gov.au), climate was assumed to be uniform across the study area.

A previously conducted soil survey (Odeh et al. 2004; J.Triantafilis, pers. comm.) provided particle size distributions for 91 soil sites in the study area (Figure 1). This soil dataset contained samples to a depth of 1.2m and provided the clay, silt and sand percentages. The data was first merged from 6 layers to 3 layers. The final merged layers were from depths 0-30cm, 30-60cm and 60-300cm. The layers created were arbitrary, but could be thought to represent a cultivation layer, subsoil and deeper subsoil. Mean profile clay contents in the area ranged from 12 to 77%, with a mean of 52 % (Table 1), typical for the heavy clays in the area. Six soil points appeared to be located in irrigation water storages, a landuse difficult to recreate using SWAP. Therefore the landuses that described the remaining 85 soils points were cotton, perennial pasture, native vegetation, and winter wheat (Table 1). Crop heights rather than crop factors were used to calculate actual evapotranspiration from ETref. Since, for native vegetation the depth of the soil profile had to be 3 m to cover the maximum root depth (Table 1), all soil profiles were extended from 1.2m to 3m. Comparing soil profiles with different depths might have given a correct estimate of deep drainage below the rootzone, but would not have been comparable in terms of potential recharge to the groundwater (Jolly and Cook 2002).

Table 1. Statistical summary of the soil properties in the study area indicating mean, maximum, minimum and standard deviation % of the clay silt and sand contents in the samples, and the resulting hydraulic conductivity based on the pedotransfer functions in Neurotheta.

	Clay (%)	Silt (%)	Sand (%)	Calculated Hydraulic conductivity (mm/hr)
Mean	51.5	27.2	20.9	119.9†
Maximum	77.7	58.6	78.3	2030.5
Minimunn	12.3	8.3	2.3	4.3
Standard deviation	11.7	9.7	11.7	0.66†

^† mean calculated using log-transformed values and back transformed, standard deviation based on log-transformed values

To calculate hydraulic properties, pedotransfer functions (PTFs) were applied to the data using the program Neurotheta (Minasny and McBratney, 2002). Neurotheta uses a calibrated neural network to predict parameters of the van Genuchten function to describe the relationships between water content and pressure head and hydraulic conductivity, from soil particle size distributions (Minasny and McBratney 2002). In this case the only data available were the clay, silt and sand contents (Table 1). Neurotheta is additionally able to calculate 50 samples out of the distribution of hydraulic properties associated with the specific soil properties, and this incorporates the uncertainty in estimating the average hydraulic conductivity. The values generated by Neurotheta sometimes contained unrealistic estimates for water retention parameters, and these parameters were removed prior to the simulations. The number of unrealistic estimates increased when the uncertainty in the estimation of the average hydraulic properties increased. On average more than 40 values for the hydraulic properties were generated for each soil point for Monte Carlo simulations. Irrigation was stipulated for cotton only, and was based on 6 irrigations of 100 mm (Table 2). Ponding of 10cm was only allowed at the irrigated sites. Preferential flow due to soil cracks was not employed in this study, due to lack of information on the shrinkage characteristic of the soil. Most of the other inputs followed standard suggestions for SWAP as outlined in the user manual (van Dam et al. 1997).

Table 2. Description of Crop inputs for the SWAP (Soil-Water-Atmosphere-Plant) model, indicating the number of points in the landscape modelled, the simulated maximum root depth, the crop calendar and whether irrigation was simulated.

The model was initially run for a simulation period of 11 years (1992 – 2002) with the average hydraulic properties for each location and the initial soil moisture conditions set at values found at-300cm to represent field capacity. However, the final soil moisture profile indicated that this overestimated the amount of moisture in the profiles and this increased deep drainage in the earlier years. The final moisture profile of these initial runs was therefore inserted as the initial soil moisture profile for subsequent runs. The model was first run with the average hydraulic properties to derive the estimated annual deep drainage at each point. Statistical summaries, such as means and standard deviations were calculated after discarding the first year of data (1992). Subsequently, to include the uncertainty in the estimation of the soil hydraulic properties, the model was run for all values of the hydraulic properties, resulting in more than 40 simulations per soil point. From these results, statistical summaries were calculated after discarding the first year of data (1992). The probability of a certain level of deep drainage occurring in a single year was calculated by summing the daily values for all simulations at a point and calculating the fraction of years producing deep drainage greater than either 100 or 5 mm/year. The 5 mm/year value was based on suggested deep drainage under natural vegetation (Walker et al. 2002), while the 100 mm/year can be seen as a possible environmental threshold above which DD is assumed to be environmentally harmful.

Interpolation

The annual values and calculated probabilities were extrapolated over the whole area by regression kriging using the different bands in the LANDSAT-7 image (Odeh et al. 1995). This involved using stepwise regression to develop regression equations between the different model output values and the spatial coordinates, LANDSAT-7 bands and, in our case, land use. Land use is of course not independent of the LANDSAT-7 bands, but since field based land use classification was also incorporated, it included additional information. The LANDSAT-7 bands were multiplied by the regression coefficients and variograms were fitted to the regression errors using an exponential model and a maximum range of 20 km in VESPER (Minasny et al. 2002). The regression errors were then kriged using VESPER on the 50×50 m grid and summed with the calculated values from the regressions (Odeh et al. 1995). Maps of these probabilities were generated and displayed using ARCGIS.

Results and Discussion

Land use was fairly well predicted using the supervised classification. Misclassification amounted to about 15 of the 85 locations, which were mainly related to the wheat and pasture land uses because the image was taken in January. At this time, wheat land was fallow, while unirrigated pastures were dry. This meant these land uses had similar spectral signatures. In contrast, natural vegetation and cotton had very distinct spectral signatures and were therefore fairly well predicted. There were some small areas planted under sorghum, but none of the soil points were located in these areas and this land use was thus not used in the simulations.

Table 3. Mean, standard deviation, median, maximum and minimum values from simulations using the average soil hydraulic properties and the Monte Carlo simulations. Note the increase in the coefficient of variations if the uncertainty in determining the hydraulic properties is included.

All values in mm/year	Wheat	Cotton	Pasture	Native Vegetation	Wheat	Cotton	Pasture	Native Vegetation
	Average hydraulic properties				Monte Carlo simulations
Mean deep drainage	61.1	138.1	15.2	2.6	52.6	164.5	8.3	16.6
Median deep drainage	53.4	135.0	11.6	1.0	11.84	111.28	0.01	4.31
Maximum	273.2	373.7	60.4	11.2	557.5	5324.5	387.0	447.3
Minimum	0.2	67.5	-0.3†	-0.001	-12.7	-80.0	-102.9	-96.4
Coefficient of variation (%)	132.3	32.0	119.1	150	151.7	129.7	273.4	228.3

^† negative deep drainage values indicate water uptake from the soil below 3 meters.

Deep drainage (DD) under all land uses was highly episodic and somewhat correlated to rainfall events (e.g. Asseng et al. 2001; Walker et al. 2002). Annual DD values based on simulations using the average hydraulic properties ranged from 373.7 mm/year for a cotton location to an uptake of 0.3 mm/year for a pasture site (Table 3). Average annual DD values indicated that the highest levels occurred under irrigated cotton resulting in a high probability of exceeding any of the cut-off values in any given year (Figure 3). The overall range of DD values was well within the range of values found in earlier research (i.e. Walker et al. 2002; Silburn and Montgomery 2001). In contrast, annual DD values under pasture and natural vegetation were low and had strongly skewed distributions, resulting in large coefficients of variation and differences between mean and median values. This meant that the probability of exceeding 100 mm/year or 5 mm/year DD in any give year increased in the order of pasture, natural vegetation, wheat and cotton (Figure 2). In contrast, DD under cotton always exceeded 5 mm/year and also exceeded 100mm/year in 75% of the time in the 10 years simulated.

A)

B)

Figure 2. Distribution and probability graphs for the occurrence of deep drainage under the different land uses for the simulations using average hydraulic properties (A) and the Monte Carlo simulations (B) which included the uncertainty in determining the hydraulic properties. Indicated are the distributions (blue bars), cumulative probaility of exceedance (solid line) and the 100 mm/year (dashed line) and 5 mm/year (dotted line) deep drainage values.

The variability in DD increased when the uncertainty in the determination of the hydraulic properties was included (Monte Carlo simulations, MC). Interestingly, the mean DD values for cotton and natural vegetation increased, while the mean values for wheat and pasture decreased. All distributions of DD values were highly skewed, even more so than if the average hydraulic properties were used. This was indicated by large coefficients of variation and large relative differences between mean and median. The maximum and minimum values presented in Table 3 for the MC simulations were not considered to be realistic values observed in the field below the rootzone. For example the maximum simulated DD under irrigated cotton was larger than the amount of water applied in irrigation and rainfall in any given year. However this value resulted from a build up of soil moisture in two dry years before (possibly balanced by a negative DD in a year) followed by a release of water in a wet year. In addition, Silburn and Montgomery (2001) quoted lysimeter studies in southern Queensland, which reported DD of 1010 mm/year under furrow irrigated cotton, indicating that very high values occur in some years. This extreme variability highlights the need to run long term simulations benchmarked by field studies (Silburn et al. 2004). The probability of exceeding 5 mm/year DD increased for natural vegetation, but exceedance probabilities decreased for all other land uses, due to the increase in the skew of the distribution. Since the probability of exceeding a certain DD value equated to risk, we calculated the risk of exceeding 100 mm/year DD and 5 mm/year DD. In general these results indicated that including the variability in hydraulic properties had a substantial effect on the average deep drainage simulated and the resulting DD risk. In other studies (i.e Asseng et al. 2001; Bui et al. 1996), such uncertainty has not been included, which probably led to an underestimation of the average deep drainage and an overestimation of the risk. On the other hand, this study included a much shorter temporal variability compared to the Asseng et al. (2001) study, which could have led to similar over- or underestimations.

The correlations between the DD risk based on average hydraulic properties and the ancillary landscape variables (LANDSAT, spatial coordinates and land use) were reasonably strong (Table 4). However the regressions were dominated by land use, indicating the strong influence that this factor had on the DD risk. The regressions between the predicted MC DD risks and the ancillary landscape variables all had fairly low correlations (Table 4). For both MC DD risk cut-off values, eastings, land use and the shortwave-infrared channel were significant predictors. Eastings would be an expected predictor, since many of the environmental signatures in this region (whether from soil, climate or vegetation) display an East-West trend, due to the strong East-West climate trend. The shortwave-infrared channel was strongly related to the vigour of the vegetation growth (Huete et al, 1994) and the significance of this band indicated the correlation between deep drainage and land use. The lower correlations for the MC DD risks regressions were caused by the increased variability introduced by the uncertainty in hydraulic properties, which was more related to soil properties. Introducing further variability due to rainfall and/or land use would decrease the correlation between a snapshot (the LANDSAT image) and the dynamic variables (DD risk) even more. This also highlighted the dynamic nature of DD risk in contrast to hazard.

The predicted variograms for the regression errors indicated only a slight spatial trend, meaning that there was little spatial structure in this data. This might indicate that, as DD risk was mainly correlated to land use, and land use did not exhibit a spatial trend, the DD risks also exhibited little spatial trend. The resulting risk maps clearly indicated the high risk of DD in the cotton growing locations (compare Figure 3 & Figure 1), particularly in the map based on the simulations with the average hydraulic properties. There was some difference in the magnitude of the DD risk depending on the underlying soil type in the MC DD risk map, and there appeared to be an East-West trend in the map, with the western side of the area having a lower risk. This was considered reasonable because the soils generally contained more clay from East to West. It is clear that including the uncertainty arising from the estimation of the hydraulic properties changed the risk map and de-emphasised the effect of land use. The map can mainly be used to identify areas that have a high risk of DD under the current land use, such as the dark area in the southeastern corner. Another approach would be to simulate the same land use for all locations and as such indicate the risk of DD purely based on differences in soil type (Asseng et al. 2001; Triantafilis et al. 2003). However, this would not reflect the current risk, but some hypothetical risk if all land use would have been changed to, for example, wheat production or irrigated cotton. While these types of risk analysis can be useful for assessing land use change, such information can also be gained from a hazard map (i.e. Bui et al. 1996), and the additional computational effort would not be warranted.

Table 4. Overview of the results of the regression equations to predict the probability of exceeding an environmental cut off value from ancillary information.

Variable	Significant predictors	r2
Average Hydraulic properties
Risk DD > 100 mm/year	Band 2 – Band 5 and land use	0.54 (p<0.001)
Risk DD > 5 mm/year	Band 1 & 2, Band 4 & 5 and land use	0.63 (p<0.001)
Monte Carlo simulations
Risk DD > 100 mm/year	East, Band 5 and land use	0.23 (p<0.001)
Risk DD > 5 mm/year	East, North, Band 4, 5, and 6 and land use	0.31 (p<0.001)

Figure 3. DD risk map of exceeding 100 mm/year in any given year from the simulations based on the average hydraulic properties (A) and the Monte Carlo simulations (B) Darker colours indicating higher risk and this appears aligned with cotton production in (A), but less so in (B). Legend indicates the probability of exceedance, or risk.

The risk analysis presented here was useful in that it allowed identification of current land use and landscape combinations that are under high risk of DD. These areas can be targeted for management or land use change. In addition, it is interesting to note that on a landscape scale the overall risk of DD >100 mm/year is low (in the order of 0.2), which might indicate that the overall landscape deep drainage is thus in the order of 20 mm/year. In terms of groundwater management and the risk of rising water tables, such an overall landscape value might help drive larger landscape scale models. The method developed here can thus also be used to assess catchment level DD for catchment management purposes.

Inclusion of variability of rainfall through a stochastic rainfall simulator would be the next step. The work by Asseng et al. (1997) indicates that variability due to rainfall can have a major effect on the risk of DD. Interesting would be to test whether the uncertainty due to soil properties has a greater effect on DD risk than the uncertainty due to climate. From the inherent variability in the rainfall data it would be presumed that the uncertainty due to climate would be considerable. Another aspect not incorporated in this study is dynamic land use. Land use is considered to be static and does not change from year to year, such as under rotations or in response to available irrigation water (e.g. Dudley and Hearn 1993).

Conclusions

A combination of a one dimensional mechanistic model and regression kriging proved an effective method for constructing a DD risk map for the Border Rivers area of NSW and Queensland. Including the uncertainty of estimating the soil hydraulic properties increased the variability in DD estimates, but, in general, decreased the risk of exceeding 100 mm/year DD under all land uses.

Acknowledgements

This work was performed as part of the Northern Murray Darling – Water Balance Group efforts. The authors would like to acknowledge the Australian Cotton CRC for an honours scholarship for Claire Glendenning, the Border Rivers catchment management association and Queensland department of NRME for assistance in data collection and Drs B Minasny and J Triantafilis for assistance and data during this study.

References

Asseng S, Fillery I, Dunin F, Keating B, Meinke H (2001) Potential Deep Drainage under wheat crops in a Mediterranean climate I. Temporal and spatial variability. Australian Journal of Agricultural Resources 52, 45-56.

Bui EN, Smettem KRJ, Moran CJ, Williams J (1996) Use of soil survey information to assess regional salinisation risk using Geographical Information Systems. Journal of Environmental Quality 25, 433-439.

Collins A, Donaldson S (2002) DRAFT Salinity Hazard Mapping in the Barwon Region. DLWC Tamworth, NSW.

Dudley NJ, Hearn AB (1993) Systems modelling to integrate river valley water supply and irrigation decision making under uncertainty. Agricultural systems 42, 3-23

Gilfedder M, Walker G (2001) Dryland Salinity Risk, a review of assessment methods. Natural Resource Management 4, 3-9.

Huete AR, Post DF, Jackson RD (1984) Soil spectral effects on 4-space vegetation discrimination. Remote Sensing of Environment 15, 155-165.

Jolly ID, Cook PG (2002) Time lags in salinity control: controlling recharge and halting watertable rise. Natural Resource Management 5, 16-21.

Minasny B, McBratney AB 2002. The neuro-m method for fitting neural network parametric pedotransfer functions. Soil Science Society of America Journal 66, 352-361

Minasy B, McBratney AB, Whelan BM (2002) VESPER version 1.6. Australian Centre for Precision Agriculture, McMillan Building A05, the University of Sydney, NSW 2006.

Odeh I, McBratney A, Chittleborough D (1995) Further results on prediction of soil properties from terrain attributes: heterotopic cokriging and regression-kriging. Geoderma 67,215-226

Odeh IOA, Cattle SR, McBratney, AB (2004) Soil Database Assistant Version 1.6: Summary and Users’ Guides. CDROM Publication, Australian Cotton Cooperative Research Centre, The University of Sydney, Sydney NSW.

Silburn M, Montgomery J (2001) Deep Drainage Under Irrigated Cotton In Australia – A Review (in progress). In: Cotton Consultants Assoc Meeting; 21-22 June: Dalby

Silburn M, Vervoort RW, Schick N (2004) Report on Northern Murray-Darling Water Balance Workshop 2. ‘Deep drainage – so what?’ ‘Where is it going and what is it going to do and when?’ Cotton CRC occasional paper.

Triantafilis J, Huckel A, Odeh I (2003) Field-scale assessment of deep drainage risk. Irrigation Science 21: 183-192.

Triantafilis J, Odeh I, Short M, Ahmed M (2002) Mapping Deep Drainage Risk in the Lower Gwydir and Macquarie Valleys. In: Proceedings of The Irrigation Association of Australia, National Conference and Exhibition. Darling Harbour, Sydney May 21-23.

Van Dam JC (2000) Field-scale water flow and solute transport. SWAP model concepts, parameter estimation and case studies. Koninklijke Bibliotheek: Den Haag; 167pp

Van Dam JC, Huygen J, Wesseling JG, Feddes RA, Kabat P, van Walsum PEV (1997) Simulation of transport processes in the soil-water-air-plant environment. SWAP user’s manual. DLO-Winand Staring Centre, Wageningen, The Netherlands

Vervoort RW, Silburn M, Kirby M (2004) Near surface water balance in the Northern Murray-Darling Basin. Water Science and Technology 48, 207-214

Walker GR, Zhang L (2002) Plot-scale models and their application to recharge studies. In: Volume 10 Studies in Catchment hydrology. The basis of recharge and discharge. Zhang L, Walker G (eds), CSIRO Publishing: Victoria; 29pp

Walker G, Zhang L, Ellis T, Hatton T, Petheram C (2002) Estimating impacts of changed land use on recharge: review of modelling and other approaches appropriate for management of dryland salinity. Hydrogeology Journal 10, 68–90