Remote Sensing and GIS Technologies for Denudation Estimation in a Siwalik Watershed of Nepal

Remote Sensing and GIS Technologies for Denudation Estimation in a Siwalik Watershed of Nepal

Kioshi Honda**, Lal Samarakon*, akichika Ishibashi*, yasushi Mabuchi*, Shigechika Miyajima***
*Nippon Koei Co., Ltd., Research & Development Center,
2304 Takasaki, Kukizaki, Japan
**STAR/SERD, Asian Institute of Technology, Bangkok, Thailand
***Water induced Disaster Prevention Technical Centere, Katmandu, Nepal

Abstract
Wide spread and intensive soil erosin has resulted gradual land degradationin most parts of the mountains ofNepal threatening the prosperity of Nepalese society and the natural environment. The rate of denudation is highly varying with natural factors together with the ever increase demand for forest resources and utilsable land for the growing population. Comprehensive studies extending from watershed basis to regional level have to be carried out for introducing appropriate counter measures in reducing devastating disasters.

Present study is an attempt to use remote sensing to identify the land degradation in Ratu watershed in the Central Siwalik area, and establish a method to estimate the rate of denudation in GIS perspective, Landsat data procured for a period of 20 years from 1973 to 93 were analyzed for the change of forest cover in the waterhed and topographical parameters were used in a model to estimate the probable annual soil loss. Model was further improved for flood event soil yield estimation using factors that have identified in the field as main causes for intensive erosion during a heavy rainy season.

1. Introduction
Nepal is a mountainous country with an area of 147.181 square kilometers. Two third of the total area of the country is occupied by hills and mountains. Population stands at 18.8 million (1992), with a growth rate of 2.1%. It is estimated that 75% of the total energy comes from forest, (WECS, 1992). This utilization practices lead to high degree of degradation of the land and loss of natural forest resources. Exact reasons and the actual rate of land degradation is not properly investigate for the entire country, but most of the literature published blaming the farmer for destruction of forest for energy subsidy and mountain farming. The detail examination of sedimentation in watershed basis may be a better soluition to understand the cause and effect of the flood and soil erosion process, but the limited resources may hinder this efforts.

Lack of knowledge of natural event, their occurrence, recurrence time could hinder proper planning or disaster prevention, and mitigation. It was though satellite remote sensing data could offset the information gap as they can provide timely, and repetitive data of aparticular phenomena without physically visiting the area of interest. In this study, it was decide to investigate the applicability of satellite data in land degradation assessment, and monitoring the phenomenon. Also, a soil yield estimation method was proposed, which uses generally available information based on a GIS environement.

2. Outline of the Study Area
Ratu watershed, which is situated in the Siwalik area was selected for the present study. Ratu river originates from a Churiya hills of central Nepal at an altitude of 700 meters and drains into Terai region. The watershed lies west of Kamala river and south of sun Koshi. Ratu river crosses the East-West highway aroung 86E longitude and 27N latitude. River runs dry during the dry season and transport heavy load of sediment during the rainy season with the high water flow.

The Siwalik area is composed of young tertiary strata, and contains some of the most easily erodible lithologies including unconsolidated sand and gravel. Soft and loose sandstone, mudstones, and conglomerate are the predominant rock types in the present watershed. About 70% of the watershed area is dominated by forest with varying cover densities. Shorea Robusta, locally called Sal is the major species that is found in the area. Remaining land could be grouped into crop lands and barren land according to their land use practice.

3. Data Acquisition

3.1 Acquisition of Remote Sensing and field Verification Data
It was decided to obtain satellite remote sensing data over a long time span to estimate the changes,

Specifically to observe the forest cover and its change that could be used as a key to land degradation. Also field investigation was planed to identify the physical nature of the soil erosion. It transport and sedimentation to evaluate the volume of sediment production that could use to compare satellite data evaluated soil production. A major flood event was considered to estimate the volume of sedimentation by field observations. The most recent flood event, 1993 July was selected for this as it could give the opportunity to collect the information on the extent and nature of the flood by visiting the field where no other means of information is available. Referring to available satellite data sources, data shown in the Table 1 was selected for the study. Further, during the selection of data, special attention was given to choose data that are acquired in the same season of the year to minimize irradiance variations on the surface due to positional difference of the sun. The LISS sensor data was acquired irradiance variations on the surface due to positional differences of the sun. The LISS sensor data was acquired for investigation of its capability in replacing the largely used Landsat sensor for land information extraction.

