Modelling the influence of irrigated terraces on the hydrological response of a small basin

Environmental Modelling & Software 16 (2001) 31–36 www.elsevier.com/locate/envsoft

Modelling the influence of irrigated terraces on the hydrological response of a small basin I. Gatot a, P. Perez b

b,*

, J. Duchesne

c

a CSAR, Jl H. Juanda no. 98, Bogor 16123, Indonesia CIRAD, CRES, Australian National University, Canberra, ACT 0200, Australia c ENSA, Lab. PSNGR, 65 route de St Brieuc, 35042 Rennes Cedex, France

Received 21 February 2000; received in revised form 8 May 2000; accepted 10 August 2000

Abstract A geomorphologic-unit-hydrograph based model (H2U) was used to simulate the discharge of the Kali Garang river (Central Java, Indonesia). Even if the peak discharge and concentration time fitted well the observed data, the base time was systematically underestimated. Assuming that the terraced areas delayed part of the streamflow, we first developed a hydraulic model, based on the Cascaded Reservoirs Law, to simulate the outflow of a system of n terraces. This model was calibrated with observed data collected in a three-terrace-system monitoring. Afterwards, the digitized hydrographic network was divided into 2 components: terraced and non-terraced areas. Within the latter, the H2U model simulated directly the transfer to the outlet and it was assumed that the given peak discharge was used to calibrate the global Runoff Coefficient (kr) of the catchment. Concerning the terraced areas, the delayed outflow computed by the hydraulic model was then routed to the H2U model. The total discharge was given by the sum of the convolutions. The method was applied on the Banyumanik sub-basin (75 km2). The results fit well the observed discharge, even the descending limb. But some more assumptions are needed concerning the way farmers are managing the water within the terraces or the way to assess the average values of the parameters for the entire basin.  2001 Elsevier Science Ltd. All rights reserved. Keywords: Hydrology; Indonesia; Unit hydrograph; Terrace; Linear reservoir; Simulation

1. Introduction Hydrological models using the concept of Instantaneous Unit Hydrograph (IUH) have to be calibrated with already existing rainfall and discharge records (Chow, 1964). This is a major constraint for the extension of flooding early warning systems. In Indonesia, many concerned basins have not been monitored yet or present poorly reliable data sets. In order to predict discharge from ungaged basins, attempts have been carried out in order to link the hydrological response to the characteristics of the watershed. Rodriguez Iturbe and Valdes (1979) proposed to use some geomorphological parameters to describe the instantaneous hydrograph function. These parameters

* Corresponding author. Tel.: +61-2-6249-2955; fax: +61-2-62798395. E-mail address: [email protected] (P. Perez).

mostly consist of the mainstream length value L⍀, where ⍀ corresponds to Strahler’s order, of the stream number ratio (RB), of the stream length ratio (RL) and of the stream area ratio (RA). According to these parameters, the peak of discharge (Qp) and the time to peak (tp) can be derived even to simulate IUH. Unfortunately, the authors did not find a global analytical solution applicable in the product of convolution (Gupta et al., 1980; Rodriguez Iturbe et al., 1982). More recently, Duchesne et al. (1997) proposed a deterministic model, called H2U (Universal Unit Hydrograph). The model is based on fractal characteristics of most of the draining networks. The H2U model requires two parameters: n, the order of the basin according to Strahler’s classification and L, the stream average length. The two parameters are extracted from a digitised map of the hydrological network. This model was satisfactorily tested with several indonesian basins. But in some cases, repeatedly underestimations of the recession leg of the hydrograph were detected. These basins were

1364-8152/01/$ - see front matter  2001 Elsevier Science Ltd. All rights reserved. PII: S 1 3 6 4 - 8 1 5 2 ( 0 0 ) 0 0 0 6 1 - X

32

I. Gatot et al. / Environmental Modelling & Software 16 (2001) 31–36

characterized by large areas of irrigated terraces (Gatot Sumarjo et al., 1997). Cascaded terraces can be described as series of reservoirs interconnected by spillways. Even if locals use many different shapes and types of material to build them, we define here such a spillway as a pipe made of PVC and characterized by its diameter (di). Each terrace (Ti) is described by its area (Ai) and its water level (hi(t)). To predict the output discharge (Qn) from the last terrace (Tn), we need to model the individual discharge of each terrace with time (Qi(t)), knowing the above parameters and the initial conditions (h0, Q0), according to the Cascaded Reservoirs theory developed by Dooge (1959) and then discussed by Chow (1964) and Huggins and Burney (1982). The authors have attempted to simulate the hydraulic behaviour of these terraces. After calibration of the hydraulic parameters it was decided to couple this model with H2U. First, this paper briefly presents theoretical aspects of the two models. Then, the basin characteristics are described and the assumptions taken for the simulations presented. Finally, the results of the simulations are given and discussed.

