Abstract

Bed topography in river bends is highly non-uniform as a result of the spiral motion of fluid and sediment transports related to channel curvature. To grasp a full understanding of geomorphology and hydrology in natural river bends, detailed bed topography data are necessary, but are usually not of high enough quality and so require further interpolation for sophisticated studies. In this paper, an algorithm is proposed that is particularly suited to bathymetry interpolation in rivers with apparent bends. The thalweg and the two banks are used as geographical features to ensure that a concave cross-sectional bed-form can be found in bend geometry, while linear interpolations are conducted in the in accordance with secondary and main stream currents, respectively. In comparison with conventional spatial interpolation methods, the proposed algorithm is validated to ensure better performance in generating smooth and accurate bed topography in channel bends, which results in better predictions of river stage by 2D hydrodynamic simulation in practical field tests.

INTRODUCTION

The formation and the evolution of bed topography are very complicated in natural rivers. In bend sections, bed topographies are more non-uniform than those in straight sections due to the increase of flow complexity. When river flow enters a channel bend, the centrifugal force drives the near-surface flow currents toward the outer bank, which results in a superelevation of water surface and a difference of water pressure between the outer and the inner banks. It drives the near-bed flows from the outer bank toward the inner one and causes the bed topography to be eroded at the outer bank and deposited at the inner one. To distinguish from the main stream currents along the channel direction, this cross-sectional flow motion is called the secondary currents, and their combination exhibits a 3D spiral pattern of flow motion in channel bends (Thompson 1876; Rozovskii 1961; Yen 1965; Yen 1967; de Vriend 1981). To fully analyze the interactions between flow and sediment in channel bends, a number of 2D/3D models for fluvial research have been developed in recent years (Leschziner & Rodi 1979; Wu et al. 2000; Duan & Nanda 2006; Abad et al. 2008; Jang et al. 2011). Along with the improvement of numerical models, the demand for high-resolution bed topography has increased as well.

In the field, SONAR (Sound Navigation And Ranging; Smoot & Novak 1969; Simpson & Oltmann 1993) and GPR (Ground-Penetrating Radar; Spicer et al. 1997; Costa et al. 2000) are two of the most commonly used techniques to retrieve bed elevations below water surface. The measurements of bed topography in natural rivers are usually conducted at discrete points across a cross-section and at discrete cross-sections along a channel. However, limited by budgets and manpower, the point intervals in the longitudinal direction (typically about 500 meters along the river) are often larger than those in the transverse direction (typically within 30 meters across the river). Hence, further bathymetry interpolation on the observed dataset is required to obtain continuous bed topography for sophisticated flow and sediment analysis in river channels.

Bathymetry interpolation determines the elevation at a given cell by weightedly averaging the elevations of surrounding observation points according to their relative distances. The weighted method can be geometric, e.g. IDW (inverse distance weight), which weights the values of neighboring points linearly by their distances to the interested cell; or geostatistical, e.g., Kriging, which weights the values of measured points by their overall spatial arrangement described by statistical models. In general, the distance-based weighted methods perform well when the measured data are sufficient and isotropic around the interpolated cell. However, due to the anisotropy and inadequacy of bed data in natural river bends, bias in bathymetry interpolation can be significant without extra considerations.

Many researchers have pointed out the importance of dealing with the anisotropy in river bathymetry interpolation (Carter & Shankar 1997; Burroughes et al. 2001; Andes & Cox 2017). A practical method to account for the anisotropy associated with river direction is to conduct bathymetry interpolation in a flow-oriented curvilinear coordinate system (Goff & Nordfjord 2004; Merwade et al. 2006, 2008a, 2008b). Schäppi et al. (2010) further demonstrated that the error of bathymetry interpolation can also be well confined in a Cartesian coordinate system by implementing linear interpolation in the directions lateral and longitudinal to a river channel, separately, between several manually defined breaklines. Caviedes-Voullième et al. (2014) used an algorithm similar to Schäppi et al. (2010) but replaced the breaklines with cubic Hermite splines parallel to the main flow direction, called the thalweg trajectory interpolation. Chen & Liu (2017) proposed a linear interpolation method similar to Caviedes-Voullième et al. (2014), with the utilization of three breaklines to distinguish the main channel. In fact, these variations of breakline number can be attributed to the difference in the geometry of experiment channels, e.g., the river channel used by Schäppi et al. (2010) is relatively straight while those used by Caviedes-Voullième et al. (2014) and Chen & Liu (2017) are highly sinuous.

