Downscaling Land Surface Parameters for Global Soil Erosion Estimation Using no Ancillary Data

Xiaoyang Zhang, Nick Drake and John Wainwright
Department of Geography, King's College London, Strand, London WC2R 2LS, United Kingdom
Email: Xiaoyang.Zhang@kcl.ac.uk

Abstract

The assessment of soil erosion at the hillslope or cathment scale has been successful using models based on detailed data measured at such scales. One such model is that of Thornes (1985) which has been parameterised at the plot and field scale and employs soil erodibility, overland flow, slope and vegetation cover in order to estimate rainsplash erosion. There is, however, a need to estimate soil erosion at larger scales but at the continental and global scale only coarse resolution data are available. Investigations into the sensitivity of this soil erosion model to the scale of some of its parameters show increased scale leads to decreased predicted erosion. Investigating the effect of vegetation cover pattern and spatial resolution on predicted soil erosion showed that when cover is uniformly distributed the minimum possible erosion rate is produced and predicted erosion is stable at all the scales, however, when cover is heterogeneous erosion is reduced with an increase in pixel size. Furthermore, erosion is reduced exponentially by two orders of magnitude (from 0.03 to 0.0007 mm.month-1) as the DEM from which slope is calculated is reduced in spatial resolution from 30 arc seconds to 10 degrees.

As vegetation cover and slope scale non-linearly in the soil erosion model methods that overcome this need to be developed. To do this we have taken the approach of downscaling these parameters to the plot scale (30m) based only on coarse resolution data (as this is the only data available in many areas of the globe). To overcome the reduction in slope as resolution decreases we have developed a method that uses the fractal properties of topography to estimate slope at a specified scale using a global scale DEM. To overcome the problem of a reduction in the spatial variability of the vegetation cover with decreasing resolution we use the block variance (in an 8x8 pixel window) and a Polya probability mixture mass function to estimate the subpixel frequency distribution of vegetation cover at a specified finer resolution. These methods have been validated using multiple scale datasets and produced low RMS errors at all scales investigated. Finally, monthly soil erosion was calculated at the global scale using these techniques and validated using available data.

1. Introduction

Soil erosion is recognised as a major problem arising from agricultural intensification, land degradation and possibly global climatic change. However, the extent of the problem is hard to quantify as field measurements of erosion are rare, time consuming, and are usually only acquired over restricted time and spatial scales. To date global studies of erosion has involved looking at sediment yield, questionnaire surveys or modelling. Walling and Webb (1983) generated global maps of sediment yield data, however, problems associated with coupling and the variable sediment dynamics of the fluvial system mean that this data is not directly comparable to hill slope erosion. Middleton and Thomas (1997) used questionnaire surveys to gain an insight into erosion severity but this approach only provides qualitative information and suffers from bias on the part of the people being questioned. Modelling can provide a quantitative and consistent approach to estimating erosion. Using remote sensing and GIS to parameterise such models allows them to be applied over global regional and local scales. Furthermore, prediction of the effects of climate change on erosion requires a modelling approach. However, the changing scale of model parameters can have a significant, but poorly understood, impact on the apparent variability of environmental quantities and the complex interaction between spatial scale and spatial variability provides a substantial obstacle to large scale environmental modelling.

The objective of this study is to investigate and develop methods for applying the following soil erosion developed at the plot scale to the global scale:

............. (1)

where E is erosion (mm/day), K is a soil erodibility coefficient, OF is overland flow (mm/day), S is the slope and Vc the vegetation cover (%). Since the model was developed from the results of plot scale experiments(1-30m) (Thornes 1985), it cannot be directly applied at regional or global scales. In order to overcome the scaling problem we have developed methods of downscaling topography and vegetation cover from the global scale (10 arc minutes) to the plot scale (30 m). These methods are developed and tested using multi-scale databases and then used to implement a global scale erosion model.

2. Calculation of Parameters in the Soil Erosion Model

