10 1 1 512 1747 PDF

Remote Sensing and hydrological modeling
of the Upper Bhima catchment
Draft
W.W. Immerzeel
A. Gaur
P. Droogers
Research paper no. 3

P.O. Box 35
2600 AA Delft
The Netherlands
Kluyverweg 1
2629 HS Delft
The Netherlands
Tel. +31 (0)15 278 80 25
Fax +31 (0)15 262 30 96
Email [email protected]
www.nivr.nl
NATIONAL USER SUPPORT PROGRAMME (NUSP)

2001-2005
http://www.ao-go.nivr.nl
The National User Support Programme 2001-2005 (NUSP) is executed by the Netherlands
Agency for Aerospace Programmes (NIVR) and the Space Research Organization of the
Netherlands (SRON). The NUSP is financed from the national space budget. The NUSP
subsidy arrangement contributes to the development of new applications and policy-
supporting research, institutional use and use by private companies.
The objectives of the NUSP are:
To support those in the Netherlands, who are users of information from existing and
future European and non-European earth observation systems in the development of
new applications for scientific research, industrial and policy research and operational
use;
To stimulate the (inter)national service market based on space-based derived

operational geo-information products by means of strengthening the position of the
Dutch private service sector;
To assist in the development of a national Geo-spatial data and information

infrastructure, in association with European and non-European infrastructures, based
on Dutch user needs;
To supply information to the general public on national and international space-based

geo-information applications, new developments and scientific research results.
Remote Sensing and hydrological
modeling of the Upper Bhima
catchment
Authors: W. Immerzeel (FutureWater)

A. Gaur (IWMI)
P. Droogers (FutureWater)
This research was funded by: NIVR/SRON/NWO User Support Programme 2
Generaal Foulkesweg 28 6703 BS Wageningen Tel: 0317-460050

E-mail: [email protected] Website: www.futurewater.nl
Remote Sensing and hydrological modelling of the Upper Bhima catchment June 2006
2/69 FutureWater Science for Solutions

June 2006 Remote Sensing and hydrological modelling of the Upper Bhima catchment
Preface
This report contributes to a research project undertaken by FutureWater entitled: Remotely Sensed
based hydrological model calibration for basin scale water resources planning: embedding case for
Krishna Basin, India (GO-2005/025). This project is financially supported by NIVR (Nederlands
Instituut voor Vliegtuigontwikkeling en Ruimtevaart) in the context of Tijdelijke subsidieregeling
Nationaal Programma Gebruikers Onderstuining (GO).
The report describes the use of Remote Sensing in hydrological modeling of the Upper Bhima
catchment.
FutureWater
Generaal Foulkesweg 28
6703 BS Wageningen
tel: 0317 460050
email: [email protected]
FutureWater Science for Solutions 3/69


Table of contents
1 INTRODUCTION 7
2 SOIL AND WATER ASSESSMENT TOOL 9

2.1 Background 9
2.2 Evaporation 10
2.3 Crop growth 11
2.4 Irrigation 11
2.5 Groundwater 12
2.6 Reservoirs 13
3 LAND USE CLASSIFICATION 15

3.1 Introduction 15
3.2 Remote Sensing and land use classification 16
3.3 MODIS 18
3.4 Classification of the Upper Bhima sub basin 19
4 EVAPOTRANSPIRATION MAPPING 23
4.1 SEBAL 23
4.2 Reference evapotranspiration 28
4.3 Potential Evapotranspiration 29
4.4 Actual Evapotranspiration 30
4.5 Biomass production 32
5 PRECIPITATION MAPPING USING REMOTE SENSING 33

5.1 Introduction 33
5.2 Tropical Rainfall Monitoring Mission 34
5.3 Downscaling and calibration of TRMM data 35
6 MODELING THE UPPER BHIMA CATCHMENT 39

6.1 Topography 39
6.2 Soils 40
6.3 Hydrological Response Units 42
6.4 Meteorology 43
6.4.1 Solar radiation 43
6.4.2 Temperature 44
6.4.3 Relative humidity 45
6.4.4 Wind speed 45
6.4.5 Precipitation 45
6.5 Management practices 46
6.6 Reservoirs 48
7 SENSITIVITY ANALYSIS 51
7.1 Exploration of reference evapotranspiration 51
7.2 Differences in acual evapotranspiration 53

8 CALIBRATION 57
8.1 Results 57
8.2 Discussion and conclusions 61
9 RESULTS 64
9.1 Water balances 64
9.2 Biomass and crop production 64
9.3 Scenarios 64
10 CONCLUSION 66
11 REFERENCES 68

1 Introduction


2 Soil and Water Assessment Tool
2.1 Background
SWAT1 was developed primarily by the United States Department of Agriculture (USDA) to predict the
impact of land management practices on water, sediment and agricultural chemical yields in large
complex watersheds with varying soils, land use and management conditions over long periods of
time. The SWAT model has been extensively used, is in the public domain and can be considered as
becoming the de-facto standard in spatial decision support systems.
SWAT represents all the components of the hydrological cycle including: rainfall, snow, snow-cover and
snow-melt, interception storage, surface runoff, up to 10 soil storages, infiltration, evaporation,
evapotranspiration, lateral flow, percolation, pond and reservoir water balances, shallow and deep
aquifers, channel routing. It also includes irrigation from rivers, shallow and deep groundwater stores,
ponds/reservoirs and rivers, transmission losses and irrigation onto the soil surface. It includes
sediment production based on a modified version of the Universal Loss Equation and routing of
sediments in river channels. SWAT tracks the movement and transformation of several forms of
nitrogen and phosphorus in the watershed. It also tracks the movement and decay of pesticides. All
include channel routing components and carrying of pollutants by sediments. SWAT has a modular set-
up and it goes beyond the scope of this report to get into detail on each of these modules, but
reference is made to the theoretical documentation (Neitsch et al, 2001).
For modeling purposes, a watershed may be partitioned into a number of sub-watersheds or sub-
basins. The use of sub-basins in a simulation is particularly beneficial when different areas of the
watershed are dominated by land uses or soils dissimilar enough in properties to impact hydrology. By
partitioning the watershed into sub-basins, the user is able to reference different areas of the
watershed to one another spatially. Input information for each sub-basin is grouped or organized into
the following categories: climate; hydrologic response units or HRUs; ponds/wetlands; groundwater;
and the main channel, or reach, draining the sub-basin. Hydrologic response units are lumped land
areas within the sub-basin that are comprised of unique land cover, soil, and management
combinations.
No matter what type of problem studied with SWAT, the water balance is the driving force behind
everything that happens in the watershed. To accurately predict the movement of pesticides,
sediments or nutrients, the hydrologic cycle as simulated by the model must conform to what is
happening in the watershed. Simulation of the hydrology of a watershed can be separated into two
major divisions. The first division is the land phase of the hydrologic cycle, depicted in Figure 1. The
land phase of the hydrologic cycle controls the amount of water, sediment, nutrient and pesticide
loadings to the main channel in each sub-basin.
1
http://www.brc.tamus.edu/swat/index.html

Figure 1: Water balance SWAT.
The second division is the water or routing phase of the hydrologic cycle which can be defined as the
movement of water, sediments, etc. through the channel network of the watershed to the outlet.
Once SWAT determines the loadings of water, sediment, nutrients and pesticides to the main channel,
the loadings are routed through the stream network of the watershed using a command structure. In
addition to keeping track of mass flow in the channel, SWAT models the transformation of chemicals in
the stream and streambed.
2.2 Evaporation
Evapotranspiration is a collective term for all processes by which water in the liquid or solid phase at or
near the earth's surface becomes atmospheric water vapor. Evapotranspiration includes evaporation
from rivers and lakes, bare soil, and vegetative surfaces; evaporation from within the leaves of plants
(transpiration); and sublimation from ice and snow surfaces. The model computes evaporation from
soils and plants separately as described by Ritchie (1972). Potential soil water evaporation is estimated
as a function of potential evapotranspiration and leaf area index (area of plant leaves relative to the
area of the HRU). Actual soil water evaporation is estimated by using exponential functions of soil
depth and water content. Plant transpiration is simulated as a linear function of potential
evapotranspiration and leaf area index. Potential evapotranspiration is the rate at which
evapotranspiration would occur from a large area completely and uniformly covered with growing
vegetation which has access to an unlimited supply of soil water. This rate is assumed to be unaffected
by micro-climatic processes such as advection or heat-storage effects. The model offers three options

for estimating potential evapotranspiration: Hargreaves (Hargreaves et al., 1985), Priestley-Taylor

(Priestley and Taylor, 1972), and Penman-Monteith (Monteith, 1965).
2.3 Crop growth
For each day of simulation, potential plant growth, i.e. plant growth under ideal growing conditions is
calculated. Ideal growing conditions consist of adequate water and nutrient supply and a favorable
climate. The biomass production functions are to a large extend similar to SEBAL; First the Absorbed
Photosynthetical Radiation (APAR) is computed from intercepted solar radiation as a function of LAI,
followed by a Light Use Efficiency (LUE) that is in SWAT essentially a function of carbon dioxide
concentrations and vapor pressure deficits. Whilst LAI is simulated in SWAT, the fractional
Photosynthetical Active Radiation (fPAR) is measured in SEBAL. The crop yield is computed as the
harvestable fraction of the accumulated biomass production across the growing season.
Figure 2: Parameterization of crop production processes for estimating crop yiel
2.4 Irrigation
Irrigation may be scheduled manually or applied automatically by the model as response to a water
deficit in the soil. Irrigation water applied to a sub basin is obtained from one of five types of water
sources: a reach, a reservoir, a shallow aquifer, a deep aquifer, or a source outside the watershed. For
this study it is assumed that all water originated from the shallow aquifer. If automatic irrigation is
applied all soil layers are filled up to field capacity. If manually scheduling is used all scheduled water is
applied and potential excess water percolates to the shallow aquifer. In this study irrigation water is
applied automatically based on a predefined water stress criterion per sub basin. Water stress is 0.0
under optimal water conditions and approaches 1.0 as the soil water conditions vary from the optimal.
Water stress is simulated by comparing actual and potential plant transpiration:
Et ,act wactualup
wstrs = 1 =1
Et Et

