Sediment Transport Modeling Using GIS in Bagmati Basin, Nepal

1. Introduction

Soils play an important role in the maintenance of global food supplies, an ever increasingly important role as total population expands. The first ‘Global Assessment of Human-Induced Soil degradation’ (GLASOD) was published in 1990 and estimated that 1.97 billion hectares, equivalent to an area of 15% of total land cover, suffered degradation from the mid 1940’s up to 1990. The more recent GLASOD (Global Assessment of Soil Degradation) survey has indicated more than 10⁹ ha of the land surface of the world are currently experiencing serious soil degradation as a result of water erosion. For total suspended sediment yield from the land to the oceans, values closer to 15-20 x ${10}^9$ tons year ${}^{-1}$ have been most frequently cited, whereas average global specific total sediment load is approximately 140-188 tons km^-2 year^-1 [23]. On a global scale, the loss of 75 billion tons of soil costs the world about US$400 billion/year (at US$3/ton of soil for nutrients and US$2/ton of soil for water), or approximately US$70/person/year [13].

Erosion prediction is the most widely used and most effective tool for soil conservation planning and design. Because it is impossible to monitor the influence of every farm and ranch management practice in all ecosystems under all weather conditions, erosion predictions are used to rank alternative practices with regard to their likely impact on erosion. These erosion predictions are thus an essential part of soil conservation programs. Assessment of soil erosion as to how fast soil is being eroded is helpful in planning conservation work. Estimates of the rate of soil loss may then be compared with what is considered acceptable and the effects of different conservation strategies can be determined. Modeling can be an effective method of predicting soil loss under a wide range of conditions as it can provide a quantitative and consistent approach to estimating soil erosion and sediment transport. Using remote sensing and GIS to parameterize such models allows them to be applied over local, regional and global scales.

Two main types of model: empirically based and process based are available for predicting soil erosion and sediment transport. Empirically based technology means regression or lumped mathematical models, which were developed using the experimental data of plot studies on erosion by water. Zingg [29] and Musgrave [18] equations are examples of initial steps towards the empirical soil erosion models. Universal Soil Loss Equation, USLE [26], later revised as Revised USLE or RUSLE [20] is one such model developed in the USA with more than 10,000 plot years of research data and experience of soil scientists. It is the most widely used model for soil erosion estimation because of the simplicity. It is based on the set of mathematical equations that estimate average annual soil loss from inter-rill and rill erosion. In addition, the equation combines interrelated physical and management parameters such as soil type, rainfall pattern, and topography that influence the rate of erosion. Erosion Productivity Impact calculator (EPIC) model [25], which was developed to assess the effect of soil erosion on soil productivity, also uses USLE and Modified USLE (MUSLE) model [24] to simulate erosion process. Chemical, Runoff and Erosion from Agricultural Management Systems (CREAMS) model [11], Agriculture Non-point Source Pollution model (AGNPS) model [28], and Soil and Water Assessment Tool (SWAT) model [1] are the examples of hybrid models which are based on USLE/MUSLE/RUSLE for the erosion estimation but use the sediment transport approach on the basis of continuity equation for sediment yield estimation.

Physically or process based models are intended to represent the essential mechanisms controlling erosion and sediment transport process. These models are the synthesis of individual component that affect the erosion and transport process. Aerial Non-point Source Watershed Environmental Response Simulation (ANSWERS) model [2], Kinematic Runoff and Erosion model (KINEROS) model [27], European Soil Erosion model (EUROSEM) model [17], and Water Erosion Prediction Program (WEPP) [19] are examples of process based models. Although physically based models try to emulate the physical processes involved in soil erosion and sediment transport, the weakness of these models is numerous parameters they need for calibration and also suffer from the problem of equifinality [3].

The overall aim of the study is the modeling of soil erosion and transport processes in distributed manner so that erosion, deposition and sediment yield can be computed and verified with the observations in data limited conditions. To achieve this objective, an empirical model was framed within Geographic Information System (GIS) to predict soil erosion in distributed manner. Then, the sediment delivery approach is used to predict sediment yield in this study. For the empirical approach, the revised form of the USLE model, RUSLE, is used to predict erosion potential on a cell-by-cell basis in conjunction with SEDD model to determine the catchment sediment yield by using the concept of sediment delivery ratio [7].