In this study, we focus on reconstructing the bed topography in channels with apparent sinuous bends from the viewpoint of cross-sectional flow/sediment transportation. The algorithm by Chen & Liu (2017) is modified by adopting the thalweg as the only breakline to feature a typical concave bed-form between two river banks. The procedure of interpolation is described in detail and the Zeng-Wen River located in southern Taiwan is selected as the study area to validate the proposed algorithm through comparison between the interpolated and the measured bed elevations. The interpolated bed topographies obtained by the proposed algorithm are then put into a 2D flow model to test the correctness of simulated river stage with two historical rainfall events that occurred in 2012 and 2013. Two other conventional spatial interpolation methods, IDW and ordinary Kriging, are also included for comparison in terms of morphological and hydrodynamic performance.

STUDY AREA

The Zeng-Wen River originates from Mt Alishan in central Taiwan, flowing westward into the Taiwan Strait, with a total length of 138.5 km, a basin area of 1,177 km2, and an average slope of 1/200. Historically, the Zeng-Wen River was a frequently meandering river before the completion of embankment in 1938. Every couple of years, the bed topography is resurveyed and updated by the Water Resource Agency in Taiwan. In the present study, the downstream segment of the Zeng-Wen River from Er-Xi Bridge to the estuary is adopted for research. According to the survey in 2012, there are 12,093 data points of bed elevation spread among 128 cross-sections in the study segment, comprising two bends separated by three water-level stations located at Er-Xi Bridge, Xin-Zhong Station, and Ma-Shan Bridge, as shown in Figure 1.

Figure 1

The study segment of Zeng-Wen River.

Figure 1

The study segment of Zeng-Wen River.

BEND-ORIENTED BATHYMETRY INTERPOLATION (BBI)

In Chen & Liu (2017), three breaklines, including the thalweg and two relatively large bed elevations that delineate the main channel, are adopted to describe a cross-sectional bed-form between two banks. Their methodology includes four steps: the first step is to determine the point with the lowest elevation and two other points with maximum transverse bed slopes, by which each cross-section is divided into four parts including left bank, main channel in left, main channel in right, and right bank; the second step is to redistribute the number of cross-sectional points based on these four parts and to conduct linear interpolation in the transverse direction; the third step is to generate extra cross-sections between each pair of cross-sections to make sure each part at each cross-section has the same number of points; the final step conducts linear interpolation in the longitudinal direction to acquire the bed elevation between cross-sections. The uncertainty of Chen & Liu's (2017) method lies in the determination of the two points with maximum transverse bed slopes, which vibrates a lot due to the insufficiency of survey data.

To cope with it, a BBI algorithm is proposed, by which the thalweg is adopted as the only breakline and the interpolation is implemented separately in the two divisions bounded by the thalweg and the two banks along the channel direction. This is because, for a river channel with typical bends, the transverse sediment transport driven by secondary currents will ensure a concave shape of cross-sectional bed-form, which can be well described by the three deflections of thalweg and two banks. The BBI adopts linear interpolation because it is easily implemented and has been tested and found to be more robust than spline or cubic spline methods in preventing overshooting for elevation interpolation (Schäppi et al. 2010). Meanwhile, as the thalweg symbolizes the locations where transverse bed slope deflects, the bed topography is smooth and continuous on each side of the thalweg so can be interpolated by linear algorithm.

The BBI includes three steps: (1) determine the thalweg and river banks based on the measured bed elevation data; (2) conduct transverse interpolation at all measured cross-sections in the divisions bounded by the thalweg and the two banks; (3) conduct longitudinal interpolation in the divisions between two adjacent measured cross-sections. These three steps are described in detail in the following sections.

Determination of thalweg and banks

By definition in fluvial geomorphology, thalweg is a line joining the lowest bed elevations along a river that defines the deepest channel in a river valley. At each cross-section in a channel bend the thalweg represents the point where the largest bed erosion occurs under the effects of secondary currents, and the bed elevation basically increases from the thalweg to the two banks, as shown in the schematic diagram of Figure 2. In the Cartesian coordinate system, the coordinates for the thalweg are determined by the data points with the lowest cross-sectional bed elevations, while the coordinates for the right bank and the left bank are determined by the points with the highest bed elevations on the right-bank and left-bank sides of the thalweg facing river downstream, respectively; where i is the number of cross-section. As the embankment of the Zeng-Wen River was built high enough to defend from flooding under a 100-year recurrence period, the overflow of river water occurred rarely and the coordinates of the two banks can be simply determined by the dike locations. In addition to field survey, the thalweg and two banks can also be delineated from the interpretation of satellite images, such as Figure 1, where thalweg and banks can be clearly observed. In the case of natural areas with lots of vegetation in the floodplain, the river banks can be delineated by: (1) identifying the boundary of riparian vegetation from the interpretation of multispectral information in satellite images (Lin et al. 2013); (2) using penetrating LiDAR or laser to delineate the bank where abrupt discrepancy of ground elevation occurs (Scherer et al. 2012); (3) manual measurement of ground elevation across the floodplain by handheld GPS devices (Wang 2011).