The model was implemented on a lat/long projection for a region covering Europe, Asia, Arabia and northern Africa on a 10 arc minute grid. The parameters of the erosion model were investigated based on available data in this area.

2.1 Vegetation cover

NDVI is currently the only globally available remote sensing estimate of vegetation cover. Cover was calculated from NDVI using the following regression relationship derived by collating the literature on the relationship between AVHRR NDVI and Vc (Drake et al. 1995, Figure 1):

According to the maximum monthly global normalised difference vegetation index (NDVI) data with 10' resolution from 1985 to 1990, which is contained in the Global Ecosystems Database (GED) CD-ROM (Kineman and Ohrenschall, 1992), the monthly vegetation cover was estimated from monthly NDVI based on the regression equation.

Figure 1. Scattergram showing the relationship between percent vegetation cover and NDVI from different studies.

2.2 Slope

There are numerous global DEMs available at different scales from which slope can be calculated using a GIS (e.g. 10', 5' and 30''). Using Arc/Info software, the slope used here is calculated from the finest resolution DEM, the 30'' DEM (http://edcwww.cr.usgs.gov/landdaac/), and then the slope is averaged to 10'.

2.3 Erodibility

Erodibility K was calculated from the FAO soil map of the world by employing the relationship between soil texture and K outlined by Olsen (1981). This relationship assumes soils have a 2% organic matter content, fine granular structure, and moderate permeability. Hence the erodibility was calculated in terms of the available soil data in GED.

2.4 Overland flow

Overland flow is the main source of energy causing water soil erosion. Because of the complex interaction among precipitation, evapotranspiration, infiltration and overland flow, a number of models have been developed for computing overland flow in a range of conditions. Based on the SCS model (1972), the Carson and Kirkby (1972) model and the Philip (1957) infiltration equation, six models were developed to estimate monthly overland flow at the global scale and implemented GED data (Zhang et al. 1996). When tested against rainfall simulation measurements and large scale estimates of runoff coefficients the Carson and Kirkby (1972) model performed best. The model assumes that the overland flow occurs when the daily rainfall total exceeds a critical value which represents the soil water storage, and that the frequency distribution of rainfall in excess of the critical value is exponential. If it is assumed that the daily rainfall amounts approximate an exponential frequency distribution within each month in a year from the long term point of view, then:

where J(i) is the rain day frequency density, ri(i) is the daily rainfall (mm), N(i) is number of rain days per month, ro(i) is mean rain intensity per rain day each month (mm.day-1), P(i) is total monthly precipitation (mm), Q(i) is total initial soil water (mm), rc is water storage capacity (mm), OF(i) is monthly overland flow (mm) and i is the i-th month (from 1 to 12). Q(i) was updated each month using the Brutsaert and Stricker (1979) evapotranspiration model.

3. Downscaling Slope

The problem with calculating slope from global scale DEMs is that it is progressively underestimated at coarser scales. In order to investigate the sensitivity of the soil-erosion model to slope scales, the 30'' DEM was degraded to a set of lower resolution data. The slope values were calculated from these DEMs and then reduced to 10' by averaging. Erosion was calculated using average annual values of OF, Vc and K (also at 10' resolution). Thus the only changing parameter is the spatial resolution of the slope. Average slope is reduced from 3.93% at 30'' pixel size to 0.39% at 30' pixel size in the Eurasia and northern African. Predicted erosion is reduced exponentially by two orders of magnitude (from 0.03 to 0.0007 mm.month-1) as slope is reduced by these scaling effects (Zhang et al. 1997) (Figure 2).

To overcome this problem we developed a way of calculating slope independently of scale using a fractal approach. Most landscapes can be considered to be fractal and it can be shown that the percentage slope (S) is related any scale d by the equation (Zhang et al, submitted):

.......... (6)