where wstrs is the water stress for a given day, Et is the maximum plant transpiration on a given day
(mm H2O), Et,act is the actual amount of transpiration on a given day (mm H2O) and wactualup is the
total plant water uptake for the day (mm H2O). The water stress criterion is used in calibrating
simulated Et,act with the measured SEBAL Et,act.
2.5 Groundwater
Recharge to unconfined aquifers occurs via percolation of excessively wet root zones. Recharge to
confined aquifers by percolation from the surface occurs only at the upstream end of the confined
aquifer, where the geologic formation containing the aquifer is exposed at the earths surface, flow is
not confined, and a water table is present. River courses and irrigation canals are connected to the
groundwater system, and surface water groundwater interactions are taken care for.
After water is infiltrated into the soil, it can basically leave again the ground as lateral flow from the
upper soil layer which mimics a 2D flow domain in the unsaturated zone or from return flow that
leaves the shallow aquifer and drains into a nearby river. The remaining part of the soil moisture can
feed into the deep aquifer, from where it can be pumped back by means of artificial extraction. The
total return flow thus consists of surface runoff, lateral outflow from root zone and aquifer drainage to
river.
Figure 3: Schematic diagram of the partitioning of infiltration into sub-surface water fluxes after water
uptake by roots have taken place
SWAT simulates two aquifers in each subbasin. The shallow aquifer is an unconfined aquifer that
contributes to flow in the main channel or reach of the subbasin. The deep aquifer is a confined
aquifer. Water that enters the deep aquifer is assumed to contribute to streamflow somewhere outside
of the watershed (Arnold et al., 1993). The effects of groundwater extractions on baseflow (Qgw),
defined as the contribution of the shallow aquifer to stream flow, is of specific relevance in this study.

Figure 4: Schematic representation of shallow and deep aquifers in SWAT (Neitsch et al, 2001)
Base flow calculations are based on a combination of Hooghoudt (1940) and Smedema and Rycroft
(1983) according to
[ ] [
Q gw,i = Q gw,i 1 exp gw t + wrchrg , sh (1 exp gw t ])
where Qgw,i is the groundwater flow into the main channel on day i (mm H2O), Qgw,i-1 is the
groundwater flow into the main channel on day i-1 (mm H2O), gw is the baseflow recession constant,
t is the time step (1 day), wrchrg,sh is the amount of recharge entering the shallow aquifer on day i
(mm H2O). To enable the simulation of the effect of groundwater extractions a component was added
to the model which assumes a minimal baseflow defined as a percentage of the actual amount of
water stored in the aquifer. For calculations of water table fluctuations a specific yield is assumed of
0.15 m3/m3 is assumed.
2.6 Reservoirs
Reservoirs are located within a sub basin off the main channel. Water flowing into these water bodies
must originate from the sub basin in which the water body is located. Reservoirs are located on the
main channel network. They receive water from all sub basins upstream of the water body. A
schematic representation of reservoirs in SWAT is shown in Figure 5.

Figure 5: Schematic representation reservoirs in SWAT (Neitsch et al, 2001)
The water balance for reservoirs includes inflow, outflow, rainfall on the surface, evaporation, seepage
from the reservoir bottom and diversions.
The model offers three alternatives for estimating outflow from the reservoir. The first option allows
the user to input measured outflow. The second option, designed for small, uncontrolled reservoirs,
requires the users to specify a water release rate. When the reservoir volume exceeds the principle
storage, the extra water is released at the specified rate. Volume exceeding the emergency spillway is
released within one day. The third option, designed for larger, managed reservoirs, has the user
specify monthly target volumes for the reservoir.

3 Land use classification
3.1 Introduction
Based on a study for the entire Krishna basin a land use map was available at the start of the project.
This land use map (Figure 6) was derived from a temporal series of MODIS2 Normalized Difference
Vegetation Index (NDVI) data. The first results of the evaporation mapping revealed however that the
land use map was not accurate enough, specifically in the irrigated areas and the forested areas in the
Western Ghats. The Upper Bhima sub-basin has a number of specific characteristics which caused
errors in the land use classification.
Figure 6 Land use map based on Krishna basin classification
First of all the acreage of irrigated agriculture seems to be underestimated. The command areas of the
Khadakwasla, Ujani and Yadgaon and Bhatghar reservoirs are mainly cultivated with sugarcane. Sugar
cane is a C4 crop, produces large amounts of biomass with yields up 70 t/ha. The growing season of
sugarcane has length of 11 months and sugar cane can thus be grown year around. These unique
properties might have led to an inaccurate representation of the sugarcane areas. Secondly there
seems to be large areas in and around the Western Ghats which are classified as rain fed agriculture
but in practice they should be classified as low density forests.
2
http://modis.gsfc.nasa.gov/

Because of these errors a new land use classification has been performed for the Upper Bhima sub
basin. This chapter describes the methodology used and the results of the land use classification.
3.2 Remote Sensing and land use classification
The application of Remote Sensing is a common technique in land use and land cover (LULC) mapping.
Traditionally multi-spectral data are used to classify pixels in the image as a certain LULC unit. The
basic assumption being that each LULC unit exhibits unique spectral behavior. A LULC has unique
reflectance and emittance properties at different location along the electromagnetic spectrum. This is
illustrated in Figure 7. It shows the spectra of water, bare soil and vegetation and the bands of the
Landsat TM satellite. The combination of digital numbers (DN) for each band for each pixel allows the
classification of the pixel to a certain land use class. This technique is summarized as spectral pattern
recognition.
Figure 7: Spectral behaviour of major land uses and the spectral bands of the Landsat satellite
Multi-spectral image classification is usually performed using one or a limited number of images. In
agricultural areas however there may be distinct spectral and spatial changes during the growing
season which necessitates the use of multi-temporal imagery. A single image at the onset of the
growing would prohibit the distinction between bare soil and a crop. Problems may also arise when
dual-cropping patterns are prevalent. Winter wheat in fall may show a similar spectrum as a pasture in
spring. If imagery of multiple dates are used the chances on erroneous classifications are much less
(Immerzeel and Droogers, 2005, Bastiaanssen et al., 2006)
There is a distinction between supervised and unsupervised classification techniques (Lillesand and
Kiefer, 2000).
In a supervised classification the spectral signatures of training areas with known LULC are specified to
the classification algorithm. After the training stage the actual classification process can take place.
During the classification process each pixel is compared to the specified spectral signatures and
labelled to the corresponding LULC class. There are four popular classification methods which are
commonly used:

The Box classifier method

The Minimum Distance classifier
The Minimum Mahalanobis Distance classifier
The Maximum Likelihood classifier
The Box classifier method is based on the distances towards class means and the standard deviation
per band of each class. Multi-dimensional boxes are drawn around class means based on the standard
deviation of each class. The user can insert a multiplication factor (usually > 1) to make all boxes a bit
wider. If the spectral values of a pixel to be classified fall inside a box, then that class name is
assigned. If the spectral values of a pixel fall within two boxes, the class name of the box with the
standard deviation is assigned. If the spectral values of a pixel do not fall within a box, the undefined
value is assigned. The default multiplication factor is the square root of 3 for 3 bands.
The Minimum Distance classifier is based on the Euclidean distances towards class means only. For the
spectral values of a pixel to be classified, the distances towards the class means are calculated. If the
shortest (Euclidean) distance to a class mean is smaller than the user-defined threshold, then this class
name is assigned to the output pixel. Else the undefined value is assigned.
For the spectral values of a pixel to be classified with the Minimum Mahalanobis Distance classifier, the
distances towards the class means are calculated as Mahalanobis distance, which depends on the
distances towards class means and the variance-covariance matrix of each class. The class name with
the shortest Mahalanobis distance is assigned, if this distance is smaller than the user-defined
threshold value. Else, the undefined value is assigned.
The Maximum Likelihood classification assumes that spectral values of training pixels are statistically
distributed according to a multi-variate normal probability density function. For each set of spectral
input values, the distance is calculated towards each of the classes is calculated using Mahalanobis
distance. Another factor is added to compensate for within class variability. The class name with the
shortest distance is assigned, if this distance is smaller than the user-defined threshold value. Else, the
undefined value is assigned.
In an unsupervised classification the image data are first of classified by aggregating them into spatial
clusters based on the statistical properties of the pixel values (e.g. average and variation). The number
of clusters is usually chosen larger then the actual number of LULC classes which are expected. After
the clustering the LULC identity is assigned to each land use class based on ground reference data,
which is usually not an easy task. In case of the unsupervised classification the classification step is
followed by the training step. The clusters are natural groupings of pixels with similar spectral
behavior. There are numerous clustering algorithms that can determine these natural groupings. One
common form of clustering is called the K-means approach. The user identifies the number of clusters
to be located in the image. The algorithm then randomly locates the number of specified cluster
centers in the multidimensional feature space. Each pixel is then assigned to the cluster whose
arbitrary mean vector is closest. After all pixels are classified revised mean vectors for each cluster is
calculated and the process is iteratively repeated until no significant changes in the location of the
mean vectors occur. The convergence criterion used to stop the iteration is specified by the user.
Other commonly used methods are based on textural information of the image. They use information
from neighboring pixels to assign classes to pixels and are referred to as contextual classifiers.

Spectral rationing is another important technique that can support image classifications. It refers to an
enhancement technique by combining pixel values from different bands. The most commonly used
ratio techniques are ratios to enhance the vegetation cover such as the Normalized Difference
Vegetation Index (NDVI). The NDVI is defined as:
NIR RED
NDVI = Eq. 1
NIR + RED
Where NIR is the reflectance in the near infra-red part of the spectrum (wavelength = 0.7-1.1 m) and
RED is the reflectance in the red part of the spectrum (wavelength = 0.6-0.7 m). Healthy vegetation
free of environmental stress reflects well in the NIR part of the spectrum and reflects poorly in the red
part of the spectrum. The NDVI amplifies this effect.
In this study a combination of spectral rationing and multi-temporal unsupervised classification is

applied to make a land use classification of the Upper Bhima sub basin. A time series of 16 NDVI
images acquired with the Moderate Resolution Imaging Spectroradiometer (MODIS) is used to perform
an unsupervised classification.
3.3 MODIS3
MODIS is a key instrument aboard the Terra (EOS AM) and Aqua (EOS PM) satellites. Terra's orbit
around the Earth is timed so that it passes from north to south across the equator in the morning,
while Aqua passes south to north over the equator in the afternoon. Terra was launched on December
18 1999 and Aqua was launched on May 4 2002. These data will improve the understanding of global
dynamics and processes occurring on the land, in the oceans, and in the lower atmosphere. MODIS is
playing a vital role in the development of validated, global, interactive Earth system models able to
predict global change accurately enough to assist policy makers in making sound decisions concerning
the protection of our environment.
The MODIS instrument provides high radiometric sensitivity (12 bit) in 36 spectral bands ranging in
wavelength from 0.4 m to 14.4 m. The responses are custom tailored to the individual needs of the
user community and provide exceptionally low out-of-band response. Two bands are imaged at a
nominal resolution of 250 m at nadir, with five bands at 500 m, and the remaining 29 bands at 1 km. A
55-degree scanning pattern at the EOS orbit of 705 km achieves a 2,330-km swath and provides
global coverage every one to two days.
In this study the MOD02 product acquired with the Aqua satellite is used. The Level 1B data set
contains calibrated and geolocated at-aperture radiances for 36 bands generated from MODIS Level 1A
sensor counts (MOD 01). The radiances are in W/(m2 m sr). The NDVI is calculated using band 1
(620-670 nm) and band 2 (841-876 nm) which both have a spatial resolution of 250 meter:
band 2 band1
NDVI MODIS =
band 2 + band 1
3
http://modis.gsfc.nasa.gov/index.php