Figure 2

Schematic diagram for determining thalweg and banks.

Figure 2

Schematic diagram for determining thalweg and banks.

Transverse interpolation

In the transverse direction, interpolations are conducted only at measured cross-sections. To ensure a concave shape of bed topography with the lowest elevation at the thalweg for each cross-section, the interpolation is implemented independently on the right-bank and the left-bank sides. Meanwhile, at all cross-sections, the interpolated points on the same side of the channel are designed to be the same in number with uniform divisions. Thus, a one-on-one relation of interpolated points between different cross-sections can then be established in the longitudinal direction. The relation between the measured and interpolated points is schematically displayed in Figure 3, in which represent the coordinates for the points to be determined while represent the coordinates for the measured points. At each cross-section i, the can be determined linearly as below:  
formula
(1)
where j is the number of interpolated points at each cross-section; C symbolizes the spatial coordinate for the interpolated points (X or Y); and c represents the coordinate for the measured points (x or y). From Equation (1), one can see that for , for , and for due to the overlay of the measured and the interpolated points at the thalweg and two banks. It should be noted that the numbers of R and L must be large enough to capture all measured points. An overall check on the distance between every two adjacent measured points for all cross-sections shows that the minimum measured distance is 2 m. Thus, to make sure the minimum sampling division is smaller than 2 m, the number of cross-sectional divisions is required in this study.
Figure 3

Schematic diagram for transverse interpolation.

Figure 3

Schematic diagram for transverse interpolation.

Once the coordinates are determined, the bed elevation can be interpolated linearly as below:  
formula
where  
formula
(2)
in which and are the coordinates and bed elevations for the measured points at cross-section i closest to on the opposite sides; is the distance between and ; is the distance between and .

Longitudinal interpolation

In the longitudinal direction, the coordinates and bed elevations of the interpolated points between two adjacent measured cross-sections are calculated linearly as below:  
formula
(3)
in which i = 0, … , n, representing the number of measured cross-sections along the river channel; k = 0, … , m, representing the number of interpolated cross-sections between two measured cross-sections i and . Although the amount of k can differ from section to section, a constant value is assigned herein because the distances between measured cross-sections are roughly the same. This deliberate arrangement of symbolism in Equation (3) ensures that equals and when and , respectively. In total, there will be a mesh system containing grids, in which the mesh points have to be interpolated. The grid system for longitudinal interpolation is shown in Figure 4 and the projective correspondence between the grid mesh and the interpolated bed topography is schematically displayed in Figure 5. It is seen that, with the BBI algorithm, a concave and smooth cross-sectional bed-form can be obtained throughout a channel bend.
Figure 4

Schematic diagram for longitudinal interpolation.

Figure 4

Schematic diagram for longitudinal interpolation.

Figure 5

Spatial correspondence between grid mesh and interpolated bed topography.

Figure 5

Spatial correspondence between grid mesh and interpolated bed topography.

VALIDATION

The performance of the proposed algorithm, BBI, is tested separately in aspects of geomorphology and hydrology. Two conventional interpolation methods, IDW and ordinary Kriging, are included for comparison by the discrepancy between the observations and predictions in terms of bed elevation and river stage. The bed elevations measured in 2012 and the river stages observed during two historical rainfall events in 2012 and 2013 are selected for case studies. The IDW and ordinary Kriging are chosen because they are two of the most used interpolation algorithms without considering anisotropy. The simplicity of the two algorithms echoes the BBI algorithm in this study, which uses linear interpolation instead of anisotropic interpolation across and along the river direction.

Bed elevation accuracy