Table 1 Information on the satellite data obtained for the study

Sensor	Date	Sun Angle	Sun Azimuth	Processing Level
MSS	1973.03.14	Not Known	Not Known	System Corrected
MSS	1977.03.20	Not Known	Not Known	System Corrected
TM	1993.03.16	47°	127°	System Corrected
TM	1995.03.22	46°	127°	System Corrected
LISS-II	1995.03.15	47°	127°	System Corrected

Referring to the availability of satellite data, and considering major flood events in the past, aerial photographs acquired in 1979 and 1992 were obtain for the study. The scale of the photographs are 1:40,000, and decide to photograph the Ratu watershed using a helicopter. A normal camera that uses visible light with 50mm focal length was used to obtain vertical photographs. Further, some of the sub-streams were surveyed for recognize the sediment formation, distribution and cludes to differentiate the age of the deposits.

3.2 Geometric Correction of Satellite Data
The spatial resolution of the acquired satellite data are 30m, 57x79m and 34m for TM, MSS and LISS, respectively. In addition, their spatial orientation also dissimilar. Therefore, the procured satellite data were brought into a common map projection (UTM) by constructing mapping function through identifying control points on 1:25,000 topographical maps. First order transformation functions were established for all the datasets, and nearest neighbor resampling was used for re-mapping the satellite data into UTM projection.

3.3 Radiometric Rectification of Satellite Data
Digital comparison of multi-sensor data requires further adjustment as the observations are made in sensor specified discrete spectral bands. The amount of energy received at the sensor from a particular earth feature is a function of received energy, reflectance, atmospheric propagation, sensor sensitivity, and the spectral bandwidth. It we consider a homogeneous land surface and similar atmospheric condition, the amount of received bandwidth. It we consider a homogeneous land surface and similar atmospheric condition, the amount of received energy is a function of bandwidth of the spectral channel. In the present study, the satellite data of three sensors adjusted for the differences in the sensor bandwidths considering linear spectral reflectance within each respective Landsat TM spectral bands. Further, histogram equalization was carried out to compensate for different atmospheric conditions that would have prevailed at the time of satellite passes having 1993 landsat TM as reference dataset.

3.4 Database Creation
The conventional maps, topographical and vegetation were digitized in establishing a GIS database in the UTM map projection for land degradation and soil erosion analysis. Further, the aerial photographs and helicopter photographs were scanned, rectified and registered into UTM projection. Thus the completed GIS database for the present study included multi-sensor temporal satellite data, aerial photographs, elevation and land cover information.

4 Method of Analysis

4.1 Land Degradation and satellite Data
As the satellite data are obtained for the dry season, classification of forest and other lands was straightforward as the crop lands are uncultivated. This can easily be accomplished by satellite data derived vegetation indexes. Further, the changes in the forest cover in its density and any deforestation can also be

Estimated by vegetation indices. Most of the vegetation indices are based on empirical evidences, and nearly all of commonly used vegetation indices are concerned only with red and near-infrared regions of the electromagnetic spectrum. NDVI based on the equation 1 was used with multi-sensor digital data that were adjusted for bandwidth differences and compensated for atmospheric variations at the time of satellite passes.

1992 aerial photographs were taken as ground truth information to classify 1992 Landsat TM data. Also, these photographic information were further compared with other dates of satellite data as the observation conditions were re-established to 1992 satellite pass. Aerial photographs covering the study area were classified, and digitally compared with NDVI images of fours dates for evaluation of forest land changes in the watershed.