3.4 Classification of the Upper Bhima sub basin
In this study a time series of 16 cloud free MODIS AQUA NDVI images is used from October 2004 to
May 2005. Due to the monsoon cloud free images from June to September could not be obtained. The
dates and time of acquisition of the images used are shown in Table 1.
Image # Date Acquisition time (GMT)

Hour Minute
1 16/10/2004 8 10
2 21/10/2004 8 30
3 04/11/2004 8 40
4 22/11/2004 8 30
5 03/12/2004 8 10
6 17/12/2004 8 20
7 07/01/2005 8 40
8 18/01/2005 8 20
9 12/02/2005 8 15
10 19/02/2005 8 20
11 14/03/2005 8 30
12 23/03/2005 8 20
13 10/04/2005 8 10
14 19/04/2005 8 5
15 14/05/2005 7 55
16 17/05/2005 8 25
Table 1: Acquisition dates and time of MODIS imagery. Local time in the Upper Bhima catchment is
GMT +5:30
The 16 NDVI images which were used in the classification are shown in Figure 8. The NDVI peak is
clearly just after the monsoon in October and for the major part of the sub basin the NDVI gradually
decreases except for the forested parts in the Western Ghats and the irrigated areas downstream of
the main reservoirs (Khadakwasla, Ujani and Yadgaon and Bhatghar). Just before the onset of the
monsoon by the end of May the NDVI is lowest in the basin.

Figure 8: 16 MODIS NDVI image used in the land use classification
The 16 NDVI images were stacked into one image and an unsupervised classification using the K-
means clustering was applied. Initially 15 clusters were chosen and a 0.95 convergence criterion was
applied. The 15 classes were grouped into 7 major land use classes based on a field visit and a visual
inspection using the Google Earth application4, the old land use map and an additional decision rule
using the slope derived from the SRTM DEM. The additional decision rule was required to distinguish
between supplemental irrigated agriculture and low density forests and between irrigated agriculture
and high density forests. Areas with slopes over 4% were classified as low density and high density
forest respectively. The specifics of each class are shown in Table 2.
Class name Cluster Slope Area (ha) Area %

# criterion
Water 1 N/A 38950 0.9
Rangelands 2-4 N/A 1228175 26.8
Rainfed agriculture 5-8 N/A 1536194 33.6
Supplemental 9-11 < 4% 613575 13.4
irrigated agriculture (11)
Irrigated agriculture 12-15 < 4% 910775 19.9
(15)
Forest high density 11 > 4% 143363 3.1
Forest low density 15 > 4% 103369 2.3
4
http://earth.google.com

Table 2: Specifics of the land use classes
Rain fed agriculture covers the largest area in the basin (33.6%), but there are also considerable
irrigated areas. The supplemental irrigated agriculture class covers around 13.4% of the basin area.
These are areas which do not obtain their irrigation water from the large scale reservoirs but from
groundwater or small tanks. The large command areas downstream of the major reservoirs cover a
total area of 19.9%. The final LULC map is shown in Figure 9.
Figure 9: LULC map Upper Bhima
Table 3 shows a comparison between the acreages of each class in the old and the new classification.
The differences are considerable. The area under rain fed agriculture has been reduced considerably
(from 52.8% to 33.6%), while the area under irrigated agriculture has increased (4.4% to 19.9%).
The forested areas in the Western Ghats are now also well depicted in the classification.
Class name Area % Area %

(old) (new)
Water 0.7 0.9
Rangelands 29.6 26.8
Rainfed agriculture 52.8 33.6
Supplemental 11.9 13.4
irrigated agriculture
Irrigated agriculture 4.4 19.9
Forest high density 0.5 3.1
Forest low density 0.0 2.3
Table 3: comparison between the two land use classifications


4 Evapotranspiration mapping
4.1 SEBAL
SEBAL provides a way to estimate and monitor actual ET with spatial sensitivity, without having to
have soil moisture, land use and vegetation conditions. SEBAL solves the surface energy balance for
heterogeneous terrain on the basis of reflected solar radiation and emitted thermal radiation (surface
temperature). The actual ET (ETact) fluxes from SEBAL reflect the effects of various natural factors
that influence ET, such as moisture availability, presence of pests and disease, salinity, and other
factors. The standard ET equations are designed to compute potential ET, or the level of ET that would
occur under optimal or pristine, although sometimes general corrections are applied for conditions
where water is limiting limitations by using a reduction coefficient (ETact= ETpot, e.g. Budyko, 1969).
SEBAL is one of the first mathematical procedures that can operationally estimate spatially distributed
ETact from field to river basin scale over an unlimited array of land use types, including desert soil,
open water bodies, sparse natural vegetation, rainfed crops, irrigated crops, etc. The SEBAL model
solves the energy balance for every every individual MODIS and Landsat pixel, thereby providing the
spatial sensitivity. Satellite images need to be cloud-free to be processed for energy balance purposes.
The three primary bio-physical inputs from MODIS and Landsat images into SEBAL are (i) surface
temperature, (ii) surface albedo and (iii) Normalized Difference Vegetation Index (NDVI) (see Figure
10). All of these parameters are measured directly or derived from measurements recorded by
satellite-based sensors. In addition to that, a water mask is created. The water mask is meant for the
assignment of particular values that are applicable to water only: emissivity, surface resistance, and
soil heat flux/net radiation fraction. The latter fraction is relevant because the equations for soil heat
flux for land and water are completely different. An existing generalized land use map is necessary to
assign vegetation heights for the computation of the surface roughness for all pixels. This vegetation
height is only considered for the turbulent parameters (i.e. momentum flux) and not for anything else.
The inputs to SEBAL consists of (i) satellite multi-spectral radiances, (ii) routine weather data and (iii)
DEM.

Figure 10: Data flow and key steps for the determination of spatially distributed ET fluxes according to
the SEBAL method
SEBAL uses a set of algorithms to solve the energy balance at the earths surface. The instantaneous
ET flux is calculated for each pixel within a remotely sensed image as a 'residual' of the surface energy
budget equation:
E = Rn - G - H (1)
where E is the latent heat flux (W/m2) (which can be equated to ET), Rn is the net radiation flux at
the surface (W/m2), G is the soil heat flux (W/m2), and H is the sensible heat flux to the air (W/m2).
Figure 11: Main terms of the Surface Energy Balance

Rn represents the actual radiant energy available at the surface. It is computed by subtracting all
outgoing radiant fluxes from all incoming radiant fluxes. This is further specified in the surface
radiation balance equation:
Rn = RS - RS + RL - RL - (1 - 0) RL (2)
where RS is the incoming short-wave radiation (W/m2), is the surface albedo (dimensionless), is
the incoming long wave radiation (W/m2), RL is the outgoing long wave radiation (W/m2), and o is
the surface thermal emissivity (dimensionless).
In Eq. (2), the net short-wave radiation (RS - RS) that remains available at the surface is a
function of the surface albedo (). The broad band surface albedo is derived from the narrow band
spectral reflectances () measured by each satellite band. The incoming short-wave radiation (RS) is
computed using the solar constant, the solar incidence angle, a relative earth-sun distance, and a
computed broad band atmospheric transmissivity. In this application of SEBAL, atmospheric
transmissivity was obtained from a few selected automatic weather stations in Nevada that sustained
the data integrity check. The incoming long wave radiation (RL) was computed using a modified
Stefan-Boltzmann equation with an apparent emissivity that is coupled to the shortwave atmospheric
transmissivity and a measured air temperature. Outgoing long wave radiation (RL) is computed using
the Stefan-Boltzmann equation with a calculated surface emissivity and surface temperature. Surface
temperature is derived from the narrow-band satellite measurements of thermal infrared radiation.
The challenge in energy balance modeling is to partition net radiation into sensible (H) and latent heat
fluxes (E). To guide this partitioning process, extreme values of H are defined for hot and cold
pixels selected by the operator. The hot pixel has a high surface temperature associated with the
absence of evaporative cooling. Characteristics that qualify pixels for consideration as hot pixels are
as follows: (i) their surface temperatures occur near the upper end of the frequency distribution of all
pixels in the image, (ii) they have relatively sparce vegetative cover (as indicated by NDVI) and (iii)
they appear to consist of essentially bare land on the false color composite.
In contrast, characteristics that qualify pixels for consideration as cold pixels are as follows: (i) their
surface temperatures occur near the lower end of the frequency distribution, (ii) they appear to be
within open water surfaces or well irrigated fields and (iii) they have a low reflectance in the visible
part of the spectrum. For each image processed, the operator must decide which pixels to select as
the hot and cold pixels.
At the hot pixel it is assumed that ET is zero; thus, HRn-G. At the cold pixel, it is assumed that
sensible heat is very small or zero, because all or most of the net available energy is used for ET; thus,
ERn-G. The sensible heat does not necessarily have to be exactly zero, but it should be small.
Select of the hot and cold pixels for each image constrains the range of sensible heat within realistic
bounds, essentially providing a self-calibrating feature, with the H fluxes for all other pixels lying
between. Interpolation between these two bounds is done according to a linear function of surface
temperature where the a and b coefficients are obtained from linking H to the surface temperature at
the hot and cold pixels:

H = {a cp u* (a + b T0)} / {ln(z2/z1) - h (z2,L) + h (z1,L)} (3)
Where a (kg/m3) is the moist air density, cp (J/kg/K) is the specific heat at constant pressure, u*
(m/s) is the friction velocity, z (m) is the near-surface reference height to which the a,b coefficients as
well as the constant H flux applies, h (-) is the stability correction for heat transport and L (m) is the
Monin Obukhov length. Eq. (3) requires single layer wind speed observations and surface roughness
for the determination of friction velocity u* to be known:
u* = {u(z2) k} / {ln(z2/z0m) - m (z2,L) + m (zom,L)} (4)
where z0m (m) is the surface roughness length for momentum transport and k (-) the von Karmans
constant. SEBAL uses an iterative process to correct for atmospheric instability caused by buoyancy
effects of surface heating. The surface roughness z0m is computed from estimates of vegetation
height based on a simple land use mask and the actual LAI following the simple roughness model of
Raupach (1994) and suggestions of model performance tests published by Verhoef et al. (1997).
In Eq. (1), the soil heat flux (G) and sensible heat flux (H) are subtracted from the net radiation flux at
the surface (Rn) to compute the "residual" energy available for evapotranspiration (E). Soil heat flux
is empirically calculated as a G/Rn fraction using instantaneous surface temperature and an
approximation for near-surface soil temperature at sunrise. A light interception reduction function is
applied for the presence of foliage using NDVI as an indicator for canopies.
The E time integration in SEBAL is split into two steps. The first step is to convert the instantaneous
latent heat flux (E) into daily E24 values by holding the evaporative fraction constant. The
evaporative fraction EF is:
EF = E/ (Rn - G) (-) (5)
Field measurements under various environmental conditions indicate that EF is nearly constant with
time during the diurnal cycle. Thus, because EF ~ EF24, the 24 hour latent heat flux can be
approximated by:
E24 = EF (Rn24 G24) (W/m2) (6)
where (-) is an advection enhancement parameter. It is assumed that the evaporative fraction EF
specified in Eq. (5) remains quasi-constant during daytime hours and that variations can be related to
advection . Experimental work has demonstrated that this holds true for environmental conditions
where soil moisture does not significantly change (e.g. Shuttleworth et al., 1989; Brutsaert and Sugita,
1992; Nicols and Cuenca, 1993; Kustas et al., 1994; Crago, 1996; Farah, 2001). This needs to be true
for single days only to get an appropriate value for E24 during satellite overpass days.
Advection occurs when hot desert air is blown over irrigated areas, typically in the afternoon, by winds
driven by differences in air densities. The regional scale atmospheric circulation adds extra energy to
moist vegetation. This extra energy can increase ET to the extent that it exceeds the net available
energy (Rn-G0), which from an energy balance point of view implies that the sensible heat flux is
directed towards the surface (H is negative).

