3.1. Bank Retreat: The Middle Rio Grande
The first model application site is located at River Mile (RM) 111, downstream of the San Acacia Diversion Dam, on the middle Rio Grande, New Mexico. A series of historical aerial photographs were available between 2002 and 2012 and they delineated the bank retreat changes over the years. A five-year period, October 2005 to September 2010, is simulated for a combined vertical and lateral modeling. The river bathymetry in October 2005 is reconstructed from the cross section surveys conducted in July–August 2005 by Tetra Tech., Inc., the dry areas for bank top and floodplains are obtained from the 2003 DTM and 2012 LiDAR data. The SMS software is used for bathymetry development.
The practical stream modeling procedure is established in the process of this study and it is described. The modeling consists of two stages: hydraulic analysis and bank erosion modeling. Hydraulic analysis is a necessary step to check and verify the model development and setup and a key model input, the Manning’s roughness coefficient, is calibrated if necessary. The calibration is done by matching the model results with the measured stages. The flow analysis results are used as the initial condition for the subsequent bank erosion modeling. Hydraulic analysis also helps guide a proper selection of the model domain and the development of an appropriate 2D mesh. In general, the hydraulic analysis is carried out in the following steps: (a) solution domain selection; (b) mesh generation; (c) topography and flow roughness representations on the mesh; (d) model calibration if the Manning’s coefficient is unknown; and (e) model application. Bank erosion modeling is carried out after the hydraulic analysis. It follows the following step: (a) bed gradation representation on the mesh; (b) key mode input parameters; (c) model calibration; and (d) model application. It is recommended to use a very coarse mesh for initial model setup and tests. Once verified with the model, a final mesh is used for project applications. The above procedure is illustrated below with the Rio Grande modeling.
The model domain is dictated by the area of interest for the study questions and often limited by the availability of bathymetric and terrain data. The model domain for the present study is chosen as shown in
Figure 3a. The domain covers cross sections SA-1246 in the upstream to SA-1262 in the downstream and has a reach longitudinal length of 2900 m. The 2D mesh generated for the final simulation consists of 11,722 mixed quadrilateral and triangular cells (11,584 nodes) as shown in
Figure 3b. The terrain data reconstructed from the survey data, representing the October 2005 condition, is interpolated onto the 2D mesh using SMS. The 2D represented terrain is displayed in
Figure 3c.
The flow roughness and initial surface and subsurface sediment gradations are found to be the two major model inputs with the vertical model. In this study, the Manning’s coefficient with a value of 0.026 is used in the main channel and bare bars based on the previous studies of the same river [
5]—no calibration is needed. For the vertical erosion modeling, the initial surface and subsurface sediment gradations are specified in 13 zones corresponding to 13 survey points shown in
Figure 4a. Zones are delineated using the nearest-distance criterion and each zone has the same gradation as the measured data at the corresponding point. The bed gradation was surveyed by Bauer [
30] during the periods of 15–20 June 2006 and 17–19 July 2006. The measured gradations at the survey points are plotted in
Figure 4b,c. The bed gradation is one of the most important inputs for the vertical erosion simulation based on the previous experience in modeling many other streams.
A time-series flow hydrograph, from July 2005 to the end of 2010, is imposed as the upstream boundary condition based on the data at the USGS San Acacia Floodway station (#08354900) (
Figure 5a). The upstream sediment supply is based on the capacity equation of Parker [
31] as no measured sediment flux data is available. Discharges less than 28.32 cubic meter per second (m
3/s) are excluded from the modeling per suggestions by Lai [
32,
33]. Low flows do not change the channel morphology and their exclusion can speed up the simulation significantly. At the downstream boundary, the normal depth boundary condition is applied.
A total of seven sediment size classes are used, from 0.002 to 125 mm, which represent the size range of the bed material load (
Table 1). The vertical erosion model uses the Parker equation [
31] as the sediment transport capacity. The bedload adaptation length is chosen to be 51.82 m (about the average channel width of the study reach) and the active layer thickness is 30 mm which is about twice the mean particle diameter in the main channel. These parameters are recommended as the standard inputs by SRH-2D.
The final set of needed model inputs are associated with the bank properties for the lateral bank erosion modeling. For the study case, a section of the outer bank, on the river right, is simulated (see
Figure 5b) as this section experienced bank erosion. Initial bank toe and top lines are delineated and represented by two mesh lines. The polygon formed by them represents the bank zone to be simulated for lateral retreat. The bank zone consists of quadrilateral cells and has a total of five lateral nodes. Longitudinally, fifteen bank profiles are used; the toe and top nodes of each bank are displayed in
Figure 5b as filled dots. At each bank, a set of bank erosion specific parameters are needed as inputs; they are summarized below:
constant time step of 1 h is used for the bank erosion model which is compared with 5 s time step for the vertical model.
Bank material properties are needed as inputs: Bank porosity = 0.4, critical shear stress = 0.0 Pa; and volumetric compositions in terms of the seven size classes in
Table 1 = 13%, 29%, 29%, 29%, 0%, 0%, and 0%.
The erodibility coefficient is calibrated and the final values are: 2.0 × 10−5 for banks 2 to 10, 5.0 × 10−6 for banks 11 to 12, and 2.0 × 10−6 from banks 13 to 15.
The combined vertical and lateral simulation is carried out from October 2005 to September 2010. First, a constant flow discharge of 28.32 m3/s is imposed for the flow modeling (without invoking the sediment and bank erosion models). The flow results are used as the initial condition. It is also used as the cutoff discharge for the combined vertical and lateral modeling.
The predicted flow velocity and shear stress distributions are shown in
Figure 6. The results show that flow in the upstream section above the bend has higher velocity and shear stress than the section downstream due to narrower cross sections upstream. The outer bank of the RM 111 bend is subject only to moderate fluvial loading with a shear stress about 1.0 Pa. This is a bit surprising considering that this bend has been subject to bank erosion over the study period. At the study site, the hydraulic flow results may indicate possible bank erosion mechanisms. First, the outer bank must be made of weaker materials than those in the upstream; it can be owing to geological reasons and/or vegetation covers. A field trip by the author in February 2013 confirmed the above conjecture. It was found that banks at RM 111 composed of very fine, unconsolidated, easily erodible materials while the upstream section was covered mostly by mature bushes. Bank retreat occurred mostly along the bare section of the bend. Second, non-fluvial processes such as seepage and/or piping may contribute to the bank erosion but it cannot be confirmed during this study. If this is the case, rainfall events would also influence the erosion process, and the fluvial stress based models alone are insufficient to capture all processes. At best, the erodibility coefficient may be used to take the non-fluvial processes into consideration, but only to a limited extent.
The combined vertical and lateral simulation is carried out next and results are presented. Among several model inputs, the erodibility coefficient is found to be the major calibration parameter for the simulation. The predicted 2010 bankline, shown as red dots in
Figure 7 representing bank toe and top nodes, is compared with the data (blue solid line). In the figure, the January 2006 and summer 2010 historical aerial photographs are also displayed for comparison. It is seen that the calibrated model does a reasonably good job in capturing the bank retreat over the five-year simulation period. Up to 23.8 m of bank retreat distance is predicted at several banks and it is in agreement with the historical data.
It is noted that the calibrated erodibility coefficient varies along the outer bank with its value reduced towards the downstream. For example, the erodibility value is 2.0 × 10
−5 for banks 2 to 10, 5.0 × 10
−6 for banks 11 to 12, and 2.0 × 10
−6 from banks 13 to 15. The reduction along banks 11–15 may be justified as the banks in this section were reinforced by the presence of bushes as shown in
Figure 7a. The presence of bushes with roots deep in the bank was also confirmed during the field trip in February 2013 by the author.
The five-year modeling at this site shows that the present combined vertical and lateral model is applicable to practical streams. The model is robust and simple to apply for engineering applications. The moving mesh approach works well for the case and the model results compare well with the available data. Despite the success, the simulation reveals also areas that need future attention. First, field investigation is recommended for the bank erosion site so that the dominant bank erosion mechanisms may be identified. The present model is valid primarily for fluvial processes. If non-fluvial processes are present, the model would be less accurate, although the erodibility coefficient may be adjusted to take those effects into account to a certain extent. In addition, the terrain data in the main channel of the retreating banks are too scarce (
Figure 3a). The October 2005 terrain reconstructed from the cross section data is probably inaccurate. High uncertainty in terrain may cast doubts on the reasonableness of the calibrated erodibility coefficient. Erodibility coefficient is a physical parameter; however, it is possible that the calibrated value can contain non-physical components contributed by the terrain data uncertainty. A new bathymetric survey, with denser and more accurate measuring equipment, is recommended.
3.2. Bank Retreat: The Choshui River
The present model is next applied to simulate the vertical and lateral morphological processes of a Chosui River reach in Taiwan. Choshui River is the longest river in Taiwan with about 187 km in length, draining from an area of 3157 square kilometers watershed. Details about the river hydrology and geomorphic features may be found in Lai et al. [
15].
Following the modeling procedure described for the Rio Grande case, the model domain is selected first, as shown in
Figure 8a. The study reach is about 16.7 km long and an average of 1.9 km in width. It is located near the river mouth with an average slope of 1:900. The upstream boundary is located 2.9 km upstream of the Ziqiang bridge and the downstream boundary is 1.8 km downstream of the Xibin bridge. The study reach is challenging to simulate as it is of the multi-channel braided type and subject to significant morphological changes. In addition, both moving mesh and fixed mesh approaches are used for the modeling to evaluate the pros and cons of each. With the model domain determined, a 2D mesh is generated consisting of 16,698 mixed quadrilateral and triangular mesh cells. This initial mesh is used by both the moving mesh and fixed mesh approaches.
A three-year simulation is carried out from July 2004 to August 2007. The initial river bathymetry was based on the field measured bed elevation in April 2004 (black lines in
Figure 8a) and the final bathymetry used by the model is shown in
Figure 8b. Interpolation from the field data to the model is performed with the SMS software. The key model inputs are as follows. The flow resistance is represented by the Manning’s coefficient (
n) and a constant value of
n = 0.027, estimated using the field data, is adopted. The initial bed gradation is based on the measurement carried out in April 2004. The data showed that the bed had a medium diameter of 1.44~1.78 mm during the simulation period. Five data points within the study area are plotted in
Figure 9a. Five sediment size classes are used to represent the range of sediments in the river: (1) 0.0625–0.125 mm; (2) 0.125–0.25 mm; (3) 0.25–0.5 mm; (4) 0.5–2.0 mm; and (5) 2.0–32 mm. The Engelund–Hansen sediment capacity equation [
34] is used to compute the erosional rate potential as this equation is suitable for sandy rivers. At the upstream boundary, the hourly recorded discharge data are used as the boundary condition (
Figure 9b) while the sediment supply rate is estimated to be 75% of the capacity computed by the Engelund–Hansen equation. At the downstream boundary, the stage-discharge rating curve developed from the measured data is used.
Model inputs specific to the lateral bank erosion model are presented next. A section of the left bank is selected for bank retreat modeling along with vertical processes in the channels. The bank zone, represented by the 2D mesh, is shown in
Figure 10 in which bank toe and top bank lines are also marked. A total of 21 banks, shown as black circles, are used for bank retreat modeling. The bank erosion parameters are as follows: (a) a constant time step of one hour is used with the bank module (the time step for the vertical model is 5 s); and (b) all banks have a porosity of 0.3 (measured), the critical shear stress of 0.5 Pa (estimated from field data), erodibility of 3.0 × 10
−5 (calibrated), the angle of repose of 30° (estimated), and the volumetric sediment composition of 0.5, 0.35, 0.15, 0.0, and 0.0 for the five size classes (measured), respectively. The critical shear stress of 0.5 Pa is estimated based on the 0.045 critical Shields parameter and 0.65 mm diameter. Again, this modeling shows that the erodibility coefficient is the only major parameter that should be calibrated using the measured bank retreat data. For the present simulation, the same erodibility is used for all banks as the bank properties do not vary much from each other based on field observations. During the continuous dynamic simulation, the materials eroded from the bank is deposited to the nearby channels.
The combined vertical and lateral simulation is carried out from July 2004 to August 2007 (3-year) using both the moving mesh and fixed mesh approaches. A comparison of the predicted and measured net erosion and deposition depth over the three-year period is shown in
Figure 11 for both mesh approaches. Large mismatch is observed upstream of the Ziqiang bridge near the upstream boundary. It is expected as the upstream boundary is located within a bend and results near the boundary are less accurate. Immediately downstream of the Ziqiang bridge, the model correctly predicts the “straightening” trend of the channel. Both immediate bars downstream of the bridge are eroded but the corresponding opposite banks of the bars experience aggradation. The model predicts the observed erosion trend in the area marked as “Left_Upstream Zone” in
Figure 11c. The predicted erosion depth, however, is much less than the measured data and the extent is also not well predicted. The mismatch may be attributed to unaccounted flows: the area was impacted by the tributary flow nearby which is not incorporated by the numerical model due to lack of data. Bank erosion marked as “Right-Middle Zone” is not predicted by the model at all. Attempts have been made to simulate this zone by activating the lateral model along the banks. Such efforts are unsuccessful as it is found that there is very few water actually flowing towards this area. It is conjectured that the bank erosion in this zone is due probably to undocumented events. Possible causes include: (a) local peculiar topographic features that caused water directed towards the area but not captured by the model; (b) impact of the tributary flow on the opposite side of the bank which is not included by the model; and (c) bank failure due to non-fluvial processes. The bank marked as “Left-Downstream Zone” has experienced large and rapid bank retreat during the simulation period. Both moving and fixed mesh approaches encounter problems. Both models predict reasonable bank erosion but retreat distance is under-predicted. With the moving mesh approach, further increase of the erodibility coefficient to match the measured bank retreat distance is not possible; the model would diverge as the mesh representing the eroding bank has been distorted too much. With the fixed mesh approach, the simulated bank retreat distance is highly dependent on the mesh density. With the current mesh used, it is insufficient to represent each bank well leading to gross under-prediction of the bank erosion. The needed increase in mesh density results in prohibitively long computing time and it is abandoned during the study. Overall, the model predicts the overall erosion and deposition patterns, but the agreement is more qualitative than quantitative.
This application case demonstrates that the present vertical and lateral model can predict the gross erosion and deposition patterns and bank erosion processes although the river is subject to rapid and significant bank retreat. Agreement with the measured data is good quantitatively in some areas but only qualitatively in others. The uncertainty in bed gradation spatial distribution and lack of such data are partially responsible for the less satisfactory model predictions. In the current study, the bed gradation survey was carried out mostly at exposed bars which tend to underestimate the sediment size in the main channel. The modeling shows that even current state-of-the-art models are not yet able to make a quantitative prediction of bank erosions for braided or morphologically dynamic river systems such as the Chosui River. For such complex rivers, accurate measured data of bed gradation distribution is important. Without such data, the model may only be used for qualitative prediction of the future morphological trend and identification of potential erosion and deposition trouble spots.