Effect of Environmental Changes on Future Hydrology of the Northern Taiwan using Remote Sensing

(1)

航測及遙測學刊第十五卷第 3 期第 243-262 頁民國 99 年 9 月 243

Journal of Photogrammetry and Remote Sensing Volume 15, No.3,September 2010, pp.243-262

1 Postdoctoral Fellow, Research Center for Environmental Changes, Academia Sinica Received Date: Mar. 30, 2010

2 Professor, Department of Landscape Architecture, Chinese Culture University Revised Date: Apr. 29, 2010

3 Professor, School of Forestry and Resource Conservation, National Taiwan University Accepted Date: May. 10, 2010

4 Professor, School of Tourism, Ming-Chuan University

*Corresponding Author, Phone: 886-2-28610511 ext.41512, E-mail: zqq@faculty.pccu.edu.tw

Effect of Environmental Changes on Future Hydrology of the Northern Taiwan using Remote Sensing

Chih-Da Wu¹ Chi-Chuan Cheng ^2* Hann-Chung Lo ³ Yeong-Kuan Chen ⁴

ABSTRACT

This study aims to improve on stream flow simulation using the Generalized Watershed Loading Functions (GWLF) model by including remote sensing techniques to estimate the cover coefficient (CV), as well as integrating the SEBAL model, the CGCM1 model, and the Markov model to predict land-use and ET changes. Moreover, the results were adopted to assess the effect of environmental changes on future hydrology or the north of Taiwan. The processes include land-use classification using hybrid approach and Landsat-5 TM images, a comparison of stream flow simulations using the GWLF model with two CV values derived from remote sensing and traditional methods, and finally the prediction of future land-use and CV parameters for assessing the effect of land-use change and ET change. The results indicated that the study area was classified into seven land-use types with 89.09%

classification accuracy. The stream flows simulated by two estimated CVs were different, and the simulated stream flows using the remote sensing approach presented more accurate hydrological characteristics than a traditional approach. In addition, the consideration of land-use change and ET change indeed affected the predicted stream flows under climate change conditions. The results of the hydrology analysis based on the SRES scenarios of CGCM1 model predicted that the river flows of north Taiwan will become greater due to the effects of climate change, land-use change and ET change.

Therefore, the results obtained from this study can be extrapolated to the future studies of global environmental change and water resource management.

Keywords: Evapotranspiration cover coefficient, Environmental change, Stream flow simulation, GWLF model, Remote Sensing

1. INTRODUCTION

Considerable attention has been given for investigating the interactions among climate change, human activities, and nature ecosystems around the world (Arnell and Reynard 1996;

Manning and Nobre 2001). To the hydrology system, precipitation and evapotranspiration (ET) are the main driving forces for ground hydrology cycles (Tung and Haith 1995), but quantity and distribution of these two factors have been changed due to the increase of average temperature of the earth. This phenomenon also affects other hydrological components such as infiltration, and stream flow, moreover, to disturb the available water resources for human.

If the current trend does not change, the impact of global warming on future climatic condition and hydrology process would become a major concern (Wu and Haith 1993). Estimates of global warming are usually based on the application of General Circulation Models (GCMs), which attempt to predict the impact of increased atmospheric CO2 concentrations on climatic variables. Furthermore, in order to assess the sensitivity of hydrological regimes to the climatic changes associated with global warming, previous studies often rely on the GCMs coupled with a stream flow model, and have used these models to predict the impact of climate change on hydrologic characteristics in different areas (Bhakdisongkhram et al. 2007;Yu et al. 2002). In recent years, the Generalized

(2)

244 Journal of Photogrammetry and Remote Sensing Volume 15, No.3, September 2010

Watershed Loading Functions (GWLF) model has been widely applied to estimate human, natural, and climate effects on hydrology systems because parameters in the GWLF model can be adjusted according to the land-use types, soil characteristics and climate conditions of a watershed (Haith and Shoemaker, 1987; Haith et al., 1992; Tung and Haith 1995; Tung 2001;

Tien 2003). During the GWLF simulated procedures, the amount of ET of a watershed was calculated using an evapotranspiration cover coefficient (CV). The CV for each land-use type is defined as a ratio of actual ET to potential ET.

However, the estimation of actual ET and CV on a large spatial scale is problematic. Traditionally, CV can be determined from the published seasonal values based on crop types such as those given in the user’s manual of the GWLF model (Haith et al., 1992). This approach often requires estimates of crop development (e.g., planting dates, time to maturity, etc.) which may not be available. Moreover, a single set of consistent values is seldom available for all of a watershed’s land-use. A cursory setting of CV would influence the accuracy of stream flow simulations (Haith et al., 1992; Davis and Sorensen, 1969).

The increasing availability of remote sensing technology now produces satellite images that can easily and effectively provide large-scale spatial and temporal surface information. For hydrology studies, ET can be computed without quantifying other complex hydrological processes through remote sensing techniques (Morse et al., 2000). Thus, previous researches (Menenti and Choudhury, 1993;

Laymonet al., 1998; Mauser and Schadlich, 1998; Morse et al., 2000; Chen et al., 2006) adopted remotely sensed data to calculate the energy balance parameters such as surface temperature, net radiance, sensible heat flux, soil heat flux, and then estimated the ET according to these parameters. If the CV and ET estimated from remote sensing techniques could be integrated into impact of climate change on hydrology, it will improve the predictability of climate change effects on water flow changes.

In addition to the climatic factors, land-use change also strongly affects the hydrology cycle.

In the GWLF model, parameters such as curve number (CN) and CV are related to the land-use status of a watershed. Most previous studies assumed that catchment land-use remains consistent over the periods of time (Arnell and

Reynard 1996). This assumption may reduce the accuracy of model prediction. Many existing spatial simulation models have been applied in various fields (Muller and Middleton 1994;

Turner 1993). In those models, the Markov model is widely used on the projection of land-use change. In a Markov model, the change in an area is summarized by a series of transition probabilities from one state to another over a specified period of time. These probabilities can be subsequently used to predict the land-use properties at specific future time points (Burham 1973). Many researchers have applied the Markov model to monitor the land-use and landscape change (Cheng et al. 2005; Lindsay and Dunn 1979; Muller and Middleton 1994;

Turner 1993). However, few studies integrate Markov predictions into hydrological assessments under changing climate conditions.