Within SEBAL, advection is incorporated into the parameter that varies with the weather conditions
and the moisture of the underlying soils. The impact of weather conditions on advection is expressed
by the changes of the evaporative fraction of the reference crop EFgrass between satellite overpass
and the 24-h counterpart governed by day and night weather conditions. The influence of advection
decreases non-linearly with soil moisture content, which implies that only wet surface are exposed to
advection.
Depending on the time scale chosen, different time integrations of (Rn- G0) need to be obtained. For
the daily time scale, ET24 is formulated as:
86400 10 3
ET24 = E24 (mm d-1) (7)
w
where Rn24 (W m-2) is the 24-h averaged net radiation, (J kg-1) is the latent heat of vaporization,
w (kg m-3) is the density of water.
The second step is the conversion from a daily latent heat flux into 14-day values, which was achieved
by application of the Penman-Monteith equation:
E24PM = (sa Rn24 + acp e/ra) / (sa + (1 + rs/ra)) (W/m2) (8)
where sa (mbar/K) is the slope of the saturated vapor pressure curve, acp (J/m3 K) is the air heat
capacity, e (mbar) is the vapor pressure deficit, (mbar/K) is the psychrometric constant and ra
(s/m) is the aerodynamic resistance. The parameters sa, e and ra are controlled by meteorological
conditions, and Rn and rs by the hydrological conditions.
The SEBAL computations can only be executed for cloudless days. The result of E24 from Eq. (6) has
been explored to convert the Penman-Monteith equation (8) and to quantity rs inversely using E24=
EPM. The spatial distribution of rs so achieved, will consequently be used to compute E24 by means
of Eq. (8) for all days without satellite image available (e.g. Farah, 2001). The total ETact for a given
period can be derived from the shorter term average values of sa, e, ra and Rn
The standard 250 m Digital Elevation Model (DEM) has been used for the correction of air pressure and
related air density and psychrometric constant at higher elevation. The DEM is also used to correct the
absorbed solar radiation values, both for slope and aspect. Southern facing terrain due to the angle of
incidence absorb, namely, more solar radiation per unit land than the Northern facing slope.

4.2 Reference evapotranspiration
Figure 12: Total reference evapotranspiration from 1 October 2004 to 31 May 2005
120
100
80
ETref (mm)
60
40
20
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
Period
Figure 13: Bi-weekly reference evapotranspiration from 1 October 2004 to 31 May 2005

4.3 Potential Evapotranspiration
Figure 14: Total potential evapotranspiration from 1 October 2004 to 31 May 2005

4.4 Actual Evapotranspiration
Figure 15: Total actual evapotranspiration from 1 October 2004 to 31 May 2005

140 BASIN 140 WATR 140 RNGE

120 120 120
100 100
ET (mm)
ET (mm)
100
ET (mm)
80 80 80
60 60 60
40 40 40
20 20 20
0 0 0
2 4 6 8 10 12 14 16 2 4 6 8 10 12 14 16 2 4 6 8 10 12 14 16
Period Period Period
140 AGR1 140 AGR2 140 AGR3

120 120 120
100 100
ET (mm)
ET (mm)
100
ET (mm)
80 80 80
60 60 60
40 40 40
20 20 20
0 0 0
2 4 6 8 10 12 14 16 2 4 6 8 10 12 14 16 2 4 6 8 10 12 14 16
Period Period Period
140 FRST 140 FRSE
120 120 ETpot
100 ETact
ET (mm)
100
ET (mm)
80 80
60 60
40 40
20 20
0 0
2 4 6 8 10 12 14 16 2 4 6 8 10 12 14 16
Period Period
Figure 16: Bi-weekly potential and actual evapotranspiration per land use from 1 October 2004 to 31
May 2005

4.5 Biomass production
Figure 17: Total biomass production from 1 October 2004 to 31 May 2005
1600
1200
Production (kg/ha)
800
400
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
Period
Figure 18: Bi-weekly biomass production from 1 October 2004 to 31 May 2005

5 Precipitation mapping using Remote

Sensing
5.1 Introduction
A limited set of precipitation data from meteorological stations is available for the Upper Bhima
catchment. A dataset with monthly precipitation data for eight Indian Meteorological Department
(IMD) stations and for 15 additional stations from the National Data Centre (NDC) in Pune were
available. The location, source and time coverage of the data is shown in Figure 19.
Figure 19: Location, source and time overage monthly precipitation data
Since there is an extreme spatial variation in precipitation ranging from 2000 mm in the Western Ghats
to 650 mm in the most eastern part of the sub-basin it is important to have a dense network of
meteorological stations. In this study SEBAL evapotranspiration maps will be used to calibrate the
hydrological model. SEBAL data are available from October 2004 until May 2005, but precipitation data
are only available up to 2004 for a limited number of stations. For these reasons Remote Sensing
derived precipitation from the Tropical Rainfall Monitoring Mission (TRMM) is used in this study.

5.2 Tropical Rainfall Monitoring Mission5
The Tropical Rainfall Measuring Mission (TRMM) is a joint space project between NASA and the
Japanese Aerospace Exploration Agency (JAXA). TRMM is designed to measure tropical precipitation
and its variation from a low-inclination orbit combining a suite of sensors to overcome many of the
limitations of remote sensors previously used for such measurements from space. The fundamental
objective of TRMM is to understand the role of latent heat in driving the circulation of the global
atmosphere. With the inclusion of rain radar, TRMM will provide the first opportunity to estimate the
vertical profile of the latent heat that is released through condensation. The TRMM rainfall data is
particularly important for studies of the global hydrological cycle and for testing the realism of climate
models, and their ability to simulate and predict climate accurately on the seasonal time scale. Other
scientific issues such as the effects of El Nio on climate could be addressed with a reliable, extended
time series of tropical rainfall observations as well. The TRMM satellite has a latitudinal range from
50oS - 50 oN. The TRMM satellite has been launched November 27 in 1997 and the mission has
recently been extended to 2009.
The observatory for rainfall observations consists of a precipitation radar, a multi-frequency microwave
radiometer and a visible and infrared (VIS/IR) radiometer. Used in this combination in a
complementary way, these sensors, with the designed orbit of 350 km altitude and inclination of 35,
will meet the measurement accuracy for rainfall needed to fulfill the stated science objectives of the
mission. For related precipitation observations (i.e., lightning) and for additional important climate
observations, a lightning imager and an Earth radiation budget sensor have been added. A brief
description of each of the five instruments is included here:
The Precipitation Radar (PR), the first in space, will measure the 3-D rainfall distribution over
both land and ocean. More specifically, this instrument will define the layer depth of the
precipitation and provide information about the rainfall reaching the surface, the key to
determining the latent heat input to the atmosphere. A unique feature is that it will permit the
measurement of rain over land where passive microwave channels have more difficulty. The
PR is an electronically scanning radar operating at 13.8 GHz with horizontal polarization using
a 128-slotted waveguide antenna and solid state power amplifiers to develop an active
phased array. The horizontal resolution is 4.3 km at nadir, the range resolution is 250 m and
the scanning swath width is 220 km.
The multi-channel microwave radiometer, designated as the TRMM Microwave Imager (TMI),
is designed to provide information on the integrated column precipitation content, its areal
distribution, and its intensity. The horizontal resolution will allow the definition and
investigation of most rainfall types, including convective cells. The technique is best suited for
estimates over oceans, where data are needed most for climate model verification. The TMI
will operate on 5 frequencies, of which 4 will have dual polarization thus providing 9 channels
of data. The 5 frequencies are 10.65 GHz, 19.35 GHz, 22.235 GHz (single polarization), 37.0
GHz and 85.5 GHz. The horizontal resolution of the TMI will range from 5 km at 85.5 GHz to
45 km at 10.65 GHz. A scan angle of 65 will provide a swath width of 760 km.
The Visible Infrared Scanner (VIRS) will provide very high resolution information on cloud
coverage, type, and cloud top temperatures and also serve as the link between these data
5
http://trmm.gsfc.nasa.gov/

and the long and virtually continuous coverage by the geosynchronous meteorological
satellites. The VIRS is a 5 channel cross-track scanning radiometer operating at 0.63, 1.6,
3.75, 10.80, and 12.0 microns. The instrument, with a swath width of 720 km, will provide
cloud distributions by type and height and rain estimates from brightness temperatures at a
horizontal resolution of 2.1 km (nadir). Direct correlation with PR and TMI measurements will
assist in providing more accurate rain estimates from VIRS thereby enhancing its capability
for verifying or calibrating the rain estimates from operational meteorological satellites.
The Lightning Imaging Sensor (LIS) is designed to investigate the global incidence of
lightning, its relationship with the global electric circuit, and, in conjunction with the PR, TMI,
and VIRS, its correlation with rainfall. Lightning is not a necessary condition for tropical
rainfall such as warm rain or even for some convective rainfall in the tropics. However, when
updrafts and downdrafts are sufficient to cause lightning, rainfall results. LIS will be optimized
to detect the lightning location, mark the time of occurrence, and measure the radiant
energy. LIS is a calibrated optical sensor operating at 0.7774 microns and will observe the
distribution and variability of lightning over the Earth. The horizontal resolution at nadir will
be 5 km and the swath width will be 590 km.
Several orbital and gridded data products are available for download at the Goddard Distributed Active
Archice Centre (DAAC)6. For this study the 3B-43 product is used. The global rainfall algorithm (3B-43)
combines the estimates generated by combined instrument rain calibration (3B-42) and global gridded
rain gauge data from CAMS, produced by NOAA's Climate Prediction Center and/or global rain gauge
product, produced by the Global Precipitation Climatology Center (GPCC). The output is rainfall for
0.25x0.25 degree grid boxes for each month.
5.3 Downscaling and calibration of TRMM data
For the model simulation the irrigation year June 2004 to May 2005 will be used in the simulations.
SEBAL data are available from October 2004 to May 2005.
The total precipitation at the original resolution from June 2004 to May 2005 is shown on the left side
of Figure 22. The 0.25 degree resolution is still too coarse for application in a simulation model.
Therefore a downscaling procedure is developed and applied to the data. For all meteorological
stations in the sub-basin which cover the period 1998-2004 the monthly TRMM data were extracted.
For each station a regression analysis was performed on the monthly precipitation amounts.
6
http://disc.sci.gsfc.nasa.gov/data/datapool/TRMM/01_Data_Products/index.html