2. Theory 2.1. The hydrological model (H2U) The H2U model has been developed by Duchesne et al. (1997) from an initial analogy with the Maxwell’s theoretical distribution of the gas molecular velocities. The H2U model is based on the combination of two parameters: n and L. These parameters are directly extracted from the drainage network representation. The geomorphologic response (IUH) links the parameters n, L and a so-called gamma function as follows:

冉冊

n r(L)⫽ 2L

n/2

冉 冊

1 (n/2)−1 nL L exp ⫺ 2L n ⌫ 2

冉冊

冉冊

1 t k⌫(n) k

n−1

e−(1/k)

冕 t

Q(t)⫽S ieff(t)r(t⫺t)dt

(1)

(2)

The discharge is then calculated through the following

(3)

0

where Q=discharge (m3/s), t=time (s), S=basin area (m2), ieff=average rainfall intensity (m/s), r=probability density function. With the above description, we consider only the transfer function of the hydrological model. The production function, transforming rainfall into runoff, also uses a deterministic approach based on soil surface features and land use characteristics (Casenave and Valentin, 1989; Perez, 1994). As its adaptation to the Indonesian case is still in process, the authors have preferred to calibrate the runoff coefficient with the observed data for each flooding event studied. 2.2. The hydraulic model (cascaded terraces) Following Wu et al. (1997), the theoretical assumption is that cascaded irrigated terraces are analogue with theoretical linear reservoirs in series. According to this approach, the output discharge is calculated as a linear function of the water storage in the reservoir (Chow, 1964). Assuming a parallelepiped shape for the reservoir, the equation is: Q⫽kAh

(4)

where k=reservoir parameter (s⫺1), A=reservoir area (m2), h=hydraulic elevation (m). In order to limit the number of parameters, the system is considered to be constituted with identical terraces, i.e. the plot area, the outlet section and the initial water level are the same for each terrace. In this case, it has been demonstrated that (Nash, 1957): Qn⫽Q0 e−kt[1⫹kt/(1!)⫹k 2t2/(2!)⫹%

where n=Strahler’s basin order, L=mean hydraulic length, ⌫=gamma function. This analytical expression of the nugget of the Instantaneous Unit Hydrograph echoes the statement suggested by Rodriguez Iturbe et al. (1982) that the hydrological response of a basin depends only on some of the gross features, not on the details of the network geometry. There is an analogy between Eq. (1) and the one proposed by Nash (1957): u(t)⫽

product of convolution, assuming that L=tv (t=average routing time of the rain drops and v=flow mean velocity):

(5)

⫹k (n−1)t(n−1)/(n⫺1)!] where Qn=discharge from the nth terrace (m3/s), Q0=initial discharge from the first terrace (m3/s), t=time (s), k=reservoir parameter (s⫺1). Nevertheless, Eq. (5) describes the output discharge from a series of terraces continuously supplied with water, as Qn tends to Q0 when time is infinite. Besides, the proposed expression does not take into account natural factors such as rainfall, evaporation and infiltration, though they are essential in hydrological processes. Thus, for the first terrace, the global equation of continuity gives: EA dt⫹IA dt⫹Q dt⫽⫺A dh⫹iA dt

(6)

where E=evaporation from the terrace (m/s), I=infiltration in the soil (m/s), i=rainfall intensity (m/s).

I. Gatot et al. / Environmental Modelling & Software 16 (2001) 31–36

33

We assume here that evaporation, infiltration and rainfall are constant during the monitoring period, i.e. a single flooding event. In fact, we can nearly neglect the evaporation. As we are simulating the discharge from irrigated terraces, during the rainy season, we can also neglect the infiltration. The assumption concerning the rainfall intensity is more questionable, as flooding often occurs after non-uniform precipitations. But this point will be discussed later. According to Eq. (4) and neglecting evaporation and infiltration, Eq. (6) becomes: dh⫽⫺k(h⫺i/k)dt

(7)

In order to take into account the rainfall input, considering the new variable h⬘=h⫺i/k, Gatot Sumarjo (1999) has demonstrated that Eq. (5) becomes:

冋冉 冊 冉 冊 冉 冊 册

i kt i Qn⫽Ai⫹Q0 e−kt 1 1⫺ ⫹ 1⫺ kh0 1! kh0 ⫹