Figure 6 shows the comparison of the observed and the interpolated bed elevations at six random cross-sections distributed along the river segment by the three interpolation algorithms. At the outlet and entrance of the river segment, the thalwegs are roughly located at the middle of the cross-sections, as shown in Figure 6(a) and 6(f), and the three interpolation algorithms show insignificant discrepancies due to the symmetry of cross-sectional bed topographies. For the cross-sections around the bends, the thalweg in Figure 6(d) and 6(e) shift toward the right banks at the bend entrance while the thalwegs in Figure 6(b) and 6(c) shift toward the left banks at the bend outlet. This shifting of thalweg is induced by the strong secondary currents that develop and decline as bend curvature increases and decreases, respectively. Moreover, the depth of thalweg in the bend section is observed to be deeper by about 2 meters compared with that in the straight section, showing the evidence of secondary currents that increase transverse bed slope by enhancing cross-sectional erosion and sedimentation. In the bends, the bed-form generated by the BBI is much closer to the observations compared to the IDW and Kriging. In Figure 6(d) and (e), the Kriging generates significant errors in estimating bed elevations around thalwegs; while in Figure 6(c), the IDW show obvious inaccuracy in delineating the main channel. It is interesting to observe that, even for the cross-section with two bed depressions as shown in Figure 6(e), the bed elevations can be well predicted by the proposed BBI.

Figure 6

Comparison of observed and interpolated cross-sectional bed elevations.

Figure 6

Comparison of observed and interpolated cross-sectional bed elevations.

Two-dimensional bed topographies interpolated by the three algorithms are shown in Figure 7, in which the second bend of the study river segment is particularly zoomed in on for demonstration. One can see that the bed topographies along the thalweg obtained by IDW and Kriging are more rugged and discontinuous than those obtained by the BBI method. This phenomenon is attributed to the fact that the IDW and Kriging weight the neighboring measured points purely by distance to the relevant cell without considering river directions, which leads to the bias of bed prediction especially in bends where bed anisotropy increases with channel curvature. In contrast, the BBI well maintains the continuity of bed topography along the river direction even under data shortage. From the viewpoint of fluid dynamics, this discontinuity of bed topography will increase the friction between water and bed surface that further results in false estimations of flow velocity and river stage.

Figure 7

Two-dimensional bed topographies interpolated by (a) BBI, (b) IDW, and (c) Kriging.

Figure 7

Two-dimensional bed topographies interpolated by (a) BBI, (b) IDW, and (c) Kriging.

The differences between the predicted and the measured bed topographies are usually evaluated to test the accuracy of bathymetry interpolation, which, however, cannot be done when all measured data are used for interpolation without a spare for validation. In this study, an alternative validation method is adopted by removing a measured cross-section one at a time from interpolation. The interpolated results are compared with the measured data at the removed cross-section using the following indicators (Murphy 1988):  
formula
(4)
 
formula
(5)
in which is the interpolated bed elevations at point i; is the measured bed elevations at point i; is the number of points, equal to the amount of total measured points since each cross-section is sequentially removed; is the root-mean-square error of bed elevation; is the bias of bed elevation in percentage. Positive values of indicate overall overestimation of bed elevation while negative values represent overall underestimation of bed elevation.

The relative frequencies of bed elevation error by the three interpolation methods are displayed in Figure 8. It shows that the bed elevation errors by the BBI are more concentrated around zero with less extreme values compared with those by the IDW and the Kriging. It also shows that the overestimation and underestimation of the bed elevations are basically the same for individual algorithms due to the symmetry of the three lines in Figure 8. Listed in Table 1 are the error statistics by the three algorithms. The BBI has the highest accuracy with 82.58% of lying within −0.1 to 0.1 meters out of the 12,093 points. The value of by the BBI is 0.21 m, much smaller than those predicted by the IDW (1.18 m) and the Kriging (1.61 m). The values of show that the bed topography predicted by the BBI generates only 1.38% of bias related to the measurement, while the bias is six times larger by the IDW and eleven times larger by the Kriging.

Table 1

Error statistics of bed elevation interpolation

Elevation error Algorithm
 
