Parameter evaluation for soil erosion estimation on small watersheds using SWAT model

: This research was undertaken for the evaluation of soil erosion using the semi-distributed basin scale SWAT model for four subcatchments of the Dhrabi River Catchment (DRC), which is located in the Pothwar Plateau region. Two subcatchments (catchment-25 and -31) are characterized by gullies while the other two (catchment-27 and -32) are managed with terraced landuse system. The performance of the model was satisfactory with coefficient of determination ( R 2 ) = 0.67 to 0.91 and Nash-Sutcliffe efficiency (E NS ) = 0.54 to 0.85 for both surface runoff and sediment yield during the calibration (2009-2010) and validation (2011) periods. The P USLE factor was found to be the most sensitive parameter during model calibration. It was observed that all of the rainfall-runoff events occurred during the monsoon season (June to September). The estimated annual sediment loss ranged from 2.6 t/hm 2 to 31.1 t/hm 2 over the duration of the simulation period for the non-terraced catchments, in response to annual precipitation amounts that were between 194.8 mm to 579.3 mm. In contrast, the predicted annual sediment levels for the terraced catchments ranged from 0.52 t/hm 2 to 10.1 t/hm 2 due to similar precipitation amounts. The terraced catchments resulted in 4 to 5 times lower sediment yield as compared to non-terraced catchments. The results suggest that there is a huge potential for terraces to reduce soil erosion in the DRC specifically and Pothwar area generally, which have proven to be an efficient approach of establishing soil and water conservation structures in this region.