4.2 Soil Erosion Model
Soil erosion in this region could be due to tectonic activities and/or high intensity continuous rainfall, (lves, 1989, Summerfield, 1991). The tectonic activities are yet to define precisely when quantifying erodibility, and these activities are rather a continental process and may not be evaluated in micro level watershed analysis. On the other hand, water erosion is more regional and can be considered in micro scale, but lack of field measurements on precipitation and runoff could hinder the application of most of the empirical soil erosion models. Current soil erosion models, Universal Soil Los Equation and, most of the ther models that are commonly used are developed by various authors on extensive field observations. Further, most of these equation requisite long term rainfall records, and measured denudation rates under different land use conditions.

Estimation of soil erosion in the Ratu watershed, utilizing a model that incorporate rainfall data is questionable due to lack of field observations. The other major factor that decide the rate of soil erosion in this watershed is the surface topography, hence a model that relate the surface topography was considered for erosion potential estimation. Honda, 1993 proposed a method for soil erosion in mountainous region by incorporating surface gradient. He has demonstrated the accuracy applying in a well monitored watershed in Japan. The annul average denudation E can be estimated as;

Where, E30 - rate of denudation at a slope of 30o; S-gradient of the pixel under consideration; S30 - tan (30o)

The denudation factor has to be identified for every spatial location of the watershed or, has to be measured for well distributed points to apply the model as it is. This is not realistic even in a well monitored suggested. As no field measurements are carried out, denudation of dominant land cover types are extracted from documented materials, and a relationship with the respective NDVI values are established. Given the NDVI value for any location, the rate of denudation is estimated based on the relationship. As the relationship between denudating and NDVI plays a critical role in this method, attempt was made to identify the potential of NDVI in representing the state of the forest cover in the watershed.

4.3 Erosion Characteristics of Storm Event
I was observed in the field that the landslides along the banks of streams, undercutting, and earth topples are some of the factors that could highly contribute to sediment transport during a storm event. This phenomenon was specially visible in the sub-streams, (JICA, 19960. Therefore, it was decided to include the distribution of sub-streams or a factor that represents the valleys in the flood event soil erosion model.

5 Results and Discussion

5.1 land Degradation Assessment by Vegetation Index
The actual size of the bareland within the forest cover is almost smaller than the spatial resolution of the sensors, therefore, photographs were overlaid and compared to investigate the possibility of identifying them using satellite data. Degraded mountain ridges, and soil deposited tributaries are very prominent in 1992 aerial photographs and helicopter photographs. These photographs were classified assigning a value zero for forest covered pixels, and 100 for non-forested pixels. Subsequently, original pixels were resampled into 100x100 meter to represent the average of original pixels fell in each new pixel. Thus the digital value of each newly created pixel represented the bareland or forest area percentage within the pixel concerned. Accordingly, the Landsat TM NDVI image of 1993 was resampled into same pixel size generating average NDVI values for each of the 100

Meter pixels. The relationship of NDVI and the bareland ratio was compared to identify the potential of NDVI in recognizing the change of forest density. Categorized bareland percentages into 10% intervals, and the corresponding averages of NDVI values are show in Table 2.

Table 2 Distribution of photographs derived vareland percentage and 1993 values

Bareland%	0	1-10	10-20	20-30	30-40	40-50	50-60	60-70	70-80	80-90	90-99	100
NDVI-1993	125.37	124.44	122.77	120.91	120.04	118.9	114.65	113.1	107.07	104.58	100.93	98.71

This table shows that the NDVI decreases with the increase of bareland percentage. The highest NDVI for this area was 125, and lowest was 99 for complete forest and 100% bare, respectively. This inverse relationship shows that NDVI can be used to represent the land degradation in the Ratu watershed.

Similarly, NDVI values of 1973, 77,95 and 95-LISS data were compared with the bareland percentages. Distribution of the NDVI values and the aerial photographs established 1992 bareland percentages are shown in Figure 1. The 1979-MSS and 95-LISS data are shown in break lines as they show some incompatibility with the other three dates data. 1977 MSS data shows higher NDVI values than others, and the deviation is irregular. It was found the there was a thin cloud cover in this dataset, and thid might have formed higher level of scattering that was not presented in the other datasets. The LISS dataset also shows different trend of NDVI. The reason for this is not clear, and required to study the internal calibration procedures of Landsat and IRS sensors. These two dates data, 1977 MSS and 1995 LISS were excluded in the proceeding analysis for soil erosion estimation.

Figure 1 Distribution of bareland percentages and NDVI values
Further, the helicopter photographs incorporated in the GIS database were compared with NDVI image of 1993 to investigate the degree of accuracy in identifying the denuded ridges and silted tributaries in the watershed. It was observed that most of the degraded redges, and relatively larger sub-streams can clearly be the densely forested areas. This could be due to the spatial resolution limitation of the sensor and mixed land cover that could present in a single pixel of TM data. As a whole, the comparison of aerial photographs and helicopter photographs with satellite data showed that the remotely sensed data can satisfactorily be used in continuous nature of NDVI facilitate the interpretation of forest density more realistically than classifying the forest cover into discrete density classes.

5.2 Temporal Changes of Average Annual Soil Yield
Estimation of soil yield in the Ratu watershed was carried out based on the equation 2. the rate of denudation required in this was acquired from published literature. Table 3 shows denudation rate in the Himalayan region and in a regulated watershed in Japan.

Rate of erosion for East Nepal, 780-3680 tons/km2/year represents about 0.4 mm/year, and 18.5 mm/year, respectively assuming the density of the sediment is 2 g/cm3. These figures are for extreme conditions

In the Eastern Siwalik area.

Table 3 Land denudatinfor some Nepal and Japan

Description of the site	Rate (tons/km2/year)	Reference
SiwalikL East Nepal, S-aspect, sandstone foothills, landuse
from forest to grazing lands	780-3,680 (0.4-1.8mm/year)	Chatra 1976
Siwalil: Far West Nepal, S-aspect,sandstone foothills
Degraded Forest	2,000 (1 mm/year)	Laban, 1978
Degraded forest, gullied land	4,000 (2mm/year)	Laban, 1978
Severely degraded, heavily grazed	20,00 (10mm/year)	Sakya
Middle MountainsL Katmandu valley,steep slopes 800	(0.4mm/year)	Laban, 1978)
Japan: Asio region n Tochigi Prefecture for 30o slope
Bare lands	20 mm/year	Honda, 1993
Grasslands	1 mm/year	Honda, 1993
Forest	0.1 mm/year	Honda, 1993