2. Methodology

A very popular empirical model, known as USLE is used to estimate soil erosion in this study. Then, sediment delivery approach is used to estimate the sediment yield which a part of eroded sediment that appears at watershed outlet. Empirical methods such as the USLE have been found to produce realistic estimates of surface erosion (and also sediment yield) over areas of small size [26, 10 4]. Sediment delivery distributed (SEDD) model couples USLE with a spatial disaggretion criterion of sediment delivery processes. The revised form of USLE, commonly known as RUSLE, is expressed as:

$$A=R*K*L*S*C*P$$ E1

Where, A = average annual soil loss predicted (ton ha^-1), R = rainfall runoff erosivity factor (MJ mm ha^-1 hr ^-1), K = soil erodibility factor (ton ha hr MJ^-1 ha^-1 mm^-1), L = slope length factor, S = slope steepness factor, C = cover management factor and P = support practice factor.

The value of RUSLE factors are computed using the following methods as described in the Agricultural Handbook 703 [20].

$$R=\frac{1}{n}{\displaystyle \sum _{i=1}^n\left({\displaystyle \sum _{j=1}^m{E}_j{\left(I_{30}\right)}_j}\right)}$$ E2

Where, n = total number of years, m = total number of rainfall storms in i ^th year, I ₃₀ = maximum 30 minutes intensity (mm hr ^-1), E _j = total kinetic energy (MJ ha^-1) of j ^th storm of i ^th year and it is given as:

$$E_j={\displaystyle \sum _{i=1}^p{e}_k*d_k}$$ E3

Where, p = total number of divisions of j ^th storm of i ^th year, d _k = rainfall depth of k ^th division of the storm (mm), e _k = kinetic energy (MJ ha^-1 mm^-1) of k ^th division of the storm and is given as [20]:

$$e_k=0.29(1-0.72e^{(-0.05i_k)})$$ E4

Where, i _k = intensity of rainfall of k ^th division of the storm (mm hr^-1)

If λ is the horizontal projection of the slope length (in meter), then L factor is given as:

$$L={\left(\frac{\lambda }{22.1}\right)}^m$$ E5

Where, λ = contributing slope length (in meter), m = variable slope length exponent.

The slope-length exponent ‘m’ is related to the ratio β of rill erosion (caused by flow) to interrill erosion (principally caused by raindrop impact) by the following equation:

$$m=\beta /(1+\beta )$$ E6

For moderately susceptible soil in both rill and inter-rill erosion, McCool et al. (1989) suggested the equation:

$$\beta =\frac{(Sin\theta /0.0896)}{3.0{(Sin\theta )}^{0.8}+0.56}$$ E7

Where, θ = slope angle (degrees).

The slope steepness factor S is evaluated from the following equations (McCool et al., 1987):

$$\begin{array}{l}S=10.8sin\theta +0.03 fors9\%\\ S=16.8sin\theta -0.50fors\ge 9\%\end{array}$$ E8

Where, s = slope in percentage.

C and P factors are assigned to different grid according to land cover while K factor is estimated using the soil data.

In a catchment, not all eroded soil reaches the catchment outlet but a part of the soil eroded in an overland region gets deposited within the catchment. The values of ratio of sediment yield to total surface erosion, which is termed as sediment delivery ratio (D _R ), for an area are found to be affected by catchment physiography, sediment sources, transport system, texture of eroded material, land cover etc. [23]. However, variables such as catchment area, land slope and land cover have been mainly used as parameters in empirical equations for D _R [9, 12].

Ferro & Minacapilli [5] and Ferro [1997] hypothesized that D _R in grid cells is a strong function of the travel time of overland flow within the cell. The travel time is strongly dependent on the topographic and land cover characteristics of an area and therefore its relationship with D _R is justified. Based on their studies on probability distribution of travel time, the following relationship was assumed herein for a grid cell lying in an overland region of a catchment:

$$D_R=\exp (-\gamma t_i)$$ E9