where A is the fractal constant and D the fractal dimension of topography. To test this method we calculated the local (3x3 window) ( and D from the 1 km DEM and then calculated slope for a d of 200 m and 30 m, compared them to two DEMs of this resolution and found RMS errors of 3.92% at 30m and 8.63% at 200m. Slope has good statistical scaling properties in most areas although scaled slopes are affected where elevations change rapidly. This occurs in about 4% to 6% of the area when the pixel size of the DEMs is reduced by more than 32 times (Zhang et al, submitted).

4. Downscaling Vegetation Cover

It is therefore necessary to consider the effects of variations in the fine scale spatial distribution of vegetation cover for accurate erosion modelling at regional and global scales. Methods of increasing the spatial resolution of remotly sensed imagery are one approach (Atkinson 1997, Ramstein and Raffey 1989). These techniques generally employ some form of pixel splitting by assuming or measuring the spatial dependence at a finer scale. Such methods are hard to apply at the global scale because of the large variations in spatial dependence from region to region. To bypass this problem we developed an different approach that uses the frequency distribution of vegetation cover within a large pixel to describe the spatial variation of vegetation at a finer scale.

4.1 Developing a Multi-Scale Database of Vegetation Cover

In order to investigate the change in the frequency distribution of vegetation cover with scale and to develop methods for predicting these changes we have developed a multi-scale database of vegetation cover images which were derived from ground measurements of cover and mixture modelling of aerial photography (0.55m), TM (30m) and AVHRR (1.1 km) imagery in the Guadelintin Basin, southern Spain acquired between the 5th and 25th of April 1996.

We have used both NDVI and linear mixture modelling to estimate vegetation cover. When mixture modelling was applied TM and AVHRR imagery pure pixels of some endmembers did not exist. To overcome this the radiance of endmembers was determined by using the inverse Linear inequality constrained mixture modelling (LICMM) to solve for the endmembers R such that:

where F a vector of known proportions derived from ground measurements or estimates from higher resolution imagery than that being mixture modellied, Y is a pixel spectra, R is a matrix of unknown endmember spectra, and V and H are constraints provided by a knowledge of the spectral reflectance properties of the materials that comprise by the image. For example the spectral reflectance of one material might be known to be larger or less than that of another for a specified wavelength.

Once the endmembers are known the linear inequality constrained mixture modelling (LICMM) was used to calculate the proportions by minimising the error on:

Correlations between ground measurements of vegetation cover and the air photo estimates were high (r square =0.86), as were correlations between the air photos and TM (r square =0.9) and TM and AVHRR (r square =0.79).

4.2 The Poyla Mass Function

Analysis of this multi-scale data set showed that the histogram of the very high resolution (0.55m) vegetation cover images goes through large but predictable changes as their spatial resolution is decreased and can be simulated quite effectively using a Poyla mass function to describe the general shape of the vegetation frequency distribution at a specified scale (Figure 4).

The Polya distribution is a mixture distribution of a beta and binomial distribution and is represented as:

where n is the number of events, alpha and beta are the parameters defined by variance and mean value, x is a random variable between 0 and n. The frequency distribution shown in Figure 3A is not regular since there are peaks due to the large amount of bare soil pixels. When the simulated distribution was compared with that measured from the image the statistics are that r-square is 0.83, RMS error is 0.0027 and chi-square is 2.38(103 (Figure 4A). The statistics of the regression and the RMS error indicates that the simulated proportion is matched very well with the measured values, however, the chi-square value reveals the big bias associated with the sharp peaks in the histogram at low vegetation cover. When the image pixel size was enlarged 10 times the frequency distribution of vegetation cover changed considerably, notably in that the pixels with bare soil disappear. The predicted distribution matches the real distribution quite well (Figure 4B) and the statistics are improved since r-square is 0.85, RMS error is 0.0046 and chi-square is 82.9.

By analysing the other images we have shown that the Polya mass function can successfully simulate the frequency distribution of numerous spatial distributions and scales of vegetation cover, the relationship between measured and predicted frequency of vegetation cover is significant and the regression slope is around 1. However, the chi-square is sometimes very greater than the critical value because of peaks in the histogram at the extreme vegetation cover levels are not predicted very well. These groups decrease with reduced spatial resolution but the fit is not always better at lower resolution.

4.3 Scaling the Polya mass function

The results of changing the spatial resolution also show that though the mean is stable across scales the variance is reduced at increasingly smaller scales. Therefore some methods of predicting this reduction in variance are needed in order to employ the Poyla function to predict the frequency distribution at the fine scales from coarse resolution data. The subimage variance, which was estimated by dividing the air photo estimates of vegetation cover into 119 subimages with 7x10 pixels, shows a consistant pattern in variance reduction with increasing pixel size (Figure 5). The curve shapes in different subareas are similar and are a function of resolution. To assess the relationship between variance and spatial resolution, the numerical models outlined in Table 1 were applied. The correlation coefficients between variance and scale in all the subimages show that Var=a2+b2ln(d)is the best function to describe the decline of subimage variance with increasing spatial resolution, however, the subimage variance in every subimage is significantly correlated with the pixel size in all four models (Table 1).

Table 1. Correlation coefficients between subimage variance and cell size in all the subareas.Var represents the subimage variance, d is the pixel size, and a and b are constants.

Thus subimage variance can be predicted across scales using a numerical model but this requires a reduction in spatial resolution so that a and b can be calculated for each subimage. If a and b can be calculated using a small subimage size then the less spatial resolution will be lost. The minimum subimage size was investigated by gradually reducing it, predicting the variance, and comparing it to the actual variance. It was found that it is just about possible to estimate the variance at a specified scale d using a 5x5 subimage but the results using a 8x8 subimage are statistically better (Figure 6).

4.4 Testing the Polya function at the different scales

In order to test how the technique works at progressively coarser scales, all images were divided into four validation regions where the relatively coarser image is composed of 8x8 pixels. Then the frequency distribution of vegetation cover at 0.55 m resolution derived from TM was validated using the air photos and the AVHRR validated using TM. The results of the former comparison shows that large areas of bare soil cause problems as do those histograms with multiple modes, however, there are significant relationships between predicted and measured vegetation frequency in all four subareas (Table 2A). This is also the case when AVHRR was used to predict the frequency distribution at 30 m and was validated by TM (Table 2B). Another possible source of error is the different average vegetation cover estimates derived from the different sensors.

Table 2. Statistics of vegetation cover frequency at: A) 0.55 m scale between that predicted from TM and measured from aerial photography, B) 30 m scale prediction from AVHRR and measurement from TM.

