Improving river stage forecast by bed reconstruction in sinuous bends

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.

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 km², 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.

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.

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).

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.