Where, t _i = travel time (hr) of overland flow from the i ^th overland grid to the nearest channel grid down the drainage path and γ = coefficient considered as constant for a given catchment.

The travel time for grids located in a flow path to the nearest channel can be estimated if the lengths and velocities for the flow paths are known. The direction of flow from one cell to a neighboring cell is often ascertained by using an eight direction pour point algorithm in grid-based GIS analysis. Once the pour point algorithm identifies the flow direction in each cell, a cell-to-cell flow path is determined to the nearest stream channel and thus to the catchment outlet. If the flow path from cell i to the nearest channel cell traverses m cells and the flow length of the i ^th cell is l _i (which can be equal to the length of a square side or to a diagonal depending on the direction of flow in the i ^th cell) and the velocity of flow in cell i is v _i , the travel time t _i from cell i to the nearest channel can be estimated by summing up the time through each of the m cells located in that flow path:

$$t_i={\displaystyle \sum _{i=1}^m\frac{l_i}{v_i}}$$ E10

In this study, the method of determination of the overland flow velocity proposed by the US Soil Conservation Service was chosen due to its simplicity and the availability of the information required (SCS, 1975). The flow velocity (v _i ) is considered to be a function of the land surface slope and the land cover characteristics:

$$v_i=a_i*S_i{}^b$$ E11

Where, b = a numerical constant equal to 0.5 [22, 5], S _i = slope of the i ^th cell and a _i = a coefficient related to land use [8]. Introducing equations (10) and (11) into equation (9) gives

$$D_R=\exp (-\gamma {\displaystyle \sum _{i=1}^m\frac{l_i}{a_i{S}_i{}^{0.5}}})$$ E12

It should be noted that l _i / S _i ^0.5 is the definition of travel time used by Ferro & Minacapilli [5]. Values of the coefficient a _i for different land uses were adopted from [8] and are presented in Table 4.

If S _E is the amount of soil erosion produced within the i ^th cell of the catchment estimated using equation (1), then the sediment yield for the catchment, S _y , was obtained as follows:

$$S_y={\displaystyle \sum _{i=1}^n{D}_R*S_E}$$ E13

Where, n = the total number of cells over the catchment and the term D _R = the fraction of S _E that ultimately reaches the nearest channel. Since the D _R of a cell is hypothesized as a function of travel time to the nearest channel, it implies that the gross erosion in that cell multiplied by the D _R value of the cell becomes the sediment yield contribution of that cell to the nearest stream channel. The D _R values for the cells marked as channel cells are assumed to be unity.

3. Study Area

The study area selected for this study is Bagmati Basin, Nepal. The basin is chosen because of its bio-climatic diversity due to elevation differences from valley floors to mountain summits, and related land use changes having influence on soil erosion, which is considered typical for the Middle Mountains of Nepal. Bagmati is the draining river from the Kathmandu city which is the capital of Nepal. The Bagmati basin covers an area of 3,500 km² in total and drains out of Nepal across the Indian State Bihar to reach the Ganges. The watershed with the elevation ranging from 57 m to 2,913 m is situated at latitude of 26° 30’ to 28°N and longitude 85° to 86°E. The watershed can be divided into three main areas: the upper, middle and the lower Bagmati watershed areas (BWA). The Upper Bagmati Watershed Area covers the whole of the Kathmandu valley including its source at Shivapuri. From the Chovar gorge, the river flows into the Middle Bagmati watershed Area across the Mahabharat and Siwalik ranges. The catchment area of upper and middle Bagmati basin is about 2,800 km². The terrain of the upper and middle BWA is rugged and comprised of several steep mountains except Kathmandu valley. The area of upper and middle Bagmati basin draining to Karmaiya is considered in the study on the basis of data availability.

Figure 1.

Bagmati basin, Nepal