In order to evaluate the soil yield, the denudation rate of each pixel has to be estimated using reference information given in Table 3, or otherwise. The values obtained by comparing aerial photographs in the Table 2 represent the average NDVI for forest density percentages in the Ratu watershed. The highest value of NDVI, 125 was for 100% forest cover in this area, and the lowest value represented bareland. In establishing denudation rate in the Rate watershed, three critical soil loss rates and corresponding NDVI values were recognized. Publish documents identify the Siwalik area as severely degraded, hence the least degraded areas in the Siwalik area was assigned a denudation rate of 0.4 mm/year (800 tons/km2/year), which is similar to Middle Mountains. The average NDVI of the Middle Mountain forest coverage was 144, and this value was used in conjunction with the denudation rate of this area. Accordingly, the 100% forest coverage within the watershed, which is an average for 100x100 meter pixel unit, was assigned a denudation rate of 2.0 mm/year (4000 tons/km2/year). Finally, bareland was assigned 20mm/year which was used for a Japanese watershed. All of these denudation rates are assumed to be for land surface with 30o slope angle. In order to estimate the soil yield for the whole watershed a relationship was considered for denudation and NDVI, (Honda, 1994). He has demonstrated that NDVI is related to common logarithm of the denudation rate. The interpretation of denudation with respect to NDVI for each pixel in their original resolution was carried out as shown in Figure 2.

Figure 2 NDVI and Log E₃₀ relationship
This relationship gives the rate of denudation for any pixel in the watershed with respect to 1993 land cover condition. The other factor that is required for soil yield estimation is the slope gradient, and it was calculated from the digital elevation data incorporated in the GIS database. Considering no appreciable change of the topography during the 22 year period from 1973 to 95, proposed equation for soil erodibility was used for estimation of the denudation rate and the volume of the soil production based on NDVI derived forest cover densities. The estimated values are tabulated in Table X.