Introduction
Globally, problems related to soil and water resources are a growing concern.It has been estimated that between 549-1094 million hm 2 of land are affected by wind and water erosion, respectively [1,2] .The land degradation due to soil erosion is a severe problem which threaten the soil resources and agricultural productivity [3,4] .Borrelli et al. [5] predicted that the highest rates of soil erosion in Southeast Asia, Sub-Saharan Africa and South America occur in the least developed economies.According to Food and Agriculture Organization (FAO) [6] , South Asia countries are losing at least 10 billion USD per year due to agricultural production losses occurring from land degradation.Annually, the global productivity loss ranges from 13 billion to 28 billion USD for dryland (rainfed cropland) production [7] .More than 10 billion hm 2 of land worldwide are experiencing severe soil degradation due to water erosion, based on a global survey [8] .Current water erosion rates are accelerating at a second order of magnitude on arable land which creates a major imbalance between soil formation and soil loss by a factor 10-100 [9,10] .The formation of new soil is a notoriously slow, gradual and continuous process as evidenced by the formation of one-inch of new topsoil, which requires 100-1000 years [11] .However, the rate varies widely, depending on climate, time, parent material, topography and living organisms.Globally, soil is being lost from land areas 10-40 times faster than the rate of renewal and annually about 10 million hm 2 of cropland is lost due to soil erosion [12] .
The majority of Pakistan is located in arid to semi-arid climatic zones, which accounts for a total geographical area of 80 Mhm 2 .About 16 million hm 2 of this overall area is vulnerable to soil erosion, and 11 Mhm 2 that is particularly exposed to water erosion [13] .Several factors are involved in accelerating soil erosion such as urbanization, deforestation, overgrazing, improper tillage practices, leaving the land fallow resulting in low organic matter, land-tenure system, small and fragmented land holdings, and overall poverty [14] .Land degradation is one of the most important issues in the rainfed area of the Pothwar Plateau, where 1.21 million hm 2 out of 2.2 million hm 2 , is affected by gully erosion and only 0.61 million hm 2 is cultivated [15] .The major reasons for extensive soil erosion in the region are uneven topography, steep slopes, erratic and high-intensity rainfall and absence of appropriate management practices.It was estimated that 150 t/hm 2 to 165 t/hm 2 of soil is eroded annually [14] .This rise in soil erosion has endangered soil and water conservation structures, increased the loss of fertile soil and vegetation, resulted in reservoir depletion, and surface and groundwater contamination.
Annual soil loss over different landuses in the Soan River catchment of the Pothwar region ranges from 18.70 t/hm 2 to 63.5 t/hm 2 using Morgan approach integrated with GIS and RS and it was estimated that the rate of soil erosion mainly depends on the nature of vegetation cover, overland flow and texture of the soil [16] .It was estimated that 75% area of the Dhrabi River Catchment (DRC) is impacted by soil erosion with mean rates of 82 t/hm 2 using the Revised Universal Soil Loss Equation (RUSLE) [17] and Water Erosion Prediction Project (WEPP) [18,19] models [20] .This results in a high variability of soil fertility and productivity within the area, which diminishes the storage and filtering function of the soil [20] .Soil erosion is a three step process involving detachment, transportation and deposition which causes onsite as well as offsite problems.Onsite effects are the removal of organic matter and soil nutrients which reduce agriculture production.
Offsite problems are often more severe include river silting, impaired water quality of reservoirs, reduced reservoir storage, and exacerbation of floods and landslides [21] .Due to these effects the soil and water resources are under threat and the productivity of land is decreasing which ultimately leads to a reduction in agricultural production.The accelerated soil erosion is a serious agro-environmental threat to food security and agriculture sustainability worldwide [3,22,23] .
Topography, landuse/land cover (LULC), soil type, soil structure and climatic conditions are major related factors that influence soil erosion.
Assessment of soil erosion at the catchment scale is a difficult task and is most realistically performed using available soil erosion modeling techniques including: (1) Universal Soil Loss Equation (USLE) [24] , (2) physically-based models such as WEPP, (3) and a combination of empirical and physically based methods that are used in the Soil and Water Assessment Tool (SWAT) [25] .These tools can be very much helpful in implementing the planning and management of soil and water conservation projects.
SWAT exemplifies a compromise between empirical and physical algorithms; e.g., a modified version of USLE (MUSLE) [26] that is used to simulate water erosion.Furthermore, it is considered a more suitable tool for agricultural management practices in watersheds, compared with other models [27] .It was developed in the early 1990s to assist water resource managers in assessing the impact of management and climate on water supplies and non-point source pollution in watersheds and large river basins [28] , and can be used in small agricultural watersheds to simulate water and soil loss [29][30][31][32][33] .SWAT is a watershed-scale ecohydrological model that has been tested for a wide variety of watershed scales and environmental conditions worldwide [34][35][36][37][38][39] , and has been applied for an extensive range of climate change, landuse change, Best Management Practice (BMP) assessments and other scenario analyses.For example, Mosbahi et al. [40] used SWAT for soil erosion risk assessment to adopt the appropriate management intervention and recommended it for prioritization of vulnerable areas in semi-arid catchments.At present, roughly 3500 SWAT-related studies that were published between 1993 and 2018 have been identified in the literature [41] .
Samad et al. [42] applied SWAT model for the assessment of sediment yield in Rawal Dam catchment of Pothwar region.The effectiveness of the model was evaluated with Nash and Sutcliffe coefficient (0.79).The hydrological modeling of Simly Dam watershed located in Soan River basin Pothwar region was studied in 2015 using GIS and SWAT.The estimated water balance results revealed that the inflow was successfully reproduced with Coefficient of Determination (R 2 ) of 0.75 and the author recommended that SWAT can be used efficiently in semi-arid regions to support water management policies [43] .To the best of our knowledge, the application of SWAT reported here is the first time that it has been applied for DRC subcatchments.Thus the specific objectives of this study were to: (1) test SWAT for four small subcatchments that are representative of typical DRC conditions, and (2) to evaluate soil erosion related parameters used in SWAT in the context of subcatchments that have been treated with terraces versus other subcatchments that do not have terraces.

