Abstract
We propose a simple computer based phenomenological model to simulate the characteristics of river flow velocity in large scale. We use shuttle radar tomography mission based digital elevation model in grid form to define the terrain of catchment area. The model relies on mass-momentum conservation law and modified equation of motion of falling body in inclined plane. We assume inelastic collision occurs at every junction of two river branches to describe the dynamics of merged flow velocity.
1. Introduction
The complexity of river network has been studied intensively over the past decades in order to describe its pattern formation and dynamics [1]. The evolution of river network pattern has been investigated based on various approaches such as continuum [2–4], stochastic [5], and cellular [6] approaches. On the other hand, one of the most important dynamical characteristics of a river is its flow velocity which corresponds to the water discharge along the river lines. Many different empirical models have been proposed to describe the associated dynamics in various scales and for various purposes as well [7–12]. Recently, Gauckler-Manning-Strickler (GMS) formula based empirical model has been implemented to describe the global river flow in certain area of Europe [9].
In the meantime, one of the theoretical descriptions of dynamics of open/close river channel is given by the so-called Saint Venant equations (SVEs) in one or two dimensions. These spatiotemporal equations are developed on the basis of mass-momentum conservation law in the form of coupled continuity and momentum differential equations of average water discharge and cross-sectional river flow channel functions where both need initial and boundary conditions [7, 13–19]. Recently, one-dimensional SVE has been adopted in a coupled hydrodynamic analysis model to investigate the large scale water discharge dynamics at basin of Yangtze River, China, in the presence of complex river-lake interaction [7].
In general, the river flow dynamics is affected by friction, roughness, slope, and number of rocks. Usually, in SVEs all these parameters are incorporated in the GMS formula which is found empirically to approximate the average flow velocity in short distance [17]. However, solving SVEs to simulate the river flow dynamics in a large scale such as the whole area of catchment area is not an easy task. It can be done by solving either the two-dimensional SVE or the one-dimensional SVE incorporated by junction model with internal boundary conditions [7, 19, 20]. It is likely that, for large scale cases, the calculation to solve the related SVEs based on rigorous semianalytical or numerical methods can be very difficult and probably time consuming. One of the difficult challenges to solve these equations numerically for the corresponding cases is how to determine the corresponding initial conditions [7].
In this report, we propose a simple approach to develop a computer based discrete phenomenological model of river flow velocity. We still rely on mass-momentum conservation law by assuming that an inelastic collision occurs at a river junction. Instead of using GMS formula, to accommodate friction, roughness, and number of rocks in single parameter, here we introduce the notion of effective gravitational acceleration, whereas the flow velocity follows the equation of motion of falling body in inclined plane.
The procedure that we followed to develop the model includes the following: (i) determination of river catchment area in the form of a computational window grid, (ii) development of algorithm for flow velocity model, and (iii) examining the model characteristics with respect to the variations of initial water discharges at headwaters.
To our best knowledge, the proposed model has never been reported elsewhere. We compare our model with field observation data in order to justify the results. It should be emphasized from the beginning that this simple model is an attempt to determine characteristics of river flow velocity and it is not intended to replace SVEs. In principle, however, this model can also provide initial condition of the whole river network which is needed to solve numerically the corresponding time-dependent SVEs.
We considered the catchment area of Ciliwung River (CR) situated in West Java Province, Indonesia, for simulating the model. To define the corresponding terrain we used digital elevation model (DEM) in grid form which is constructed based on free shuttle radar topography mission (SRTM) data provided by USGS [20]. Each grid in DEM represents the average height within the corresponding terrain. We only take into account the surface runoff and exclude the subsurface runoff as well as landscape characteristics of surrounding river lines such as forest canopy, vegetation, land use, and sediment transport.
We organize this report as follows: in Section 2 we discuss the associated DEM of CR catchment area and then this is followed by the explanation of flow velocity model in Section 3. The results of our simulation and its discussion are given in Section 4. We conclude our report by a summary in Section 5. Sections 2–4 discuss consecutively the abovementioned procedure that we followed.
2. DEM of Ciliwung River Catchment Area
The CR is one of the most important rivers in Java Island, Indonesia. It spans over 120 km from Mount Pangrango in West Java to Java Sea, north to the capital city Jakarta. Its catchment area is 149 km2. In recent years, during the rainy season, this river contributes around 30% of flood in Jakarta.
To plot the CR catchment area we used DEM with m2 grid resolution and in order to get a more detailed terrain which still can be handled by modest computational facility, the associated DEM was corrected by m2 resolution DEM as well as data from field survey. Using GIS software this modified DEM was converted into numbers that can be read by Excel software. Shown in Figure 1(a) is the plot of m2 DEM of the corresponding CR river network pattern. It is obvious that water will flow to the lower land due to gravitational force, such that we can use this fact as a rule in our computer code to determine the corresponding river pattern based on the given DEM data after identifying the positions of all headwaters. Therefore, it is clear that the boundaries of all river branches which define a two-dimensional computational window are automatically provided by DEM, where each grid represents the corresponding computational mesh. In Figure 1(b) we give an example of specific small area of the corresponding computational window grid (see figure caption) which is highlighted with grey color grid. The number inside each grid represents DEM value in meter unit at specific position.
(a)
(b)
Based on our modified DEM, we identified that there are 126 headwaters. According to Strahler classification of river branches, all these headwaters converge into order 4 river channel at Katulampa area with channel dimension of ~70 m width and 0.3–0.5 m water level. Similar to the most of tropical rivers, catchment upper area of CR is also filled up by rocks with relatively large size.
3. Flow Velocity Model
We consider a discrete model where the water discharge is represented by the corresponding two-dimensional grid in the modified DEM in which its position is denoted by indices . After determining the river network pattern, we model the corresponding flow velocity by assuming that its magnitude at particular position of water discharge along a river channel is given by the following equation: where and denote the flow velocity of water discharge at higher and lower positions, respectively, while is the height difference, where the value of is provided by DEM. The symbol represents the effective gravitational acceleration which is defined to accommodate the real condition of river channels and denotes water level from river bed at the position . Flow velocities in (1) can be interpreted as the average of velocity in the associated grid. In the meantime, the width and water level of each river branch are assumed to be uniform.
In general, we assume that the parameter depends on position and water level which can be determined empirically through field observation at particular time. The assumption of water level dependency is taken to accommodate the fact that water flow is more efficient as water level increases. But for the sake of simplicity, here we assume that is uniform along all river channels and it will be determined later numerically. Comparing with the flow velocity expression in GMS formula, this parameter plays similar role with the ratio where is hydraulic radius of a river channel and represents roughness coefficient [9].
Next, to model the merging of flows from two different branches at a river junction we assume that the collision of water at each junction is inelastic satisfying the mass-momentum conservation law in the following discrete numerical scheme: Here, and are the water discharge and momentum, respectively, as illustrated in Figure 2. Consider and , where and denote the water level from river bed and width of the river channel, respectively, while and are, respectively, the flow velocity vector and its magnitude at the corresponding position. Here, we assume that the channel cross-section is rectangular. Based on assumption that the collision is inelastic we can find from (1) and (2) the expression of flow velocity as follows: where is the angle between and . Note that after collision we only consider the related flow velocity magnitude, while its direction follows the shape of the river channel. To ensure that the collision is inelastic, it is important to emphasize that our computer code is set to run the water discharge with slower velocity first, which will be collided by the faster one.
It is worth to note that this junction model in principle can be extended to include more than two merging branches, namely, by adding another discharge and momentum terms on the left hand side of (2) and (3), respectively. Furthermore, this model can also be modified to accommodate the river branch bifurcation cases.
An example of a more complicated junction model has been proposed previously in [21]. In that model, flow from different river channel does not experience inelastic collision. Instead, by defining control volume to each channel, the water from both channels is assumed to flow side by side in the subsequent channel, where the momentum change of each flow is determined by hydrostatic pressure and pressure that related to the change of its control volume as well as interfacial shear force between them. Another junction model is given in [22] which also considers the application of continuity and mass-momentum conservation law incorporated with momentum and energy correction parameters. The implementations of these models are likely not easy in large scale cases and local boundary conditions are needed. Comparing with these two models, it is clear that our junction model is far simpler in its implementation, but it is probably less accurate indeed.
In addition, we also assume that after collision the width and water level will change according to the following chosen relations: for and for . The symbols and are small and dimensionless real parameters. These assumptions are taken to model the width and water level of river channel that resulted from the merging of two different branches. We consider the width and water level as two independent parameters since both are naturally predetermined. Here, these parameters are assumed to be uniform and they will be determined numerically. In general, however, one might consider that their values should not be the same for different junction and can be determined empirically. Indeed, one can choose another type of relation to model the width and water level . However, based on our numerical experiment which will be presented in the next section, the relations given by (5) and (6) so far are the most applicable ones to meet the corresponding river parameters and certainly this assumption must be further clarified by field observation which is beyond the scope of this report.
To this end, it is well known that based on the following definition of Froude number for a channel with rectangular cross-section where is the actual gravitational acceleration, one can classify the flow regimes into subcritical , critical , and supercritical flow [23]. It will be discussed on the next section that the model seems to accommodate only the case of subcritical flow.
4. Results of Model Simulation and Discussion
To simulate our model, we use the water discharge data sample found from our field surveys during dry season in August, 2013. Based on those data samples we consider dividing the 126 headwaters into five regions and choose Katulampa as the simulation point as shown in Figure 1(b). The corresponding data of headwaters and Katulampa point were considered as the associated initial and final conditions, respectively.
The initial flow velocity, width, and water level from field survey in each region are assumed to be uniform as given in Table 1. Prior to discussing the characteristics of our model with respect to the initial flow velocity and water level variations, we first simulate the flow velocity, water level, and width at Katulampa with initial values given in Table 1. Our simulation yields the following results: m/s, m, and m. On the other hand, from field survey, we observed the following data: m/s, m, and m. To meet the simulation results that match Katulampa data, we adjusted by trial and error the parameters in (1), (5), and (6) and found the following values: m/s2, , and , whereas in our computer code we run the initial water discharge from region 1 first and then from regions 2, 3, 4, and 5, respectively, since . It is worth to note that by assuming m/s2 all these velocities at headwaters and Katulampa are related to subcritical flows according to Froude number given by (7).
The simulation results of flow velocity, water level, and width of the whole catchment area of CR are given in Figure 3. Comparing to our data sample, the simulated flow velocity width and water level do not represent most of the actual conditions although at Katulampa these quantities are in good agreement. This can be explained as a consequence of not applying actual data of the whole headwaters to our model. One might expect that the use of actual data of all headwaters and considering position-dependent , , and will give fairly accurate results with tolerable deviations. As shown in Figure 3(a), it is found that at every river branch junction, the flow velocity of merged water discharge is getting slower due to inelastic collision given by (3), which corresponds to kinetic energy losses phenomenon.
(a)
(b)
(c)
Next, to simulate the characteristics of flow velocity at Katulampa, we vary the initial flow velocity at headwater of particular region while keeping the velocity of other regions the same. Given in Table 2 are five different variations of the first one, denoted by letters A, B, C, D, and E, with the corresponding results depicted in Figure 4. In variation A, we vary the initial flow velocity m/s. It is found that the flow velocity at Katulampa exhibits nonlinear variation as shown in Figure 4(a). However, it should be noted that the change is relatively small. In contrast, as shown in Figure 4(b), we observe for variation B, where m/s variation is considered, that the flow velocity at Katulampa does not vary at all. On the other hand, the other different characteristics are found for variations C and D, with the corresponding variations given as follows: m/s and m/s, respectively. As shown in Figures 4(c) and 4(d), both variations exhibit opposite characteristics, where in variation C the flow velocity is increasing almost linearly, while for variation D the flow velocity is decreasing linearly. In the meantime, the variation E given in Figure 4(e) is similar to variation A where its variation is nonlinear with respect to the linear variation of initial flow velocity m/s.
(a)
(b)
(c)
(d)
(e)
For comparison, we also consider another five variations of the second one, denoted by letters F, G, H, I, and J, with the corresponding variations of initial flow velocity of each region given in Table 3. Here, we exchange the initial flow velocity between region 4 and region 5 such that . The results are depicted in Figure 5. Comparing variations A, B, and C and variations F, G, and H of Figures 4(a)–4(c) and Figures 5(a)–5(c), respectively, we see that similar characteristics are exhibited for variations of regions 1, 2, and 3, respectively, but with different flow velocity magnitudes. Interestingly, this similarity is not the case for variations I and J given in Figures 5(d) and 5(e), respectively, where both variations demonstrate flat characteristics. These results are in contrast to the previous variations D and E given in Figures 4(d) and 4(e).
(a)
(b)
(c)
(d)
(e)
To complete our investigation on the flow velocity characteristics at Katulampa with respect to initial flow velocity variation, we also simulate the model by uniformly increasing the value of initial flow velocity of each region. The results are given in Figure 6 showing nonlinear variation of velocity at Katulampa with respect to the uniform increment of all initial flow velocities by .
Clearly, since it is considered that the water level does not vary and that the related velocity at Katulampa tends to grow exponentially, then the corresponding Froude number at that point, which is given by (7), can probably be of . Therefore, based on this fact it is safe to say that our model most likely cannot describe the supercritical flow in a fairly accurate way.
All these results obviously indicate that the influence of initial flow velocity variation at headwaters will lead to specific characteristics of flow velocity at lower land. However, the associated characteristics cannot be easily explained due to complex interaction under inelastic collision among the corresponding river branches. This situation can be exemplified by the nonlinear characteristic of variation E which is clearly different compared to variation I with flat characteristic, although both variations have the same initial flow velocity variation on the related region (0.5270–14 m/s) and the location of the corresponding headwaters is relatively close to Katulampa as shown in Figure 1.
We have discussed previously in Section 3 that in general the parameter depends on water level. Therefore, in addition to the aforementioned variations, we have also investigated the characteristics of flow velocity and water level at Katulampa by increasing linearly the water level of each headwater as follows: with m, where is an integer. In the meantime, the effective gravitational acceleration is considered to increase linearly based on the following equation: which is taken under the assumption that higher water level will flow more easily. Here is given in Table 1. The result is shown in Figure 7. It is found that the flow velocity variation at Katulampa increases nonlinearly with tendency of saturation while the corresponding water level increases linearly.
(a)
(b)
Finally, it should be admitted that our simple model must be further examined by field observation which should be conducted in relatively long time period. However, this model can be used to determine and predict qualitatively the characteristics of river flow velocity in large scale. In addition, for a relatively small grid DEM, in principle, the calculated velocity, water level, and width can be used as initial condition of time-dependent SVEs.
5. Summary
We have discussed our simple computer based model for simulating the river flow velocity. We consider Ciliwung River as the area of simulation where the corresponding digital elevation model (DEM) is based on SRTM data from USGS with m2 resolution and modified by m2 as well as data from field survey. The flow velocity is described by equation of motion of falling body in inclined plane due to gravitational force. We define the effective gravitational acceleration parameter to accommodate the real condition of river channel. At the junction of two river branches we assume that the merged flow velocity is determined by inelastic collision. It is found that our model can determine the flow velocity at Katulampa in good agreement with observation data. The variations of initial flow velocity of certain region in the Ciliwung River headwaters show various characteristics that we suspect to occur due to complex interaction among the river branches. This simple model can be used to determine and predict the characteristics of river flow velocity at subcritical level as well as providing an initial condition for time-dependent Saint Venant equations.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
This work was funded by BOPTN Research Grant of year 2013 from Directorate of Higher Education, Ministry of National Education and Culture, Republic of Indonesia, under Contract no. 198/IT3.41.2/L2/SPK/2013.