5. Results

For calculating soil erosion at the global scale accurately, the Thornes model need scaling up. Therefore, the slope S is scaled down to a plot scale using the equation (6). According to the subpixel frequency distribution of vegetation cover scaled by using Polya function, the vegetation cover is divided into 101 levels from 0% to 100%. The frequency is expressed as fi in each level. Based on the partitioning method the erosion model was scaled up to:

When the model was run in its unscaled form (equation 1) in qualitative terms the spatial distribution of erosion seems reasonable, however, the rates of erosion are very low in most areas (Table 3). After running the scaled model (equation 10), the monthly soil erosion rate is much higher though this is obscured somewhat by the loss in spatial resolution caused by the Polya function (Table 3, Figure 7). The spatial pattern of soil erosion suggests that high soil erosion can be accounted for in terms of the steep unstable terrain, the high erodible soil, high monthly precipitation, and vegetation removal by human activity or seasonal factors. Conversely, the low soil erosion reflects the lower relief, the greater density of the vegetation canopy, and the areas of low precipitation. There is lot of erosion in southern east Asia, India, along Himalaya and Alps, and west Africa. On the contrary, there is no water erosion in desert areas. The seasonal distribution is that erosion rate is high from Spring to Autumn in southern east China, Summer and Autumn in India, from April to October in west Africa. Relatively high erosion rates occur in western Europe during the winter months while very low soil loss occurs in summer.

Table 3. Annual soil erosion (mm/y) in typical environments from models. E1 is soil erosion with no scaling; E2 scaling all the parameters.