Study area description
This study was conducted for four DRC subcatchments which are referred to as subcatchment-25, subcatchment-27, subcatchment-31 and subcatchment-32 (Figure 1).The DRC drains an area of 196 km² between latitudes 32°42ʹ36″N to 32°55ʹ48″N and longitudes 72°35ʹ24″E to 72°48ʹ36″E in District Chakwal, Pothwar, Pakistan.Precipitation is the main source of freshwater in the DRC.The undulating and uneven topography has deep to shallow gullies, large to small terraces and low to medium hills between elevations of 465-919 masl.The slope steepness varies from 2% in areas characterized by relatively flat plain conditions to >30% along the hillsides.The dominant soil type is a sandy loam that has low (<1%) organic matter.Generally, the climate is hot in the summer season and cold during the winter.The summer season extends from April to September, with the highest temperatures occurring during June and July (30°C-35°C).The winter season spans the months of October to March, with the coldest temperatures occurring in December and January (0°C-5°C).About 65%-70% of the annual precipitation occurs during the monsoon season (July to September) while the remaining 30%-35% of the annual precipitation occurs during the winter season (December to March).
The average annual precipitation is about 630 mm [15] .
The major landuse classifications of this area are: Agricultural Land (22%; 43 km 2 ), Barren Land with Shrubs and Bushes (32%; 62 km 2 ), Fallow/Range Land with Range Grasses (33%; 65 km 2 ), Residential Areas (4%; 9 km 2 ), Water Bodies (3%; 7.0 km 2 ) and Forests (6%; 11 km 2 ).The location of the DRC in Pakistan and the locations of the four DRC subcatchments are shown in Figure 1.In the Pothwar Plateau region, the agriculture fields are not flat and the crops are grown on terraces that consist of wide and deep gullies.The field terraces are situated at different elevation levels as shown in Appendix A (Figure A1).The terrace systems typically fail due to breaching of field embankments/bunds when intense precipitation events occur (Figure A2).Locally loose stone structures have thus been installed in clusters in the upper, middle, and lower parts of terraced catchments to help mitigate this problem (Figures A3-A5).
Subcatchments-25 and -31 are characterized by eroded gullies while terraces have been installed on subcatchments-27 and -32 (Figure 2).The gullies are visible with deep and wide beds while the terraces are installed in incised gullies with an average vertical interval of about 0. A physical topographical survey of the subcatchments was conducted using a Global Positioning System (GPS; eTrex Venture® HCx with accuracy <10 m, 95% typical).A DEM of each subcatchment was generated as a function of point source elevation data in GIS and the Inverse Distance Weighted (IDW) method [44,45] as shown in Figure 3.The IDW was used because it is popular deterministic method, and has been used at varying spatial and temporal scales because of simplicity.The landuse maps were generated using Google Earth [46] as shown in Figure 4.The landuse classifications of each subcatchment with the corresponding percentages of landuse area are given in Table 1.The barren land with shrubs and bushes are used for grazing in all four subcatchments, while the fallow/range land are converted and cultivated during the cropping season in catchment-27 and -32.The salient features of the study subcatchments are discussed below, which are fully representative of the DRC area and have well defined boundaries.
Subcatchment-25: The catchment contains a deep gully having a wide gully bed.The upslope portion of the gully has bushes, scrub trees and range grasses.The main crop grown during the winter season on the agricultural area is wheat.Subcatchment-25 has a total area of 2.0 hm 2 with an average slope of 10.5% and elevation range of 527.1 m to 539.8 m.
Subcatchment-27: This subcatchment has an area of 3.0 hm 2 , with a total of 11 terraces that have been installed across a gully incision.The average vertical interval is about 0.5 m between the terraces.Arable crops such as winter wheat are grown on the terraces.The subcatchment elevation ranges between 528.3 m to 540 m with an average slope of 5.8%.
Subcatchment-31: The subcatchment is characterized by a gully incision with grass growing on the gully slopes.The area is 1.5 hm 2 with an average slope of 10%.The elevation is between 509.15 m to 520.67 m.
Subatchment-32: This subcatchment has an area of 3.3 hm 2 and 7.6% average slope.The incised gully bed has been modified with 13 terraces, with an average vertical interval of about 0.5 m between the terraces, which are used for growing arable crops such as sorghum and millet mixed fodder and wheat during the winter season.The elevation ranges from 513.22 to 518.58 masl.