The climate of the Bagmati watershed can be subdivided into three altitude/climate zones. These are: (a) Subtropical sub humid zone below 1,000 m: the southern most parts of the Bagmati watershed area including the Siwaliks region lie in this zone, (b) Warm temperate humid zone between 1,000-2,000 m: a large part (more than 60%) of the BWA lies in warm temperate humid zone between 1,000 – 2,000m altitudes and (c) Cool temperate humid zone between 2,000-3,000 m: only a small portion (about 5%) of the Bagmati watershed falls above 2,000 m. The annual average rainfall in the watershed is about 1,800 mm and it produces 1,400 mm of runoff per year on average, which accounts for about 75% of annual average rainfall. In the basin, steep slope in mountainous area and land use change are the major factors of soil erosion, which is considered typical for the Middle Mountains of Nepal. Total population in the catchment is about 1.5 millions. Figure 1 shows the map of the catchment along with streams and tributaries.

S N	Data type	Stations	Location		Duration	Remarks
Lat. (N)	Long.(E)
1	Rainfall	Daman	27° 36’	85° 05’	1987-97	Daily
		Hetauda	27° 25’	85° 03’	1987-97	Daily
		Godavari	27° 35’	85° 05’	1987-97	Daily
		Airport	27° 42’	85° 22’	1990-97; 1993-97	Daily; Hourly
		Nagarkot	27° 42’	85° 31’	1990-97	Daily
		Sindhuligadhi	27° 17’	85° 58’	1990-97	Daily
		Karmaiya	27° 07’	85° 28’	1990-97	Daily
2	Sediment	Karmaiya	27° 07’	85° 28’	1990-91,93,95-97	Daily

Table 1.

Description of hydrologic data set

S N	Description	Scale or; grid resolution	Source	Remarks
1	DEM	90 m (SRTM DEM)	USGS	Raster
2	Landuse	1:25,000	BIWMP	Vector
3	Soil	1:25000	BIWMP	Vector
4	Watershed boundary	1:25,000	BIWMP	Vector
5	River network	1:25000	BIWMP	Vector

Table 2.

Description of spatial data set