Table 4 soil yield and denudation estimated for the Ratu watershed

	1973	1993	1995
Soil Yield m³	271.110	321,156	292,096
Rate of denudation mm/year	3.29	3.89	3.55

Reliable date for comparison with satellite data estimation was not available. The estimated denudation rates are very much comparable with published materials presented earlier justifying the satisfactory application of present model in soil yield prediction in Ratu watershed.

5.3 Soil Yield in a Storm Event
The foregoing section evaluated the annual average soil yield on the basis of land use and the surface topography. This model is not suitable for estimation of soil production during a cloud bust incidence that could bring about 500 mm of rainfall within few hours. The excessive rainfall could accelerate the surface erosion, but it was identified that the landslide along the banks of the streams, undercutting, or earth topple are some of the main factors highly contribute to soil erosion during 1993 flood event. Also, it was found the areas around the sub-streams are very much vulnerable to soil erosion during a flood event due to these factors. Attempt was made to modify the equation used in the previous section to use in storm event.

In establishing a storm factor for the present mode, it was considered the field observation are true and correct. Therefore, the soil yield estimated for the 1993 storm event for each sub-watershed was considered as the actual amount of sediment produced, (UICA, 1996). Hence, equating these with the estimations of the model that is adjusted for a storm event prodeuced a storm event factor (a) of the Raut watershed for a cloud burst even similar to that of 1993. The calculation was carried out as shown below:

Here (v)_i is the estimated production of watershed i by field investigation
a_i Storm event factor for watershed i
å_i^alpixels Estimation of soil production of watershed i using annual average model

The average of ai for the sub-watersheds was 11.2 with a standard deviation of 6, Very high deviation were found in sub-watersheds with larger extents. This only relates field observation with the model estimation therefore, further analysis was carried out to study the relationship of ai with topographical characteristics of subwatersheds.

Soil production carried out by field observation was based on stream length, size of landslide scars, depth of deposited materials, vertical drop of valleys, and stream gradient, (UICA, 1996). Therefore , the stream length could be considered as a suitable factor to represent the soil production during a storm event. A strem length factor (Si) of each watershed was calculated as shown in equation 4, and compared with respective ai to identify the possibility of replacing the ai with the stream length factor.

S_i = (L_i x L_i) / A_i (4)
Whee L_i is stream length and A_i is the area of a sub-watershed i

The average value of the above calculation was 9.2. This figure is comparable with the factor estimated for the 1993 storm event (11.2), hence a stream length factor as shown in equation 4 could replace the ai in the equation 3 in estimating the soil yield for a storm event for the Ratu watershed. In the present calculation the main stream were excluded as it was observed that the explained soil lumps are mainly distributed in the tributaries with relatively steep riverbanks.

6. Conclusion
In concluding it can be said that satellite data can satisfactorily be used in land degradation estimations. Multi-temporal data from the same satellite can use through graphical correction for spectral bandwidth and atmospheric dissimilarities. Further analysis are warranted for comparison of multi-sensor spectral data. The model proposed to average annual soil yield produced land denudation rate that is very much similar to the values published by intensive field work. This validates the potential of this model in estimating soil yield in a watershed where the accurate field measurements are scarce to find. Modification of the model for a storm event was suggested incorporating stream densities in the watershed. It is important to say that the present model was used in a watershed in the Siwalik area, and more case studies are required to identify it as a general model.

References

Honada, K., Murai, S., and Shibasaki, R.: Prediction of vegetation restoration by erosion control works in Asio copper mine Japan, 13th IGRSS, 1993
lves, J.D. and Messerli, B.L The Himalayan Dilemman, Routledge London and New York, 1989
JICA; Research Report on the Investigation of Landslides and Soil Erosion in Nepal using Remote sensing Technology, March 1996
Summerfield, M.A.: Global Geomorphlogy, Longman Scientific & Technical, 1991