Compared with the traditional stream flow simulation, which calculates a CV using the published reference values and without involving the effect of land-use change, our present efforts presents an approach to estimate the CV by Markov model and surface energy balance algorithm for land (SEBAL;

Bastiaanssen et al. 1998) model and considering of future land-use status and ET change. Our study was motivated by three questions: Is the accuracy of stream flow simulation improved by using the CV estimated from remote sensing? Is the integration of SEBAL model and Markov model a feasible scheme to predict the future land-use and ET parameters for stream flow simulations? Does the consideration of land-use change and ET change affect the results of hydrology assessment under climate change conditions in north Taiwan?

2. MATERIALS AND METHODS

2.1. Study Area and Materials

The study area selected for the empirical analysis was the north of Taiwan (Fig. 1). There were four mainstreams (Dan-Shui river, Kee-Lung river, Xin-Dian river, and Da-Han creek) in this area and divided the region into seven watersheds. The area covers 734589.7 ha and involves five counties (Kee-Lung, I-Lan, Taipei, Tao-Yuan, and Hsin-Chu). Several industrial or scientific centers in this area demand water resource. Therefore, it is important to realize how the hydrology system

(3)

would be Fiv study, L digital te inventor flow obs TM ima Novemb adopted.

spectral the midd the blue, IR; band above si Band 6 i The exte the raw maps. Se provided bureau. N generate Adminis regarded

Effect of Env

e changed in ve kinds of d

andsat-5 TM errain mode ry data in 199

served recor ages from d ber 25, 199 . Landsat-5 bands rangi dle infrared , green, and r d 5 and band ix bands ha is the therma ent of the st images and econdly, DTM d by the ae

National land ed by the

stration, Min d as the par

Fig. 1 T

F

Chih-Da Wu, vironmental Cha

n the future.

data were co M images of el (DTM), n 95, daily wea rds. Firstly, different date 5; January TM images ing from the (mid-IR). B red bands; ba d 7 are the m

s a 30m gro al IR with 12 tudy area w d used to ge

M with 40 m erial survey d-use invent e Departm nistry of th ragon for ev

The study are

Fig. 2 Locatio

Chi-Chuan Ch anges on Future

onsulted for 1995 and 20 ational land ather and str three Lands es (July 20 4, 2002) w s contain se e visible blu Bands 1 to 3 and 4 is the n mid-IR. Each ound resolut 20 m resolut as clipped f enerate land m resolution office, fore ory data in 1 ment of L he Interior valuating im

ea: the north

on of the sele

heng, Hann-Ch e Hydrology of

this 002, d-use ream sat-5 and were even ue to

are near h of tion.

tion.

from d-use was estry 1995 Land are mage

clas sele clim pre met Cyu and wat area sam the flow flow Alt stat sim Sha met thes hyd loca loca stat

hern Taiwan c

ected meteor

hung Lo, Yeon f the Northern T

ssification. T ected as the matic condi cipitation d teorological u-Chih, Sin- d Sih-Yuan tershed, whi a in norther mple for mod Cyu-Chih m w observati w station we hough the tions is less milar elevati ang-Guei-Sh teorological se two drological a ation of stu ation of the tions.

covered by L

rological and

ng-Kuan Chen Taiwan using Re

The period fro e baseline, r

tions”. Dai ata were co

stations. The -Wu, Cyue-C climatic stat ch occupies rn Taiwan, w del validation meteorologic

ions at Sh ere collected

distance be then one km on (Cyu-Ch an-Ciao stati

stream flow stations is analysis. Fi

udy area. F selected met

Landsat-5 TM

d flow station

n emote Sensing

rom 1995 to 2 representing ily temperat ollected fro ese are: Fu-G Ci, Mei-Hua

tions. The D s the largest

was assigne n. Climatic r cal station an hang-Guei-S d from 1995

etween two m, both stati hih station:

ion: 70.9m).

w data collec s compara igure 1 sh Figure 2 sh

teorological

M image

ns

245

2002 was

“current ture and

m seven Guei-Jiao, a, Yi-Lan Dan-Shui

drainage ed as the records at nd stream han-Ciao to 2002.

o gauged ions have : 90.0m;

Thus the cted from

ble for hows the

hows the and flow

(4)

2.2 Methods

2.2.1 Land-use classification using hybrid algorithm

Hybrid classification is an integrated algorithm which combines the advantages from both traditional supervised and unsupervised approaches (Lillesand and Kiefer 2000; Lo and Choi 2004). The analytical procedures included 4 steps: (1) Eight blocks were firstly selected from the study image according to ground land-use information. Each block contained at least three to four kinds of land-use types; (2) In the unsupervised stage, the selected blocks were clustered into spectral subclasses by unsupervised classification and then merged or deleted subclass signatures as appropriate based on transformed divergence (TD), see equation 1.

TD ranged from 0 to 2000. If two classes can be separated completely, then the TD approaches 2000; (3) Spectral signatures obtained from each block were then combined and integrated into a single spectral signature; (4) In the supervised stage: supervised classification method with the single spectral signature were finally applied to generate the land-use map of the study area.

 

  



i j i j i j ^T



i i

m m m m Cov Cov

tr

Cov Cov

Cov Cov tr D