Station Slope Intercept R2

Ambavade-2 0.26 21.94 0.29
Askheda 0.80 18.15 0.80
Barur (Takali) 0.59 16.34 0.51
Chaskman 1.06 11.79 0.81
Kashti 0.40 18.49 0.37
Khamgaon 0.86 18.29 0.65
Late 0.47 2.65 0.50
Pargaon 0.84 20.15 0.64
Rosa (Kolgaon) 0.71 15.34 0.57
Sakhar 1.09 0.84 0.82
Shirur 0.32 12.69 0.45
Sidhewadi 0.58 19.03 0.51
Wegre (Muthe) 0.70 37.58 0.68
2
Table 4: Slope, intercepts and R for the 13 selected meteorological stations
Table 4 shows the slopes, intercepts and R2 values of the regression analysis. The R2 values vary from
0.29 (Ambavade-2) and 0.82. For Askheda and Sakhar the scatter plots are shown in Figure 20.
800 Ashkeda (r2 = 0.80) 1000 Sakhar (r2 = 0.82)
800
600
P station (mm)
P station (mm)
600
400
400
200
200
0 0
0 200 400 600 800 0 200 400 600

P TRMM (mm) P TRMM (mm)
Figure 20: Scatterplots for monthly TRMM precipitation and observed precipitation from 1998-2004 for
two meteorological stations
The next step in the downscaling procedure is the spatial interpolation of the slopes and intercepts of
the meteorological stations (Figure 21). To avoid erroneous precipitation amounts, especially in the dry
seasons, only those station with an R2 > 0.50 and an intercept < 20 are used in the interpolation (7
stations). The slopes range from 0.57-1.10 and the intercepts range from 0-22 mm/month.

Figure 21: Slope and intercept used to scale TRMM derived precipitation
The final step in the procedure is to convert the centers of the TRMM grid boxes to points and to
interpolated these points to a raster of a higher resolution (0.005 degrees by 0.005 degrees). A spline
tension interpolator was used with weight 1 and using 4 surrounding points (Franke, 1982). The result
of this downscaling procedure is shown in Figure 22.
Figure 22: Total precipitation from June 2004 to May 2005. The scaled TRMM data is shown in the left
figure and the downscaled TRMM data is shown in the right figure.
This downscaling procedure is repeated for each month between June 2004 and May 2005. Figure 23
shows the monthly precipitation in the basin. Precipitation is clearly concentrated in the monsoon
months from June to September with the peak in July 2005. There is a clear decreasing trend from the
western to the eastern part of the sub basin.

Figure 23: Monthly downscaled TRMM derived precipitation in the Upper Bhima sub basin

6 Modeling the Upper Bhima catchment

The SWAT model of the Upper Bhima catchment is created in several consecutive steps, which are
described in this chapter. Based on the topography the catchment is subdivided into smaller units
(sub- basins)
6.1 Topography
A Remote Sensing derived digital elevation model (DEM) acquired with Shuttle Radar Topographic
Mission (SRTM)7 is used to partition the catchments in sub-basins and rivers. GIS based hydrological
surface function are used in an automated procedure of the SWAT GIS interface (AVSWAT). The
original DEM has a spatial resolution of 90m, but to reduce calculation times the DEM is resampled to a
spatial resolution of 250 m.
Figure 24: Digital elevation model, sub-basins, and major rivers
Figure 24 show the DEM of the catchment. The elevation ranges from 414m in the east of the
catchment to 1458m in the western Ghats. The elevation distribution of the entire catchment is shown
in Figure 25. It shows that nearly 95% of the area of the catchment is below 800m.
7
http://www2.jpl.nasa.gov/srtm/

100
Area (%) 80
60
40
20
0
0 500 1000 1500
Elevation (m)
Figure 25: Elevation distribution in the catchment
The catchment has an area of 4.567.781 ha and a total of 115 sub-basins are delineated with a
minimum area of 20.000 ha (Figure 24).
6.2 Soils
The FAO UNESCO digital soil map of the world is used in this study (FAO, 1995). The Digital Soil Map
of the World (DSMW) CD-ROM is based on the FAO/UNESCO Soil Map of the World. The FAO-Unesco
Soil Map of the World was published between 1974 and 1978 at 1:5.000.000 scale. The soil map of the
world was prepared on the base of the topographic map series of the American Geographical Society
of New York . The legend of the original soil map of the World (FAO, 1974) comprises an estimated
4930 different map units, which consist of soil units or associations of soil units. When a map unit is
not homogeneous, it is composed of a dominant soil and component soils. The latter are: associated
soils, covering at least 20 % of the area; and inclusions, important soils which cover less than 20 % of
the area. The legend of the soil map of the world (1974) comprises 106 soil units (from Af to Zt),
grouped in 26 major soil groupings. The soil map for the upper Bhima catchment is shown in Figure
26.
The majority of the alluvial plains in the catchment consist of chromic vertisols. Vertisols are churning
heavy clay soils with a high proportion of swelling 2:1 lattice clays. These soils form deep wide cracks
from the surface downward when they dry out, which happens in most years. Dry Vertisols have a
very hard consistence; wet Vertisols are (very) plastic and sticky. It is generally true that Vertisols are
friable only over a narrow moisture range but their physical properties are greatly influenced by soluble
salts and/or adsorbed sodium. Infiltration of water in dry (cracked) Vertisols with surface mulch or a
fine tilth is initially rapid. However, once the surface soil is thoroughly wetted and cracks have closed,
the rate of water infiltration becomes almost zero. (The very process of well/shrink implies that pores
are discontinuous and non-permanent.) If, at this stage, the rains continue (or irrigation is prolonged),
Vertisols flood readily. The highest infiltration rates are measured on Vertisols that have a considerable
shrink/swell capacity, but maintain a relatively fine class of structure. Not only the cracks transmit
water from the (first) rains but also the open spaces between slickensided ped surfaces that developed
as the peds shrunk. Data on the water holding capacity of Vertisols vary widely, which may be