Data required and collection
Three types of data were required for modeling sediment yield: (1) spatial/raster data including DEM, masked DEM, landuse, soil and slope data, and (2) daily meteorological data including precipitation, temperature (maximum and minimum) in a lookup table, and (3) observed runoff and sediment data.
Meteorological, measured runoff and measured sediment data were obtained from the Soil and Water Conservation Research Institute (SAWCRI), District Chakwal [47] .The subcatchments were monitored for three years during the 2009 to 2011 time period.SAWCRI installed an automatic weather station and several recording rain gauges to collect metrological data (Figure 1).An experimental setup was constructed at the outlet of each subcatchment for the measurement of runoff and sediment yield as shown in Figure 5.The runoff discharge measurement was performed using a sharp crested rectangular weir and an automatic water-level recorder was installed to measure the runoff depth.A settling basin was used for sediment collection.After the runoff event, the trapped sediment in the settling basin was collected, air dried and weighed.Further description of the runoff and sediment yield measurement processes are given in Iqbal et al. [48] Figure 5 Experimental setup for runoff and sediment yield

SWAT Model description
SWAT is a comprehensive, semi distributed, physically-based basin scale hydrological model used to predict the impact of landuse and agricultural management practices on water and sediment yield in large as well as small watersheds over long durations of time [25,49] .
The ArcGIS SWAT (ArcSWAT) interface [50,51] uses geographic information systems (GIS) spatial algorithms to spatially link multiple model input data, such as catchment topography (DEM), soil, land use, land management and climatic data.The delineation of a catchment into sub-catchments in ArcSWAT is performed on the basis of drainage area and topography.
Typically, each sub-catchment is then further sub-divided into hydrological response units (HRUs) based on spatial uniformity in landuse, soil type and slope.The model computes surface runoff for each hydrologic response unit (HRU) by using the modified USDA-SCS curve number method [52] or Green and Ampt infiltration method [53] .Sediment yield from each HRU is estimated with the Modified Universal Soil Loss Equation (MUSLE) [26] as given in the following mass balance equation: 0.56 .
11.8( ) where, S.Y = Sediment yield, t/hm 2 ; Q surf = Surface runoff, mm/hm 2 ; q peak = Peak runoff rate, m 3 /s; area hru = Area of hydrological response unit, hm 2 ; K USLE = USLE Soil erodibility factor; C USLE = USLE cover and management factor; P USLE = USLE support practice factor; LS USLE = Slope length and slope steepness factor; CFRG = Coarse fragment factor.A detailed description of these and other components are available in the SWAT theoretical documentation [49] .

SWAT Model setup and simulation
During the watershed delineation process in ArcSWAT, each of the four subcatchments was maintained as a single subbasin, because of the relatively small size of each subcatchment (from 1.5 hm 2 to 3.3 hm 2 ).Then, each of four subcatchments was discretized into HRUs, resulting in the following number of HRUs for subcatchments-25, -27, -31 and -32, respectively: 8, 12, 6 and 14.An appropriate database of subcatchment parameters and a comprehensive topographic report of the catchment are generated upon successful execution of the terrain processing module within the ArcSWAT interface.The respective look-up tables, land use and soil maps were supplied to the model for reclassification according to SWAT coding conventions.
Further, the subcatchment was classified into different slope categories and HRUs with unique land cover, soil and slope classes by overlaying all three maps with the HRU threshold set at 0% to allow all of the original landuse, soil and slope data to be included in the creation of HRUs.This was necessary to prevent the loss of any of the input data, which occurs when the HRU development thresholds are set at a non-zero percentage [54] , and because SWAT pollution outputs are sensitive to the resolution of HRU and subcatchment delineations [34,54] .
The weather station location table, and tables of daily precipitation and temperatures (maximum and minimum) data, were incorporated into the model setup and linked with the required input files.The surface runoff and sediment yield predicted by the un-calibrated model were initially compared versus different rainfall and sediment loss events on a daily basis using default parameter values.The general algorithm used for the sediment yield simulation for the four subcatchments is given in Figure 6.The initial un-calibrated simulation was saved as Sim1 which was further modified for model calibration.

