My Account: Log In | Join | Renew
1st Page

Journal of Environmental Quality - Article



This article in JEQ

  1. Vol. 39 No. 1, p. 26-34
    Received: Jan 27, 2009
    Published: Jan, 2010

    * Corresponding author(s):


A Methodology to Estimate the Future Extent of Dryland Salinity in the Southwest of Western Australia

  1. Peter Caccetta *a,
  2. Robert Dunnea,
  3. Richard Georgeb and
  4. Don McFarlanec
  1. a CSIRO Mathematical and Information Sciences, 65 Brockway Rd., Floreat, Western Australia
    b Dep. of Agriculture and Food, Bunbury, Western Australia
    c CSIRO Land and Water, Floreat, Western Australia


In the southwestern agricultural region of Western Australia, the clearing of the original perennial vegetation for annual vegetation-based dryland agriculture has lead to rising saline groundwater levels. This has had effects such as reduced productivity of agricultural land, death of native vegetation, reduced stream water quality and infrastructure damage. These effects have been observed at many locations within the 18 million ha of cleared land. This has lead to efforts to quantify, in a spatially explicit way, the historical and likely future extent of the area affected, with the view to informing management decisions. This study was conducted to determine whether the likely future extent of the area affected by dryland salinity could be estimated by means of developing spatially explicit maps for use in management and planning. We derived catchment-related variables from digital elevation models and perennial vegetation presence/absence maps. We then used these variables to predict the salinity hazard extent by applying a combination of decision tree classification and morphological image processing algorithms. Sufficient objective data such as groundwater depth, its rate of rise, and its concentration of dissolved salts were generally not available, so we used regional expert opinion (derived from the limited existing studies on salinity hazard extent) as training and validation data. We obtained an 87% agreement in the salinity hazard extent estimated by this method compared with the validation data, and conclude that the maps are sufficient for planning. We estimate that the salinity hazard extent is 29.7% of the agricultural land.

    DEM, digital elevation model TM, thematic mapper

In the southwest of Western Australia (Fig. 1 ), the clearing of native perennial vegetation for annual vegetation-based agriculture has lead to rising saline groundwater levels and the subsequent loss of previously productive land to dryland salinity. The 18 million ha of agricultural land is predominantly used for the growing of cereal grains including wheat (Triticum aestivum L.), barley (Hordeum vulgare L.), and lupins (Lupinus spp.), with sheep (Ovis aries) and to a lesser extent, cattle (Bos taurus) grazing on annual-based legume and grass pastures. The region has a semiarid climate comprised of a cool, wet period between the months of May and August, and a hot dry period between December and March. Average rainfall across the broadacre agricultural areas ranges from 650 mm per annum to the west and 280 mm per annum inland.

Fig. 1.

The extent of the Land Monitor project– the southwest agricultural region of Western Australia. Standard Landsat TM scene extents are shown.


It has long been documented that the clearing of native perennial vegetation within this environment causes saline groundwater levels to rise, resulting in increased concentrations of dissolved salts in runoff and the loss of previously productive agricultural land (McFarlane et al., 1992). In this paper we concentrate on the latter issue, and in particular the problem of estimating the future likely extent of salt-affected land as it relates to agricultural production. For the purposes of this paper we define a soil as saline if it contains sufficient soluble salts to limit plant growth under either waterlogged or nonwaterlogged conditions, and/or when salinity affects soil structural properties such that its productive use is detrimentally affected. We define “salinity hazard” as any area of land within wheatbelt valleys that have the potential to become salt affected (shallow watertable with rising trend) and “salinity risk” as a specific area that has a quantified consequence, likelihood, and timescale to impact (Spies and Woodgate, 2005).

Ferdowsian et al. (1996) reviewed the numerous methods that had been used in Western Australia to estimate the extent of salinity and its likely long-term impact at equilibrium. They concluded that because of the large physical extent of the region, and that each region is not uniformly affected, it is of interest to have spatially explicit estimates in the form of maps to aid regional and farm-scale planning.