attributed to the complex pore space dynamics. Water is adsobed at the clay surfaces and retained
between crystal lattice layers. By and large, Vertisols are soils with good water holding properties.
However, a large proportion of all water in Vertisols, and notably the water held between the basic
crystal units, is not available to plants. Investigations in the Sudan Gezira have shown that the soil
moisture content midway between large cracks changes very little, if at all, when the clay plain is
flooded for several days or even several weeks. The soil's moisture content decreases gradually from
more than 50 percent in the upper 20 cm layer to 30 percent at 50 cm depth. Deeper than 100 cm,
the soil moisture content remains almost invariant throughout the year.
Figure 26: FAO/UNESCO soil map of the Upper Bhima catchment
The upland areas mainly have a luvisol soil type. This soil type is common in upland areas with a
distinct wet and dry season. Luvisols have favourable physical properties; they have granular or crumb
surface soils that are porous and well aerated. The `available' moisture storage capacity is highest in
the argic horizon (15 to 20 volume percent). The argic horizon has a stable blocky structure but
surface soils with a high silt content may be sensitive to slaking and erosion. Most Luvisols are well
drained but Luvisols in depression areas with shallow groundwater may develop gleyic soil properties
in and below the argic horizon. Stagnic properties are found where a dense illuviation horizon obstructs
downward percolation and the surface soil becomes saturated with water for extended periods of time.
The physical properties needed for the SWAT model are also derived from the FAO soil map of the
world. These are summarized in Table 5.

Table 5: Soil physical properties
Code Description Depth AWC Clay Silt Sand Bulk Organic

density carbon
content
cm mm/mm % % % g/cm3 %
Be66-2c Eutric Cambisols 1280 0.21 22 35 43 1.21 0.79
Bv12-3b Vertic Cambisols 1430 0.10 50 21 30 1.59 0.77
Hh11-2bc Haplic Phaeozems 1368 0.15 28 24 48 1.51 0.52
I-Hh Lithosols 950 0.17 21 21 58 1.45 0.56
Lc75-1b Chromic Luvisols 1430 0.14 29 21 50 1.53 0.60
Nd51-2b Dystric Nitosols 2250 0.19 34 18 48 1.19 0.74
Vc43-3ab Chromic Vertisols 2150 0.20 32 16 52 1.20 0.52
Vc44-3a Chromic Vertisols 1430 0.07 54 23 23 1.70 0.59
Vc45-3a Chromic Vertisols 1550 0.07 58 19 23 1.71 0.41
Vp42-3a Pellic Vertisols 1550 0.08 57 23 21 1.65 0.51
6.3 Hydrological Response Units
Hydrological Response Units (HRUs) are the smallest unit of calculation in the SWAT model and are
defined as unique combinations of soil and land use. Before determining the HRUs a preprocessing
step is taken. However an area threshold is specified for both soils and land use. This threshold
determines whether a certain soil or land use class is used as an HRU based on the proportion of the
total sub basin area the soil class or land use type has. The thresholds used are 10% and 5 % for soil
and land use respectively. The HRUs are created using an automated procedure in the AVSWAT
interface. A total of 768 HRU are generated. A map with the HRUs is shown in Figure 27.
Figure 27: Hydrological Response Units in the catchment

6.4 Meteorology
The SWAT model requires the following daily meteorological parameters:
Solar radiation
Minimum and maximum temperature
Relative humidity
Wind speed
Precipitation
Different sources of information have been used to derive the daily weather data from 1 June 2004 to
31 May 2005. For all parameters data are used for two different meteorological stations; one in Pune
and one in Sholapur (Figure 28). Similar to the procedure in SEBAL: the data of the two stations have
been averaged. For precipitation a virtual meteorological station is created for each of the 112
subbasins.
Figure 28: Location of the meteorological stations
6.4.1 Solar radiation
The solar radiation is defined as the amount of radiation reaching a horizontal plan on the earths
surface. First the extraterrestrial radiation at the top of the atmosphere was calculated, which is then
corrected for scattering, reflection and adsorption by atmospheric gases, clouds and dust. The
extraterrestrial radiation is depending on the latitude and the Julian day in the year. The amount of
solar radiation (Qs (MJ m-2 day-1)) reaching the surface of the earth is calculated using the Angstrm
formula:
n
Q s = a s + bs Q a
N
Where as(-) is a regression constant, expressing the fraction of extraterrestrial radiation reaching the
earth on overcast days, as+bs is the fraction of extraterrestrial radiation reaching the earths surface on

clear days, n/N (-) is the relative sunshine duration and Qa (MJ m-2 day-1) is the extraterrestrial
radiation. Values of 0.25 and 0.50 are used for as and bs respectively.
Solar radiation is calculated for both Pune and Sholapur (Figure 29). The relative sunshine duration is
derived from the IWMI climate atlas8.
26
24
Qs (MJ m-2 day-1)
22
20
18
16
14
31/7/04 29/9/04 28/11/04 27/1/05 28/3/05 27/5/05
Date
Figure 29: Average daily solar radiation Pune and Sholapur
6.4.2 Temperature
Public domain data9 are used for daily maximum and minimum temperature (Figure 30).
50
40
30
T (oC)
20
10
0
31/7/04 29/9/04 28/11/04 27/1/05 28/3/05 27/5/05
Date
Figure 30: Average daily maxium and minium temperature
8
http://www.iwmi.cgiar.org/WAtlas/atlas.htm
9
http://www.wunderground.com

6.4.3 Relative humidity
9
Public domain data are used for daily relative humidity.
0.8
0.6
RH (-)
0.4
0.2
0
31/7/04 29/9/04 28/11/04 27/1/05 28/3/05 27/5/05
Date
Figure 31: Average daily relative humidity
6.4.4 Wind speed
9
Public domain data are used for wind speed.
6
Windspeed (m s-1)
0
6/1/04 7/31/04 9/29/04 11/28/04 1/27/05 3/28/05 5/27/05
Date
Figure 32: Average daily windspeed
6.4.5 Precipitation
Because of the spatial variation in precipitation a different approach has been adopted. The TRMM
precipitation rasters as described in Chapter 5 are used to generated monthly precipitation for each
sub basin from October 2004 to September. At the center of each sub basin a virtual meteorological
station is generated. The monthly precipitation data are converted to daily values using the following
procedure:

A random day in the month is selected.

A random amount of precipitation is generated. In the monsoon months (June, July, August,
September) this amount is between 30 and 60 mm day-1. In all other months this amount is
between 10 and 20 mm day-1.
This process is repeated until the total amount of monthly precipitation is distributed.
An example of the daily precipitation of two sub basins is given in Figure 33. Sub basin 57 is located in
the Western ghats in the south western part of the catchment and sub basin 61 is located in the
eastern part of the catchment.
60 60
Sub 57 (1744 mm) Sub 61 (618 mm)
P (mm day-1)
P (mm day-1)
40 40
20 20
0 0
6/1/04 7/31/04 9/29/04 11/28/04 1/27/05 3/28/05 5/27/05 6/1/04 7/31/04 9/29/04 11/28/04 1/27/05 3/28/05 5/27/05
Date Date
Figure 33: Daily precipitation for sub basins 57 and 61
6.5 Management practices
Quantifying the impact of land management and land use on water supply and quality is a primary
focus of environmental modeling. SWAT allows very detailed management information to be
incorporated into a simulation. In the model for the Upper Bhima catchment the following operations
are used:
The plant operation initiates plant growth. This operation can be used to designate the time
of planting for agricultural crops or the initiation of plant growth in the spring for a land cover
that requires several years to reach maturity (forests, orchards, etc.). The plant operation will
be performed by SWAT only when no land cover is growing in an HRU. Before planting a new
land cover, the previous land cover must be removed with a kill operation or a harvest and kill
operation. Information required in the plant operation includes the timing of the operation
(month and day or fraction of base zero potential heat units), the total number of heat units
required for the land cover to reach maturity, and the specific land cover to be simulated in
the HRU.
The harvest and kill operation stops plant growth in the HRU. The fraction of biomass
specified in the land covers harvest index (in the plant growth database) is removed from the
HRU as yield. The remaining fraction of plant biomass is converted to residue on the soil
surface.

The fertilizer operation applies fertilizer or manure to the soil. Information required in the
fertilizer operation includes the timing of the operation (month and day or fraction of plant
potential heat units), the type of fertilizer/manure applied, the amount of fertilizer/manure
applied, and the depth distribution of fertilizer application.
Application of irrigation water. Information is required on the amount of irrigation water to be
applied, the salt content of the irrigation and the timing of the operation.
The agricultural land use classes need to be parameterized for different crops and management
activities such as fertilization and irrigation. As described in chapter 3 three agricultural land use
classes have been defined, which have been parameterized as follows:
Table 6: management operations rainfed agriculture (AGR1)
Date Crop Operation Remarks

Jun-10 Sorghum Planting of crop
Jun-22 Sorghum Fertilizer application 300 kg/ha Urea
Sep-30 Sorghum Harvest
Table 7: management operations supplemental irrigated agriculture (AGR2)

Jun-10 Sorghum Planting of crop
Jun-10 Sorghum Fertilization 300 kg/ha Urea
Aug-01 Sorghum Fertilization 300 kg/ha Urea
Sep-29 Sorghum Harvest
Oct-15 Winter wheat Planting of crop
Oct-17 Winter wheat Fertilization 300 kg/ha Urea
Nov-01 Winter wheat Fertilization 300 kg/ha Urea
Nov-10 Winter wheat Irrigation 40 mm
Nov-20 Winter wheat Irrigation 40 mm
Dec-15 Winter wheat Irrigation 100 mm
Jan-15 Winter wheat Irrigation 100 mm
Feb-10 Winter wheat Fertilization 300 kg/ha Urea
Table 8: management operations intensive irrigated agriculture (AGR3)


Jun-02 Sugarcane Planting of crop
Jun-02 Sugarcane Fertilization 300 kg/ha Urea
Jun-15 Sugarcane Irrigation 60 mm
Jul-15 Sugarcane Irrigation 40 mm
Aug-15 Sugarcane Irrigation 50 mm
Sep-01 Sugarcane Fertilization 300 kg/ha Urea
Sep-15 Sugarcane Irrigation 20 mm
Oct-15 Sugarcane Irrigation 50 mm
Nov-15 Sugarcane Irrigation 40 mm
Dec-15 Sugarcane Irrigation 40 mm
Jan-01 Sugarcane Fertilization 300 kg/ha Urea
Jan-15 Sugarcane Irrigation 40 mm
Feb-15 Sugarcane Irrigation 50 mm
Mar-15 Sugarcane Irrigation 100 mm
Apr-15 Sugarcane Irrigation 100 mm
May-15 Sugarcane Irrigation 100 mm
May-30 Sugarcane Harvest
6.6 Reservoirs
The upper Bhima catchment is a complex catchment and stream flow is for a large part determined by
reservoirs in the upstream area of the catchment which receives most precipitation. Data on reservoirs
is scarce and sensitive given the water disputes between the different states. A total of 13 reservoirs
have been modeled in SWAT. The actual reservoirs and the locations of the reservoirs in the model are
shown in Figure 34. SWAT locates the reservoirs for routing purposes on the sub basin boundary. In
the model a user specified release rate is used to model outflow out of the reservoir. Besides monthly
releases the following inputs are required for each reservoir:
RES_ESA: Surface area of the reservoir when filled to the emergency spillway (ha)
RES_PSA: Surface area of the reservoir when filled to the principal spillway (ha)
RES_EVOL: Volume of water held in the reservoir when filled to the emergency spillway (104
m3 H2O)
RES_PVOL: Volume of water held in the reservoir when filled to the principal spillway (104 m3
H2O)
RES_VOL: Initial Volume of water held in the reservoir when filled at the beginning of the
simulation (104 m3 H2O)

Figure 34: Location of reservoirs
For modeling purposes the reservoirs are grouped together into two reservoirs; one for the Ujani
branch in the middle and one for the Bhatghar branch in the south. It is assumed that all areas with
AGR3 land use receive their irrigation water directly from the reservoirs. The irrigation year 2004-2005
is an average year and it is assumed that the irrigation amounts specified in Table 8 are actually
required. The required stream flow to fulfill these irrigation requirements for both branches are shown
in Figure 38.
300
Ujani
250 Bhatghar
200
Qirr (m3/s)
150
100
50
0
Fe ary
ne
ry
em t
Au y
ec er
ch
ov er
O er
ay
Ja r
ril
s
be
l
Ju
Se gu
ua
b
Ap
ob
b
Ju
M
ar
nu
em
em
br
M
ct
pt
Month
Figure 35: Stream requirements to fullfill irrigation requirements

The difference between the amount of water that flows into the reservoirs and the amount that is
required for irrigation is released as stream flow to downstream basins. The distribution across the
different months is based on the historical monthly stream flow distribution presented in.
0.50
0.45
0.40
0.35
Fraction
0.30
0.25
0.20
0.15
0.10
0.05
0.00
Fe ary
e
ry
em t
Au y
ec er
ch
e r
O er
ay
Ja er
ril
us
N obe
l
n
Ju
ua
D mb
Ap
b
b
Ju
M
ar
nu
g
em
br
M
ct
ov
pt
Se
Month
Figure 36: Monthly fractions of annual streamflow

7 Sensitivity analysis
Major objective of the project is to use the spatial patterns of ET in the auto calibration of the SWAT
model. The model independent parameter estimation package PEST will be used for this purpose
(Doherty, 2004).
The purpose of PEST (which is an acronym for Parameter ESTimation) is to assist in data
interpretation, model calibration and predictive analysis. Where model parameters and/or excitations
need to be adjusted until model-generated numbers fit a set of observations as closely as possible
then, provided certain continuity conditions are met, PEST should be able to do the job. PEST will
adjust model parameters and/or excitations until the fit between model outputs and laboratory or field
observations is optimised in the weighted least squares sense. Where parameter values inferred
through this process are nonunique, PEST will analyse the repercussions of this nonuniqueness on
predictions made by the model.
The first step in setting up PEST is however to define a logical set of parameters and observation
which are to be used in the optimization. Using a spatial set-up results in a number of new challenges.
The SEBAL analysis provides an eight month time series with a biweekly time step of both reference
and actual ET. SWAT provides a 12 month monthly time series of reference ET (ETref) and actual ET
(ETact) per HRU. The parameter optimization is performed on actual ET. The SEBAL data are first
aggregated to monthly data and summarized per HRU in order to be able to directly compare model
outputs to the SEBAL analysis. Before setting up PEST there are two topics which need to be
addressed:
There is a need to explore whether the reference ET generated by SEBAL and SWAT is similar
in time and space and if not a simple manual calibration will need to be performed. Since ETref
is the basis for the calculation of ETact and it is very important that the observed and modeled
ETref time series per subbasin are similar before embarking on any optimization effort with
PEST.
Secondly we need to know if there are significant relations ships between the difference in
observed ETact and modeled ETact (ETact) and time, land use type, soil type, water
management and topography. Before selecting the parameters we would like to optimize we
need to know which part of the system can explain the ETact in both time and space.
These two topics are described in this chapter.
7.1 Exploration of reference evapotranspiration
The reference evapotranspiration is only affected by climatic parameters and can be computed by
weather variables only and is valid for the hypothetical reference crop under well watered conditions
(Allen et al., 1998). The meteorological data for Pune and Sholapur are used for this purpose both in
SWAT and in SEBAL. Conceptually ETref in SWAT is calculated per sub-basin, because it only depends
on the meteorological data and not on soil and land use characteristics.

There are differences in space and time in ETref between SEBAL and SWAT mainly caused by spatial
altitude dependent operations which SEBAL performs on important parameters. The DEM is used to
calculated and correct air pressure and density and thus the psychometric constant. The DEM is also
used to correct the absorbed solar radiation values, both for slope and aspect. Southern facing terrain
due to the angle of incidence, absorb more solar radiation per unit land than the Northern facing slope.
Another cause of this difference is the fact that SEBAL uses grass (0.12m height) as a reference crop
and SWAT uses alfalfa (0.40m height).
35
30
25
20
ETref
15
10
0
Oct-04 Nov-04 Dec-04 Jan-05 Feb-05 Mar-05 Apr-05 May-05
Figure 37: Monthly difference in ETref between calibrated and uncalibrated model
Figure 37 shows the monthly difference in ETref for the entire Upper Bhima catchment. Summed over a
year these differnce may accumulate to over 150 mm annualy. To overcome these differences a simple
calibration procedure has been applied and the amount of incident solar raditation has been correccted
per subbasin per month to calibrate the ETref. Assuming a lineair relation between solar raditaion and
ETref the correction factors have been derived iteratively. The largest correction is required in the
month March (average = 6.1% ) and the smallest in May (average = 1.1%). Figure 38 shows the
scatterplot of the uncalibrated and calibrated annual ETref per subbasin. Calibrated SWAT ETref exhibits
perfect correlation with SEBAl ETref.

1250
1200
1150
SWAT ETref
1100
1050
SWAT base
SWAT calibrated
1000
1200 1205 1210 1215 1220 1225 1230 1235 1240
SEBAL ETref
Figure 38: uncalibrated and calibrated sum of ETref from October 2004 to May 2005 per SWAT
subbasin
7.2 Differences in acual evapotranspiration
After calibration of ETref it is visually verified through a set of box whisker plots whether ETact (ETact
SEBAL ETact SWAT) is explained by different land use, different soil type, month, or precipitation
zone. This analysis leads to the identification of a number of PEST optimisation runs.
Figure 39 shows box-whisker plots of monthly ETact per land use, soil type, month and precipitation
class respectively. The order of magnitude of SWAT ETact is similar to SEBAL, ETact is on average
slightly positive indicating that ETact measured by SEBAL is slightly higher than SWAT simulated
values. The distribution of ETact resembles a normal distribution in most cases, but ranges between
the first and third quartile are considerable and vary between land uses, soil types, months and
precipitation classes.

150 150
100 100
.ETact (mm)
.ETact (mm)
50 50
0 0
-50 -50
-100 -100
-150 -150
-3a
R2
R3
-3b
h
-1b
b
SE
GE
-3a
Vp 3a
R1
ST
-2c
-2b
-2b
I-H
-3a
-
FR
AG
AG
AG
FR
42
44
45
66
12
75
RN
51
43
11
Be
Vc
Vc
Lc
Bv
Nd
Vc
Hh
150
150
100
100
.ETact (mm)
50
.ETact (mm)
50
0
0
-50
-50
-100
-100
-150
-150
dec
feb
y
oct
r
apr
jan
ma
ma
no
p4
p5
p1
p3
p2
Figure 39: Box whisker plots of monthly ETact (SEBAL SWAT) per land use class (top left), soil class
(top right), month (bottom left) and precipitation class (bottom right; p1 = 0-800 mm yr-1, p2 = 800-
1100 mm yr-1, p3= 1100-1400 mm yr-1, p4=1400-1700 mm yr-1, p5=1700-2000 mm yr-1 )
There are clear differences between land uses, soil type and month. The differences between the
precipitation classes are less pronounced. This led to the formulation of the following optimisation
runs:
Available water capacity (AWC). The AWC is defined as the difference between the field capacity of the
soil and the permanent wilting point. It is defined per soil layer per soil type and determines to a large
extent the water buffer capacity of the soil. Ten different soil types with two layers each lead to 20
different parameters to be optimized.
Maximum potential leaf area index (BLAI). The LAI is the leaf area divided by the land area. The BLAI
is one of six parameters that determine leaf area development of a crop in SWAT and determines the
maximum threshold. BLAI is specified per land use type, with the exclusion of water surfaces this lead
to 6 parameters to be optimized (Error! Reference source not found.).
Monthly rainfall increment (RFINC). The RFINC is specified per month and per sub basin and is defined
as the relative monthly adaptation in rainfall. The assumption is made that the spatial distribution of
the TRMM derived precipitation is correct, however that for specific months relative corrections are
allowed within clear boundaries. This leads to one variable to be optimized per month.
Groundwater revap coefficient (REVAP1 and REVAP2). In SWAT water may conceptually move from
the shallow aquifer into the overlying unsaturated zone. In dry periods water is evaporated from the
capillary fringe is evaporated and this process is referred to as groundwater revap and is in SWAT

quantified by the revap coefficient (r) multiplied by ETref. Two optimisation runs are designed. For
the first optimisation run one r for each land use, except water, is defined (REVAP1, 6 variables). For
the second optimisation run it is assumed that r varies per land use and per elevation zone. Four
different elevation zones are defined (0-500 m, 500-600 m, 600-700m and >700m). In combination
with land use this results in 21 unique r resulting from unique combinations of elevation zone and
land use class (REVAP2, 21 variables)


8 Calibration
8.1 Results
All optimisation runs calculate the based on monthly data at sub basin level. In total there are 920
observations (8 months times 115 sub basins). Two final optimisation runs are performed based on
the results of the individual runs (COM1 and COM2). COM1 combines AWC, RFINC and GWREVAP and
COM2 combines AWC, RFINC and GWREVAP2.
Table 9: Results of different optimisations runs. #var and #obs are the number of variables and
observations used in the optimisations. The RMSE () is the Root Mean Square Error or objective
function, is the average of the residuals and # model calls is the number of model calls required to
reach the optimisation results.
PEST run Variable # var # obs RMSE () r2 (mm) # model calls
BASE - 0 920 5.29E+05 0.40 5.20 -

AWC Available water content 20 920 4.49E+05 0.49 5.90 664
BLAI Maximum plant leaf area index 6 920 5.19E+05 0.41 5.50 68
RFINC Monthly rainfall increment 12 920 2.54E+05 0.70 0.80 324
GWREVAP Groundwater revap coefficient 6 920 2.78E+05 0.68 1.60 55
GWREVAP2 Groundwater revap coefficient 21 920 2.66E+05 0.70 1.70 792
COM1 Available water content, monthly 53 920 1.63E+05 0.81 0.50 2987
rainfall increment, groundwater
revap coefficient
COM2 Available water content, monthly 38 920 1.77E+05 0.79 -0.45 1610
rainfall increment, groundwater
revap coefficient
Table 9 shows the results for the different optimisation runs. It shows that the monthly rainfall
increments and the groundwater revap coefficient have the largest effect on the RMSE. Variation of
BLAI has the least effect and only increased r2 by 0.01, which is negligible. The best results are
achieved by combining optimisation runs AWC, RFINC, GWREVAP2. The RMSE is reduced by 69% and
the r2 is 0.81. This means that 81% of the variation in monthly ETact for all 115 sub basins is
explained by the optimised model. COM2 results do not deviate significantly from COM1 results.
However the number of variables used in the optimisation, 38 versus 53, results in a much lower
number of required model calls (1610).
There are several ways to evaluate the reliability of any optimisation of a distributed hydrological
model across time and space. Figure 40 the scatter plots for COM1 between SEBAL ETact and SWAT
ETact on catchment, sub basin and HRU level respectively. It also shows the individual monthly data
and the eight month sum of ETact. The figure shows that the goodness of fit decreases with spatial
and temporal detail. The r2 of monthly catchment ETact for example is as high as 0.90, while at HRU
level the r2 is only 0.35. In time we see similar patterns. The r2 at monthly sub basin level is 0.81
while if the eight month sum is analyzed the r2 increase to 0.92.

Monthy ETact Sum ETact
250 1600
r2 = 0.90 r2 = 1.00
ET SWAT (mm)
ET SWAT (mm)
200 1200
Catchment
150
800
100
50 400
0 0
0 50 100 150 200 250 0 400 800 1200 1600

ET SEBAL (mm) ET SEBAL (mm)
250 r2 = 0.81 1600 r2 = 0.92

ET SWAT (mm)
ET SWAT (mm)
200 1200
Sub basin
150
800
100
50 400
0 0
0 50 100 150 200 250 0 400 800 1200 1600

250 r2 = 0.35 1600 r2 = 0.39

ET SWAT (mm)
ET SWAT (mm)
200 1200
150
HRU
800
100
50 400
0 0
0 50 100 150 200 250 0 400 800 1200 1600

Figure 40: Scatter plots of SEBAL and SWAT ETact. Monthly data are shown on the left side graphs and
the eight month sum is on the right side of the graphs. Spatial detail increases from top to bottom and
ranges from catchment, sub basin to HRU level respectively. SWAT results relate to COM2
optimisation.
Figure 41 shows the map with the eight month ETact sum for SWAT and SEBAL at sub basin and at
HRU level. At sub basin level the spatial patterns between SWAT and SEBAL are highly consistent. At
HRU level there are considerable differences. The general spatial patterns are well depicted however
some HRUs within a sub basin evaporate more than measured with SEBAL and some less, however
aggregated over the entire sub basin these differences are levelled out.

Figure 41: Eight month sum of ETact for SWAT and SEBAL on sub basin and HRU level respectively.
SWAT results relate to COM2 optimisation.
We have shown that it is possible to accurately predict monthly ETact at sub basin level the results are
verified by cross checking the discharges in the river system. No measured discharged are available for
the simulation period (June 2004 May 2005). There is however a dataset available from 1970 to
1996. The discharges of the Sina tributary are used, because the Sina catchment is less disturbed by
large reservoirs and man made structures than the Bhima tributary.

Table 10 shows the observed and simulated discharges. The simulated discharges in 2004-2005 are
well within one standard deviation of the average measured discharges between 1970 and 1996. It
should be noted though hat the coefficient of variation in the observed discharges are large ranging
from 68% in August to 170% in March.

Table 10: Discharges Sina branch. Qswat is the modelled dischare from june 2004 to may 2005. Qobs is
the average observed discharge from 1970-1996 and CVobs is the coefficient of variation in the
observed discharges.
3 3
Month QSWAT (m /s) Qobs (m /s) CVobs (%)
june 4 90 161
july 240 540 96
august 366 740 68
september 456 562 81
october 78 297 114
november 47 61 117
december 45 28 112
january 16 15 98
february 9 8 97
march 2 5 170
april 0 3 149
may 0 6 149
8.2 Discussion and conclusions
The spatial distributed hydrological model SWAT can be successfully calibrated using the GML
algorithm on the basis of Remotely Sensed derived evapotranspiration in a data scarce area. The best
results were obtained by optimising a combination of soil, meteorological and groundwater related
parameters on an eight month time series of sub basin actual evapotranspiration. Optimising a total of
53 variables using 920 monthly observations increased the r2 from 0.40 to 0.81. A validation with
measured discharges reveals that the modelled discharges are well with the one standard deviation of
the observed data.
Separate optimisation runs revealed that ETact is more sensitive to the groundwater and
meteorological parameters than the soil and land use parameters. On sub basin level the ETact shows
least response to the land cover dependent maximum leaf area index.
A first important step when calibration on ETact is to ensure that ETref is optimised first. We show that
a near perfect match between observed and measured ETref by small adaptations (average = 3.5%)
of monthly radiation.
It is also concluded that at the HRU level more work is required to fine-tune the calibration procedure.
The calibration is only reliable at the spatial and temporal scale on which the observations, used in the
optimisation, are based. Future work should focus on calibration strategy that incorporates HRU level
ETact observations and discharges at a high temporal resolution in the objective function.
Recently interest in using simulation models in ungauged or sparsely gauged basins has increased
leading to some concerted actions. The most relevant is the Prediction in Ungauged Basin (PUB)
initiative; an International Association for Hydrological Sciences (IAHS) initiative for the decade of
2003-2012, aimed at uncertainty reduction in hydrological practice (Sivapalani et al., 2003). PUB

focuses the development of new predictive approaches that are based on "understanding" of
hydrological functioning at multiple space-time scales. This study provides an ET based innovative
approach at different temporal and spatial scale that fits well into PUBS science program.
Differences in space and time in ETref between SEBAL and SWAT are caused by spatial altitude
dependent operations which SEBAL performs on important parameters. The DEM is used to calculated
and correct air pressure and density and thus the psychometric constant. The DEM is used to correct
the absorbed solar radiation values, both for slope and aspect. Southern facing terrain due to the angle
of incidence, absorb more solar radiation per unit land than the Northern facing slope. Another cause
of this difference is the fact that SEBAL uses grass as a reference crop and SWAT uses alfalfa.
Traditional calibration on a limited number of discharge stations lumps all hydrological processes
together and chances on the equifinality problem are much larger. In this study we show that using
spatially distributed observations with a monthly temporal resolution provide much better results. The
success of the approach lays in the spatial and temporal isolation of the calibration problem at hand.
Information content of a time series of discharges at the outlet of a catchment is simply insufficient to
attribute deviations between observation and simulation to specific processes at a specific location at a
specific point in time. Although we have found very good results at the sub basin level on a monthly
time step, more work is required to increase the reliability of the results at HRU level and eventually
with a daily time step. Promising in this respect could be the use of a combined objective function
including both spatial evapotranspiration and discharges at a high temporal resolution.
One of the variables used in the optimisation is a monthly rainfall increment. It is generally not
common practice to vary model excitations in a calibration procedure, however we believe that in this
case it is legitimate within tight bounds, since precipitation is acquired with TRMM and subject to larger
errors than station data.
In this study the GML algorithm was used in the optimisation and there are no indications that
optimisations did not converge into a global minima. The main cause is probably that a time series of
spatial ET exhibits more linear behaviour than discharge at a limited number of locations. Moreover it
has been shown that global search algorithms require much more function calls to identify the global
minimum. The GML algorithm is much more efficient in this respect (Skahill and Doherty, 2006). Given
the fact that the best optimisation in this study required 2987 model calls this is a critical requirement.

FutureWater
Science for Solutions 63/69
9 Results
9.1 Water balances
9.2 Biomass and crop production
9.3 Scenarios

FutureWater
10 Conclusion

FutureWater
11 References
Allen, R., Pereira, L.A., Raes, D., Smith, M., 1998, Crop evapotranspiration; guidelines for computing
crop water requirements, FAO Irrigation and Drainage Paper No. 56, FAO, Rome
Bastiaanssen, W.M.A., Klaasse, A., Zwart, S., Immerzeel, W. and Droogers, P., 2006, The hydrological
flow path and options for sustainable water resources management in the overexploited Rio Bravo
Basin: A preliminary analysis from remote sensing and hydrological modeling, WaterWatch report,
Wageningen
Brutsaert, W. and M. Sugita, 1992. Application of self preservation in the diurnal evolution of the
surface
energy balance budget to determine daily evaporation, J. Geophys. Res. 97(D17): 18377-18362
Crago, R.D., 1996. Conservation and variability of the evaporative fraction during the daytime, J. of
Hydr. 180: 173-194
Doherty, J., 2004, PEST; Model-Independent Parameter Estimation; User Manual: 5th edition.
Watermark Numerical Computing
Farah, H.O., 2001. Estimation of regional evaporation under different weather conditions from satellite
and meteorological data, Ph.D. thesis, Department of Agrohydrology, Wageningen University, The
Netherlands: 170 pp.
FAO, 1974. FAO-Unesco Soil Map of the World, 1:5.000.000. Unesco, Paris.
FAO, 1995. FAO-Unesco digital Soil Map of the World and derived soil properties, 1:5.000.000. Unesco,
Paris.
Franke R. 1982. Smooth Interpolation of Scattered Data by Local Thin Plate Splines. Comp. & Maths.
with Appls. 8: 237-281.
Hargreaves, G.L., G.H. Hargreaves, and J.P. Riley. 1985. Agricultural benefits for Senegal River Basin.
J. Irrig. and Drain. Engr. 111(2):113-124.
Hooghoudt, S.B. 1940. Bijdrage tot de kennis van enige natuurkundige grootheden van de grond.
Versl. Landbouwkd. Onderz. 46: 515-707.
Immerzeel, W.W. and Droogers, P., 2005, Exploring evaporation reduction in the Hai basin; analysis
using the SWAT model, FutureWater report, Wageningen
Lillesand, T.M. and Kiefer, R.W., 2000, Remote Sensing and Image Interpretation, 4th edition, John
Wiley & Sons, New York

Monteith, J.L. 1965. Evaporation and the environment. p. 205-234. In The state and movement of
water in living organisms. 19th Symposia of the Society for Experimental Biology. Cambridge Univ.
Press, London, U.K.
Neitsch, S.L., Arnold, J.G., Kiniri, J.R., Williams, J.R., 2001, Soil and Water Assessment Tool Theoretical
Documentation Version 2000. Texas: Grassland, Soil and Water research Laboratory and Blackland
Research Center.
Nicols, W.E. and R.H. Cuenca, 1993. Evaluation of the evaporative fraction for the parameterization of
the surface energy balance, Water Resources Research 29(11): 3681-3690 NOAA National Weather
Service, West Gulf River Forecats Center.
Priestley, C.H.B. and R.J. Taylor. 1972. On the assessment of surface heat flux and evaporation using
large-scale parameters. Mon. Weather Rev. 100:81-92.
Ritchie, J.T. 1972. A model for predicting evaporation from a row crop with incomplete cover. Water
Resour. Res. 8:1204-1213.
Shuttleworth, W.J., R.J. Gurney, A.Y. Hsu and J.P. Ormsby, 1989. FIFE: the variation in energy
partitioning at surface flux sites, IAHS Publ. 186: 67-74
Sivapalani M, Takeuchi K, Frank SW, Gupta VK, Karambiris H, Lakshmi V, Liang X, McDonell JJ,
Mendiondo, EM, OConell PE, Oki T, Pomeroy JW, Schertzer D, Uhlenbrook S, Zehe E. 2003. IAHS
Decade on Predictions in Ungauged Basins (PUB), 20032012: Shaping an exciting future for the
hydrological sciences. Hydrological Sciences 48: 857-880
Skahill BE, Doherty J. 2006. Efficient accommodation of local minima in watershed model calibration.
Journal of Hydrology (in press).
Smedema, L.K. and D.W. Rycroft. 1983. Land drainageplanning and design of agricultural drainage
systems, Cornell University Press, Ithica, N.Y.
Tucker, C.J., 1979, Red and photographic infrared linear combination for monitoring vegetation.
Remote Sensing of Environment, 8, 127-150

10 1 1 512 1747 PDF

Uploaded by

Copyright:

Available Formats

10 1 1 512 1747 PDF

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

10 1 1 512 1747 PDF

Uploaded by

Copyright:

Available Formats

Remote Sensing and hydrological modeling

of the Upper Bhima catchment

Research paper no. 3

NATIONAL USER SUPPORT PROGRAMME (NUSP)

The objectives of the NUSP are:

To stimulate the (inter)national service market based on space-based derived

To assist in the development of a national Geo-spatial data and information

To supply information to the general public on national and international space-based

Authors: W. Immerzeel (FutureWater)

This research was funded by: NIVR/SRON/NWO User Support Programme 2

Generaal Foulkesweg 28 6703 BS Wageningen Tel: 0317-460050

2/69 FutureWater Science for Solutions

FutureWater Science for Solutions 3/69

4/69 FutureWater Science for Solutions

2 SOIL AND WATER ASSESSMENT TOOL 9

3 LAND USE CLASSIFICATION 15

5 PRECIPITATION MAPPING USING REMOTE SENSING 33

6 MODELING THE UPPER BHIMA CATCHMENT 39

FutureWater Science for Solutions 5/69

6/69 FutureWater Science for Solutions

FutureWater Science for Solutions 7/69

8/69 FutureWater Science for Solutions

2 Soil and Water Assessment Tool

FutureWater Science for Solutions 9/69

Figure 1: Water balance SWAT.

10/69 FutureWater Science for Solutions

for estimating potential evapotranspiration: Hargreaves (Hargreaves et al., 1985), Priestley-Taylor

2.3 Crop growth

Figure 2: Parameterization of crop production processes for estimating crop yiel

FutureWater Science for Solutions 11/69

12/69 FutureWater Science for Solutions

FutureWater Science for Solutions 13/69

Figure 5: Schematic representation reservoirs in SWAT (Neitsch et al, 2001)

14/69 FutureWater Science for Solutions

3 Land use classification

Figure 6 Land use map based on Krishna basin classification

FutureWater Science for Solutions 15/69

3.2 Remote Sensing and land use classification

16/69 FutureWater Science for Solutions

The Box classifier method

FutureWater Science for Solutions 17/69

In this study a combination of spectral rationing and multi-temporal unsupervised classification is

18/69 FutureWater Science for Solutions

3.4 Classification of the Upper Bhima sub basin

Image # Date Acquisition time (GMT)

FutureWater Science for Solutions 19/69

Figure 8: 16 MODIS NDVI image used in the land use classification

Class name Cluster Slope Area (ha) Area %

20/69 FutureWater Science for Solutions

Table 2: Specifics of the land use classes

Figure 9: LULC map Upper Bhima

Class name Area % Area %

FutureWater Science for Solutions 21/69

22/69 FutureWater Science for Solutions

FutureWater Science for Solutions 23/69

Figure 11: Main terms of the Surface Energy Balance

24/69 FutureWater Science for Solutions

FutureWater Science for Solutions 25/69

H = {a cp u* (a + b T0)} / {ln(z2/z1) - h (z2,L) + h (z1,L)} (3)