CC BY
299
УДК 631.4
DOI: 10.19047/0136-1694-2019-100-36-52 Cite this article as:
Papaiordanidis S., Gitas I.Z., Katagis T., Soil erosion prediction using the Revised Universal Soil Loss Equation (RUSLE) in Google Earth Engine (GEE) cloud-based platform, Dokuchaev Soil Bulletin, 2019, V. 100, pp. 36-52, DOI: 10.19047/0136-1694-2019-100-36-52
Soil erosion prediction using the Revised Universal Soil Loss Equation (RUSLE) in Google Earth
Engine (GEE) cloud-based platform © 2019 г. S. Papaiordanidis*, I. Z. Gitas**, T. Katagis***
Laboratory of Forest Management and Remote Sensing, Aristotle University of Thessaloniki, B Bld. 59 Mouschounti Str., Foinikas, 55134 Thessaloniki, Greece, https://orcid.org/0000-0002-3925-2550, e-mail: [email protected], https://orcid.org/0000-0003-0056-5629, e-mail: [email protected], https://orcid.org/0000-0002-1322-7699, e-mail: [email protected], Received 17.10.2019, Revised 18.10.2019, Accepted 18.12.2019
Abstract: High-quality soils are an important resource affecting the quality of life of human societies, as well as terrestrial ecosystems in general. Thus, soil erosion and soil loss are a serious issue that should be managed, in order to conserve both artificial and natural ecosystems. Predicting soil erosion has been a challenge for many years. Traditional field measurements are accurate, but they cannot be applied to large areas easily because of their high cost in time and resources. The last decade, satellite remote sensing and predictive models have been widely used by scientists to predict soil erosion in large areas with cost-efficient methods and techniques. One of those techniques is the Revised Universal Soil Loss Equation (RUSLE). RUSLE uses satellite imagery, as well as precipitation and soil data from other sources to predict the soil erosion per hectare in tons, in a given instant of time. Data acquisition for these data-demanding methods has always been a problem, especially for scientists working with large and diverse datasets. Newly emerged online technologies like Google Earth Engine (GEE) have given access to petabytes of data on demand, alongside high processing power to process them. In this paper we investigated seasonal spatiotemporal changes of soil erosion with the use of RUSLE implemented within GEE, for Pindos mountain range in Greece. In addition, we estimated the correlation between the seasonal
components of RUSLE (precipitation and vegetation) and mean RUSLE values.
Keywords: soil erosion prediction, RUSLE, Google Earth Engine, Pindos mountain range.
INTRODUCTION
Soil erosion is one of the many natural processes that take place in ecosystems; however, accelerated soil erosion has a negative impact on agriculture and silviculture, hydrological systems, land degradation, and loss of non-renewable soil resources (Lal, 1998; Morgan, 2009). In this scope, the estimation of soil erosion, as well as the temporal and spatial distribution of the process, is of great importance in order to prevent soil degradation and sustain high-quality soils. Estimating soil erosion, however, is a difficult task due the many impacting factors, such as climate, terrain, soil, vegetation, and land cover (Lu et al., 2004).
Traditional field measurements of soil erosion, despite being accurate and reliable, are very expensive and time consuming (Castillo et al., 2012), thus many scientists turned to predictive models that use satellite data to calculate soil erosion (Wischmeier and Smith, 1978; Lane et al., 2003; Pandey et al., 2007; Rahman et al., 2009). One of these methods is the Universal Soil Loss Equation (USLE), and its descendant, the Revised Universal Soil Loss Equation (RUSLE) (Renard et al., 1991). This equation has been widely used to predict soil erosion in many different ecosystems (Millward and Mersey, 1999; Angima et al., 2003; Fernandez et al., 2003). RUSLE uses multispectral satellite images, as well as satellite-acquire elevation models of the terrain, along with precipitation and soil data (Renard et al., 1997).
Gathering and preprocessing these large amounts of data can often become difficult, and quite time consuming. Emerging and modern technologies offer now new possibilities regarding data processing. More specifically, cloud-based services are widely used by scientists to acquire, analyze and process satellite data on the fly. Such popular platforms are AppEEARS (Application for Extracting and Exploring Analysis Ready Samples), GFW (Global Forest Watch) and GEE (Google Earth Engine). GEE can be accessed either from its online integrated development environment (IDE), or using the Application Program
Interface (API) that is being provided. GEE grants access to high performance processing power using cloud-based technologies, as well as; to very large amounts of data which are stored in cloud-based databases. The data in GEE's databases come from many different sources, including satellite images, geospatial datasets, meteorological data, land cover/land use maps, topographic data, and even social and economic data (Gorelick et al. 2017).
The main goal of this study was to predict soil erosion in the mountain range of Pindos during a full seasonal cycle. The specific objectives are:
• To estimate the monthly soil erosion by employing the Revised Universal Soil Loss Equation (RUSLE) with the Google Earth Engine (GEE) cloud-based platform.
• To assess the correlation between precipitation and soil erosion predictions.
• To assess the correlation between vegetation (both type and density) and soil erosion predictions.
MATERIALS AND METHODS
The study area is located at the centermost part of Greece, along the mountain range of Pindus, extending from 38°49'51.16" to 39°41'23.04" North and from 21°3'26.74" to 22°14'34.50" East (Fig. 1). The altitude ranges from sea level to 2300 meters, with slopes up to 45 degrees, and the surface area is 12,431.25 km2. Concerning the climate of the study area, temperature varies with elevation, with mean monthly temperatures ranging from 0.9 °C to 21.4 °C, and annual rainfall from 1,000 to 1,800 mm. The vegetation consists of black pine (Pinus nigra) and common beech (Fagus sylvatica) in the middle altitudes (1,000 m - 1,600 m), with Balkan pines (Pinus heldreichii Christ.) covering the higher altitudes (> 1,600 m) (Touchan et al. 2012).
Data. The satellite data that we used in this study were four images that were derived from the Sentinel-2 MultiSpectral Instrument (MSI), and were acquired in January 29, April 9, July 3, and October 26, all during 2018. Sentinel-2 is a high-spatial resolution (10 m), multi-spectral constellation, used for monitoring of vegetation, soil and water cover.
Fig. 1. Study area; mountain range of Pindos, Greece.
The data were level 1 products, which means that they were already geometrically and radiometrically corrected (Drusch et al., 2012). The four images that were used represent the seasons in order to assess how soil erosion changes throughout the year. A map of broad-leaf and coniferous trees was used as well, provided by the Copernicus Land service. The precipitation data that were used were in the form of R factor maps, acquired from European Soil Data Centre (ESDAC). The R factor, also known as Rainfall Erosivity factor, is the average (in the present case) monthly sum of the kinetic energy products of each storm (Renard et al., 1997a). The R factor maps were created using the best available datasets in Europe, namely the Rainfall Erosivity Database on the European Scale (REDES). The R-factor values were normalized to temporal resolutions of 30 minutes using linear regression (Panagos et al., 2015a; Ballabio et al., 2017). Similarly, LS factor maps where obtained from ESDAC. The LS factor, or the combined slope length and slope angle factor, has the greatest influence on soil erosion at the European scale and describes the effect of topography on soil erosion (Panagos et al., 2015b). The LS-calculation was performed using the original equation proposed by Desmet and Govers (1996) and
implemented using the System for Automated Geoscientific Analyses (SAGA), which incorporates a multiple flow algorithm and contributes to a precise estimation of flow accumulation (Desmet and Govers, 1996). The DEM used to calculate the LS factor map is a new highresolution (25 m) DEM of the European Union (EU-DEM) which is a hybrid created from Shuttle Radar Topography Mission (SRTM) and Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Model (GDEM). The EU-DEM is a product of the Copernicus program, and it is statistically validated with an overall vertical accuracy of 2.9 m root mean square error. Pindos mountain range has one of the greatest LS factor values in Europe (Panagos et al., 2012). Finally, a K factor map was acquired from ESDAC as well. The K-factor expresses the susceptibility of a soil to erosion, and it is depended upon some specific parameters such as organic matter content, soil texture, soil structure and permeability (Renard et al., 1997b). The K factor map was created by field measurements recorded during the Land Use/Cover Area frame Survey (LUCAS) soil survey in 2009. The equation used to create the K factor map includes five soil parameters (texture, organic matter, coarse fragments, structure, and permeability) (Wischmeier and Smith, 1978) and can be seen bellow:
K = [(2.1 x 104M114 (12 - OM) +3.25 (s - 2) + 2.5 (p - 3))/100] x 0.1317 (1).
Where:
M - Textural factor with M = (msilt + mvfs) • (100 - mc); mc (%) - Clay fraction content (< 0.002 mm); msilt (%) - Silt fraction content (0.002-0.05 mm); mvfs (%) - Very fine sand fraction content (0.05-0.1 mm); OM (%) - Organic matter content;
s - Soil structure class (s = 1: very fine granular, s = 2: fine granular, s = 3, medium or coarse granular, s = 4: blocky, platy or massive);
p - Permeability class (p = 1: very rapid, ..., p = 6: very slow).
Soil erosion modeling. In this paper the RUSLE equation was used to model soil erosion. RUSLE predicts soil erosion by using the following formula:
A = R x K x LS x C x P (2).
Where:
A - Computed spatial and temporal average soil loss per unit area, expressed in the time unit of R;
R - Rainfall erosivity factor;
K - Soil erodibility factor;
LS - Combined effect of slope length (L) and slope steepness
factor (S);
C - Cover-management factor;
P - Practice factor.
The maps representing the R, K, and LS factors were uploaded as assets in Google Earth Engine. The R factor maps included four images, one for each of the months January, April, July, and October. The Sentinel-2 collection was filtered using GEE's Integrated Development Environment (IDE) in JavaScript programming language, in order to acquire the least cloudy image for each month. The NDVI index was then calculated, and the C factor was computed using the following formula (Van der Knijff et al., 1999) (3):
C = exp[-a( NDVI )] (3).
Where a equals 2 and P equals 1, which are unitless parameters that determine the shape of the curve relating to NDVI and the C-factor (Mallinis et al., 2016).
Since the R factor data had the lowest spatial resolution, all the other maps were scaled down to fit the 500 m resolution, along with the DEM that was later used to calculate slopes in the study area. The P factor represents the effects of practices, such as direction of tilling in fields, on the reduction of soil erosion. Since no such data are available and they cannot be derived from satellite imagery, the P factor was set to 1. The final erosion maps were downloaded from Google Earth Engine along with the slope map. The potential correlations between RUSLE variable and variables such as vegetation and precipitation were investigated using the Pearson's correlation coefficient. Pearson's
correlation coefficient (r) is a measure of the linear association of two variables (Doyle 2011).
1
0.8 0.6 0.4 0.2 0
0 0.2 0.4 0.6 0.8 1
NDVI
Fig. 2: The curve that represents the relation between C factor and NDVI.
RESULTS AND DISCUSSION
It should be mentioned that the use of the GEE platform reduced by a lot the time required to access, download, and preprocess images. In case of data unavailability in GEE's databases, uploading data in the cloud and combining them with already available data can be a very user-friendly and time-efficient process. Since GEE's IDE allows the user to write their own programming scripts, automation of processes can be performed, thus reducing computation time and effort.
Initially, the erosion maps that were produced had a pixel size of 500 meters (0.25 km2), and in each of those pixels a value was assigned signifying the predicted erosion in that area in tons per hectare per month. The predicted erosion values were then classified into five classes and were represented in the table (Table 1) and corresponding the-
matic maps (Fig. 3) for each of the four months. As shown in the results below, RUSLE appears to vary through the year. This can be expected since two of the factors (R and C factors) have a highly seasonal character (Fig. 2). It appears the months with the highest erosion predictions are October and January with corresponding mean erosions 3.25 and 2.88 tons per hectare per month and total erosion of the study area 162,027 and 143,010 tons per hectare per month (Table 1, Fig. 3-5).
Table 1. Mean RUSLE and the sum for each month
RUSLE (ton/ha/month)
October January April July
< 5 49,815 36,221 32,767 11,804
5-15 60,876 45,213 19,942 3,992
15-25 26,809 25,298 2,999 514
25-35 10,852 16,449 820 82
> 35 13,675 19,829 974 82
Mean 3.25 2.88 1.16 0.33
Sum 162,027 143,010 57,502 16,474
Fig. 3. Values of RUSLE for the months October and January classified in 5 classes.
Ш '•'"
N
тйг. L - .t, ; - ,
% vNV '
Fig. 4. Values of RUSLE for the months April and July classified in 5 classes.
Sum RUSLE (ton/hectare/month)
180 000 160 000 140 000 120 000 100 000 80 000 60 000 40 000 20 000 0
October January April July
Fig. 5. Sum of RUSLE at different seasons of the year.
The correlation between the seasonal factors mentioned above (R and C factors) and the erosion predictions was investigated, and the following results were produced (Tables 2 and 3).
Table 2. RUSLE values at different precipitation classes
RUSLE (Mean) (ton/ha/ month)
Precipitation (mm/month) Autumn Winter Spring Summer
0 - 30 1.2699 Q.1Q11 1.4913 Q.3Q51
30 - 80 1.7397 Q.6785 1.2634 Q.68Q1
80 -130 5.5519 5.7765 2.9181 1.8491
130 -180 4.9851 4.3Q97 6.2453 -
> 200 1Q.Q173 4.48Q5 - -
Pearson's correlation coefficient 0.1981 0.2943 0.1993 0.2323
It appears that the higher the precipitation, the higher the predicted erosion for that time, almost in a linear manner. Although, the Pearson's correlation coefficient was calculated and found to be lower than expected (r: 0.1981-0.2943), which indicates a weak correlation between precipitation and the value of RUSLE.
Table 3. RUSLE values at different NDVI classes
RUSLE (Mean) (ton/ha/month)
NDVI Autumn Winter Spring Summer
<0 10.1027 23.1567 6.6812 1.0206
0 - 0.3 5.1609 12.4148 3.9266 1.7364
0.3 - 0.6 4.1757 2.6822 1.5712 0.8577
>0.6 1.1864 0.4635 0.4635 0.2032
Pearson's correlation coefficient -0.2229 -0.4816 -0.4565 -0.32984
On the other hand, NDVI and RUSLE values seem to have a negative correlation, meaning the higher the NDVI values, the lower the RUSLE values. This finding agrees with other studies that have presented a negative correlation between vegetation and soil erosion exists (Mohammad and Adam 2010). The Pearson's correlation coefficient ranges from -0.22 (weak correlation) to -0.48 (medium correlation). There is a case to be made concerning the species of vegetation of the study area. As mentioned earlier, one of the dominant species of the lower altitude areas is the common beech, which is a deciduous tree. The NDVI index is calculated based on the spectral values of foliage and their property of reflecting near-infrared light while absorbing red light. Since part of the trees are deciduous, lower NDVI values are expected at the winter and autumn seasons. From the NDVI index the C factor map was produced, and it was supposed to represent the ability of an area to resist erosion based on "how much" healthy vegetation was present. The visual difference in NDVI values can be seen in figure 6 (Fig. 6).
Fig. 6. NDVI index in four seasons.
As other studies have shown, the ability of trees to reduce erosion risk is not due to their foliage only, but also due to their root system and the leaf litter on the ground (Wilcox et al. 2003). That suggests that the decrease in NDVI values and therefore the decrease in C factor values are not proportionate to the decrease of erosion resistance in reality. To test the impact of deciduous trees on the predictions of soil erosion by RUSLE we investigated the mean soil erosion prediction on areas with broadleaf (Fagus sylvatica) and coniferous trees across the year (Table 4, Fig. 5).
No significant differences in mean soil erosion predictions were observed according to what type of vegetation was present (broadleaf or coniferous). This in turn suggests that deciduous trees don't cause RUSLE to overestimate.
Fig. 7. Broadleaved and Coniferous forests (provided from Copernicus Land Monitoring Service) on a Sentintel-2 NDVI image.
Table 4. Mean RUSLE values at broadleaf and coniferous forest areas across the year
Mean RUSLE/Species Broadleaved Coniferous
October 3.21 3.48
January 2.78 2.36
April 1.57 1.27
July 0.36 0.39
CONCLUSIONS
The main goal of this study was to predict soil erosion in the mountain range of Pindos during a full seasonal cycle. Resources from GEE and ESDAC were used to produce the predicted soil erosion maps by calculating RUSLE for January, April, July, and October 2018. After comparing the different results both by their sum and mean erosion per hectare per month, the seasonal factors of RUSLE were investigated. From the analysis the following can be concluded:
• the employment of RUSLE in GEE results in the production of soil erosion prediction maps in a time-efficient manner within in a user-friendly environment,
• the RUSLE model predicted higher erosion risk in October and January and lower in April and July,
• erosion values have a low positive correlation with precipitation and medium negative correlation with NDVI values,
• the type of vegetation (deciduous or evergreen) did not cause mean RUSLE values to vary in different seasons.
It would be interesting to investigate in future studies, the predicted erosion using RUSLE in full year cycles for many consecutive years, using GEE cloud-based platform and Landsat satellite imagery.
REFERENCES
1. Angima S., Stott D., O'neill M., Ong C., Weesies G., Soil erosion prediction using RUSLE for central Kenyan highland conditions, Agriculture, ecosystems & environment, 2003, Vol. 97, No. 1-3, pp. 295-308.
2. Ballabio C., Borrelli P., Spinoni J., Meusburger K., Michaelides S., Beguería S., Klik A., Petan S., Janecek M., Olsen. P., Mapping monthly rainfall erosivity in Europe, Science of the Total Environment, 2017, Vol. 579, pp. 1298-1315.
3. Castillo C., Pérez R., James M.R., Quinton J., Taguas E.V., Gómez J.A., Comparing the accuracy of several field methods for measuring gully erosion, Soil Science Society of America Journal, 2012, Vol. 76, No. 4, pp. 1319-1332.
4. Desmet P., Govers G., A GIS procedure for automatically calculating the USLE LS factor on topographically complex landscape units, Journal of soil and Water Conservation, 1996, Vol. 51, No. 5, pp. 427-433.
5. Doyle C., A dictionary of marketing, Oxford University Press.
6. Drusch M., Del Bello U., Carlier S., Colin O., Fernandez V., Gascon F., Hoersch B., Isola C., Laberinti P., Martimort P., Sentinel-2: ESA's optical high-resolution mission for GMES operational services, Remote sensing of Environment, 2012, Vol. 120, pp. 25-36.
7. Fernandez C., Wu J., McCool D., Stöckle C., Estimating water erosion and sediment yield with GIS, RUSLE, and SEDD, Journal of soil and Water Conservation, 2003, Vol. 58, No. 3, pp. 128-136.
8. Gorelick N., Hancher M., Dixon M., Ilyushchenko S., Thau D., Moore R., Google Earth Engine: Planetary-scale geospatial analysis for everyone, Remote sensing of Environment, 2017, Vol. 202, pp. 18-27.
9. Lal R., Soil erosion impact on agronomic productivity and environment quality, Critical reviews in plant sciences, 1998, Vol. 17, No. 4, pp. 319-464.
10. Lane S.N., Westaway R.M., Murray Hicks D., Estimation of erosion and deposition volumes in a large, gravel-bed, braided river using synoptic remote sensing, Earth Surface Processes and Landforms: The Journal of the British Geomorphological Research Group, 2003, Vol. 28, No. 3, pp. 249-271.
11. Lehner B., Verdin K., Jarvis A., New global hydrography derived from spaceborne elevation data, Eos, Transactions American Geophysical Union, 2008, Vol. 89, No. 10, pp. 93-94.
12. Lu D., Li G., Valladares G.S., Batistella M., Mapping soil erosion risk in Rondonia, Brazilian Amazonia: using RUSLE, remote sensing and GIS, Land degradation & development, 2004, Vol. 15, No. 5, pp. 499-512.
13. Mallinis G., Gitas I.Z., Tasionas G., Maris F., Multitemporal monitoring of land degradation risk Due to soil loss in a fire-prone Mediterranean landscape using multi-decadal Landsat imagery, Water resources management, 2016, Vol. 30, No. 3, pp. 1255-1269.
14. Millward A.A., Mersey J.E., Adapting the RUSLE to model soil erosion potential in a mountainous tropical watershed, Catena, 1999, Vol. 38, No. 2, pp. 109-129.
15. Mohammad A.G., Adam M.A., The impact of vegetative cover type on runoff and soil erosion under different land uses, Catena, 2010, Vol. 81, No. 2, pp. 97-103.
16. Morgan R.P.C., Soil erosion and conservation, John Wiley & Sons.
17. Panagos P., Ballabio C., Borrelli P., Meusburger K., Klik A., Rousseva S., Tadic M.P., Michaelides S., Hrabalikova M., Olsen P., Rainfall erosivity in Europe, Science of the Total Environment, 2015a, Vol. 511, pp. 801-814.
18. Panagos P., Borrelli P., Meusburger K., A new European slope length and steepness factor (LS-Factor) for modeling soil erosion by water, Geosciences, 2015b, Vol. 5, No. 2, pp. 117-126.
19. Panagos P., Van Liedekerke M., Jones A., Montanarella L., European Soil Data Centre: Response to European policy support and public data requirements, Land use policy, 2012, Vol. 29, No. 2, pp. 329-338.
20. Pandey A., Chowdary V., Mal B., Identification of critical erosion prone areas in the small agricultural watershed using USLE, GIS and remote sensing, Water resources management, 2007, Vol. 21, No. 4, pp: 729-746.
21. Rahman M.R., Shi Z., Chongfa C., Soil erosion hazard evaluation - an integrated use of remote sensing, GIS and statistical approaches with biophysical parameters towards management strategies, Ecological Modelling, 2009, Vol. 220, No. 13-14, pp. 1724-1734.
22. Renard K.G., Foster G.R., Weesies G., McCool D., Yoder D., Predicting soil erosion by water: a guide to conservation planning with the Revised Universal Soil Loss Equation (RUSLE), United States Department of Agriculture Washington, DC.
23. Renard K.G., Foster G.R., Weesies G.A., Porter J.P., RUSLE: Revised universal soil loss equation, Journal of soil and Water Conservation, 1991, Vol. 46, No. 1, pp. 30-33.
24. Touchan R., Baisan C., Mitsopoulos I.D., Dimitrakopoulos A.P., Fire history in European black pine (Pinus nigra Arn.) forests of the Valia Kalda, Pindus mountains, Greece, Tree-Ring Research, 2012, Vol. 68, No. 1, pp. 4551.
25. Van der Knijff J, Jones R., Montanarella L., Soil erosion risk assessment in Italy, Citeseer, 1999.
26. Wilcox B.P., Breshears D.D., Allen C.D., Ecohydrology of a resource-conserving semiarid woodland: Effects of scale and disturbance, Ecological Monographs, 2003, Vol. 73, No. 2, pp. 223-239.
27. Wischmeier W.H., Smith D.D., Predicting rainfall erosion losses-a guide to conservation planning, Predicting rainfall erosion losses-a guide to conservation planning, 1978, Iss. 537, 58 p.