The approach described in this paper used training and validation data derived from expert knowledge, variables derived from high resolution digital elevation models (DEMs), and remnant vegetation presence/absence maps, to predict the extent of future salinity. This approach has been developed progressively. Early work by Caccetta (Caccetta, 1997) derived information from relatively low resolution DEMs, summarizing the information as “landform” classes, and then demonstrated that the landform classes associated with lower relative positions in the terrain were more likely to be saline than other classes. Caccetta et al. (1995) demonstrated how catchment parameters including those relating to the amount of clearing could be derived and subsequently used to form predictions of the future extent of salinity. Evans et al. (1995) considered the use of a decision tree classifier to produce maps of future likely extent using these DEM-derived variables. However, early results were limited by the availability of only low resolution DEMs, which were derived from existing contour data having elevations sampled at intervals varying from 5 to 20 m.

Success of this and related research (Furby et al., 1995) lead to the formation of the “Land Monitor” project (Caccetta et al., 2000) to map and monitor the extent of perennial vegetation, extent of salinity, and the prediction of the future extent of salinity. Noting the limitations of existing DEMs, the project produced high resolution elevation models (raster DEMs with 10 m horizontal sampling and approximately 1 m accuracy in elevation). Greater local accuracies (as compared with expert-derived local estimates) in predictions were obtained (Evans and Caccetta, 2000), though the method did not extrapolate well as the decision tree classifier was applied on a pixel-by-pixel basis. A limitation of the approach, from a planning perspective, was that it did not allow for scenario modeling of alternate land uses and crops. Noting this, an attempt was made to estimate the likely future groundwater depth using physical-based simulation models incorporating estimates of groundwater recharge and groundwater depths as estimated from bore hole observations (Campbell et al., 2000). This proved unsuccessful for a number of physical (geological complexities at various scales) and statistical (insufficient borehole records for estimation of hydraulic head surfaces, groundwater trends and nonrandom placement of bores in the landscape) reasons.

In this paper we describe the methodology, analyses, and algorithms used to produce a broad-scale forecast of salinity hazards in the southwest agricultural region of Western Australia using ground data and DEM-derived variables.

Materials and Methods

The main physically-based assumption used in our analysis is that, locally low-lying areas are more prone to dryland salinity than higher areas. The second assumption is that groundwater levels are correlated with surface topography, such that given similar geological characteristics, the resultant landforms and mechanisms resulting in the extent of future salinity in one region will be similar to that in other regions having similar watershed morphological and land-clearing characteristics. Given these assumptions, the strategy in creating an estimate of salinity hazard was to:

  1. Stratify the region into hydrogeological zones;

  2. Acquire data including raster DEMs, remnant vegetation extent maps, geographic stratification data, and ground-truth data (described below) for training and validation;

  3. Identify large true depressions such as lakes, particularly in areas having poorly defined surface flow and false depressions (artifacts created from DEM processing). Landsat data acquired in wet seasons were used for this purpose;

  4. Derive watershed morphological variables from the DEM;

  5. Derive watershed-clearing variables from the DEM and the remnant vegetation extent map;

  6. Use the information from catchments where hydrologists have made future extent prediction (based on numeric modeling where possible) as training data to build a classifier using the variables in 2, 3, 4, and 5 as the predictors. This was a three-step process where (i) the variables were used to predict the position of the channel head defining the boundary between future saline (later termed “valley hazard” by George et al. [2005]) and nonsaline land; (ii) flow path operations were used to assign downstream flow paths of these heads to also be classed as having a salinity hazard; and (iii) the areas adjacent to, and within a relative elevation range of the salinity hazard flow paths, were assigned as having salinity hazard. Typically the relative elevation range was up to 2 m, roughly representing the minimum rooting depth of the majority of plants grown in these areas, and that these plants could potentially be accessing water tables having a salinity hazard.

  7. Validate the predictions. Where possible, the accuracy of the salinity hazard was determined by calculating the accuracy of the map as compared to ground data not used in training the classifier.

The following subsections describe each of these points in turn.


Geographic Stratification

The “Characteristic soils map of southwestern Australia” (Department of Agriculture and Food, 2002) is based on a hierarchical classification. This Soil Landscape Zone (SLZ) level information was used as the basis for stratification. From the map, a stratification of amalgamated soil map classes based on an expert hydrologist's perceived salinity hazard was conceived (Fig. 2 ).

Fig. 2.

The Soil Landscape Zones and the 11 training areas for the agricultural region of Western Australia.