SWAT calibration and validation
The determination of the most sensitive parameters for model calibration and validation was considered as the first step in this study [55] .The user determines which variables to adjust based on sensitivity analyses, literature data, local expertise and other relevant sources of information [56] .Sensitivity analysis is the rate of change in model output values with respect to changes in model input parameters.For determining the most sensitive parameter for model calibration the sensitivity analysis was performed in the ArcSWAT interface using five parameters for sediment yield [49] : P USLE , C USLE , K USLE , linear parameter for calculating the maximum amount of sediment that can be re-entrained during channel sediment routing (SPCON), and the exponent parameter for calculating sediment re-entrained in channel sediment routing (SPEXP).The second step is the calibration process.The model calibration was performed with adjustments of key model parameters until the model output matched as closely as possible with the observed data.An iterative approach is usually used for manual calibration involving the following steps: (1) perform the simulation; (2) compare observed and simulated values; (3) assess if reasonable results were obtained; (4) if not, further adjust the input parameters based on expert judgment and other guidance within reasonable parameter value ranges; and (5) repeat the process until it is determined that the best results have been obtained [55] .The manual calibration for 2009 to 2010 was performed using this five-step iterative approach and the procedure given in Figure 7 for the four subcatchments.
The final step is the model validation which was done by comparing the model results with observed data using climate data inputs for 2011.No additional parameter adjustments were performed during the validation period.Figure 7 SWAT manual calibration flowchart used for surface runoff and sediment yield (from Engel et al. [57] ; adapted from Santhi et al. [58] )

Model Performance Evaluation
Agreement between observed and simulated values are commonly assessed by using efficiency criteria such as the R 2 , Nash Sutcliffe Modeling Efficiency (E NS ) and Index of Agreement (d) [59] .The SWAT calibration and validation results were evaluated using the E NS and R 2 , which are the most widely used statistical methods as previously reported by Gassman et al. [34,35] and Bressiani et al. [39] The model performance was judged based on performance evaluation criteria (PEC) for catchment-scale models suggested by Moriasi et al. [60] their criteria state that satisfactory performance occurs on a daily, monthly and annual basis for a flow simulation if R 2 >0.6 and E NS >0.5, and for a sediment prediction if R 2 >0.45 and E NS >0.45.

Nash-Sutcliffe coefficient (E NS )
The E NS coefficient [59,61] is a normalized statistic that indicates the efficiency of a model by relating the goodness-of-fit of the model's predictions to the variance of the observed data.The E NS provides a measure how well the simulated results match the observed data along a 1:1 line: where, E NS is the Nash-Sutcliffe coefficient; Q oi is the observed data; Q si is the simulated data and 0 Q is the average of observed data.The E NS ranges from -∞ to 1. E NS = 1 indicates a perfect match between the simulated and observed data, E NS = 0 indicates that mean of the observed data is as accurate the simulated data, and E NS < 0 occurs when the observed mean is a more accurate predictor than the simulated data.E NS is sensitive to extreme values due to squared differences [59] .

R 2
The coefficient of determination is computed as: where, R 2 is the coefficient of determination; Q si is the simulated data; s Q is the average of simulated data; Q oi is the observed data, and 0 Q is the average of observed data.The R 2 is a measure of strength of linear correlation between observed and predicted outcomes by the model [62] .The R 2 statistics ranges from 0 to 1. R 2 = 0 means there is no correlation while a value of 1 indicates that there is a perfect correlation between the observed and simulated data.The R 2 statistics provides an estimate of how well the variance of observed values are replicated by the model predictions [59] .