Hydrologic data (rainfall, evaporation, suspended sediment concentration) for the basin are obtained from Department of Hydrology and meteorology (DHM). Digital Elevation Model (DEM) data, in 90 m resolution, was obtained from obtained from United States Geological Survey (USGS) (available at: http://srtm.usgs.gov). STRM DEM provides comprehensive and consistent global coverage of topographically derived data sets, including streams, drainage basins and ancillary layers. Other spatial data set such as: soil, land use, basin boundary, river network are obtained from Bagmati Integrated Watershed Management Programme (BIWMP). The details of hydrologic data are provided in Table 1 while Table 2 contains the details about spatial data set.

As observed in the DEM of the watershed (Figure 2), the elevation varies significantly from as low as 137 m to as high as 2913 m from mean sea level. Lower part of the watershed is relatively flat compared to the upper and middle part. Kathmandu, the capital of Nepal lies in the upper part of the watershed. One third of the watershed is relatively flat as 34% of the watershed area has slope in the range of 0 - 10%. About 50% of the area has mild slope ranging from 10 - 30%. Remaining 15% watershed contains high slope with slope value more than 30%.

Figure 2.

Digital Elevation Model (DEM) of Bagmati basin (SRTM DEM, 90 m resolution)

The land use in the watershed is observed to be mixed type. Cultivated land is major land use pattern in the upper part of the watershed while in middle and lower part of the watershed, forest area is seen to be dominant land use type. Majority of built-up area falls on the upper part of watershed, which represents Kathmandu. The land use pattern in the watershed is presented in Figure 4. 5. More than half of the watershed area (58%) is covered by forest. Cultivated land accounts for 38% of the area of the watershed while nearly 4% of the land in the watershed is barren. The land use distribution in the watershed is presented in Figure 3. The most extensive soils in the area are Dystrochrepts, Hapludalfs and Haplumbrepts, which occupy most of the hilly and mountaineous land. The texture of these soils is sandy/loamy in nature that varies from sandy clay to loam. The Dystrochrepts are also the most important soils in the inner Terai valleys. Soil type Rhododtalfs is commonly found in the gently undulating slopes and restricted to scattered, quasi-subtropical areas in the lower Hiamlayas. These soils are prone to severe soil erosion. The soil in the south face on the low altitude Mahabharat range is Dystrochrepts and Hapludalfs. These soils are mostly cultivated. The Haplaquepts are the dominant soils in the Terai plain as well as on paddy fields in hilly areas and elsewhere. Major soil types in the mountainous lands are Haplumbrepts and Dystrochrepts. Loamy soil texture is dominant in the watershed as demonstrated in Figure 4.

Figure 3.

Land use distribution in Bagmati basin

Figure 4.

Soil map of Bagmati basin

Data Preparation And Simulation

Revised Universal Soil Loss Equation is one of the simplified models, which predicts soil erosion from hillslopes. The factors such as rainfall runoff erosivity factor (R) associated with the model represent the effects of climatic parameters in soil erosion while soil erodibility factor (K) represents the nature of the soil, its characteristics and influence in soil erosion. Topography and land use practices are other major factors incorporated in the model to account their effects in soil erosion.

Out of seven rainfall stations in Bagmati basin, one station measures hourly rainfall while remaining six other stations measures daily rainfall. So, rainfall data from these seven stations are analyzed to find the correlation in the rainfall pattern. The analysis of daily, monthly and annual rainfall trends of these stations showed that the trend was similar for all these stations. This helped to in disaggregating daily rainfall data into hourly data for the remaining six stations. For the basin, rainfall erosivity index “R” value was computed for monthly basis R value for was computed using equations (2), (3) and (4) since sediment yield information was available on monthly basis. Soil erodibility (K) factor values were assigned on grid by grid basis on the basis of soil texture [21] of the basin and assigned K values are presented in Table 3. Topographical parameters (L, S) were extracted from 90 m resolution SRTM digital elevation model (DEM) obtained from USGS (http://srtm.usgs.gov). Equation (5) was used for L factor calculation while S factor was computed using equation (8) for each cell. Similarly, C value, which depends on land use, was obtained from different literature [21], [16]. The C values used assigned for different land use in the basin are tabulated in Table 4. In case of P factor, the value is taken 0.5 for agricultural land and for rest of the land use; P value is assigned to be 1.

Textural Class	Organic matter content (%)
0.5	2	4
Fine sand	0.0211	0.0184	0.0132
Very fine sand	0.0553	0.0474	0.0369
Loamy sand	0.0158	0.0132	0.0105
Loamy very fine sand	0.0580	0.0501	0.0395
Sandy loam	0.0356	0.0316	0.0250
Very fine sandy loam	0.0619	0.0540	0.0435
Silt loam	0.0632	0.0553	0.0435
Clay loam	0.0369	0.0329	0.0277
Silty clay loam	0.0487	0.0422	0.0343
Silty clay	0.0329	0.0303	0.0250

Table 3.

Soil Erodibility factor by soil texture in SI unit (ton ha hr MJ^-1 ha^-1 mm^-1)

S N	Land Use	C value basis	C Value	a value
1	Cultivation	Crops, disturbed land	0.4000	1.55
2	Water body	Depositional sinks	0.0001	3.08
3	Forest	Forest	0.0020	0.76
14	Barren land	Fallow	1.0000	3.08
15	Built-up area	Paved, occasional construction	0.0005	6.19

Table 4.

Cover management factor (C) on the basis of land use

Results And Discussion

"Once RUSLE parameters for Bagmati basin was computed following the procedure outlined alobe, sediment delivery ratio (SDR) map for Bagmati basin computed using Equation (12)". The SDR map for the basin is presented in Figure 5 below. It is observed that flat areas around the south and north parts of the watershed has low sediment delivery ratios while the hilly areas within the watershed had higher values for sediment delivery ratio. This finding is consistent with the fact that steep areas are supposed to have higher sediment delivery ratio compared to flat areas. In terms of watershed management perspective, the areas with higher values of SDR should be given higher priority compared to areas with lower SDR values for implementation of erosion control measures in this watershed.

Figure 5.

Sediment Delivery Ratio (SDR) map for Bagmati Basin

The sediment yield data are available for only few months of the year for Bagmati basin. So, it was not possible to analyze the long term sediment yield value and thus, monthly computation is carried out. Soil erosion map and SDR map was used to compute the sediment yield value at the watershed outlet. The Observed monthly sediment yield was compared with the computed as seen in Figure 6 below. The simulated result using this approach is fairly consistent with the observed data although this methodology slightly overpredicted sediment yield for the most of the observed months. There can be several reasons which can lead to overestimation of sediment yield values. For example, only one rainfall station had hourly measurement while remaining stations recorded daily values. It was assumed that the rainfall pattern over the watershed was similar. If rainfall data with finer temporal resolution were available for all the stations, the computed of R value would have been more reliable.

Figure 6.

Comparison of simulated and observed sediment yield

The comparison of observed and computed sediment yield also indicate that great care is required in the selection of input values for the rainfall (R) and soil erodibility (K) factors. The USLE model was developed from the data suing the experiments that were carried out on a standard plot of 22.1 m length of uniform 9% slope. So, USLE-based performance can expected to be better for finer (for example 30 m) DEM resolution. Earlier studies have demonstrated that DEM resolutions can affect the outcome of RUSLE based simulations and better agreement can be obtained using fine DEM resolution [4]. Similarly, RUSLE results may be improved if more detailed soil, land use/cover data are available.

The model prediction may have been improved if γ coefficient was calibrated using the measured sediment yield values at mean annual scale for SDR computation. During SDR calculation, the sensitivity analysis of the parameter γ showed that the computed Sy was not very sensitive to γ in equation (12). The variation of γ value by 15 times (from 0.1 to 1.5) changed the Sy value only 10%. Since large variation in γ affected Sy insignificantly during sensitivity analysis, γ value was taken as 1 in the computation for simplicity. The sensitivity analysis has supported the findings of Jain & Kothyari [10] where they had reported that Sy was not very sensitive to γ in their study.

Conclusion

Soil erosion is a natural process. Modeling a natural process using mathematical simulation involves use of complex relationships. The number of factors associated with such complex process imposes their effect in various degrees. It is, thus, essential to consider only those factors, which are likely to have dominant effects in the process while carrying out mathematical simulation. This simplifies the process and is acceptable in most cases. Universal Soil Loss Equation (USLE) (and its revised form, RUSLE) is one of such simulation model, which predicts soil erosion from hillslopes. The factors such as rainfall runoff erosivity factor (R) associated with the model represent the effects of climatic parameters in soil erosion while soil erodibility factor (K) represents the nature of the soil, its characteristics and influence in soil erosion. Topography and land use practices are other major factors incorporated in the model to account their effects in soil erosion.

This study is an attempt to estimate soil erosion and sediment yield at Bagmati River basin using existing conceptual methods and GIS. This methodology can be used for the identification of sediment source areas and prediction of sediment yield at a catchment scale with available optimum data sets. ArcGIS was used for discretizing the catchment into grid cells of different resolutions. Grid cell slope, drainage direction and catchment boundary were generated from DEM using pour point method. The DEM was further analyzed to classify the grid cells into overland flow and channel region by using channel initiation threshold area approach. After preparing different USLE parameter layers, the gross surface erosion map was computed. The sediment delivery ratio of overland flow cell was assumed to be a function of the travel time of overland flow from given cell to the nearest downstream channel cell. For channel cells, the sediment delivery ratio was assumed to be unity. The computed and observed values were observed to have some discrepancy for monthly sediment yield. The variation is resulted by the few assumptions made during the analysis. In the study, computation of soil erodibility value (K) was based on soil texture only. Similarly, constant cover management factor (C) values were used instead of time varying because of the lack of series of land-use map for different years. Use of finer resolution DEM can also improve the estimation of slope length (L) and Slope steepness (S) factor. Improved results can be expected if these enhancements are incorporated. The proposed modeling framework is simple and can be a useful tool in conservation planning with reasonable reliability at data scarce areas.

Acknowledgments

The author would like to thank Department of Hydrology and Meteorology, Nepal and Bagmati Integrated Watershed Management Programme for providing the data used in this study. The present work benefited from the input of Dr. Dushmanta Dutta (CSIRO, Australia), who provided valuable comments and assistance to the undertaking of the research summarised here.