Ground-truth Data

Sufficient objective data such as groundwater depth, its rate of rise, and its concentration of dissolved salts were not available (Campbell et al., 2000), so we used a combination of field survey and regional expert knowledge derived from the limited number of previous hydrogeological investigations as salinity hazard extent ground truth.

The ground truth consisted of expert hydrological opinion variously captured as: (i) estimated salinity hazard maps showing the predicted extent of “hazard” areas based on detailed hydrogeological investigations of particular regions (Type 1); and (ii) feedback on interim produced hazard maps from hydrologists expert in each region.

An example of the training data of Type 1 is depicted in Fig. 3 , and the spatial location of all the regions used for collecting ground-truth data is depicted in Fig. 2

Fig. 3.

Example of the fours areas of training data derived from Broomehill. The gray areas within the four training areas were considered by an expert hydrologist to have a lower salinity hazard on the basis of the current extent and forecast risk. Conversely, the white areas were considered to have a higher likelihood of becoming saline.


Remnant Vegetation Map

The Land Monitor remnant vegetation map from the year 1996 was used. The remnant vegetation map was derived from the analysis of time-series remotely sensed Landsat TM data.

Digital Elevation Model Data and Derived Variables

The Land Monitor DEM is a raster DEM produced from soft-photogrammetry techniques applied to 1:40,000 scanned aerial photography, producing a DEM with vertical accuracy of the order of 1 to 2 m and produced as a raster file having 10 m pixel size. The original 10 m data had some artifacts in it, and it was subsequently necessary to smooth the data. As it was desirable that the Land Monitor salinity hazard maps could be used concurrently with the salinity and vegetation monitoring results (and indeed the vegetation results are used to form some of the predictor variables), the DEM was resampled to 25 m corresponding to the chosen resolution of the resampled Landsat satellite data used in the project.

From the DEM a number of variables were derived, including: FlSlp, slope in the direction of flow; Plan, plan curvature; Profile, profile curvature; Tangent, tangent curvature; AvHgt, average upslope relative elevation; AvPlan, average upslope plan curvature; AvProf, average upslope profile curvature; AvTan, average upslope tangent curvature; AvFlSlp, average upslope slope; UpArea, upslope area; UpClr, upslope cleared area; PerUpClr, percentage upslope cleared area.

The derivation of these variables is described in Appendix A

Analysis of Data to Form Areas of Salinity Hazard

Based on the hydrologists' understanding of the processes of salinization for the region, the SLZ regions in Fig. 2 were amalgamated into the six regions shown in Fig. 4 , with each region treated independently in the derivation of the salinity hazard maps. Regions 1 and 6 have poorly-defined surface water channels and many lakes; Regions 3 and 5 have incised, well-defined flow paths with many such paths per square kilometer; and Regions 2 and 4 are characterized by broad, flat valley palaeo-drainage systems.

Fig. 4.

Six Soil Landscape Zones identified for independent analysis and salinity hazard model development.


For areas with defined surface flow features (i.e., Regions 2, 3, 4, and 5), the steps used for predicting salinity hazard areas with shallow groundwater levels were as follows: (i) the ground data, provided by the hydrologist, was digitized to produce an image indicating areas at hazard/not at hazard; and (ii) from this training data the flow paths were modeled as being at hazard or not at hazard. The modeling of the flow paths involved: (i) extracting the flow paths from the DEM; (ii) cleaning up the flow paths in the instances where they had isolated segments; and (iii) using the digitized ground data and the DEM-derived variables to produce a decision tree for each zone. The decision tree produces a hazard map showing which flow paths were predicted to have a potential of becoming saline in the future.

A decision tree methodology (Atkinson and Therneau, 1997) was used for several reasons. For any given pixel there is a high correlation between the DEM-derived variables. Decision tree predictions are resistant to problems caused by correlations between the variables and are resistant to the inclusion of additional uninformative variables. These variables are simply ignored in the model building process. In addition, decision trees are not influenced by the scale of the variables. That is, if the variables are rescaled, the model will remain the same but on a new scale.

In addition, and perhaps more importantly, a decision tree can be converted into a series of “if …, then …” statements. This allows for an automated conversion of the decision tree into a set of rules that can be easily applied using software tools designed for spatial data.