Results and discussion
The calibration of SWAT was successfully performed for surface runoff and sediment yield (Figures 8 to 11) based on literature guidance and the calibration techniques described in the SWAT user manual.The key parameter adjustments that were performed for calibrating the surface runoff included the runoff curve number (CN2 = 65), average slope length (SLSUBBSN = 60) and average slope steepness (HRU_SLP = 0.016).The default and final values of the parameters that were adjusted during the sediment yield calibration, and their relative ranking in terms of sensitivity, are given in Table 2.It was observed that the values of the parameters selected during the model calibration process were in the range of default values.The P USLE factor was found to be the most sensitive as compared to the other parameters.The value of soil erosion parameters used during calibration were similar to those recommended by Klik et al. [20] , who used the RUSLE and WEEP models to estimate the average annual soil loss in the DRC.The K, LS, and C factors values were comparable according to soil type, topography and vegetation cover (Table 2), respectively.

Table 2 Soil erosion parameter ranking used for model
calibration [63] Value The statistical evaluation (R 2 , E NS ) of the model performance for surface runoff and sediment yield is given in Table 3.The resulting statistics for the four small DRC subcatchments all met the satisfactory criteria, and several of the statistics exceeded the good or very good criteria, as suggested by Moriasi et al. [60] Furthermore, the high R 2 values indicate that there is a strong correlation between the observed and simulated surface runoff and sediment yield levels, while the E NS values show that the observed versus predicted runoff and sediment yield plots fit the 1:1 line well.The erratic and intensive nature of rainfall has a profound impact on the generation of several peak runoff events during the monsoon (rainy) season which leads to severe soil erosion on steep sloped areas of the DRC.It was observed that 63% of annual rainfall occurred in the monsoon season (June to September).In SWAT, the erosive impact of rainfall is generally estimated in term of peak runoff generation, so the results obtained during calibration and validation are represented in Figures 8-11 for surface runoff and sediment yield for all four subcatchments.The analysis was performed for each total rainfall event and the respective total surface runoff and sediment yield generated by each event.
The overall statistical results indicated that the performance of SWAT was satisfactory for all of the subcatchments and showed the simulated values generally matched the corresponding observed values well.
However, model adequacy should be further evaluated by how well the model captures high and low rainfall events, specifically regarding the replication of fluctuations in the resulting hydrographs and sediment yields.The graphical results (Figures 8-11) revealed that SWAT was able to satisfactorily reproduce most of the low flow and sediment yield events (due to low rainfall events) although some relatively low sediment yields were considerably overpredicted; e.g., sediment yield events on August 7, 2011 (Figure 8) and December 8, 2011 (Figures 9 and  11).In contrast, it was also found that SWAT typically underestimated or overestimated high flow and sediment yield events, in response to high rainfall events.For example: (1) a maximum intensity rainstorm on July 29, 2010 resulted in overestimation of surface runoff and sediment yields for all four subcatchments, and (2) another maximum intensity rainstorm on July 29, 2009 which resulted in underestimations for subcatchment-25 (Figure 8) versus overestimations for the other subcatchments (Figures 9-11).These discrepancies may occur due to inaccuracies in observed climate, runoff and sediment data, such as some of the rainfall events not being measured properly which in turn leads to underestimations or overestimations of runoff peaks.Another possible reason could be related to short, rapid rainfall events, which can lead to an overestimation discrepancy because small subcatchments have low times of concentration and thus a low capacity to minimize peak runoff.Also the curve number (CN) technique cannot accurately predict runoff for days that experience several storms.The underestimation and/or overestimation of sediment yield was also observed during high intensity rainstorm events which may be due to uncertainties in runoff simulation measurements as well as uncertainties in model parameterization.This may be also due to observed data used for model calibration and validation.Relatively short term events with several storms having high intensity may be not captured well by sampling of sediment data, including inaccurately high loads measured during short term events which leads to an overestimation in sediment yield.