Kriging IDW BBI (L = R = 100) BBI (L = R = 250) 
Points with (%) 54.57 69.80 78.18 82.58 
(m1.61 1.18 0.71 0.21 
(%) 17.15 9.41 4.59 1.38 
Elevation error Algorithm
 
Kriging IDW BBI (L = R = 100) BBI (L = R = 250) 
Points with (%) 54.57 69.80 78.18 82.58 
(m1.61 1.18 0.71 0.21 
(%) 17.15 9.41 4.59 1.38 
Figure 8

Frequency distribution of bed elevation error.

Figure 8

Frequency distribution of bed elevation error.

To further test the validity of BBI, sensitivity analysis has been conducted with respect to the density of sampling points inside a cross-section and to the number of measured cross-sections as well. The bed elevations interpolated with a lower density of cross-sectional points are compared to those with in Table 1 and Figure 8. Table 1 shows that the values of and decrease by 0.5 m and 3.21%, respectively, when the number of sampling points increases from 100 to 250. With respect to the cross-section number, the performances of bed interpolation are tested under three scenarios: (1) keep all measured cross-sections in the bend; (2) remove one every two measured cross-sections in the bend; (3) remove all central measured cross-sections in the bend. The topographies and the errors of bed interpolation under these three scenarios are displayed in Figure 9 and summarized in Table 2, respectively. Figure 9 shows that, for scenarios (2) and (3) above, the thalweg and two banks delineated from satellite images coincide well with the measurement, whereas some details of bed topography are inevitably missing due to the lack of data. Table 2 shows that the errors are the highest when removing the central cross-sections, second highest when removing one every two cross-sections, and the smallest when keeping all cross-sections. However, even in the worst case, there are still 60% of the sampling points having elevation error smaller than 0.1 m, with and . Overall, the sensitivity analysis shows that the BBI algorithm is adequate for bed reconstruction in bends under insufficient measured data.

Table 2

Sensitivity analysis of BBI in channel bend

Elevation error Scenario
 
Keep all cross-sections Remove one every two cross-sections Remove central cross-sections 
Points with (%) 81.27 74.79 60.09 
(m0.29 0.89 1.43 
(%) 1.89 7.99 13.33 
Elevation error Scenario
 
Keep all cross-sections Remove one every two cross-sections Remove central cross-sections 
Points with (%) 81.27 74.79 60.09 
(m0.29 0.89 1.43 
(%) 1.89 7.99 13.33 
Figure 9

Bed topographies in the bend interpolated by BBI by (a) keeping all measured cross-section, (b) removing one every two measured cross-sections, and (c) removing all central measured cross-sections.

Figure 9

Bed topographies in the bend interpolated by BBI by (a) keeping all measured cross-section, (b) removing one every two measured cross-sections, and (c) removing all central measured cross-sections.

River stage simulation accuracy

As bed and flow affect each other, the robustness of bathymetry interpolation can be indirectly verified by the comparison between the observed and the simulated river stages under different interpolated bed topographies. In order to more precisely capture the influence of bed variation on the flow field in channel bends, a 2D unstructured-grid hydrodynamic model is used in this study for river stage simulation. The model solves the Saint-Venant equations (Barré de Saint-Venant 1871) (also known as shallow water equations) using Galerkin's finite-element scheme (Ciarlet 1978) in dealing with weighted residuals. The governing equations of the 2D flow model with hydrostatic form and Boussinesq approximation in Cartesian coordinate system are given as:  
formula
(6)
 
formula
(7)
 
formula
(8)
in which is the free water surface elevation; and are the horizontal velocities in x, y direction, respectively; g is the acceleration due to gravity; is the earth's tidal potential; is the effective earth elasticity factor; is the density of water; is the atmospheric pressure at the free water surface; and are wind shear stresses in x, y direction, respectively; and are bottom shear stresses in x, y direction, respectively; and , where is the bed elevation. The formulas of bottom shear stress can be calculated as below:  
formula
(9)
 
formula
(10)
where is the bottom drag coefficient which can be parameterized using the Manning formulation:  
formula
(11)
where is Manning's roughness coefficient. This approach provides a depth-dependent bottom drag coefficient that gives larger values for shallow-water areas and smaller values for deep-water areas.

Since the bed data in 2012 are adopted for interpolation, two rainfall events close to the survey time, Rainstorm 0612 (June 10–16 in 2012) and Typhoon Kongrey (August 28–September 1 in 2013), are selected for case studies. In both cases, river overflow did not occur. The hourly discharges at Er-Xi Bridge recorded during the two events (peak discharge 3,680 cm for Rainstorm 0612; peak discharge 6,854 cm for Typhoon Kongrey) are input as the upstream conditions for hydrodynamic simulation. Figure 10(a) and 10(b) show the comparison of the observed and the simulated river stages during Rainstorm 0612 at Ma-Shan Bridge and Xin-Zhong Station, respectively, in which the circles represent the observation and the lines are the river stages simulated under the bed topographies interpolated by the BBI, the IDW, and the Kriging algorithms. At Ma-Shan Bridge or Xin-Zhong station, the discrepancies between the simulated and observed river stages increase sequentially for BBI, IDW, and Kriging, which coincides with the performance of bed interpolation. Figure 11(a) and 11(b) show the comparison of the observed and the simulated river stages at Ma-Shan Bridge and Xin-Zhong station during Typhoon Kongrey, respectively. Again, the agreement between the simulated and the observed river stages for BBI is much better than the other two, especially around the peak.

Figure 10

Simulated river stages for Rainstorm 0612 in 2012 at (a) Ma-Shan Bridge and (b) Xin-Zhong Station.

Figure 10

Simulated river stages for Rainstorm 0612 in 2012 at (a) Ma-Shan Bridge and (b) Xin-Zhong Station.

Figure 11

Simulated river stages for Typhoon Kongrey in 2013 at (a) Ma-Shan Bridge and (b) Xin-Zhong Station.

Figure 11

Simulated river stages for Typhoon Kongrey in 2013 at (a) Ma-Shan Bridge and (b) Xin-Zhong Station.

The accuracy of the simulated river stage can be evaluated by the indicators similar to those for bed topography, and they are shown as below:  
formula
(12)
 
formula
(13)
where is the observed river stage at time stage i at a specific water-level station; is the simulated river stage at time stage i at the observed point; and is the number of hours in which river stages are output for validation. It should be noted that the definition of herein is different from that of in Equations (4) and (5); the former denotes the number of temporal time stages while the latter denotes the number of spatial points.

The error statistics of river stage simulation for the three interpolated algorithms are summarized in Table 3. In the 0612 event, the equals 3.1 m and the equals 20% for BBI on average, much smaller than those for IDW (; on average) and Kriging (; on average). Notably, the errors in simulated river stage become larger at the downstream station (Ma-Shan) in comparison with those at the upstream station (Xin-Zhong), implying that the hydraulic error that accumulates is attributed to the inaccuracy of bed topography as flow travels downstream. This phenomenon is reasonable due to the fact that fluid plays the role of a sponge that absorbs topographical errors proportional to the area it sweeps. However, for the BBI, the increase of (about 0.17 m) is the smallest among the three algorithms. Similar patterns can be found in the high discharge event of Kongrey, but the errors in river stage simulation generally reduce due to the fact that the influence of bed condition degenerates as fluid momentum increases. Again, in the Kongrey event, the river stage simulation errors for BBI are still the smallest (with and ), showing that accurate bathymetry interpolation can greatly benefit the river stage forecast at high or low discharges.

Table 3

Error statistics of river stage simulation

River stage error Rainstorm 0612 (2012)
 
Typhoon Kongrey (2013)
 
Kriging IDW BBI Kriging IDW BBI 
Ma-Shan  (m5.63 3.88 3.20 4.55 4.07 2.49 
(%) 91.72 47.39 25.54 45.75 17.59 5.32 
Xin-Zhong  (m4.28 3.32 3.07 3.82 3.06 2.46 
(%) 32.33 18.47 15.33 25.52 13.33 2.36 
River stage error Rainstorm 0612 (2012)
 
Typhoon Kongrey (2013)
 
Kriging IDW BBI Kriging IDW BBI 
Ma-Shan  (m5.63 3.88 3.20 4.55 4.07 2.49 
(%) 91.72 47.39 25.54 45.75 17.59 5.32 
Xin-Zhong  (m4.28 3.32 3.07 3.82 3.06 2.46 
(%) 32.33 18.47 15.33 25.52 13.33 2.36 

In order to discover the influence of bed topography on two-dimensional flow field, the simulated peak water depths over the channel during Rainstorm 0612 and Typhoon Kongrey are displayed in Figures 12 and 13, respectively. In the low discharge event of Rainstorm 0612, the simulation by the BBI shows that the water is mainly concentrated around the thalweg with water depth gradually decreasing to zero towards the banks. In the cases of IDW and Kriging, the water surface becomes wavier due to the rugged bed-form, which disperses extra water from thalweg towards the two banks inconsistent with typical circular motion in natural river bends. In the case of high discharge event of Kongrey, a similar phenomenon can be found but the water level is boosted more abruptly. On the right-bank side, the water depth simulated under IDW and Kriging reaches up to 10 meters which could greatly increase false flood alarms.

Figure 12

Spatial distribution of simulated water depth for Rainstorm 0612 in 2012 using the topographies generated by (a) BBI, (b) IDW, and (c) Kriging.

Figure 12

Spatial distribution of simulated water depth for Rainstorm 0612 in 2012 using the topographies generated by (a) BBI, (b) IDW, and (c) Kriging.

Figure 13

Spatial distribution of simulated water depth for Typhoon Kongrey in 2013 using the topographies generated by (a) BBI, (b) IDW, and (c) Kriging.

Figure 13

Spatial distribution of simulated water depth for Typhoon Kongrey in 2013 using the topographies generated by (a) BBI, (b) IDW, and (c) Kriging.

CONCLUSIONS

Bed topography in natural river bends is highly non-uniform and usually lacks detailed measured data due to the limitation of budgets and manpower. To realize the variation of bed topography in areas without sufficient data, proper interpolation of bed topography is crucial for sophisticated fluvial studies. However, without considering the anisotropy of bed topography resulting from bend effects, bed interpolation error may be significant. In this study, a BBI algorithm is proposed to consider the bend effects by interpolating river bathymetry in the directions across and along the channel direction, separately. Meanwhile, the thalweg and two banks are delineated as the boundaries in the process of interpolation to ensure the continuity and smoothness of bed topography. To examine the interpolation algorithm, a segment of Zeng-Wen River in central Taiwan is selected as the study area, and two algorithms without considering bend effects, IDW and ordinary Kriging, are included for comparison.

Results show that the bed elevations interpolated by the BBI agree much better with the measured data, having errors and deviations clearly smaller than those predicted by the IDW and the Kriging. Moreover, the bed topography generated by the BBI has a smoother and more realistic pattern in response to the continuity of flow and sediment transportation. To demonstrate the influence of bed interpolation on river hydrology, a 2D hydrodynamic model is introduced to simulate the river stages during two rainfall events with low and high discharges in 2012 and 2013, respectively. In both events, the errors in simulated river stage under the bed topography generated by BBI are the smallest among the three algorithms. In the case of BBI, the water in bends is conveyed more smoothly via the thalweg channel, which helps maintain the structure of secondary currents and the reasonability of water surface near banks. Simulated results also indicate that hydraulic errors can be accumulated due to the inaccuracy of bed topography as river water travels downstream. The success of BBI implies that the poor performance of IDW and original Kriging does not come from the algorithm simplicity but from the lack of consideration of bend influence. This coincides with the findings of Merwade et al. (2006), who showed that anisotropic Kriging requires additional transformation according to river orientation before performing better than ordinary Kriging.

Overall, the study has demonstrated that the BBI is an accurate and easily applied algorithm for bed topography interpolation when measured data are limited in natural river bends. The utilization of a global Cartesian coordinate system in both bathymetry interpolation and hydrodynamic simulation is very helpful to reduce the efforts in coordinate transformation in association with river sinuosity. Nevertheless, the proposed algorithm should be more applicable for channels with sinuous bend geometry that generate circular motions strong enough to shape up a concave cross-sectional bed-form. For a relatively straight channel, bed topography is shaped mainly by longitudinal fluid forces that can result in many bars and thalwegs at the same time staggered across a section. In this case, more complex interpolation algorithms in previous research should be considered. In the future, sensitivity tests should be carried out to study the performance of different algorithms related to river sinuosity and data density.

ACKNOWLEDGEMENTS

The authors would like to express their sincere gratitude to the National Science and Technology Center for Disaster Reduction for the computer facilities and the Water Resource Agency for providing the field data.

REFERENCES

REFERENCES
Abad
,
J. D.
,
Buscaglia
,
G.
&
Garcia
,
M. H.
2008
2D stream hydrodynamic, sediment transport and bed morphology model for engineering applications
.
Hydrological Processes
22
(
10
),
1443
1459
.
Andes
,
L. C.
&
Cox
,
A. L.
2017
Rectilinear inverse distance weighting methodology for bathymetric cross-section interpolation along the Mississippi River
.
Journal of Hydrologic Engineering
22
(
7
),
04017014-1
04017014-12
.
Barré de Saint-Venant
,
A. J. C.
1871
Théorie du mouvement non permanent des eaux, avec application aux crues des rivières et a l'introduction de marées dans leurs lits (Theory of nonpermanent movement of water with application of floods of rivers and to the introduction of tides within their beds)
.
Comptes Rendus des Séances de l'Académie des Sciences
73
,
147
154
,
237–240
.
Burroughes
,
J. E.
,
George
,
K. J.
&
Abbott
,
V. J.
2001
Interpolation of hydrographic survey data
.
The Hydrographic Journal
99
,
21
29
.
Carter
,
G. S.
&
Shankar
,
U.
1997
Creating rectangular bathymetry grids for environmental numerical modeling of gravel-bed rivers
.
Applied Mathematical Modelling
20
(
11
),
699
708
.
Caviedes-Voullième
,
D.
,
Morales-Hernández
,
M.
,
López-Marijuan
,
I.
&
García-Navarro
,
P.
2014
Reconstruction of 2D river bends by appropriate interpolation of 1D cross-sectional information for flood simulation
.
Environmental Modelling & Software
61
,
206
228
.
Ciarlet
,
P. G.
1978
The Finite Element Method for Elliptic Problems
.
North-Holland
,
Amsterdam, New York, Oxford
.
Costa
,
J. E.
,
Spicer
,
K. R.
,
Cheng
,
R. T.
,
Haeni
,
F. P.
,
Melcher
,
N. B.
,
Thurman
,
E. M.
,
Plant
,
W. J.
&
Keller
,
W. C.
2000
Measuring stream discharge by non-contact methods: a proof-of-concept experiment
.
Geophysical Research Letters
27
,
553
556
.
de Vriend
,
H. J.
1981
Velocity redistribution in curved rectangular channels
.
Journal of Fluid Mechanics
107
,
423
439
.
Jang
,
J. H.
,
Ho
,
H. Y.
&
Yen
,
C. L.
2011
Effects of lifting force on bed topography and bed-surface sediment size in channel bend
.
Journal of Hydraulic Engineering
137
(
9
),
911
920
.
Leschziner
,
M. A.
&
Rodi
,
W.
1979
Calculation of strongly curved open channel flow
.
Journal of the Hydraulics Division.
105
(
10
),
1297
1314
.
Lin
,
C. H.
,
Chuang
,
H. K.
,
Lin
,
M. L.
&
Huang
,
W. C.
2013
Establishment of the watershed image classified rule-set and feasibility assessment of its application
. In:
Proceedings of the International Conference on Image Processing, Computer Vision, and Pattern Recognition (IPCV)
,
Las Vegas
.
Merwade
,
V. M.
,
Maidment
,
D. R.
&
Goff
,
J. A.
2006
Anisotropic considerations while interpolating river channel bathymetry
.
Journal of Hydrology
331
,
731
741
.
Merwade
,
V.
,
Cook
,
A.
&
Coonrod
,
J.
2008a
GIS techniques for creating river terrain models for 409 hydrodynamic modeling and flood inundation mapping
.
Environmental Modelling & Software
23
,
1300
1311
.
Merwade
,
V.
,
Olivera
,
F.
,
Arabi
,
M.
&
Edleman
,
S.
2008b
Uncertainty in flood inundation mapping: current issues and future directions
.
Journal of Hydraulic Engineering
13
,
608
620
.
Rozovskii
,
I. L.
1961
Flow of Water in Bends of Open Channels
.
The Israel Program for Scientific Translations
,
Jerusalem
.
Schäppi
,
B.
,
Perona
,
P.
,
Schneider
,
P.
&
Burlando
,
P.
2010
Integrating river cross section measurements with digital terrain models for improved flow modelling applications
.
Computers & Geosciences.
36
,
707
716
.
Scherer
,
S.
,
Rehder
,
J.
,
Achar
,
S.
,
Cover
,
H.
,
Chambers
,
A.
,
Nuske
,
S.
&
Singh
,
S.
2012
River mapping from a flying robot: state estimation, river detection, and obstacle mapping
.
Autonomous Robots
33
(
1–2
),
189
214
.
Simpson
,
M. R.
&
Oltmann
,
R. N.
1993
Discharge-measurement System Using an Acoustic Doppler Current Profiler with Applications to Large Rivers and Estuaries
.
US Government Printing Office
,
Washington, DC
.
Smoot
,
G. F.
&
Novak
,
C. E.
1969
Measurement of Discharge by the Moving-Boat Method
.
US Government Printing Office
,
Washington, DC
.
Wu
,
W.
,
Rodi
,
W.
&
Wenka
,
T.
2000
3D numerical modeling of flow and sediment transport in open channels
.
Journal of Hydraulic Engineering
126
(
1
),
4
15
.
Yen
,
B. C.
1965
Characteristics of Subcritical Flow in a Meandering Channel
.
Iowa Institute of Hydraulic Research
,
Iowa City
.
Yen
,
C. L.
1967
Bed Configuration and Characteristics of Subcritical Flow in a Meandering Channel
.
University of Iowa
,
Iowa City
.