The downside of using a decision tree approach is that the model is highly sensitive to small changes in the data. That is, a small perturbation of the data may produce a decision tree that looks very different, although the accuracy may be very similar. This instability means that we cannot recommend a definitive consensus model based on DEM variables for salinity hazard assessment and prediction.

The decision tree was applied to whole stratified zones to produce maps showing the flow paths that are predicted to be areas of salinity hazard. This was used to produce an extent map by extracting all regions within four increments of 0.5 m above the nearest “at hazard'” flow path. The appropriate height was estimated using the digitized ground data to produce the most accurate hazard estimate or prediction.

We made use of the fact that all but one (Region 6) of the Soil Landscape Zones had more than one training area. The training areas were used in pairs (train on one, test on the other) to investigate the affect of choosing a particular height above stream levels to which salinity may extend as and if groundwater levels rose.

Feedback on the interim hazard maps led to the realization that an additional secondary adjustment was necessary to account for expert feedback on the relative position of hazards along a flow line. A point on a flow line high in the landscape might have a “best match” salinity rising 0.5 m above that point (i.e., most of the forecast salinity would fall inside the 0.5-m increment), whereas at a point lower in the landscape on the same flow path, the best match might be defined by a flow path rising 1 m above the flow line. This secondary adjustment was not incorporated into the decision tree model.

Instead, a notional “scaled height” was introduced to improve the predictions. Scaled height was used to calculate the hazard so that the actual height above stream lines was higher in the broad valleys and lower in the higher ground. The values of the scaling variable used for the calculation of scaled height were chosen by making use of the experts' opinions in conjunction with the resulting agreement of the prediction against the training data. In practice the appropriate height was found to range up to 2 m vertical elevation. Again this was used to produce a salinity hazard extent map by extracting, for all regions, 0.5-m increments above the nearest “at hazard” flow path. This four category (0–0.5, 0.5–1, 1–1.5, and 1.5–2 m) map is referred to as the average height above valley floor (AHAVF) map by McFarlane et al. (2004) and valley hazard by George et al. (2005)

In areas with no defined surface flow (Regions 1 and 6), the procedure was used to identify real depressions and sinks (after pit filling the DEM) and estimate the hazard as the scaled height above these depressions.

The hazard maps were validated using ground data that had not been used in deriving the models and by inspection by a hydrologist who was expert in salinity investigations in each region. In the later case, they were often compared with outputs from hydrological models.

Results and Discussion

Using a map of the training area (such as Fig. 3) and the DEM variables, a decision tree model was fitted. Such a model (Atkinson and Therneau, 1997) provided an internal cross-validated error estimate as part of the fitting process. However, in a context like this, with spatially correlated variables and the need to extrapolate models over considerable distances, an estimate generated on the training data may be optimistic.

As a test of the ability of the model to extrapolate over large distances we considered the following three areas: (i) Date Creek, (ii) Broomehill, and (iii) Tambellup. These three areas are within SLZ Region 3. The procedure was to train on one area (e.g., Date Creek) and then evaluate the model on the other areas in turn. Figures 3 and 5 show the hydrologist's hazard map for Broomehill and the Broomehill salinity hazard estimate from the Date Creek model, respectively. As Broomehill is almost 100 km from Date Creek, it would be expected that the accuracy of the model would be significantly lower at Broomehill. However, the model had an 87% accuracy at Broomehill vs. a 86% accuracy at Date Creek, that is, this model gave essentially the same accuracy far from the training data as it did on the training data. The accuracy on the Tambellup catchment was 82%.

Fig. 5.

Predicted hazard area for Broomehill using the decision tree model trained on another catchment located in a similar Soil Landscape Zone, 100 km northwest (Date Creek catchment).


It needs to be emphasized that “accuracy” in this instance means conformity with the predictions of the expert hydrologists. These predictions are themselves based on experience and extrapolation from similar conceptual and mathematical models. This is similar to “scenario modeling”, that is, given the current understanding of the hydrologists, what is the future scenario for the spread of salinity.

The result of this exercise indicated that we could use the fitted model to extrapolate over a considerable distance within the same SLZ. Using the combined images, a map of the entire agricultural area was completed (Fig. 6 ).

