Gis Project

Large Scale River Network Routing Using RAPID Ahmad Tavakoly

Introduction Rivers play an important role in the Earth’s hydrological cycle (Figure 1), and most climate system models nowadays include continental scale river transport models (RTMs) to complete the global water balance. In the continental scale Water on land regulates heating and moistening of the atmosphere and directly has effect on Earth`s climate. Here in the US, several databases for river network are accessible and we can easily get data for all large river basins.

Figure 1: Hydrological cycle as represented in most earth system models (Ref. Larson et al, 2007).

One of the best sites is NHD1 Plus. This data base first released in 2006 and it has some water components data including: catchment characteristic, flow volume, velocity, and etc. in this web site the nation divided to 21 hydrology regions and just by click on the map for desired region all data for that region are available. Figure (2) shows region 12 based on this data set.

1 National hydrography dataset

Figure 2: Texas gulf (hydrologic region 12, ref. NHD Plus website)

Land surface model (LSMs) developed by atmosphere scientists for atmospheric model with boundary conditions (water and energy balance). River routing models in the LSMs have been classically used gridded river networks that best fit the computational domain used in LSMs (David et al. 2009). Today geographic information system (GIS) database can be used. Cedric David developed a river routing model for large-scale applications called: RAPID (Routing Application for Parallel computation of Discharge) to calculate discharge using Muskingum routing scheme. This model can be used for any river network using input data at the upstream and downstream of each river reach. My dissertation is about parallelization of St. Venant equation for a large river network. Using the GIS one can show the movement of water volume along a river direction. The purpose of my GIS project is to learn how the RAPID model works and try to run it for Texas Gulf (hydrologic region 12) or a new network. As we can see in the figure 1, the region 12 is also divided into production units: 12 a, b, c, d, e, and f. in addition, I intend to animate a map with time based on RAPID flow calculations.

Rapid Model To analyze flow variables, such as the flow depth and velocity or the flow depth and the rate of discharge, the continuity equation and the momentum or energy equation are used. However in some cases some terms of the governing equations are smaller than the other terms. These analysis procedures called approximate methods. For instance, if continuity equation solved simultaneously with a simplified form of momentum equation is called hydrologic routing. Rapid model is developed for NHDPlus river network and it can be use for parallel computing. This model uses the Muskingum method for channel routing. The backbone of RAPID is a vector-matrix version of the Muskingum method (David et al. 2009): ∴∗

( I − C1.N ) .Q (t + ∆ t )= C1 .Q (t +) e

C2 N .Q t(+ ) Q t( +) C3 .Q t ( e

)

ΜΕ Ρ Γ Ε Φ Ο Ρ Μ Α Τ 0 Where: t is time and Δt is the river routing time step. I is the identify matrix. N is the river network matrix. C1, C2 and C3 are parameter diagonal matrices and for a given reach j, they can represent as [Birkhad and James, 2002]: \*

C1 j = C2 j = C3 j =

K j X j − 0.5∆t

K j ( 1 − X j ) + 0.5∆t K j X j + 0.5∆t

K j ( 1 − X j ) + 0.5∆t

K j ( 1 − X j ) − 0.5∆t

K j ( 1 − X j ) + 0.5∆t

MERGEFORMAT0

Kj is storage constant (with dimension of a time) and Xj is a dimensionless weighting factor. Q is a vector of outflows from each reach, and Qe is a vector of lateral inflows for each reach. Lateral inflow in this model is calculated by land surface model.

Parameter estimation

Parameters K and X need to be determine in the RAPID model. For this purpose, an inverse method is used. This method optimizes the parameters so that out puts match with the observed data. Fread, (1993) said that X is between 0.1 and 0.3 in most streams. The Muskingum K parameter is estimated from Eq. (3) as follows (Fread, 1993, Tewolde and Smithers, 2006): \*

Kj =

Lj Vwj

MERGEFORMAT0

Where: Kj = storage constant, Lj = reach length, Vwj = wave celerity. The USGS IDA (http://ida.water.usgs.gov/ida/) provides 15- minutes flow measurements that can be used to determine wave celerity. For the San Antonio and Guadalupe basins data at fourteen gaging stations are obtained from IDA. Cross correlation values for each pair of stations were calculated. The maximum correlation and corresponding time lag for any of the two stations are calculated used to determine the flow wave celerity. More information about how to determine wave celerity is online available at: (https://webspace.utexas.edu/aat669/report.mht?uniq=5bqctw&xythosdownload). The celerity used here in the model simulation for Region 12 is uniformly 2.5 m/s. To avoid estimating (Kj and Xj) for all reaches, RAPID model optimizes two multiplying factors and using following equations: λk λx ,

K j = λk .

Lj

\* MERGEFORMAT0

x j = λx ⋅ 0.1

Vwj

At the end of the optimize procedure, a couple of (

,

) is determined for a

λk λx given basin in the network.

Rapid Model Applications In this study first we try to run RAPID model for Guadalupe and San Antonio River Basins (Figure 3). NHDPlus database provides the mapped streams and rivers as well as the catchments that surround them in the United States. On

this database river reach in the national network is assigned a unique integer identifier called COMID. NHDPlus catchments also have a COMID, the same COMID being used for the reach and its local contributing catchment. From the attribute table of NHD flowlines feature, the Guadalupe and San Antonio river basins have a total 5175 reaches. These reaches have an average length of 3.00 Km and average catchment size area around them is 5.11 Km2. River network and sample reach with its catchment shown in the figure 4.

San Antonio and Guadalupe Basins

µ Legend Texas Basin

0

55 110

220

330

440 Kilometers

Figure 3: Guadalupe and San Antonio basins

To run RAPID model, runoff is required and it is calculated by a land surface model called Community Noah Land Surface Model with Multi-Physics Option Noah-MP [Niu, et al., 2009]. The runoff is provided on a grid (several netCDF files, thousands of them, one for each time step of the study). This grid needs to be converted to something readable by RAPID. David do that using a "flux coupler". That's basically a program that accumulates the runoff from the land surface and dumps it to the corresponding river reach.

µ Legend NHDPlus reaches NHDPlus catchments

0 10 20

40

60

80 Kilometers

Figu re 4: NHDPlus river network and catchment for the Guadalupe and San Antonio Basins

The model is used to calculate river flow in all 5175 river reaches of the Guadalupe and San Antonio River Basins for four years (01 January 2004 – 31 December 2007). Output of the model include: volume and flow. Rapid calculates variables every three hours. The first set of values that are in the output file correspond to an average of the first three hours (between 00:00 and 03:00). The second set of values corresponds to between 03:00 and 06:00. To take daily average of model output, David developed program called “process_model_flow_map_for_Arc.f90”, in this program IS_M is a number of days. In this code, IS_M is multiplied by 8. 8 is the number of time steps per day (8 x 3 = 24 hours per day). Therefore, if we put IS_M=105, what the program (process_model_flow_map_for_Arc.f90) gives us is the average outflows of the 105th day (from 00:00 hours to 23:59 hours). This model also used for the hydrology region 12. This region includes 74615 reaches (figure 5).

Base map of Region12

µ 0

55

110

220

330

440 Kilometers

Legend USGS Gages Subbasin

Figure 5: hydrology Region 12

Results Each river named by COMID in the Rapid model and we can get the time series for every COMID which represents upstream river reach. From the model output, computed flow rate are compared with the daily data. Four stations in the San Antonio and Guadalupe river basins are selected. Figures 6 shows selected gages in these basins. As we can see two of these stations are on the main river San Antonio and Guadalupe Rivers and the others have smaller upstream basins.

SanGuad Basins with selected gages to compare reults

Guadalupe Rv nr Spring Branch, TX

b

µ

b

Legend

b

SanGuadBasins

Guadalupe Rv at Victoria, TX

b

San Antonio Rv at Goliad, TX

b

selectedgages SanGuadFlowline

San Antonio Rv nr Falls City, TX

0 12.5 25

50

75

100 Kilometers

Figure 6: location of selected stations for comparison of result.

Once the run is done, model results are compared with the observed data. Figures 7-10 show hydrographs of flow routing.

Figures 7: Hydrograph of observed and computed flow for the Guadalupe River at Victoria

Figures 8: Hydrograph of observed and computed flow for the Guadalupe Rv nr Spring Branch, TX

Figures 9: Hydrograph of observed and computed flow for the San Antonio Rv at Goliad, TX

Figures 10: Hydrograph of observed and computed flow for the San Antonio Rv nr Falls City, TX

Based on the above figures, it seems that results for gages on the main rivers going dawn along a stream have more accuracy comparing with gages on the small rivers. It is very difficult to comment on why downstream would be better than upstream. One of the potential reasons is that the land surface model that was used to create the input data for RAPID is calibrated using downstream gages. Other than that, it all due to the land surface model calculations: runoff scheme, land cover, vegetation scheme. Overall those results are very satisfactory. From the RAPID output we can also see a variation of flow over basins. Figures 11 and 12 show flow out of all reaches for San Antonio Guadalupe and Region 12 Basins respectively. Both figures represent flow routing over basins on 20 April 2004.

Flow caluculated by Rapid model for SanGuad Basins

µ

Legend

NHDFlowline_San_Guad_with_dir Qout

0.000000 - 0.448000 0.448001 - 1.510000 1.510001 - 3.070000 3.070001 - 4.950000 4.950001 - 7.940000 7.940001 - 14.500000 14.500001 - 25.700000 25.700001 - 45.300000 45.300001 - 68.900000 68.900001 - 168.000000 streamgage_San_Guad_gotQ_spa_join_2004_2007_full

0

10 20

40

60

80 Kilometers

Figures 11: Flow rate calculated for all reaches by model for San Antonio and Guadalupe Basins

Figures 11: Flow rate calculated for all reaches by model for San Antonio and Guadalupe Basins

Scalability of parallel computation In this project the model is set up for San Antonio and Guadalupe river network with 5,175 rivers and Region 12 with 74,615 reaches. The river network of the NHDPlus data set has about 3 million reaches for the United States. To simultaneously compute flow and volume of river in all reaches for such a large problem, parallel computing is needed. RAPID can be applied on several processors. In order to assess scalability of the river network model, the model is run on the Lonestar supercomputer (http://www.tacc.utexas.edu/resources/hpcsystems/#lonestar) for Region 12. Figure 12 shows performance of the model using different number of processors. In this figure two computing times are given: the wall-clock time is the time

between the start and the end time of the all computations for the last processor. The red plot in Fig. 12 is the simulation time difference, when adding one more processor. For example the red plot at number of processor=1 represents the simulation time difference when using two processor and one processor. Fig. 12 shows that simulation time decreases significantly with increasing number of processors up to 5. An increase beyond 5 processors results in small decrease in simulation time.

Figures 12: Scalability of RAPID computations

Conclusion Table 1 shows the root mean square error (RMSE) of computed flow rate for the selected gages. As this table shows the model results are satisfactory. Table 1: RMSE for the selected gaging stations Guadalupe River at Victoria

RMSE

75.08

Guadalup e Rv nr Spring Branch, TX 45. 49

San Antonio Rv nr Falls City, TX

San Antonio Rv at Goliad, TX

36.78

44.05

In this project, the scalability of RAPID model for parallel computing is also investigated for region 12. Fig. 12 shows scalability for the simulation considered.

Acknowledgement The RAPID model developed by Cedric David (2009). I would like to acknowledge him for sharing source code and his help to run model.

References: A.L. Birkhead, C.S. James, Muskingum river routing with dynamic bank storage, Journal of Hydrology 264 (2002) 113–132. David, C. H, D. R. Maidment, G.-Y. Niu, Z.-L. Yang and F. Habets., “River network routing in the Guadalupe and San Antonio River Basins”, submitted to Water Resources Research on Sept 25 2009. Fread, D. L. (1993), Flow Routing, in Handbook of Hydrology, edited by D. R. Maidment, pp. 10.17-10.18, McGraw-Hill, New York. Larson., J. W., A. P. Craig, J. B. Drake, D. J. Erickson, M. L. Branaetter, M. W. Ham, “A Massively Parallel Dynamical Core for Continental- to Global-Scale River Transport”, 2007. Niu, G. Y., et al. (2009), The Community Noah Land Surface Model with MultiPhysics Options, Journal of Geophysical Research-Atmospheres, (submitted). Tewolde and JC Smithers, Flood routing in ungauged catchments using Muskingum methods, Water SA Vol.32 (3) 2006: pp.379-388.