(8)

k 2t 2 i k (n−1)t(n−1) 1⫺ ⫹%⫹ 2! kh0 (n−1)!

3. Material and methods 3.1. The Kali Garang basin The experiment was located in the Kali Garang river basin (190 km2), Central Java, Indonesia (Fig. 1). The annual mean precipitation reaches 2500 mm in the downstream part of the basin and more than 4000 mm at the summit. Most of the precipitations occur from November to January. The highlands, between 400 and 2050 m elevation, correspond to the steep slopes of the Ungaran Mount, an ancient volcano dated from the Pleistocene. The geological substratum is mainly constituted with basaltic lava. Slopes range from 15 to 40%, even more near the summit. The basin includes more than 6000 ha of irrigated terraces, mainly located between 200 and 400 m elevation, permitting the farming successively of 3 to 4 crops a year (Perez et al., 1998). 3.2. Calibration of the hydraulic model In order to calibrate the reservoir parameter (k), three cascaded terraces were monitored. They were roughly shaped like rectangles and their areas were, respectively, 78.4 m2 (A1), 104.0 m2 (A2) and 89.6 m2 (A3). They were interconnected with PVC pipes (d=35 mm) and the initial water level in each terrace was nearly 90 mm (h0), before the simultaneous opening of the spillways. The outlet discharges were monitored at the same time, each 180 seconds. Water was collected in two-litre capacity measuring glasses. Between two discharge measurements, the water level in each terrace was also recorded.

Fig. 1. The Kali Garang basin and its main tributaries (Central Java, Indonesia).

3.3. Validation of the hydrological model The model was validated with discharge measurements from the Banyumanik sub-basin (79 km2). During 1997, this station had been equiped with a pressure head sensor. The data were recorded on a 6-minute time step. The rainfall was monitored at three locations (Banyumanik, Nyatnono and Pagersari) and the Thyessen’s method was used to compute the average precipitation in the basin. The geomorphological parameters were extracted from a digitised map, 1/25 000 scale, of the drainage network (Table 1). Several assumptions were necessary to run the simulations, coupling the two models: 3.3.1. Assumption 1 (A1) The basin is divided into terraced and non-terraced areas. The terraced area is constituted with one or several independent geographical sub-units depending on the design of the irrigated schemes. Each area controls the part of the hydrological network enclosed in its boundaries.

34

I. Gatot et al. / Environmental Modelling & Software 16 (2001) 31–36

Table 1 Geomorphologic parameters of the Kali Garang and Banyumanik basinsa Basin

Area (km2)

[n]

[L] (km)

[Lmax] (km)

Kali Garang Banyumanik Banyumanik terraced area Banyumanik non-terraced area

190 79 15 64

7 6 5 6

19.22 18.61 6.60 12.01

35.49 23.10 11.75 23.10

a

[RL] 2.68 2.51 – –

[Rc] 4.08 4.65 – –

[n] Strahler’s order, [L] stream mean length, [Lmax] stream maximum length, [RL] stream mean length ratio, [Rc] confluence ratio.

3.3.2. Assumption 2 (A2) Each Terraced area is represented by a set of three identical terraces in cascade, continuously supplied with water. Rainfall is first stored in the system and released after the opening of the terrace’s outlets. As the infiltration is neglected, when the storage capacity is full, the runoff coefficient equals 100%. When the water is released it enters the drainage network and the transfer is simulated by the H2U model with r(L). The parameters nter and Lter are extracted from the digitised map of the terraced area. 3.3.3. Assumption 3 (A3) The non-terraced area is most of the time a non-continuous addition of dependent patches, separated by irrigation schemes. The Instantaneous Unit Hydrograph is computed through r(L) with nnonter=ntotal and Lnonter=Ltotal⫺Lter. It is assumed that the observed peak discharge (Qp) is only due to the response of the nonterraced area. Thus, the runoff coefficient (kr) is optimised, for each flooding event, by comparing the simulated and observed values of Qp. 3.3.4. Assumption 4 (A4) The storage capacity (h0=50 mm) is the same for each terrace. The outlets are closed until the observed rainfall intensity falls to 5 mm/h. When the rainfall amount exceeds the storage capacity, the remainder is divided into unit volumes characterized by an average and constant intensity (i). This simplification enables the use of Eq. (8) to compute the output discharge from the terraced area.

