Forecasting the discharge into Hoa Binh reservoir by applying the connecting model marine - Imech1d

Vietnam Journal of Mechanics, VAST, Vol. 30, No. 3 (2008), pp. 149 – 157 FORECASTING THE DISCHARGE INTO HOA BINH RESERVOIR BY APPLYING THE CONNECTING MODEL MARINE - IMECH1D Nguyen Tien Cuong Institute of Mechanics, VAST, Vietnam Trinh Thu Phuong National center for Hydro Meteorological Forecasting Abstract. Da river is the biggest one of three rivers: Da , Thao and Lo which flow into Hong river, so calculating and forecasting correctly the discharge to Hoa Binh reservoir take an important role in managing the task of flood preventing in the Hong - Thai Binh river system. In this paper, a connected hydraulic-hydrological model is developed. The model was applied to forecast the discharge into Hoa Binh reservoir in the flood season 2006. The difference of this connected hydraulic-hydrological model in comparison with the other hydrological models is using the complete Saint-Venant 1D to simulate flow in rivers. Results of verification of the model show that it is good enough for practice. 1. THE COMBINED MODEL OF IMECH-1D AND MARINE The simulating hydraulic one dimension IMECH-1D programme pack was developed at Institute of Mechanics, Vietnam Academy of Science and Technology. It was close verified and applied to simulate the flood forecast for Hong-Thai Binh river system. The set of simulating hydrological programme MARINE was developed at Toulouse Institute of Fluid Mechanics, France. It was handed over to Institute of Mechanics within the framework of “FLOod Control Decision Support” (FLOCODS) project financed by European Community. The combined model of MARINE and IMECH-1D has a high ability to simulate dis- charge process of flow correctively. The Saint-Venant 1D equations [2] used in IMECH-1D are: ∂Ac ∂t + ∂Q ∂x = q, (1) ∂Q ∂t + ∂ ∂x ( β Q2 A ) + gA (∂Z ∂x + Sf ) = 0, (2) where: Q = Q(x, t) is discharge of the flow in the river section; Z = Z(x, t) is the water level in the river section; Sf is the friction slope; Ac is the cross-sectional area; q is lateral flow; β is the correcting coefficient; A is the flow area. Lateral flow q is taken from the result of MARINE. Friction slope in IMECH-1D is calculated as: Sf = n2Q |Q| A2R4/3 , (3) 150 Nguyen Tien Cuong and Trinh Thu Phuong where: R is hydraulic radius. dV dt = P + ∑ k Qk. (4) Mathematical model of a cell: V = V (Z) = V (Z(t)): volume of the cell is a function of water level Z; P : rainfall or vaporization at that cell; Qk: the water exchanged between the cell and related cells. The Preissmann sketch difference [2] is used to system of equations (1), (2), (4). After difference of equation (1), we get: 1 2∆t [ (Ac, i+1 +Ac,i)− ( ATc,i+1 + A T ci )] + 1 ∆X [ θ (Qi+1 −Qi) + (1− θ) ( QTi+1 −QTi )] , = 1 2 [ θ (Qi+1 + qi) + (1− θ) ( qTi+1 − qTi )] . (5) We write the equation (2) as following: a1 + a2 + a3 + a4 = 0, where: a1 = ∂Q ∂t ; a2 = ∂ ∂x ( β Q2 A ) ; a3 = gA ∂Z ∂x ; a4 = gASf . Using the Pressmann sketch difference we have: a1 ≈ a˜1 ≡ 1 2∆t [ (Qi+1 +Q1)− ( QTi+1 +Q T i )] , a2 ≈ a˜2 ≡ β ∆x [ θ (Q2i+1 Ai+1 − Q 2 i Ai ) + ( 1− θ )(QT 2i+1 ATi+1 − Q T 2 i ATi )] , a3 ≈ a˜3 ≡ g 2 [ θ(Ai+1 + Ai) + (1− θ)(ATi+1 + ATi ) ] 1 ∆x [ θ(Zi+1 − Zi) + (1− θ)(ZTi+1 − ZTi ) ] , a4 ≈ a˜4 ≡ g θ 2 [ (ASf )i+1 + (ASf )i ] + g (1− θ) 2 [ (ASf ) T i+1 + (ASf) T i ] . So a˜1 + a˜2 + a˜3 + a˜4 = 0. (6) From the equation (4) we get: 1 ∆t ( Vi − V Ti ) = θ ( Pi + ∑ k Qi,k ) + (1− θ) ( PTi + ∑ k QTi,k ) . (7) The linear system of equations (1) - (7) is solved by Gauss method. MARINE model simulates the process of shaping a flow by raining in the basin based on the equation of conservation of mass [5]: ∂V ∂t + u.grad(V ) = P0, (8) where: V is the volume of the liquid; U is the velocity of the flow between the cells; P0 is rainfall. Because: u.grad(V ) = div(V.u)− V.div(u). Forecasting the discharge into Hoa Binh ... 151 And for incompressible fluid, we have div(u) = 0 , using Green-Ostrogradski formula:∫∫ S div(m.u).dS = ∮ Γ m.u.n.dΓ from (8) we have: ∫∫ S ∂V ∂t .dS + ∮ Γ V.u.n.dΓ = ∫∫ S P0. (9) The velocity of the flow exchanging between the cells is calculated by [1], [2]: ‖u‖ = √pente.H 2/3 Km . (10) Because the cells are square (DEM) so putting the velocity expression (10) into the integral equation (9), we get the equation for calculating the varying of water depth: ∆H + 4∑ j=1 H 5/3 j Km . √ pente. ∆t ∆x = P0.∆t, (11) where: Pente is the slope; Km is Manning roughness coefficient; ∆x is the weight of the cell; ∆t is time increment; j is flow direction of the cell (j = 1÷ 4); H is the water depth of the cell; ∆H is the water depth changing from t to t+∆t. From the changing of water depth in each cell, we can calculate the total discharge of each cell (consisting of the discharge from rainfall, inflow and outflow) at each step equal to the varying of the volume in the cell: ∆Q = ∆H∗dx∗dx, where: dx is the cell size. For a close basin, there is only one outlet and the outflow of the basin is: q = ∆Q. For an open basin, the basin is along two banks so there are a lot of outlets. In this situation, the outflow is a sum of the exchanging discharge of the outlets: q = ∑ Q = ∑ ∆H∗dx∗dx. So, the result of MARINE is the discharge of the basins. This is the component added into lateral flow q for hydraulic model IMECH-1D. The infiltration in the MARINE can be calculated based on Green Ampt infiltration theory from continuous equation and Darcy law. Potentiality accumulating infiltration is calculated by Green - Ampt equation [3]: F (t) − ψ∆θ ln ( 1 + F (t) ψ∆θ ) = kt. (12) In which: F (t): stands for the depth of water which is infiltrated in land; ψ: capillary water column of wet surface; ∆θ = η − θi with η stands for land deficient, θi stands for land moisture; k: hydraulic conductivity. Equation (12) is a nonlinear one which can be solved by a consecutive replacement method or Newton iterative method. In case the depth of remain water ho is not ignorable it is required to replace ψ by ψ - ho before solving the problem. After finding out the 152 Nguyen Tien Cuong and Trinh Thu Phuong potential depth of accumulated infiltration F (t), velocity of potential infiltration can be defined as: f = k( ψ∆θ F + 1). (13) According to Darcy law, infiltrated output of each grid can be calculated as follow: Qinfiltrated = f∗dx∗dx Before producing remain water (t < tp), rainfall intensity is smaller than the velocity of potential infiltration and land is unsaturated. Process of water remaining starts to occur when rain intensity is higher than the intensity of potential infiltration (t = tp), and at that time, land is in saturation state. When the saturation keeps its process (t > tp), saturated areas will spread gradually into the deeper layer of land and a stream flow on the land surface running from the remaining water. Infiltration process is divided into 3 periods: - Before appearing the remaining water, all rain-water are infiltrated into land (t < tp). The velocity of potential infiltration is represented by a function of accumulated infiltration level F (t = tp). - Remain water is available when velocity of potential infiltration is smaller or equal to rainfall intensity i (t > tp). -When t = tp, the volume of accumulated infiltration at the time o remain water generation tp is solved by formula: Fp = i∗tp so f = i. Replaced in (13) making: i = k ( ψ∆θ it + 1 ) and resulting in tp = k ψ∆θ i(i− k) . (14) After, remain water is formed, the accumulated infiltration volume is calculated by the following equation: F − Fp − ψ∆ ln( ψ∆θ+ F ψ∆θ + Fp ) = k(t− tp). (15) Surface flowat grid cell P0 (Rain) Infiltration - Green Ampt Remain water layersSurface flow at grid cellConnection of grids, calculation of Q exchanged between grid cells. Process of Q~t at the exit of the valley or registered button Structure of the valley (topographic height, land, cover) Fig. 1. Diagram of MARINE model So, at each time step of infiltration module, we can get infiltrated output qinfiltrated of each containing tank. In this case of calculation of MARINE by using infiltration module Forecasting the discharge into Hoa Binh ... 153 at each time step, the total exchanging output of each compartment includes rainfall, inflow, outflow, and infiltrated water: ∆Q = ∆H∗dx∗dx− qinfiltrated. Fig. 1 shows the structure of MARINE 2. APPLICATION THE COMBINED MODEL OF MARINE AND IMECH-1D TO FORECAST FLOW INTO THE HOA BINH RESERVOIR 2.1. Data processing a. Data of basin terrain: The Da river basin on Vietnamese territory starting from the Vietnam and China border to the Hoa Binh dam in Fig. 2. It has total area of about 24.000 km2 with 36 rainfall stations and some other water and output measurement stations which are mapped in MARINE by 10 small basins connected with each other through river system (Fig. 2). The terrains used to calculate are built in a topographic map with scale of 1:50000. Land use map and flora cover map are processed and calculated according to a scale of 1:100000. Data of all basins have been examined by qualified verification problems of the model. Rainfall data taken into account is result of real measurement taken by 36 rainfall stations referring to forecasted data of numerical value model (in some cases). Rainfall distribution on the basin is based on Thiessen principle [3]. Fig. 2. Model of the Da River basin in MARINE b. Topographic data of river section starting from Vietnam-China border to the Hoa Binh dam: The whole main branches of the Da river limited from Vietnam-China bolder to the Hoa Binh dam include 121 cross-sections. In Fig. 3, one-dimensional hydraulic model IMECH-1D calculates the flow in this branch. Total lateral flow added in the middle areas of IMECH-1D is calculated by MARINE hydrological model in the whole basin and flow to river. 154 Nguyen Tien Cuong and Trinh Thu Phuong HBdamQ Upper boudary q lateral flow Cross- Sections Fig. 3. Model of the Da river in IMECH-1D 2.2. Results of pilot forecasting of flow poured into the Hoa Binh reservoir in the flood season 2006 The model MARINE-IMECH1D which is in pilot operation in National Centre for Hydro meteorological Forecasting forecasts that the flow pouring in Hoa Binh reservoir in the expected time of 48 hours to sever for the electricity generation and lowland flood prevention in the 2006 flood season. a. General assessment on flood situation of the Da river in 2006: In 2006, there were about 9 flood spells in Da river with Qmax in reservoir ≥ 3500 m3/s. From 1h/28/06/2006 to 19h/28/07/2006, there were 3 ones in succession: Qmax 1 = 6000 m 3/s at 17h/02/07/2006 Qmax 2 = 9000 m 3/s at 23h/10/07/2006 Qmax 3 = 14700 m 3/s at 12h/09/07/2006 This is the biggest output poured in Hoa Binh reservoir in this year and ranks the seventh in survey chain from 1902 until now (Qmax 1 = 22500 m 3/s in 1996; Qmax 2 = 17800 m3/s in 1945; Qmax 3 = 17200 m 3/s in 1964; Qmax 4 = 16200 m 3/s in 1971; Qmax 5 = 15800 m3/s in 1969; Qmax 6 = 15200 m 3/s in 2002). Flood season lasting from 13h/06/10 to 7h/22/10 with Qmax in Hoa Binh reservoir = 9200 m3/s 13h/12/X ranks the second in the data chain of greatest flow to Hoa Binh reservoir in October, 1902 until now (Qmax of October max = 10900 m 3/s in 1932; Qmax of October ranking the third = 9000 m3/s in 1999). b. Results forecasted by the combined model of MARINE and IMECH1D: Forecasting results of 9 flood spells in flood season 2006 are shown in Table 1. Table 1. Forecasting results for peak of a flood to Hoa Binh reservoir in flood season 2006 Qmax real Qmax forecast measurement into Time of calculation into ∆Q % error ∆τ(h) Hoa Binh appearance Hoa Binh (m3/s) (T forecast- reservoir (m3/s) reservoir (m3/s) Reality) 5500 10h/1/VII/2006 5000 -500 -9.1 14 6000 17h/2/VII/2006 5300 -700 -11.7 12 9000 23h/10/VII/2006 8900 -100 -1.1 2 14700 12h/19/VII/2006 13808 -892 -6.1 1 Forecasting the discharge into Hoa Binh ... 155 3600 23h/8/VIII/2006 3672 72 2 - 4 3800 17h/17/VIII/2006 5457 1657 43.6 17 6000 5h/19/VIII/2006 4900 -1100 -18.3 6 4000 1h/25/VIII/2006 3238 -762 -19.1 0 9200 13h/12/X/2006 8900 -300 -3.3 0 The input daily data is the combination of real rain measured by 36 rainfall stations on the basin, the discharge and water level taken by some stations, and numerical value forecasted rainfall data. It takes 30 minutes to operate the connection program MARINE- IMECH1D on computer Intel Pentium 4 3.0GHz. The results are: - Inflow poured into reservoir is forecasted in 48 hours. - Water level and discharge in some places on Da river. Hereafter are results forecasted for flood season 2006 through the connection program MARINE-IMECH1D: In Table 2, forecasting results of flood process to Hoa Binh reservoir are shown Table 2. Forecasting results of flood process to Hoa Binh reservoir (From 18/VI to 23/X/2006) Estimated time Error level Number of correct P% point/Total 6 h 18.5 %Q 97/116 82.5 12h 18.5 %Q 103/116 83.6 18h 18.5 %Q 94/116 81 24h 18.5 %Q 89/116 76.7 30h 21.5%Q 88/115 76.5 36h 24.5%Q 81/115 70.4 42h 27.5%Q 83/115 72.1 48h 30.5%Q 83/115 72.1 Forecasting result for 48 hours (in future) of discharge into Hoa Binh reservoir from 5th to 30th July 2006 in Fig. 4. Fig. 5 shows the example of forecasting process line of newscast on 18th July 2006. 156 Nguyen Tien Cuong and Trinh Thu Phuong Realmeasurement Forecast Fig. 4. Forecasting result of discharge into the Hoa Binh reservoir from 5th to 30th July 2006 HB Real measurement HB Calculation Fig. 5. Line of forecasting process of newscast on 18th July 2006 3. CONCLUSION Through the process of operating the combined model MARINE-IMECH1D for flood season 2006, it shows that: The model is rather sensitive to raining phenomenon creating flood. It will get more correct result (typically flood from 16 –22/VIII) if rainfall of the whole basin is high and steady. Result of the program will be higher than observation data (typically flood from 16 –22/VIII) if distribution of rain is unsteady. Forecasting the discharge into Hoa Binh ... 157 Forecasting result of rising flood is very good, but when flood neaps, the result is lower than observation data. The result of combined model in 2006 flood season is satisfactory. In order to get better result of combined model MARINE-IMECH1D, it is necessary to increase number of rainfall stations to approach real distribution of rain on basin, and improve quality of observation rainfall and forecast rainfall in following flood seasons. ACKNOWLEDGMENT Authors would like to thank the Board of chairman of the international project FLOCODS, the Board of chairman of the state projectKC.08-13 “Applying technology of Da river flood forecast to control Hoa Binh reservoir in flood prevention in lowlands in period of 2004-2006” for creating favorable condition to develop the model. The project is completed thanks to expense support of basic research programme. REFERENCES 1. V. Estupina-Borrell, D. Dartus, R. Ababou, Flash floodmodeling with the MARINE hydrolog- ical distributed model, Hydrology and Earth System Sciences Discussions 3 (2006)3397-3438. 2. Ven Te Chow, David R. Maidment, Larry W. Mays, Applied Hydrology, Mc Graw Ọ Hill Book Company, 1998. 3. David R. Maidment, Handbook of Hydrology, Mc Graw Ọ Hill Book company. 2000. 4. Roger J. M. De Wiest, Geohydrology, John Wiley & Son, Inc., New York, London, Sydney, 1999. 5. M. Alquier, J. Chorda, D. Dartus, V. Estupina Borrell, C. Llovel, and M. M. Maubourguet, PACTES: La chaine de prevision du Thore, Research Contract, IMFT, Toulouse, 2002. Received January 10, 2007 ỨNG DỤNG MÔ HÌNH KẾT NỐI MARINE-IMECH1D ĐỂ DỰ BÁO LƯU LƯỢNG VÀO HỒ HÒA BÌNH Sông Đà là sông lớn nhất trong ba sông Đà, Thao, Lô chảy vào sông Hồng nên việc tính toán, dự báo tốt lưu lượng vào hồ Hòa Bình có ý nghĩa rất quan trọng trong điều hành công tác phòng chống lũ trên hệ thống sông Hồng-Thái Bình. Trong bài báo này, mô hình kết nối thủy lực thủy văn được phát triển. Trong mùa lũ 2006 mô hình được ứng dụng để dự báo lưu lượng dòng chảy vào hồ Hòa Bình. Khác với các mô hình thủy văn trong mô hình kết nối thủy văn thủy lực sử dụng phương trình Saint-Venant 1 chiều đầy đủ để tính dòng chảy trong sông. Kết quả tính thử nghiệm chứng tỏ mô hình kết nối này là rất tốt.