Fig. 6.

Derived map of salinity hazard, with mapped salinity (Furby et al., 2009) overlaid, for the southwest agricultural region and (insert) for an approximately 16 km2 subregion showing the mapping detail (near Lake Magenta). Areas mapped as salt affected in 1988 are shown in orange. Areas that became salt affected by 1998 are shown in red. The salinity hazard areas are shown in blue. The greyscale background is a January 1994 Landsat TM mosaic image. The green square in the main display shows the location of the insert region.


The methodology enabled mapping of salinity hazard for almost 18 million ha of the agricultural area. Together with the estimation of current salinity (Furby et al., 2010), the datasets provided the first comprehensive estimate of the current area of saline land and future hazard. Since their estimation, the results of the Land Monitor Project have been aggregated at Soil Landscape Zone, Shire and Catchment levels by planners and land managers.

Table 1 summarizes the current extent, recent changes, and equilibrium valley hazard aggregated over a selected number of Shires, with shires being a convenient unit for planning, though having total area greater than that used primarily for agricultural production. Primary saline areas have been removed through the Land Monitor water mask. McFarlane et al. (2004) estimated that the 1996 salinity totalled 957,581 ha or 2.9% of shire land and 5.1% of agricultural land (the former includes uncleared land). In 1989, the area was 859,306 ha or 2.6% of shire land and 4.6% of agricultural land. This represented an annual increase of about 100,000 ha over the 7-yr period, or 14,000 ha per annum. Salinity hazard for the same region was estimated to be 5,464,834 ha, or 16.8% of the shire land and 29.7% of agricultural land. We note that the definition of hazard areas may include severely and moderately salt-affected areas, as well as those areas where productivity is likely to be reduced due to salinization. The latter case is generally not able to be mapped by the satellite monitoring approach (Furby et al., 2010), and thus not included in the estimates for 1989 and 1996.

View Full Table | Close Full ViewTable 1.

Shire salinity statistics (reproduced from McFarlane et al., 2004). For the sake of brevity, this table provides summary of selected shires, with those not listed explicitly combined as the “all other” category. A detailed summary by Soil Landscape Zone can also be found in George et al. (2005)

Shire Shire area 1996 Salinity Percentage of shire with salinity 1996 1989 Salinity Percentage of shire with salinity 1989 1989 to 1996 Change as percent Salinity hazard Hazard as percentage of shire
ha % ha % ha %
Boddington 191,761 2361 1.2 1195 0.6 0.61 12,767 6.7
Broomehill 117,251 4786 4.1 3902 3.3 0.75 21,395 18.2
Dumbleyung 253,993 15,206 6.0 13,505 5.3 0.67 54,169 21.3
Gnowangerup 426,640 14,009 3.3 10,578 2.5 0.80 89,850 21.1
Jerramungup 651,858 15,494 2.4 11,820 1.8 0.56 125,571 19.3
Katanning 15,1791 11,072 7.3 9911 6.5 0.77 39,306 25.9
Kent 562,415 30,683 5.5 26,275 4.7 0.78 138,454 24.6
Kulin 471,217 22,364 4.7 18,138 3.8 0.90 84,794 18.0
Lake Grace 1,038,210 79,613 7.7 72,578 7.0 0.68 273,648 26.4
Moora 376,296 28,868 7.7 23,371 6.2 1.46 81,158 21.6
Narrogin 161,723 7996 4.9 7107 4.4 0.55 45,613 28.2
Nungarin 116,236 12,913 11.1 8956 7.7 3.40 50,012 43.0
Tambellup 143,634 12,325 8.6 11,001 7.7 0.92 49,286 34.3
Wagin 194,463 12,919 6.6 11,261 5.8 0.85 56,488 29.0
West Arthur 283,085 8317 2.9 6353 2.2 0.69 72,640 25.7
Wongan-Ballidu 336,712 36,595 10.9 34,240 10.2 0.70 94,822 28.1
Woodanilling 112,882 6256 5.5 5322 4.7 0.83 36,217 32.1
All other 27,006,203 635,804 2.4 262,016 1.0 1.4 3,541,354 13.1
Totals 32,596,370 957,581 2.9 859,306 2.6 0.30 5,464,834 16.8