4. Results and discussion 4.1. Calibration of the hydraulic model The determination of k consisted in estimating the slope of the regression curve between the observed discharge and the observed water level. The fitting of the linear regression with the data set confirmed the validity of the linear reservoir law assumption. Results gave the following values: k1=1.38 10⫺4 s⫺1, k2=0.96 10⫺4 s⫺1 and k3=1.22 10⫺4 s⫺1.

Then, the model was used to simulate the discharge from the three terraces described above. As it was assumed for the mathematical development that the terraces were identical, the following parameters were taken: A=78.4 m2, d=35 mm and k=1.38 10⫺4 s⫺1, corresponding to the values of the first terrace. According to the conditions of the experiment, the discharge Q1 began to dry up at t=5760 s and Q2 at t=19,860 s. Then, it was predicted that the equation of continuity would not be verified during the total period of monitoring, due to this finite water supply. Nevertheless, Eq. (6) was used to calculate the discharge from the third terrace. As shown in Fig. 2, the simulated discharge (Q3) over-estimated the observed data during the second part of the experiment. The total volume of water lost by each terrace was also calculated. The first terrace had lost approximately 5.5 m3 (v1), the second terrace 13.0 m3 (v2) and the third one 19.5 m3 (v3). Obviously, these results are cumulative and the individual contributions of the terraces were, respectively, v1=5.5 m3, v2=7.5 m3, v3=6.5 m3. Nevertheless, these results were largely inferior to the theoretical values given by the product of the area by the initial water level: v1⬘=7.0 m3, v2⬘=10.0 m3 and v3⬘=8.5 m3. Though some water was obviously trapped in each terrace, even after the outflow stopped, it is unclear whether the kinetic of drainage was mostly influenced by surface retention and channeling processes or modification of the flow through the PVC pipes. 4.2. Validation of the hydrological model Several Horton’s parameters were extracted from the digitised map of the drainage network. The Confluence Ratio (Rc) and the Stream Mean Length Ratio (RL) gave nearly constant values between both basins as shown in Table 1. The fractal dimension of the network (D) was computed according to the expression proposed by Tarboton et al. (1992): D⫽ln(Rc)/ln(RL)

(9)

The close values of the fractal dimension, 1.43 for Kali Garang and 1.67 for Banuymanik confirmed the fractal feature of the hydrological network, permitting us

I. Gatot et al. / Environmental Modelling & Software 16 (2001) 31–36

35

Fig. 2. Evolution with time of the output discharge from the cascaded terraces. Simulation with Eq. (6) and k=1.38 10⫺4 s⫺1. Observed (dots) and simulated (lines) values.

to calculate the Probability Density Function (r(L)) of the Banyumanik basin, according to Eq. (1). Ten flooding events were extracted from the 1998 database. The selection was based upon two criteria: the peak discharge (Qp) should be higher than 11 m3/s and the rainfall pattern should be homogeneous. The latter was determined applying the Nash and Sutcliffe’s test to the Thyessen average rainfall computed for each event. We present here the results for one of the most homogeneous rainfalls, recorded on 27/03/98. The simulated discharge, coming from the non-terraced area was computed, after optimisation of Qp, with kr=8.13% (Fig. 3). The time of response of the basin was slightly underestimated. Concerning the recession limb of the curve, the model was obviously unable to simulate the observed slow recession, assumed to be generated by the terraced area. As a matter of fact, adding the contribution of the terraced area significantly improved the model fitting (Fig. 3). But the cumulative response also increased the simulated peak discharge, as the opening of all the terraces’ outlets occurred when i=5 mm/h. Nevertheless, the residual variation was still significant, but at this stage it is uncertain whether this was due to the hydraulic model assumptions or whether to more general features (rainfall spatial variability, time dependent flow velocity).

5. Conclusion The Cascaded Reservoirs theory proposed by Nash (1957), and generalized by Dooge (1959), stating that the output discharge is proportional to the water storage,

has been considered for a long time as a purely theoretical concept in hydrology. This study shows that a natural environment, modified by man, can provide linear reservoir-like entities. As a matter of fact, several aspects need further improvement. First, the continuous water supply assumption is not entirely consistent with the local irrigation device. A complementary analysis by Gatot Sumarjo (1999) has improved the results by taking into account the drying time of each experimental terrace for discretizing the global discharge (Fig. 4). As three cascaded terraces poorly represent the complexity of the actual irrigated schemes, an on-going study is aiming at optimizing the set of terrace parameters (n, A, k) according to field information. Another improvement should come from the use of a normal distribution law ruling the opening of the terrace’s outlets. Within the Banyumanik basin (79 km2), landscape features sharply differ from the upstream slopes to the downstream plateau (Perez et al., 1998). As a matter of fact, the fractal structure of the drainage network is slightly degraded as the size of the basin increases. The comparison between the theoretical and experimental (r(L))-curves show that the first one gave a peak of density located 8 km to the basin outlet, against 13 km for the latter. This difference reflected indeed the spatial organisation of the drainage network, from the branching streams running down the upper slopes to the downstream canyons. As a consequence, it was predicted that the model would anticipate the arrival of the flooding at the Banyumanik station. The same study was completed on another Kali Garang sub-basin, called Kripik basin (15 km2). As the size decreased, the spatial variability of the rainfall and of