Additional analyses of specific rainfall storm events
The detailed description of runoff generation and soil loss (sediment yield) versus the corresponding erosive events are given in Table 4 for the calibration and validation periods.It was observed that rainfall-runoff events occurred from April to September.Subcatchments-25 and -31 are natural gully systems with no engineering or vegetation protection provided by the farmers against surface runoff and soil erosion.As a result, the entire respective catchment contributed to runoff and soil erosion during most of the storms.In contrast, only the lower part of the catchments contributed to surface runoff and sediment yield in subcatchments-27 and -32 due to the terraced systems in combination with modified gentle slopes.The cumulative rainfall of 928.3 mm during the calibration period produced a total surface runoff and sediment yield of 225 mm and 44.33 t/hm 2 for subcatchment-25 (Figure 8).However, the same rainfall amount for subcatchment-27 resulted in much less total surface runoff and sediment yield of 159 mm and 8.77 t/hm 2 (Figure 9), due to the terrace system that was installed in subcatchment-27.Similarly, the total rainfall of 808.06 mm during the calibration period generated overall surface runoff and sediment yields of 159.77 mm and 51.02 t/hm 2 for subcatchment-31 (Figure 10), versus surface runoff and sediment yields of 54.58 mm runoff and 13.89 t/hm 2 for subcatchment-32 (Figure 11).During the analysis, it was observed that the maximum rainstorm events occurred in July and August which in turn resulted in peak surface runoff and maximum sediment yields.Therefore, the impact of maximum intensity rainstorms on surface runoff and sediment yield was also analyzed (Table 5).For subcatchments-25 and -27, the maximum intensity rainstorm occurred on July 29 in both 2009 and 2010.The maximum 30 minute-intensity (I30) of 84.3 mm/h on July 29, 2009 produced peak surface runoff and sediment yield amounts of: (1) 46.2 mm and 6.9 t/hm 2 for subcatchment-25, and (2) 18.36 mm and 1.1 t/hm 2 for subcatchment-27.On July 29, 2010, the maximum intensity storm of 64.1 mm/h generated 25.9 mm/31.8mm surface runoff and 7.8 t/hm 2 /1.4 t/hm 2 sediment yield for subcatchments-25/-27, respectively (Figures 8  and 9).These sediment yield results again reflect the effects of the terraces installed in subcatchment-27.The maximum intensity storm that occurred on August 12, 2011 generated the lowest total rainfall (39.6 mm), and resulted in correspondingly low peak surface runoff (7.5 mm and 3.4 mm) and sediment yields (0.6 t/hm 2 and 0.06 t/hm 2 ) for subcatchments-25 and -27 (Table 5), as compared to the impacts of the storms that occurred on July 29 during the previous two years during the calibration period.Maximum intensity rainstorm events were recorded on the same dates for subcatchments-31 and -32, but the rainfall totals and maximum 30-minute intensity levels were substantially different than those found for subcatchments-25 and -27 (Table 5).The highest peak surface runoff and sediment yield amounts that occurred in response to the storm event on August 12, 2011 for subcatchments-31/-32 resulted in surface runoff levels of 27.3 mm/15.9mm and sediment yields of 8.4 t/hm 2 /4.3 t/hm 2 , respectively.The results of the three storms reported for subcatchments-31 and -32 (Table 5) further confirm the effectiveness of terraces, which are installed in subcatchment-32, in reducing surface runoff and sediment yield.
The best fit correlations between the three components of rainfall, runoff and sediment yield for the gully subcatchments and terraced subcatchments were determined separately as shown in Figure 12 using observed and simulated results.
The correspondence of rainfall-runoff hydrographs with sediment yield is an indication of SWAT's capability to simulate the hydrological regime of the DRC subcatchments.A linear correlation was the best fit between rainfall versus both the observed and simulated runoff, resulting in respective R 2 values of 0.76 and 0.89, for the gully subcatchments.
However, a power correlation was determined between rainfall versus observed sediment yield (R 2 = 0.59) and rainfall versus simulated sediment yield (R 2 = 0.58) as shown in Figure 12.Likewise, power relationships between rainfall and runoff, and rainfall and sediment yield, were found to result in the strongest correlations for the terraced subcatchments, resulting in R 2 values from to 0.63 to 0.76.Based on this correlation technique, the basic regression equations were developed for the DRC gully subcatchments and terraced subcatchments are as follows: For  3.11 ; R 2 =0.76 (9) It was observed that rainfall storm events <10-15 mm did not produce surface runoff, and thus no soil erosion resulted for these small rainfall events.A 16 mm rainfall event produced 0.025 t/hm 2 sediment yield while rainfall amounts that ranged from 150 mm to 250 mm produced 5.0-6.0 t/hm 2 sediment yield across the subcatchments.Overall, the regression analyses indicated good correlation with field and simulated data.