D TD



 













) 2 (

1

) )(

2 ( 1

) 8 / exp(

1 2000

1 1

1

1 (1)

where, TD was transformed divergence (-);

D was divergence; Covi was covariance matrix of class i; m was mean vector of class i; and

 

A

tr was sum of the diagonal line of matrix A.

To evaluate the result of land-use classification, test areas for each cover type were selected from the generated map of November 25, 1995. All test areas were used to compare with the national land-use inventory data, and the classification accuracy was then calculated.

The same procedures were adopted to generate the land-use maps of July 20, 1995 and January 4, 2002.

2.2.2 CV Calculation using various methods

SEBAL is an image processing model that estimates actual ET by solving the terms of the surface energy balance derived from the visible,

near-IR, and thermal-IR bands of the electromagnetic spectrum (Oberg and Melesse, 2004). By the law of conservation of energy, the equation for the inputs, losses, and storage of energy for a surface on earth is expressed in Bastiaanssen et al. (1998) as:

H G R

LE



_n

 

(2)

where, LE was the latent energy consumed by ET, R_n was net radiation (W/m²); G was the soil heat flux (W/m²); H was the sensible heat flux to the air (W/m²). LE was converted into ET, expressed as a depth of water per time, by dividing by the latent heat of vaporization.

In SEBAL procedures, R_n was estimated based on the following relationship (Bastiaanssen et al., 1998; Oberg and Melesse, 2004):

) 1 ( )

1 (       

₀



_s_ _L_ _L_ _L_

n R R R R

R (3)

where, R

S↓was the incoming direct and diffuse shortwave solar radiation that reaches the surface (W/m²); α was the surface albedo, the ratio of reflected radiation to the incident shortwave radiation; R

L↓was the incoming longwave thermal radiation flux from the atmosphere (W/m²); R

L↑was the outgoing longwave thermal radiation flux emitted from the surface to the atmosphere (W/m²), ε

o was the surface emissivity, the ratio of the radiant emittance from a gray body to the emittance of a blackbody.

Soil heat flux (G) was the rate of heat storage to the ground from conduction. In the SEBAL model, an empirical relationship for G was given as:

Rn

NDVI

G

 0 . 30 ( 1  0 . 98

⁴

)

(4) where, NDVI was the normalized difference vegetation index.

Sensible heat flux (H) was the rate of heat loss to the air by convection and conduction due to a temperature difference. The calculated equation was as bellow:

(5)

Chih-Da Wu, Chi-Chuan Cheng, Hann-Chung Lo, Yeong-Kuan Chen 247 Effect of Environmental Changes on Future Hydrology of the Northern Taiwan using Remote Sensing

ah p

r dT H



C



(5)

where, ρ was the density of air (kg/m³); c

p was the specific heat of air (1004 J/kg/K); dT was the difference in temperature between the surface and the air (K); and r

ah was the aerodynamic resistance (s/m). To calculate dT, the inverse of equation (6) was considered:

p ah

C r dT H



 

( 6 )

Therefore, during the SEBAL process, the user calculated dT at two extreme “anchor pixels”

by assuming values for H at the reference pixels.

The reference pixels were carefully chosen so that at these pixels one can assume that H approximate zero at a very wet pixel (i.e., all available energy (R_n - G) is converted to ET), and that LE almost equals zero at a very dry pixel, so that H = R_n - G. These assumptions from the selected pixels provided endpoints for values and locations for H so that a relationship for dT can be established.

Once the values of H and G were calculated, the latent heat flux (LE) can be calculated from equation (2). This LE represented the instantaneous evapotranspiration at the time of the Landsat overpass. Following the computation of the evaporative fraction at each pixel of the image, one can estimate the 24-hour evapotranspiration for the day of the image by assuming that the value for the evaporative fraction (Λ) is constant over the full 24-hour period (Bastiaanssen et al. 1998). The evaporative friction is calculated for the instantaneous values in the image as:

G R

H G R

n n



 



(7)

where, the values for R_n, G, and H were instantaneous values taken from processed images. The 24 hour actual evaporation is calculated by the following equation:



) (

86400 ₂₄ ₂₄

24

G

ET   Rⁿ  ( 8 )

where, ET₂₄was daily actual evapotranspiration (cm/day); R_n24 was daily net radiation; G₂₄ was

daily soil heat flux; 86,400 was the number of seconds in a twenty-four hour period; and λ was the latent heat of vaporization (J/kg). The latent heat of vaporization allowed expression of ET24

in cm/day.

Two methods were adopted to estimate CV parameters in this study. The first method was proposed by this study, that is, to calculate the CV using remote sensing techniques (called RSCV). This procedure combines the results of land-use classification and actual ET obtained from SEBAL model, and to compute the CV for each land-use type. However, to estimate the overall CV for the study area, the percentage of land-use types was used as a weighting factor.

Due to the seasonal changes of vegetation’s structure and growing status, the values of CV should show seasonality. For the purpose of CV determination in this study, May to October was regarded as the wet season, and November to April as the dry season as suggested by Tien (2003) and Haith et al. (1992). Cover coefficients estimated from the actual ET of July 20, and November 25, 1995 were applied to represent the wet and dry season’s CV values, respectively. In addition, two cover types such as shadow and cloud within the study area had to be masked out because shadow or cloud can considerably drop the thermal band readings and cause large errors in calculating actual ET and CV (Morse et al. 2000). The equations for calculating RSCV were as follows:

t i i

PET

cvrs



ET

(9)









^m

i cvrsi wi

RSCV

1

(10)

where, cvrsi was CV of class i derived from remote sensing techniques; ETi was actual ET of class i (cm/day); RSCV was the overall CV derived from remote sensing; wi was percentage of class i; PETt was potential evapotranspiration (cm) as given by Hamon (1961):

273 021 .

0 ² ₀

 

t

t t

t T

e

PET H

(11)

where, Ht was the number of daylight hours per day during the month containing day t; e0t was the saturated water vapor pressure in millibars

(6)

on day t and T_t was the temperature on day t (℃).

Saturated water vapor pressure could be approximated as in Bosen (1960):

00136 . 0 8 . 4 8

. 1 000019 .

0

8 . 0 ) 8072 . 0 00738

. 0 [(

8639 .

0 33













t

t t

T

e (12)

The second method was the traditional approach which was based on the reference values published in the GWLF manual (Tung and Haith 1995; Davis and Sorensen 1969). The procedure was to specify the CV of each land-use according to the published reference values, and then similar to RSCV, calculate the overall CV (called REFCV). The meaning difference between cvrs derived by actual evapotranspiration and CV based on reference values is that cvrs represents the truly land-use and hydrodynamic of the objective watershed.

The other advantage of using remote sensing approach is to derive an appropriate parameter to the local area. The equation of calculating REFCV is as below:









^m

i

i w

cvref REFCV

1

(13) where, REFCV was the overall CV determined by published reference; cvrefi was CV of class i derived from GWLF manual; wi was percentage of class i.

2.2.3 Weather generation and climate change modification

The Thiessen method is a graphical method, which adjusting for the non-uniform location of gauged stations by determining their area of influence. The method is based on the construction of areas of influence centered on each point measurement, the so-called Thiessen polygons. The measurement for each point is then taken to be representative of the variable on its respective area of influence. We employed this method for the spatial integration of the seven point meteorological observations within the study area. The adjusting equations were as follows.













n

i i i

mean n i

i i mean

A RAIN RAIN

A TEMP TEMP

1

(14)

where, TEMPmean was the integrated areal temperature (℃); TEMP_iwas daily temperature of gauged station i; RAIN_meanwas the integrated precipitation (cm); RAIN_iwas daily precipitation of gauged station i; and A_i was influence coefficient of gauged station i.

A stochastic weather generation model (Tung and Haith 1995; Li et al. 2006) was used toreflect possible variations in daily temperature and precipitation heights based on the local climatestatistics from 1995~2002. Daily temperature was calculated using the first-order autoregressiveequation given in Pickering et al.

(1988) as follows.

2

1 ) 1

(   



    

 _Tm _i_ _Tm _i _T

i TEMP

TEMP (15)

where, TEMP_i was the temperature on day _i; Tm

was the mean temperature for a period ( ) (one ℃ month herein);  was the lag-one autocorrelation coefficient of temperature during the period;  was the normal sampling deviate between 0~1;

T was the standard deviation of temperature during the period.

Daily precipitation was determined by a conditional probability for the occurrence of wet and dry days based on the recorded climate sequences. Rainfall amount was generated as a logarithmic distribution function stated as below (Hong 1997):

 

^I

 

^RN

 

RAIN_i

 

_P

  ln 1 

(16)

where, RAINi was rainfall amount on day i (cm);

_P(I) was the mean rainfall on month I (cm); and RN was a random number, the range of RN was

“0＜RN＜1”.

As for the simulation of future weather data, the normalized predictor variables which were exported from the CGCM1 model (The First Version of the Canadian Global Coupled Model;

Flato et al. 2000) under SRES A2 and B2 experiments were selected to predict the possible climate change conditions. Scenario parameters acquired from the grid point N24.1°E120.0° on CGCM1 model was adopted in this study because the distance between the point and study area is the nearest.

To simulate future climatic data, we first used the weather generation model to generate a

(7)

100 yea statistic.

outputs, ( T) an△ applied sequence data was mid-term watershe weather

, ,

 

  P R

T T

t m t

m t m t

where, temperat P_t,m were precipita differenc precipita 2.2.4 St GW

His from dif GWLF m on curre records gauged simulate The regr accuracy paramete

Fig. 3

Effect of Env

ar weather re Two mon such as the nd the ratio to modify e. Finally, th s used to asse m and long ed hydrology data were st

,  





 RP

T

m m

T’_t,m and ture and pre

e the histori ation on day ce of tem ation obtaine tream flow WLF mode storic weath fferent metho model to sim

nt condition taken from station were ed accuracy

ression analy y of simulate

ers.

3 Water balan

ecord based nthly mean e difference

of precipita y the gen he simulated ess the impac g-term clim y. The modif ated as equat

&











m m t

R’_t,m were ecipitation on

ic recorded t y t; Tm and mperature ed from CGC

simulation el

her data and ods were sub mulate the stre s. The obser

the Shang- e applied to using regre ysis was used ed flows betw

nce function

on the hist ns of CGC of tempera ation (RP), w

erated wea d future wea cts of short-t mate change

fied function tions (17).

12

~ 1

12

~ 1



the modi n day t; T_t,m

temperature d RPm were

and ratio CM1 forecast n using the

d CV obtai bstituted into eam flows ba rved stream f

Guei-Shan-C investigate ession analy d to compare ween various

of the GWL

toric CM1 ature were ather ather term,

on ns of

(17)

ified and and the of ts.

ined o the ased flow Ciao the ysis.

e the s CV

mo div sha (Do disc satu uns area Ut_

St_

whe sha mo at t and surf into t (c wat disc 200 foll SF whe (cm was

LF stream flo

The stream del is show ided into allow saturat

ouglas and charge to th urated zone saturated an a as follows:

t

t R

U



1



t PC

S



1



ere Ut and allow satura istures (soil he beginning d D_t were rai

face runoff, o the deep sa cm). Daily st tershed surf charge to the 03; Tung et

low:

t t

t G Q

F

 

ere SFt was m); Gt was th

s the surface

w model (mo

m flow por wn in Fig.

lumped par ted, and de Shoemaker he stream i

. Daily wa d shallow

t Q

M

 



t t

t G D

C

 

St were t ated zone a

water in exc g of day t; R_t infall, snowm

evapotrans aturated zone ream flow is face runoff e stream (Ch al., 1999). T

the stream f e groundwat

runoff (cm).

odified from

n emote Sensing

rtion of the

3. Soil mo rameter uns eep saturate

1987) Grou is from the ater balance saturated so

t

t ET P

Q

 

the unsatura available fr cess of wilti

t, M_t, Q_t, ET_t melt water, w spiration, pe e, respectivel s given by th f and grou heng et al., 20

The equation

flows of a w ter discharge

.

m Haith et al.

249

e GWLF oisture is saturated,

d zones.

und-water shallow for the oil zones

PCt (18) (19) ated and rom soil ng point)

t, PC_t, G_t, watershed ercolation ly, on day he sum of und-water 007; Tien,

n was as

(20) watershed e (cm); Qt

1992)

(8)

Groundwater discharge is estimated by assuming shallow saturated zone as a linear reservoir using the follow equation:

t

t S

G

 r 

(21)

where, r was the recession coefficient; St was the storage of shallow saturated zone (cm).

To estimate surface runoff, the GWLF model used the curve number method to calculate runoff volumes with the considerations of land-use status and antecedent soil moisture.

The equations were as follows:

) 8 . 0 (

) 2 . 0

(

²

t t

t t t

W R

W Q R



 

(22)

4 . 2540  25



t

t CN

W CN

 100

(23)

where, Q_t was the surface runoff (cm); R_twas the daily precipitation (cm); CN was the curve number; W_t was the maximum soil runoff load (cm).

Curve number was selected as functions of land use types, soil texture and antecedent moisture. Curve number for antecedent moisture conditions 1 (driest), 2 (average) and 3 (wettest) were CN1, CN2, and CN3, respectively. CN2 could be determined by referring to the GWLF manual, and CN1 and CN3 could be calculated from CN2. The functions of CN1 and CN3 were as follows (Haith et al., 1992):

2 0059 . 0 4036 . 0 3 2

2 01334 . 0 334 . 2 1 2

CN CN CN

 

 

(24)

Calculation of ET was the most important component of this study. In the GWLF model, ET was estimated as a function of the atmospheric and surface characteristics of a watershed, as equation (25).

t

t CV PET

ET

 

(25)

where, ET_t was evapotranspiration (cm/day); CV was cover coefficient; PET_t was potential evapotranspiration (cm).

2.2.5Predictions of land-use change and future CV values

In the Markov model, it was assumed that the land-use changes of the study area could be depicted as a Markov process (Aaviksoo 1995;

Alig 1985). A transition matrix, in which the element Tij represents the amount of land-use change from type i to type j during the period from November 25, 1995 to January 4, 2002, were derived from the land-use maps. The transition probability Pij, which represents the factions of land-use changes on each land-use type, was estimated by:









^m

j ij

ij Tij T

P

1

; i = 1, 2, ... m, (26) where, Pij was the transition probability from type i to type j; Tij was the amount of land-use change from type i to type j; m was the number of land-use types.

To determine whether it is appropriate to apply the Markov model to the observed land-use changes, Goodman’s Chi-squared statistic (Goodman 1968) was used to test the null hypothesis that the land-use conditions in 1995 and 2002 were independent of each other.

If the calculated α² is larger than the reference value, we reject the hypothesis that land-use change during the period is a random procedure.

In other word, the transition process is a Markov chain procedure. We can apply the Markov model to predict the possible land-use change of the area.

 

² ²

1 1

2 ln( ) ; ( 1)







 





 

m df A

P T

m i

m j

j ij



ij (27)

where, the definitions of Tij and Pij were same as in equation (23), Aj was the fraction of land-use in each of the land-use types in 1995; df was the degree of freedom.

Assuming that the transition probabilities will remain constant in the future, the Markov model was then used to project the land-use at the next stage:

1

2 t

t P n

n

 

(28)

where, P was the transition probability matrix, nt1 and nt2 were column vectors denoting the

(9)

fractions of land-use types at t₁ and t₂.

The land-use status in 2030 of short-term, 2058 of mid-term, and 2086 of long-term future were predicted to represent the land-use conversions in the next three periods (28 years each). Furthermore, the results were integrated with the climate change data exported from the CGCM1 model to simulate the future values of CV according to equation (9) and (10). Two assumptions were made in the calculation of future CV values. The first assumption was that the land-use change of the northern Taiwan was stable. Under this assumption, the Markov model was applied to predict the percent of each land-use type using a transition probability matrix considering the land-use change from 1995 to 2002. The second assumption was that the actual ET of land-use types was invariable in the future. This simplification allows us to estimate a CV which was determined by the land-use status as shown in equation (10).

Finally, the actual ET estimated from remote sensing, the PET calculated using the future temperature, and the land-use change projected by the Markov model were integrated together to calculate the future CV value.

2.2.6 Assessment of land-use and ET changes on future stream flow simulation

A comparison of the parameters of the two methods is shown in Table 1. Future temperature and precipitation defined by SRES scenario were entered into the GWLF model to predict the stream flows based on climate change condition using various land-use and CV parameters. The objective of this process was to compare the difference of simulated flows between two approaches. The first approach is traditional approach, current land-use data was adopted and the CV parameters were acquired from the published references without considering its future change. Both of the two parameters were assumed to be consistent during the simulated procedures. In the proposed approach, future land-use data was projected using the Markov model and the future CV values were derived from the integration of SEBAL model, CGCM1 model, and Markov model. Moreover, the simulated stream flows between future and current conditions were compared for assessing the impacts of short-term, mid-term and long-term climate

change and land-use change on the hydrology of the northern Taiwan.

3. RESULTS

3.1 Land-use classification and validation

The study area was classified into seven categories using the hybrid classification method, forest, building, farmland, fallow farmland, water, cloud, and shadow. National land-use investigation data from 1995 was first combined to fit the objective categories and applied to compare with the generated maps for classification accuracy assessment. Excluding shadow and cloud, test area for each land-use type was selected and used to calculate the classified accuracy based on these examinations.

Table 2 shows the classification accuracy based on the test areas.

Table 2 implied that most errors occurred when categorizing farmland. This might be due to the fact that building and fallow farmland both have similar characteristics on spectral reflectance making the two types difficult to separate well. However, the user’s accuracy of forest, buildings, and water was 98.55%, 88.21%, and 98.01%, respectively. The overall accuracy and kappa statistic was 89.09% and 0.8525, suggesting that hybrid classification is a suitable approach to generate a land-use map. The same procedures were adopted to create the land-use maps of the northern Taiwan on July 20, 1995 and January 4, 2002, and the results are shown in Fig. 4. A further comparison of land-use status among three different dates indicated that in July 1995, forest occupied most of the study site (37.78%), then farmland (36.71%), building (10.58%), fallow farmland (7.52%), and water (7.41%) was smaller. In November 1995, forest has the largest area (34.13%), then fallow farmland (23.81%), farmland (21.67%), building (13.36%), and water (7.03%) was the smallest one. The result showed that there was a significant transition between the farmland and fallow farmland because of the cultivation and fallow time in Taiwan. As for the land-use status on January 4, 2002, forest was still the largest area (34.23%), then farmland (23.95%), fallow farmland (22.55%), building (16.11%) and water (3.16%). Obviously, the buildings have an increasing trend from 13.36% on November 25, 1995, to 16.11% on January 4, 2002.

(10)

252

Me

Trad app

Propose

Jou

Ta

ethods

ditional proach

ed approach

Table

Land t Fore Build Wate Farml

Tota Produc accura

urnal of Photog

able 1 Compa Climate c simula

CGCM1 SRES sc

2 Examinati

type Fore est 9144

ing 0 er 0

and 162

al 9306 cer's acy 98.2

Overall

a Fig. 4 Land (a: July 20,

grammetry and R

arative param change ation

model, cenario

W lan ad lan

model, cenario

Fu pr M

ion of the lan

st Buildin 4 0

9212 27 2 2472

6 11711 6% 78.66 accuracy＝8

d-use maps g 1995; b: Nov

Remote Sensin

meters of the

Land-use Without consi

nd-use chang dopted the cu

nd-use data d uture land-us rojected usin Markov mode

nd-use classi

ng Water 0 414 6539 233 7186 6% 91.00%

89.09% Ka

b generated by

vember 25, 1

ng Volume 15, N

e methods us

e status dering ge, only urrent

directly.

se status was ng the

l.

ification of N

Farmland 135 817 106 10770 11828

% 91.06%

appa statistic

the hybrid c 1995; c: Janu

No.3, September

ed in this stu E

CV values the publis without co change.

Future CV derived fr SEBAL m model, an November 25

Total U ac 9279 9 10443 8 6672 9 13637 7 40031

%

＝0.8525

c lassification uary 4, 2002)

2010

udy

ET parameter (CV value) s were acqui shed referenc onsidering it

V values wer rom the integ model, CGCM nd Markov m 5, 1995

User's ccuracy 98.55%

88.21%

98.01%

78.98%

c

)

r

ired from ces and

ts future

re

gration of M1 model.

(11)

3. 2 Cu

The and Nov To comp types, i (0.723cm water (0.267cm for Nov highest (0.185cm water (0.088cm

Tab influence stations.

integrate precipita gauged which in whole temperat potential

Effect of Env

urrent CV

e generated a vember 25, 1

pare the actu n July fore m/day), then (0.281cm/d m/day), and

vember, co value (0.3 m/day), fallo

(0.113cm/

m/day).

ble 3 was e coefficient

These coeff e the his ation observa

stations for nvolved the north Tai ture, saturate l evapotrans

Fig. 5 Act Table

V and ET

actual ET ma 995 were sh ual ET amon est had the n farmland day), fall

building (0.

nsistently f 95cm/day), ow farmland /day), a the calcu ts for seven ficients woul storical tem

ations collec r deriving a e climatic ch

iwan. The ed water vap

spiration on

tual daily eva e 3 Thiessen

Meteorolo Fu-G Cyu

Si Cy Me Y Sih

T calculati

aps from July hown as Figu

ng five land e highest v (0.530cm/d ow farml 220cm/day);

forest was then farml (0.139cm/d nd build ulated Thies n meteorolog ld be adopte mperature cted from se

weather se haracteristics mean d por pressure n July 20

a apotranspirat influence co ogical station Guei-Jiao

u-Chih in-Wu yue-Ci ei-Hua Yi-Lan

h-Yuan

ions

y 20 ure 5.

d-use value day), land

; As the land day), ding ssen gical ed to and even eries s of daily

and and

Nov the 3.

obt trad Tab (RS the REF sea from 1.24 from dry and the resu is, 1.5 valu 0.82

tion maps in oefficients of ns Thiessen

vember 25, integrated w

Table 4 an ained from t ditional appr ble 5, the ca SCV: 1.245;

value of FCV: 0.717) son or the d m the remot 45; dry seas m the traditio y season: 0.7

d to validate Dan-Shui w ults for the c the remote 18; dry seas ue than the t 26; dry seaso

b 1995 (a: Jul f seven meteo

n influence c 0.089 0.188 0.143 0.188 0.178 0.130 0.184

1995 were t weather data

nd Table 5 the proposed roaches. Com alculated CV REFCV: 0.8 the dry sea ). No matter ry season, th te sensing ap on: 0.851) w onal approac 717). To rec the followin watershed wa