The results of this work continue to be widely used in the West Australian natural resource management community. The Land Monitor “current salinity” and “valley hazard” maps are used in both planning (assessments of catchment hydrology; George et al., 2005) and public policy (e.g., Environmental Protection Authority, 2007).

A decade of use has provided feedback on the strengths and shortcomings of the approach. As such there are several recognized errors in the final results. These errors are largely associated with regional and local variation, where no attempt was made to correct them at the large scale of the project's products. The known errors are:

  1. Errors associated with stratification and associated processing: artifacts such as misfitting hazard extents at SLZ boundaries; and

  2. Errors of omission, mainly small; associated with the basic assumption (relatively low areas always most prone to salinity) and lack of appropriate spatial data (for example subsurface information), resulting in areas of salinity hazard in uplands being rejected. These omissions largely comprise of small areas (in gross area and percentage of total area) resulting from localized groundwater levels controlled by subsurface structures. This leads to an under prediction of the hazard area.

  3. Errors of commission; the mapped valley hazard includes some, often large areas of land that may have low or zero risks, such as locally elevated features (sandy lunettes with elevations within 0.5 m of saline flow paths) or ones that may be downslope of geological structures and large discharge areas (where the capillary fringe on the watertable lies sufficiently beneath the soil surface); and some soil-landform combinations that are not prone to salinity.

However at the local scale, the use of the four subclasses (AHAVF) provides some capacity for hydrologists to locally review the hazard areas and account for the over prediction noted above.

Taking these factors into account, the estimates of the extent of salinity by McFarlane et al. (2004) were updated by George et al. (2005) using the same source data. The latter authors re-assessed the areas of current salinity by adding salinity that occurred outside the DEM coverage and deleting several large errors of commission. They also revisited each hazard estimate for every Soil Landscape Zone, using the four AHAVF classes, selecting the class that best matched their revised forecast of salinity hazard. This revision took account of recent groundwater trend data and observations of the expansion of salinity between 1996 and 2005 (George et al., 2008). The resulting analyses concluded the total current salinity was 1,047,000 ha, comprising 821,000 ha of cleared farmland in 1996. The remainder was estimated to be on public land. For hazard, instead of quoting the total area as equivalent to that initially mapped, George et al. (2005) quote the revised valley hazard as likely to range between 2.8 mol L−1 ha (0–0.5 m AVAHF class) and 4.5 mol L−1 ha (0–2 m AVAHF class). They were also careful to explain that the actual area is dependant on spatial and temporal factors not able to be included in the assessment by Land Monitor (see below).

These salinity extent values are similar to farmers' observations of 933,000 ha for severely affected land (ABS 2002). However, George et al. (2005) reinforced the observations of others (Caccetta et al., 2000) that there remained some significant errors of omission (of up to 50%) in the western regions. The monitoring method (Furby et al., 2010) is known to underestimate saline areas in high rainfall areas as much of the saline land carries permanent cover where waterlogging is widespread, and overestimates salinity in some paddocks in the drier areas with consistently low productivity for reasons other than salinity.

Estimating what the equilibrium extent of salinity is not yet possible with the current methodology. Technology factors such as DEM resolution, the accuracy of the numeric predictions (used as ground truth), and others like geologic complexity, soil types, and constant changes in the water balance due to land-use and changing climate, prevent precise estimates. However, by acknowledging that the valley hazard is an estimate, and presenting it in blue (to enhance that it is a water-based hazard) and as a four-band visualization of salinity hazard, it remains a readily understood and powerful communication tool. To jointly predict location, extent, and timing, resources will be required to monitor and map the underlying drivers (rainfall, groundwater levels, and landuse) and continue acquisition of key secondary datasets. These include higher resolution DEMs, salinity extent maps, and improved rule based (spatial and temporal scales) and numerical models.

Despite these shortcomings, the Land Monitor method is the most comprehensively deployed means of mapping salinity extent and hazard at high resolutions for this region. In particular, there is now a requirement that another round of current salinity estimates be undertaken as it is 12 yr since the last analysis was done. Based on past estimates, the saline area may have increased by about 100,000 ha although dry years, in the northern parts of the agricultural region in particular, may have reduced the rate of spread. It is also time to revise the hazard maps to account for improved groundwater trend data and an increased availability of hazard estimates (e.g., new numerical models). The recent drier climate (rainfalls between 2000 and 2007 were 10–40% lower than between 1975 and 2000) has also probably reduced the salinity hazard. Quantifying the effect may aid in planning and subsequent adoption of management actions that are either inappropriate (due to diminished risk) or incapable of altering the rate of salinity spread.