Figure 7. Annual soil erosion (mm.y-1) by adding monthly soil loss derived from equation (10)

It is very difficult to validate the soil erosion rate at the global scale for a number of reasons. First, each pixel size is about 253 km2 (corresponding to the 10' pixel size) in this research, and few field measurements have been carried out in such a large scale. Secondly, the time frame field measurments are often short and do not correspond to that of the model. Thirdly, measuring soil erosion rates are fraught with methodological and practical problems and different techniques of measurement of the same erosion process give quite different results (Stocking 1996). Bearing in mind these problems some field measurements of soil-erosion rates are listed in Table 4. Though there are some discrepancies the measured soil erosion rate compare quite well with the predicted values. In east China the range from 1.48 to 7.4 mm.y-1 is within the same order of magnitude of predicted values. In west African Savannah the range between 0.37 and 3.7 mm.y-1 of measured soil erosion is very similar to the calculated value of 4.5 mm.y-1. For England and most of west Europe the values of estimated soil erosion rate are acceptable although the measured erosion rate is higher than calculated values in East Germany and Belgium (Table 4).

6. Conclusions

The non-linear scaling of soil erosion models restricts the applicability of models based on plot-scale measurements at the regional and global scales. In order to estimate the soil erosion at the global scale using a plot-scale model, it is important to calculate the parameter values at such a scale based on coarse resolution data which are the only data available at the global scale.

In order to improve the slope calculation, the fractal properties of topography, which is independent of scale, is incorporated within the slope estimates. Fractal properties of topography are analysed in relation to various slopes from catchments to global scale and a method developed to predict slope at a specified scale. Slope has good statistical scaling properties in most areas although scaled slopes are affected where elevations change rapidly. This occurs in about 4% to 6% of the area when the pixel size of the DEMs is reduced by more than 32 times.

Predicted soil erosion varies greatly with the different patterns and measurement scales of vegetation cover. Usually, soil erosion will be underestimated using coarse resolution vegetation cover data. Therefore, a mixture mass function, the Polya function, is used to simulate the subpixel frequency distribution of vegetation cover at various scales based on the knowledge of the mean value and variance in a given area. It can successfully estimate the distribution of vegetation cover in different environments.

The model was run with and without these scaling methods. The soil-erosion patterns at the global are very similar whether the parameters are scaled or not to integrate with GIS, whereas the most reasonable results are obtained when all the parameters are scaled. The results indicate that the soil erosion rate is largest in the humid subtropics and some tropical Savanna area with erosion rate mainly ranging from 1 to 15 mm.y-1. For most of semiarid areas the values usually range between 0.01 and 1 mm.y-1. Erosion is very low in the tropical rainforest and desert.

References

Atkinson, P.M. 1997. Mapping sub-pixel boundaries from remotely sensed imagines. In Zarine Kemp (eds.), Innovations in GIS 4, Taylor & Francis, London, p166-180.

Boardman, J. and Favis-Mortlock, D.T. 1993. Simple methods of characterising erosive rainfall with reference to the South Downs, southern English. In: Wicherck, S. (eds), Farm Land Erosion in Temperate Plains Environment and Hills, Elsevier, Amsterdam, pp 17-29.

Brutsaert, W. and Stricker, H. 1979. An advection-aridity approach to estimate actual regional evapotranspiration. Water Resour. Res., 15(2): 443-450.

Carson, M.A. and Kirkby, M.J. 1972. Hillslope form and process. Cambridge university press.

Craswell, E.T. 1993. The management of world soil resources for sustainable agricultural production. Pimentel, D. (eds) World Soil Erosion and Conservation. Cambridge University Press, Cambridge. pp 257-276.

Drake, N.A., Zhang, X., Berkhout, E., Bonifacio, R., Grimes, D., Wainwright, J. and Mulligan, M. 1997 (accepted). Modelling soil erosion at global and regional scales using remote sensing and GIS techniques. In: P.Atkinson (ed.), Spatial Analysis for Remote Sensing and GIS, Wiley, Chichester.

Evans, R. 1992. Assessing erosion in England and Wales. In: Proceedings 7th ISCO Conference. Sydney, vol. 1 pp 82-91.

Evans, R. 1993. Extent, frequency and rates of rilling of arable land in localities in Eland and Wales. In: Wicherck, S. (eds), Farm Land Erosion in Temperate Plains Environment and Hills, Elsevier, Amsterdam, pp 177-190.

Hurni, H. 1985. Soil Conservation Manual for Ethiopia, first draft. Ministry of agriculture, natural resources conservation and development main department, community forests and soil conservation development department. Addis Abeba, Ethiopia, p86.

Kineman, J.J. and Ohrenschall, M.A. 1992. Global Ecosystems Database, version 1.0 (on CD-ROM), NOAA, Boulder, Colorado.

Kosmas, C. et al. 1997. The effect of land use on runoff and soil erosion rates under Mediterranean conditions. Catena, 29:45-59.

Lal, R. 1993. Soil erosion and conservation in West Africa. In: Pimentel, D. (ed) World Soil Erosion and Conservation. Cambridge University Press, Cambridge, pp.7-25. Mesfin Wolde-Mariam. 1972. An Introductory Geography of Ethiopia. Berhanena Selam H.S.I. Printing Press, Addis Abeba, Ethiopia.

Olsen, G.W. 1981. Soil And The Environment: A guide to soil surveys and their application, Champman and Hall Ltd, London.

Pearce, E.A. and Smith, C.G. 1993. The World Weather Guide (third edition), Helicon, Oxford.

Philip, J.R. 1957. The theory of infiltration: 4. sorptivity and algebraic infiltration equations. Soil Sci. 84:257-267.

Pimentel, D. 1993. World Soil Erosion and Conservation. Cambridge University Press, Cambridge.

Ramstein G. and Raffy M. 1989. Analysis of the structure of radiometric remotely-sensed images. International Journal of Remote Sensing, 10(6):1049-1073.

Ryszkowski, L. 1993. Soil erosion and conservation in Poland. Pimentel, D. (eds) World Soil Erosion and Conservation. Cambridge University Press, Cambridge. pp 217-232.

Soil Conservation Service 1972. SCS National Engineering Handbook, Sec. 4, Hydrology, USDA. .

Stocking, M.A. 1996. Soil erosion. In Adams, W.M., Goudie, A.S. and Orme, A.R. (ed.) The Physical Geography of Africa, Oxford University Press, Oxford, pp326-341. Thornes, J.B. 1985. The ecology of erosion, Geography, 70: 222-234.

Walling, D.E. and Webb, B.W. 1983. Patterns of sediments yield. In Gregory, K.J. (ed.) Background to palaeohydrology, Wiley, Chichester, p69-100.

Wen, D. 1993. Soil erosion and conservation in China. Pimentel, D. (eds) World Soil Erosion and Conservation. Cambridge University Press, Cambridge. pp 63-86.

Zhang, X., Drake, N.A., Wainwright, J. and Mulligan, M. 1997 (submitted). Comparison of slope estimates from low resolution DEMs: scaling issues and a fractal method for their solution. Earth Surface Processes and Landforms.

Zhang, X., Drake, N.A., Wainwright, J. and Mulligan, M. 1997. Global Scale overland flow and soil erosion modelling using remote sensing and GIS techniques: model implementation and scaling. In Griffiths and Pearson (eds.), RSS97 Observations & Interactions- 23rd Annual Conference and Exhibition of the Remote Sensing Society, The Remote Sensing Society, Nottingham, UK. P379-384.

Zhang, X., Drake, N.A., Wainwright, J. and Mulligan, M.1996. Comparison of GIS based models for estimating monthly land surface fluxes at the global scale. In Proceedings of 1st International Conference on GeoComputation, Vol. II, 17-19 Septmenber 1996, University of Leeds, UK.