calculated CV sensing app son: 0.983) s

traditional ap on: 0.754).

y 20; b: Nov orological sta coefficients

n emote Sensing

then calcula a and shown

were the C d remote sen mparing Tab Vs of the we

842) was hig ason (RSCV r what it was he CV value pproach (we were larger th

ch (wet seaso confirm this ng stream flow

as also selec CVs were sim pproach (wet still has a la pproach (we

vember 25) ations

253

ated from as Table V values nsing and le 4 with et season gher than V: 0.851;

s the wet s derived et season:

han those on: 0.842;

situation w model, cted. The milar, that t season:

arger CV et season:

(12)

Table 4 CV obtained from the remote sensing approach

Forest Building Water Farmland Fallow farmland cvrs of July 20

(=ET/PET) 1.69 0.51 0.66 1.24 0.63

RSCV of July 20 (for wet season)

1.69× 0.3778 + 0.51 × 0.1058 + 0.66 × 0. 0741 + 1.24 × 0.3671 + 0.63 × 0.0752 = 1.245 cvrs of November 25

(=ET/PET) 1.44 0.32 0.41 0.67 0.51

RSCV of November 25 (for dry season)

1.44× 0.3589 + 0.32 × 0.1299 + 0.41 × 0. 0685 + 0.67 × 0.2239 + 0.51 × 0.2188 = 0.851

Table 5 CV obtained from the traditional approach

Forest Building Water Farmland Fallow farmland

cvref 1 0 1 1 0.3

REFCV of July 20 (for wet season)

1× 0.3778 + 0 × 0.1058 + 1 × 0. 0741 + 1 × 0.3671 + 0.3 × 0.0752 = 0.842 REFCV of November 25

(for dry season)

1× 0.3589 + 0 × 0.1299 + 1 × 0. 0685 + 1 × 0.2239 + 0.3 × 0.2188 = 0.717

3.3 Validation of stream flow simulation using the GWLF model

The Dan-Shui watershed was also used to validate the GWLF model. CV values obtained from two approaches were applied to assess their effects on stream flow simulations using the GWLF model. The basis of setting model parameters was as follows. Firstly, the values of CN₂ for five land-use types were derived from the user’s manual of the GWLF model, for example, 63 for forest, 98 for building, 98 for water, 79 for farmland and 70 for fallow farmland. The recession coefficient was adopted from the advice of Li et al. (2006) using a set value of 0.1. Daylight hours which were related to the latitude of study area were also given in the GWLF use’s manual. Finally, the CV values obtained from remote sensing (wet season: 1.518;

dry season: 0.983) and the reference manual (wet season: 0.826; dry season: 0.754) were

applied in the following hydrology simulations.

The simulated daily temperature and precipitation exported from the weather generation model based on the climate statistics of the Cyu-Chih station from 1995 to 2002 were then entered into the GWLF model with the above parameters to calculate the flow series.

The result is shown as Fig. 6. A regression analysis was used to investigate the relationship between the observed and calculated flow series.

The result shown as Table 6 indicated that the regression coefficient was 0.877 when using RSCV, and 0.853 for the use of REFCV. Such a high regression coefficient implied that, even though the simulated values are overestimates compared with the observed values, the CV value obtained from the proposed remote sensing approach could represent truer stream flow characteristics than the traditional approach.

In other words, the use of remote sensing approach to derive the parameters for the GWLF model was more suitable to simulate the hydrological flow in the northern Taiwan.

(13)

3.4 La pre

Bas 25, 199 matrix o shown statistic larger th was stat indicatin transition

For strongly period b data, to variabili future la

Fig. 6

Effect of Env

and-use an edictions

sed on the la 5 and Janu of the obse in Table 6 for the two han the refer tistically ver ng that th n during the r climate cha recommend be employed dampen th ty. In our s and-use sim

6 Observed

Fi

nd CV ch

and-use map ary 4, 2002 erved land-u 6. Goodman

periods (χ²= rence value ( ry significan he processe periods are n ange simulati

ded that at for averagin he effects o tudy, the ba mulation is a

and simulat

ig. 7 Predict

anges

ps on Novem 2, the transi

use changes n’s Chi-squa

= 76742.52) (χ²= 32.00) nt (p << 0.0 s of land not random.

ion, IPCC (2 least a 30-y ng GCM ou of inter-dec ase-line data a 7-year pe

ted stream f

tions of land

mber ition s is ared

was and 001), d-use 2007)

year utput adal a for eriod

(fro the sug per stag 205 mid The the incr 200 reac obt are CV Fro Taiw mat

flow values

d-use types

om 1995 to climate cha ggestion of iods (28 yea ge. Therefor 58, and 20 ddle-term, a e predicted r prediction, rease from 02 to 38.91%

ch 62.36% i ained from M

combined t V using equa om the result wan would tter under wh

(cm/month

in 2030, 20

2002). In or ange simulati IPCC (2007 ars, approxim e, three pred 086 represe and long-ter

results are sh , the area 13.36% in 1

% in 2030, 5 in 2086. Fur Markov mod to calculate ation (4) and t, the CV va

reveal a d hich storylin

h) of the Dan

058, and 20

n emote Sensing

rder to coinc tion, we follo 7) and adop mate to 30 ye dicted years ent the sh rm land-use hown as Fig.

of building 1995 and 14 52.13% in 2 rthermore, th del and SEBA the future v d shown as alues for the decreasing ne.

n-Shui wate

086

255

cide with owed the pted four ears) as a

in 2030, hort-term, change.

. 7. From g would 4.05% in 2058, and he results AL model values of

Table 7.

northern trend no

ershed

(14)

Table 6 Transitional pixels and probabilities from 1995 to 2002 Transition from

row to column Forest Building Water Farmland Fallow farmland Column total (1995 distribution) Forest 274

(0.823) 2

(0.006) 0 44 0.132

13

(0.039) 333

Building 0 86

(1) 0 0 0 86

Water 2

(0.033) 17 (0.283)

13 (0.217)

15

(0.250) 60

Farmland 34 (0.143)

18 (0.076)

1 (0.004)

121 (0.508)

64

(0.269) 238

Fallow farmland 14 (0.059)

33 (0.139)

7 (0.030)

80 (0.338)

103

(0.435) 237

Row total

(2002 distribution) 324 156 21 258 195 954

The numbers in parentheses indicate the transitional probability.

Table 7 Predictions of future CV values

CV 2030 2058 2086

CV for wet season

(A2 storyline) 0.997 0.896 0.813

CV for dry season

(A2 storyline) 0.719 0.636 0.564

CV for wet season

(B2 storyline) 1.000 0.902 0.813

CV for dry season

(B2 storyline) 0.723 0.643 0.573

3.5 Assessment of effects of land-use and ET changes on future hydrology

Table 8 represents the difference of the predicted stream flows between the traditional approach (REFCV parameter for current condition, without considering its future change) and proposed approach (RSCV, considering both land-use and ET change in the future). The results revealed that the simulated flows using proposed approach were lower than those using traditional approach, regardless of the monthly stream flow, mean value of monthly stream flow, and annual total volume.

To further investigate how land-use change,

ET change, and climate change affected river flows of the northern Taiwan, flow series from 1995 to 2002 was estimated to represent the current hydrological condition, and then compared with the future flows, Table 9 shows the value changes for future condition to the current condition. Positive values in the table denote an increase of flow, whereas negative values represent a decrease of flow. The results indicate that the overall values in future flow change have increasing trend. That is, the integrated effects from urban sprawl, ET decline, and atmospheric changes would lead to an increase in the rainstorm and flood occurrences of the northern Taiwan.

(15)

Table 8 Difference of stream flows (cm/month) between various approaches

A2 Storyline B2 Storyline

Month Short term

Middle term

Long term

Short term

Middle term

Long term

January 0.00 0.00 0.90 0.00 0.40 0.40

February 1.00 1.10 1.80 0.80 1.10 2.80

March -0.90 0.20 0.40 -0.90 0.00 -0.40

April 0.30 0.50 1.80 0.10 1.00 1.30

May -0.50 0.80 0.60 -0.60 0.00 1.20

June -1.60 -0.40 1.90 -1.10 0.50 1.10

July -2.00 -0.60 2.50 -2.00 -1.00 1.20

August -1.50 -0.60 0.60 -2.20 -1.00 1.00

September -1.40 -0.20 0.90 -1.70 0.10 1.40

October -2.40 -0.70 -0.90 -1.90 -1.40 0.60

November -0.70 -1.20 -2.50 -1.30 -0.80 -2.80

December -0.20 0.50 0.70 0.00 -0.40 -0.10

Total -9.90 -0.60 8.70 -10.80 -1.50 7.70

Average -0.82 -0.05 0.73 -0.90 -0.12 0.64

Table 9 Stream flow changes (cm/month) of the northern Taiwan due to land-use change, ET change and climate change

A2 Storyline B2 Storyline

Month Short term

Middle term

Long term

Short term

Middle term

Long term

January 2.6 0.8 0.9 1.3 0.1 -1.9

February 6.6 1.2 -0.9 7.2 -1.8 3.8

March 2.2 1.7 2.0 2.3 -0.7 2.6

April 1.9 -1.6 0.7 0.9 0.7 1.1

May 4.4 0.2 -6.2 -0.9 0.6 -0.4

June 0.9 -6.9 -6.8 -6.3 11.9 -0.4

July 2.2 -1.4 -6.2 0.9 -0.5 -0.6

August 14.8 3.4 6.2 9.0 5.9 8.9

September 6.5 -6.2 11.0 5.6 3.3 17.4

October -0.9 -1.7 7.9 3.5 -4.1 21.5

November -2.8 2.7 -2.1 -6.0 -1.4 -1.6

December -1.0 0.9 -5.7 -2.2 -5.3 -4.0

Total 37.4 -6.9 0.8 15.3 8.7 46.4

Average 3.1 -0.6 0.1 1.3 0.7 3.9

(16)

4. DISCUSSIONS

This study proposed the integration of the SEBAL model, the CGCM1 model, and the Markov model to estimate future land-use and ET parameters for stream flow simulations assuming a climate change scenario and then to compare the simulated results with the traditional approach. It is clear that flow values derived from the proposed method were smaller than those using the traditional approach.

Readers should be aware of some limitations regarding the numbers of CV used herein. In general, CV might be varied in different seasons due to the effects of land cover change and vegetation growth. Furthermore, cloud effects are always a serious problem in Taiwan. For the above reasons, this study was able to apply only these two clear Landsat-5 images (July 20 and November 25, 1995) for representing the CV of wet and dry seasons. This might be the reason to explain why flow simulations in specified months (such as April, July, and November) revealed opposite movement with the observational flow series no matter which kind of CV was used. If images from different months or seasons are available, certainly the monthly or seasonal RSCV value can be obtained from remote sensing technique. Thus, it is suggested that further studies can focus on the comparison of the influence of CV at different months or seasons.

In this study, the future CV values derived from the integration of the SEBAL model, the CGCM1 model, and the Markov model were based on an assumption; that is, the actual ET of each land-use type would be invariable in the future. However, ET might be affected due to the change of forest structure (species, age) and urban composition (including more concrete building or vegetation). Readers should be aware of the fact that the accuracy of future CVs can be improved if the future ET values are available.

Therefore, how to estimate the future ET values will become an important issue for further study.

In addition, The CVs of the two approaches, even though the values applied in proposed approach were decreasing due to the effects of land-use change such as urban and deforestation (see Table 7), most of them were higher than the reference CV from the traditional approach (wet season: 0.826; dry season: 0.754). Higher CV values indicate higher amounts of

evapotranspiration, which would lead to the decrease of stream flows. This explains the differences of derived flow values between the proposed and traditional approaches.

The first-order Markov chain model is commonly used in the study of land use change (Baker 1989). For example, changes in a wide variety of landscapes from predominant urban to wilderness landscapes (Bourne 1971; Marsden 1983; Hall et al. 1987). During the process, the concept of “stationary” is usually assumed although a few instances have further tested if the transition is stationary or not (Bourne 1971;

Bell 1974; Bell and Hinojosa 1977). In fact, non-stationary transition has been found in the Markov chain model of vegetation dynamics on smaller land areas (Binkley 1980; Austin and Belbin 1981; Gibson et al. 1983; Lippe et al.

1985; Rejmanek et al. 1987). But if the transition is non-stationary in reality, stationary can be assumed as a heuristic device to provide answers to “what if” kinds of questions.

However, the uncertainty of land-use simulation should be considered before applying the Markov model.

Finally, atmospheric and hydrological conditions are complex and highly changeable over time. The amounts of data are an unavoidable influence on trend analysis and the predicted changes of stream flow. Therefore, the collection of basic data records over a long period is necessary to improve the reliability of climate change simulations. This study was limited by the data acquisition and adopted only eight years records (from 1995 to 2002) of gauged weather and flow observations to assess the effect of climate change on hydrology, the result might be affected more or less.

5. CONCLUSIONS

This paper has summarized the results of an investigation into the potential effects of land-use change and ET change on future hydrology simulations using GWLF model.

Climate change scenario was based on the CGCM1 model, and differences of flow value simulations between the proposed approach and the traditional approach were assessed. Northern Taiwan was selected as the study area for the empirical analysis. The results of the study indicate that stream flow simulation using remote sensing-based CV could present truer hydrological characteristics than the traditional

(17)

approach. The integration of the SEBAL model, the CGCM1 model, and the Markov model is also a feasible scheme to estimate the future land-use status and CV values for stream flow.

The consideration of land-use change and ET change indeed affects the predicted flows. The results of the hydrology analysis based on the SRES scenarios of CGCM1 model indicated that river flows of northern Taiwan would become greater due to the effects of climate change, land-use change and ET change. Therefore, the results obtained from this study can be extended to the future studies of global environmental change and water resource management.

6. REFEREMCE

Aaviksoo K (1995) Simulation vegetation dynamics and land use in a mire landscape using a Markov model. Landscape Urban Plan 31:129-142

Alig RJ (1985) Modeling acreage changes in forest ownerships and cover types in the southeast. USDA For Serv Res Paper RM-260, Rocky Mt For & Range Exp Sta, Fort Collins, Colo

Arnell NW, Reynard NS (1996) The effect of climate change due to global warming on river flow in great Britain. Journal of Hydrology 183:397-424

Austin MB, Belbin L (1981) An analysis of succession along an environmental gradient using data from a lawn. Vegetatio 46:19-30 Baker WL (1989) A review of models of

landscape changes. Landsc Ecol 2(2):111-133

Bastiaanssen WGM, Menenti M, Feddes RA, Holtslag AAM (1998) A remote sensing surface energy balance algorithm for land (SEBAL)1: Part 1 formulation. Journal of Hydrology 212-213:198-212

Bell EJ (1974) Markov analysis of land use change: An application of stochastic processes to remotely sensed data.

Socio-Econ Plan Sci 8:311-316

Bell EJ, Hinojosa RC (1977) Markov analysis of land use change: Continuous time and stationary processes. Socio-Econ Plan Sci 11:13-17

Bhakdisongkhram T, Koottatep S, Towprayoon S (2007) A water model for water and

environmental management at Mae Moh Mine area in Thailand. Water Resour Manag 21:1535-1552

Binkley CS (1980) Is succession in hardwood forests a stationary Markov process? For Sci 26:566-570

Bosen JF (1960) A formula for approximation of saturation vapor pressure over water.

Monthly Weather Reviews 88(8): 275-276 Bourne LS (1971) Physical adjustment processes

and land use succession: A conceptual review and central city example. Econ Geogr 47:1-15

Burnham BO (1973) Markov intertemporal land use simulation model. Southern J. of Agri.

Econ 5:253-258

Chen CT, Wu ST, Chiang YF (2006) Using MODIS satellite images to estimate evapotranspiration in Taiwan. Taiwan Journal of Forest Science 21(2):249-261 Cheng CC, Wu CD, Wang SF (2005)

Application of Markov and Logit models on monitoring landscape changes. Taiwan J For Sci 20(1):29-36

Cheng CC, Wu CD, Chuang YC (2007) Influence of land-use changes and climate change on stream flow simulations: a case study of the Jiao-Long watershed. Taiwan Journal of Forest Science 22(4):483-495 Davis CV, Sorensen KE (1969) Handbook of

applied hydraulics. McGraw-Hill, New York, p1272

Douglas AH, Shoemaker IL (1987) Generalized watershed loading functions for stream flownutrients. Water Resources Bull 107:121-37

Flato GM, Boer GJ, Lee WG, McFarlane NA, Ramsden D, Reader MC, Weaver AJ (2000) The Canadian centre for climate modeling and analysis global coupled model and its Climate. Climate Dynamics 16:451-467 Gibson CWD, Guilford TC, Hambler C, Sterling

PH (1983) Transition matrix models and succession after release from grazing on Aldabra atoll. Vegetatio 52:151-159

Goodman LA (1968) The analysis of crossclassified data: Independence, quasi-independence, and interactions in contingency tables with or without missing