Appendix A

The methods described include DEM surface derivatives estimated by the use of local neighborhood finite difference approaches, watershed algorithms to estimate upslope contributing areas, and catchment surface properties including the derivatives summarized by the upslope contributing area. The upslope contributing area may be estimated by a process that simulates water flowing across the surface.

Here we employ this approach, and given the relatively flat terrain we used the method described by Quinn et al. (1991), which distributes flow from each pixel to all its down hill neighboring pixels in proportion to the slope and cross-section to that pixel. Algorithms that distribute flow among multiple downhill neighbors are sometimes referred to as a multiple flow algorithms, as opposed to algorithms that route all the flow to the pixel with the greatest drop which are referred to as single flow algorithms. Caccetta (1997) discussed the choice of this approach over a single flow approach for this application. Other variations on weightings exist (see for example the discussion by Wilson and Gallant [2000]) though are not considered here.

Typically DEMs have numerous local minima (depressions/pits), many of which are the result of errors in the DEM production process as opposed to a few that are real (e.g., lakes and soaks). The erroneous minima cause difficulty for algorithms, which simulate water flowing across the landscape, and methods for removing them need to be employed. The first step in this process requires an expert to provide information to locate all true depressions. An algorithm that eliminates all remaining false depressions can then be applied. A (raster) DEM with all erroneous depressions filled and all pixels having a defined flow direction(s) is sometimes referred to as being hydrologically sound The process of removing spurious minima is commonly referred to as pit filling.

In the following discussion we assume that the DEM has been filled and flow directions are defined. We used the algorithm described by Soille and Gratin (1994) to achieve this. As well as a DEM with pits filled and flow directions defined, the algorithm also produces an ordering of pixel locations such that the DEM locations are ordered by location from the highest to the lowest elevations. This ordering helps calculate watershed attributes described in the following section.

We now introduce some notation to aid the description of the algorithms.

We represent the raster image of elevation values comprised of l pixels by the array H : HR 2 : H = {h i ;i = 1,..,l}, and the ordering of the pixels from greatest elevation to the lowest elevation as Δ (H) = {δ j j ∈[1,..,l); j = 1,..,l;h δj+1 h δ j }. For pixel i, we represent the set of eight adjacent pixel locations as n i = {n 1,..,n 8}, and for convenience will let n i 1={n1,..,n4 } represent the four closest neighboring pixels (those horizontally and vertically adjacent), and n i 2={n5,..,n8 } the four remaining neighboring pixels (those diagonally adjacent).

For each of the neighbors n i , we calculate the value w ij R as:w ij = (h i −h j )·t ij d ij where d ij = 1 and 2 for n i 1 and n i 2, respectively, and likewise t ij =0.5 and t ij =0.354 for n i 1 and n i 2, respectively.

We define the set of neighbors with elevations greater than i as n i + = {j;∀j such that w ij <0;jn i}.

Likewise, we define the set of neighbors with elevations less than i as.n i ={j;∀ j such that  w ij >0;j ∈ n i } We define the normalized down hill flow path weights asŵ ij = w ij j;∀j∈ni +w ij Estimates of DEM surface derivatives may be estimated using local finite difference equations. Here we estimate slope, profile curvature, plan curvature, and tangential curvature as described by Gallant and Wilson (1996) and simply refer to them as Slp, Prof, Plan, and Tan, respectively, in the following and main body of text.

We define a working array as the real valued image A : AR 2 :A = {a i ;i = 1,..,l} having the same dimensions as H

To derive UpArea, we initialize a i = 1; i = 1,..,l, then calculate