Conclusions and recommendations
In this research, the SWAT model was applied on a daily time step and was used to simulate soil erosion for specific storm events for four small DRC subcatchments.Subcatchments-25 and -31 were characterized by incised gullies while subcatchments-27 and -32 consist of terraced landuse systems.The performance of the SWAT model was satisfactory with R 2 values ≥0.67 and E NS ≥0.53 for both surface runoff and sediment yield during the calibration (2009-2010) and validation (2011) periods.The P USLE factor was found to be the most sensitive factor during model calibration.It was observed that all of the rainfall-runoff events occurred during the monsoon season (June to September).
It was estimated that in 2009, the gully subcatchments produced annual sediment yield that ranged from 13.19 t/hm 2 to 20.60 t/hm 2 in response to rainfall events that were between 230.76-399.95mm and surface runoff amounts that ranged from 49.34 mm to 95.59 mm.
In comparison, the terraced subcatchments produced 1.90-3.77t/hm 2 sediment yield for the same rainfall levels and surface runoff that was measured from 10.17 mm to 54.0 mm.In 2010, the annual sediment yield for the gully subcatchments was observed to be 30.41-31.12 t/hm 2 , due to rainfall levels of 528.29-577.30mm that produced surface runoff amounts of 110.43-129.53mm.In contrast, it was observed that the terraced subcatchments yielded 0.87-10.11t/hm 2 of sediment annually, in response to rainfall amounts of 432.3-577.30mm that resulted in surface runoff that ranged from 44.4 mm to 110.43 mm.2011 was an interesting climate year of reference, because it was relatively dry year which resulted in rainfall amounts that ranged from 261.6 mm to 403.88 mm for the gully subcatchments and 194.8 mm to 403.88 mm for the terraced subcatchments.As a result, surface runoff of 28.34-82.84mm and sediment yield of 2.59-17.73t/hm 2 occurred for the gully subcatchments while the terraced catchments produced runoff of 23.61-30.47mm and corresponding sediment yield levels of 0.52-4.19t/hm 2 .It was found that the sediment yield generated by the terraced subcatchments was 4-5 times lower than the sediment yields that were exported from the gully subcatchments.Thus, there is a huge potential for terraces to reduce the soil erosion in the DRC.It is recommended that the model application should be extended for larger catchments located within the Pothwar area, to identify efficient locations for soil and water conservation structures and management practices.

Figure 1
Figure 1 Location of study catchments within the Dhrabi River Catchment (DRC), and the location of the DRC within the Pothwar area and Pakistan 5 m.Subcatchments-25 and -27 are adjacent with each other having latitude 32.8948 o to 32.8917 o and longitude 72.70 o to 72.71 o .Subcatchments-31 and -32 are also adjacent with latitude 32.9188 o to 32.9159 o and longitude 72.7109 o to 72.7100 o .The soil type is sandy loam (67%-74% sand, 14%-22% silt and 10%-14% clay) in all four subcatchments and is calcareous in nature.

Figure 6
Figure 6 General algorithm used for sediment yield simulation for the four subcatchments

Figure 12
Figure 12 Correlation between rainfall, runoff and sediment yield for the gully/ terraced subcatchments, for both the observed and simulated data Figure A1 Terraced cultivated lands in Pothwar

Figure
Figure A5 Front view of conservation structure built using local available stones