36

I. Gatot et al. / Environmental Modelling & Software 16 (2001) 31–36

Fig. 3. Observed discharge (white dots) compared with the simulated discharge coming from the non-terraced (bold line) and terraced (black line) areas. Banyumanic basin on 27–28 March 1998 (rainfall: 48.5 mm; kr: 8.13%).

Fig. 4. Evolution with time of the output discharge from the third terrace. Simulation with finite water supply. Observed (black dots) and simulated (dashed line) data.

the soil characteristics had less influence on the model’s parameters. The theoretical Probability Density Function was also very close to the observed density of the hydrological network. As a matter of fact, the results of simulation better fitted the observed hydrographs (Gatot Sumarjo, 1999).

References Casenave, A., Valentin, C., 1989. Les e´tats de surface de la zone sahelienne. Influence sur l’infiltration. Orstom Ed, Paris, France, 230 pp. Chow, V.T., 1964. Runoff. In: Chow, V.T. (Ed.), Handbook of Applied Hydrology: A Compendium of Water Resources Technology. McGraw-Hill, New York, 1418 pp. Dooge, J.C.I., 1959. A general theory of the unit hydrograph. J. Geophys. Res. 64 (1), 241–256. Duchesne, J., Cudennec, C., Corbierre, V., 1997. Relevance of the

H2U model to predict the discharge of a catchment. Water Science and Technology 36 (5), 169–175. Gatot Sumarjo, I. 1999. Modelisation de la transformation pluie-debit et de l’influence de l’amenagement de terrasses sur les crues de mousson. Doctoral thesis, ENSA, Rennes, France, 214 pp. Gatot Sumarjo, I., Duchesne, J., Perez, P., 1997. H2U, a transfer function model using fractal characteristics of the hydrographic network. In: McDonald, A.D., McAleer, M. (Eds.), Proceedings of MODSIM 97, Hobart, 8-11 December, vol. 1. MSSA, Canberra, Australia, pp. 470–478. Gupta, V.K., Waymire, E., Wang, C.T., 1980. A representation of an instantaneous unit hydrograph from geomorphology. Water Resources Research 16 (5), 855–862. Huggins, L.F., Burney, J.R., 1982. Surface runoff, storage and routing. In: Haan, C.T., Johnson, H.P., Brakensiek, D.L. (Eds.), Hydrologic Modeling of Small Watersheds, ASAE monograph, 5. ASAE, St Joseph, MO, pp. 169–225. Nash, J.E. 1957. The form of the instantaneous unit hydrograph. C.R. and Reports, Assoc. Inter. Hydrol., IUGG, Toronto, Canada, pp. 114-121. Perez, P. 1994. Gene`se du ruissellement sur les sols cultive´s du sud Saloum (Se´ne´gal). Doctoral Thesis, ENSA, Montpellier, France, 251 pp. Perez, P., Gatot Sumarjo, I., and Prasetyo, T. 1998. Kali Garang pilot project, research highlights. CIRAD and CSAR, Bogor, Indonesia, 21 pp. Rodriguez Iturbe, I., Valdes, J.B., 1979. The geomorphologic structure of hydrologic response. Water Resources Research 15 (6), 1409– 1420. Rodriguez Iturbe, I., Gonzales Sanabria, M., Marani, A., 1982. A geomorphoclimatic theory of the instantaneous unit hydrograph. Water Resources Research 18 (4), 877–886. Tarboton, D.G., Bras, R.L., Rodriguez Iturbe, I., 1992. The fractal nature of rivers networks. Water Resources Research 34 (8), 1317– 1322. Wu, R.S., Sue, W.R., Chang, J.S., 1997. A simulation model for investigating the effects of rice paddy fields on runoff system. In: McDonald, A.D., McAleer, M. (Eds.), Proceedings of MODSIM 97, Hobart, 8-11 December, vol. 1. MSSA, Canberra, Australia, pp. 422–427.