[1]as j ←|a δ j + l:∀l∈n δ j +a l ·ŵ l δ j ;j=2,..,l
The result is that UpArea = A represents an estimate, for each pixel, of the upslope watershed contributing area (in this case as measured in pixels). We also used UpArea = A to calculate the average of value A per unit contributing area, that is the quantity Ā: ĀR 2 :Ā = {ā i ;i = 1,..,l} defined as:
[2] Ā= A UpArea
Similarly, to derive Upclr, using the Land Monitor vegetation classification which maps a pixel as being woody vegetation or cleared for agriculture, we initialize for i = 1,..,l, a i = 0 or a i = 1 for land that is native vegetation or cleared for agriculture, respectively, and reperform the calculation for Eq. [1] The result is that UpClr = A represents an estimate, for each pixel, of the upslope watershed contributing area that has been cleared of perennial vegetation.

The quantity PerUpClr, is simply the per pixel ratio 100× UpClr UpArea

The quantities AvPlan, AvProf, AvTan, and AvSlp were derived by multiple applications of the algorithm, with a i; i = ,..,l initialized in each instance with the quantities Plan, Prof, Tan, and Slp, respectively, calculation of Eq. [1] performed, and then Eq. [2] calculated to get the average quantity per pixel for the upslope contributing area.

The quantity AvHgt was derived by calculating Eq. [1] with a i = h i ;i = 1,..,l, calculating Eq. [2] and then finally calculating AvHgt = ĀH

Scaled Height Above

For the quantity HgtAbv, each pixel in the image is assigned its relative elevation as compared to the nearest feature pixel. Here nearest is defined in terms of the distance as measured by overland flow.

We estimated the overland distance between two adjacent pixels l and k as:d lk * = ((h l −h k )2+d lk 2 ) Let Z:ZR 2 : Z = {z i ;i = 1,..,l} represent the array of overland flow distances, X:XR 2 :X ={x i ;i = 1,..,l} the height of feature pixels, Q : QR 2 : Q = {q i ;i = 1,..,l} the height above the nearest feature pixel, Y :Y ∈ {0,1} :Y = {y i ;i=1,..,l} an indicator variable, and S:SR 2 :S = {s i;i=1,..,l} a scaling variable. For i = 1,..l, Z, X, and Y are initialized as z i = 0, x i = h i , and y = 1 for feature pixels, and z i = ∞ (a large number > hmax), x i = ∞, and y i = 0 for nonfeature pixels. Q is initialized as q i = ∞; ∀i For i = 1,..,l, s i is a user-defined value. If ∀i, the user chooses s i = 1, then the calculations will produce the height above value. To produce the prediction image, we used different values for s i for stream segments in the upper parts of the catchment (e.g., s i = 3), as opposed to lower parts of the catchment, for example s i = 1.

Calculate, in order, the following six steps for j = l−1,..,1;(y δ j ≠ 1):

  1. z δ j ←| l:∀l∈n δ j −;zj(z l +d l δ j * )·ŵ l δ j ; (a weighted distance calculation)

  2. v l δ j = w l δ j (z l +d l δ j * ) ;∀l∈n δ j ;z l <∞ (downweight neighbor weights by overland distance)

  3. l δ j = v l δ j l:∀l∈n δ j −;zjv l δ j ; (renormalize)

  4. x δ j ←| l:∀l∈n δ j −;xl(x l + (h δ j −h l ) s δ j ) · l δ j (weighted average of downhill neighbor pixel elevations)

  5. s δ j ←| l:∀l∈n δ j −;zls l ·ŵ l δ j (weighted scaling factor)

  6. y δ j ←|1 (mark pixel as finished)


The mapping in Western Australia was performed as part of the Land Monitor project–a multi-agency state government initiative funded in part by the Natural Heritage Trust. Project partners comprised the Department of Agriculture and Food Western Australia, the West Australian Department of Land Information (now Landgate), Main Roads Western Australia, the Department of Environment and Conservation, the Department of Water and the Commonwealth Scientific and Industrial Research Organization. The results presented were produced with considerable ground-truthing effort from experts in the member agencies, in particular the Catchment Hydrology Group of the Department of Agriculture and Food.




  • All rights reserved. No part of this periodical may be reproduced or transmitted in any form or by any means, electronic or mechanical, including photocopying, recording, or any information storage and retrieval system, without permission in writing from the publisher.

Be the first to comment.

Please log in to post a comment.
*Society members, certified professionals, and authors are permitted to comment.