Revisiting Waal River Training by Historical Reconstruction

The frequency of dike breaches and flooding increased along the Waal River, a branch of the Rhine, in the early modern period. At least 29 out of 51 river dike breaches between 1757 and 1926 were provoked by ice jams (Wijbenga et al. 1993, 1994). As these ice jams formed preferentially at wider shallow sections with bars and islands, it was decided to give the river a narrow uniform width using a series of transverse groynes. For this, since 1850, the width of the Waal River was reduced to 360 m, 310 m, and finally to 260 m. To improve navigability, in the second decade of the twentieth century, the river was also straightened and dredged. Fig. 1 shows the different phases of the river training program, called Normalization.

The morphological response of the river to narrowing, straightening, and sediment extraction resulted in an incision process that reached the rate of $2\mathrm{cm}/\mathrm{year}$ (Sieben 2009). Besides floodplain draining and deterioration of riparian habitats, the lowering of the channel bed also posed problems to the hydraulic structures along the river, including the groynes. The conditions for navigation gradually deteriorated because the connection between the river and the other elements of the inland waterway network, such as canals, ship locks, and fluvial port facilities, became problematic. In addition, transverse groynes produce shallow ridges as a response to the formation of scour holes at their heads, which hampers navigation even further (Yossef 2005). Finally, a series of groynes hinders the flow at high discharges, forcing the water level to increase (Huthoff et al. 2013).

In 2013–2015, three series of transverse groynes were replaced by longitudinal training walls near the city of Tiel (the Netherlands), as an attempt to mitigate these effects (Fig. 2). The main idea of this pilot project was to create a narrow and sufficiently deep navigation channel at low flows and a two-channel system to convey water at high flows, reducing the channel incision process and lowering water levels. The new channel behind the training wall, a freely flowing stream without intensive shipping traffic and with a more natural bank, was expected to improve the river ecological conditions. However, the effectiveness of these interventions in achieving the goals has not been thoroughly investigated yet, particularly on the long-term.

Navigation is common in low-land rivers, like the Waal, especially in the final part of their course, close to the sea. The bed of these rivers often presents bars and point bars inside bends. These large sediment deposits are characterized by transverse bed slopes and produce helical flows, with the result that they locally modify the sediment transport direction. The presence of a nearby bar is thus expected to alter the distribution of sediment at the bifurcation created by a longitudinal training wall and hence to influence bifurcation stability (Bertoldi and Tubino 2007). Le et al. (2018a, b) carried out a series of experimental tests as well as numerical simulations showing that the parallel-channel system created by a longitudinal training wall tends to be unstable. The most stable configuration is obtained with equally-wide parallel channels. Both channels remain open, one deeper and one shallower, if the upstream end of the wall is located near the top of a steady alternate-bar or of a point-bar.

This work aims to assess the effectiveness of training the Waal River with a longitudinal training wall instead of groynes with a numerical model, starting from the situation of the river before the Normalization works when the river channel was almost twice as wide as now. The large natural width of the historical Waal allows splitting the river in two equally-wide channels, each of them being large enough to potentially become a deep inland water way. The other channel, expected to become much shallower and with a natural bank, seems suitable to eventually host a richer riverine ecosystem than the current groyne fields. For this reason, the channel that becomes deeper is, in this study, termed the navigation channel and the one that becomes shallower the ecological channel. The effectiveness and convenience of training the river with a longitudinal wall are assessed by comparing the long-term morphological developments obtained with the wall to the ones obtained with series of groynes along both banks and to the natural development of the river (untrained case). The focus is on channel incision, channel navigability, and flood conveyance.

The long-term morphological responses are reproduced with a two-dimensional depth-averaged (2DH) model based on the open-source Delft3D code. The numerical approach adopted in this study is supported by Le et al. (2018b), who successfully reproduced the long-term morphological developments observed in an experimental channel separated by a longitudinal training wall; by McCoy et al. (2008), who modeled a channel with a series of groynes; by Arboleda et al. (2010), who simulated floodplain sedimentation in the historical Waal River; and by Vargas-Luna et al. (2018), who successfully applied Baptist et al.’s (2003) formulation, implemented in Delft3D, for the assessment of the effects of floodplain vegetation on water levels.

The Waal River is the largest branch of the Rhine in the Netherlands, crossing the country from East to West. The study area is a 12-km long reach of the river located near Slijk-Ewijk, 6 km downstream of the city of Nijmegen (Fig. 3), before the Normalization works, when the river was twice as large as it is now (1842 AD—Latin for Anno Domini, which means in the year of the Lord, indicating how many years have passed since the birth of Jesus). The site was selected due to the availability of historical data (Schoor et al. 1999; Arboleda et al. 2010).

Data on historical bed topographies, converted from local water depths into bed levels relative to the current sea reference level are provided by Rijkswaterstaat (RWS)—Ministry of Infrastructure and Water Management of the Netherlands. An old river map of 1800 AD showing the presence of three bars [Fig. 3(a)] and a set of 17 reconstructed cross-sectional profiles (Maas et al. 1997) allow for model calibration. Another map of the same area in 1842 AD, with five reconstructed cross-sectional profiles [Fig. 3(b)], allows for model validation and provides the start conditions for the simulations of the long-term morphological developments of the river trained either with a longitudinal wall, with groynes, or without any training structures. Both historical maps fit the current situation well, showing their reliability [Fig. 3(c)]. Note that the current river channel alignment is straighter than the original one.

In the study area, the prenormalization channel width was approximately 400 m. By hypothetically dividing it in two equally-wide parts by means of a longitudinal wall, two parallel 200 m wide channels are obtained. This width appears sufficient for navigation, considering that the current river width is 230–260 m.

The historical daily discharge time series were reconstructed by Arboleda et al. (2010) based on the work of van Vuren (2005) and refer to the period 1790–1810. From these time series, the historical daily mean discharge results $\mathrm{1,550}{\mathrm{m}}^3/\mathrm{s}$, which is rather similar to the current value of $\mathrm{1,490}{\mathrm{m}}^3/\mathrm{s}$, was obtained from the time series of the period 2000–2010 (M. Schoor, personal communication, 2009–2018). The estimated historical bankfull discharge in the study area is $\mathrm{2,350}{\mathrm{m}}^3/\mathrm{s}$. The current value is approximately $300{\mathrm{m}}^3/\mathrm{s}$ higher (it depends on location) due to the channel incision process, which has followed the Normalization works.

In the early 1800s, the median grain size of bed material was 850–1,000 μm (Maas et al. 1997), slightly smaller than the present one (Frings and Kleinhans 2008; Frings et al. 2009). Arboleda et al. (2010) showed that suspended solids concentrations were similar to the actual ones, ranging between 85 and $95\mathrm{mg}/\mathrm{L}$, which is the composition of suspended sediment being silt and fine sand (Middelkoop 1997).

Three wall scenarios are analyzed. In the first two, the longitudinal wall subdivides the river in two equally-wide channels (200 m), as recommended by Le et al. (2018a, b), because this subdivision leads to a rather stable configuration in which both channels tend to remain open. Le et al. (2018b) also showed that starting the wall upstream of, or near to, a point bar top leads to sedimentation in the channel located at the bar side, whereas starting the wall much more downstream, close to or in the subsequent pool, leads to sedimentation in the channel at the opposite side. To cover these two opposite morphological developments, two cases are considered: a wall starting 500 m upstream of point bar I top (Fig. 3) and a wall starting 2,500 m downstream of the same bar top. In the third scenario, the wall, starting 500 m upstream of point bar I top, subdivides the river in two parallel channels with the eroding one having a width equal to the current navigation channel (260 m).

Two groyne-scenarios are considered: (1) a channel width between groynes heads (river-flow width) equal to half of historical Waal River width (200 m); and (2) a channel width between groynes heads equal to the current Waal River width (260 m).

The considered temporal scale corresponds to the time needed to achieve morphodynamic equilibrium in which the long-term temporal changes of two-dimensional (2D) and one-dimensional (1D) bed topography (bar configuration and longitudinal slope) become negligible.

The morphodynamic models used for the investigation are built from Arboleda et al.’s (2010) setup, who reconstructed the historical Waal River using the open-source Delft3D code. However, the present models include some important modifications, as for instance having a mobile bed instead of a fixed bed.

Delft3D solves the Reynolds equations for incompressible fluid and shallow water in three dimensions with a finite-difference scheme (Deltares 2014). The bed level changes are either derived by applying the Exner approach, i.e., by assuming immediate response of the sediment transport rate to the local hydraulic conditions, which is valid only for bed load, or by applying a sediment balance equation in which sediment deposition and entrainment rates are obtained from suspended solids concentrations, hydrodynamic parameters, and soil characteristics. In the Exner approach, the sediment transport rates are computed by means of sediment transport capacity formulas. In the sediment balance approach, the evolution of suspended solids concentration is computed with an advection-diffusion equation. In this study, sediment is assumed to be transported as a bed load. The effects of a transverse slope and spiral flow on the bed load direction are important for the simulation of the evolution of the river bed topography (Mosselman and Le 2016), in particular for the reproduction of bars (e.g., Schuurman et al. 2013), and to simulate the sediment distribution between bifurcating channels (e.g., Kleinhans et al. 2008). For the combined effect of a transverse slope and spiral flow, Struiksma et al. (1985) derived the direction ${\beta }_s$ of sediment transport (bed load) aswhere $\delta z/\delta y$ = transverse slope; $\delta z/\delta x$ = streamwise slope; ${\beta }_{\tau }$ = bed shear stress; and $f(\theta )$ = function to account for the effects of the transverse bed slope on the bed load direction in which $\theta$ represents the Shields parameter. Considering the effects of spiral flow, which become relevant if the stream lines have a sinuous path, as in channel bends and around bars, the bed shear stress (${\beta }_{\tau }$) is derived as follows (Struiksma et al. 1985)where $u$ and $v$ = streamwise and transverse velocity, respectively ($\mathrm{m}/\mathrm{s}$); $h$ = local water depth (m); $R$ = streamline radius of curvature (m); and $A$ = coefficient that weights the influence of the spiral flow, derived as follows (Struiksma et al. 1985)with $\epsilon$ ($\approx 1$) being a calibration coefficient (in this study $\epsilon =1$); $\kappa$ ($=0.4$) being the Von Kármán constant; and $g$ and $C$ being the acceleration due to gravity ($\mathrm{m}/{\mathrm{s}}^2$) and the Chézy roughness coefficient (${\mathrm{m}}^{1/2}/\mathrm{s}$), respectively.

The function $f(\theta )$ is derived by applying the formulation suggested by Koch and Flokstra (1981), extended by Talmon et al. (1995)where $D_{50}$ = median sediment size; and $A_{sh}$, $B_{sh}$, and $C_{sh}$ = calibration coefficients, with $A_{sh}$ being the most sensitive one with a value that normally falls between 0.1 and 1.0. The effects of a transverse slope on the sediment transport direction decrease if the value of $A_{sh}$ increases (Schuurman et al. 2013) and vice versa. In this study, the value of $A_{sh}$ is calibrated, whereas the other coefficients are kept constant with the values suggested by Talmon et al. (1995): $B_{sh}=0.5$ and $C_{sh}=0.3$.

In the morphodynamic model, the drying/wetting of cells is obtained based on the local water depth. The cells that become dry are removed from the computational domain and remain dry until higher discharge floods them again. In this study, the threshold value below which cells become dry is 10 cm.

In the presence of vegetation, as for instance on floodplains during high flow events, the Delft3D code corrects the local Chézy coefficient applying Baptist et al.’s (2003) approachwhere $C_b$ = bed roughness without vegetation (${\mathrm{m}}^{1/2}/\mathrm{s}$); $C_D$ = drag coefficient associated with the vegetation type, for which Vargas-Luna et al. (2015) suggest using the value of 1, corresponding to rigid cylinders; $n$ = vegetation density (${\mathrm{m}}^{-1}$), being $n=mD$, where $m$ = number of cylinders per unit area (${\mathrm{m}}^{-2}$), $D$ = cylinder diameter (m); and $h_v$ = vegetation height (m).

The following simulations are performed:

To compare the navigability and the high flow conveyance of the river in the different scenarios, the models are used again (this time with a fixed bed) to simulate the hydraulic conditions at selected values of the discharge (hydraulic models) (Castro-Bolinaga and Diplas 2014). The fixed-bed topographies are the final equilibrium ones resulting from the application of the morphodynamic model. For the Waal River, the discharge $Q_{\mathrm{low}}=680{\mathrm{m}}^3/\mathrm{s}$ is currently used as the lowest value for which it is necessary to guarantee a minimum navigation depth of 2.80 m for a channel width of 150 m. The discharge of $Q_{\mathrm{design}}=\mathrm{10,600}{\mathrm{m}}^3/\mathrm{s}$ is strongly recommended for the design of flood protection works. Another important reference discharge corresponds to the 1995 flood: $\mathrm{8,000}{\mathrm{m}}^3/\mathrm{s}$ (M. Schoor, personal communication, 2009–2018). These are the discharges considered in the hydraulic computations. The flow conveyance of the river is inferred by comparison of computed water levels at the same high discharge. Lower water levels indicate increased river conveyance whereas higher levels indicate decreased conveyance. The water levels are computed at a specific reference cross-section in the middle of the model domain.

The modeled river reach is about 12 km long and 1 km wide on average, with the width of the main channel (the historical Waal River) being approximately 400 m. The curvilinear grid follows the alignment of the historical channel with the embankments being used as land boundaries, with grid-cells of $40\times 20\mathrm{m}$ in the channel and $40\times 40\mathrm{m}$ in the floodplains. This configuration allows having a sufficient amount of data points along the cross section with a reasonable computational time. The morphological computations are not accelerated (Carraro et al. 2018).

For the calibration run, the initial main channel bed and the floodplains are horizontal in the transverse direction but sloping in the longitudinal direction, with slope $i_0=1\times {10}^{-4}$ and the floodplains being 5.5 m higher than the channel bed. For the validation run, the initial bed topography is the one derived from the 1800 map (Fig. 3). For the model simulations, the initial bed topography is the one derived from the 1842 map (Fig. 3).

The upstream hydraulic boundary condition is flow discharge. A representative yearly hydrograph [Fig. 4(b)] is derived from the mean duration curve of the period 1790–1810 [Fig. 4(a)], obtained from the data provided by Arboleda et al. (2010). The total volume of transported sediment computed using the representative yearly hydrograph [Fig. 4(b)] remains unchanged compared to the one obtained using the mean yearly duration curve [dotted red line in Fig. 4(a)]. The bed material, which is assumed uniform, has a diameter of 900 μm.

The representative yearly hydrograph is repeated every year for the entire duration of the computational period. The approach minimizes the problems related to model instability because it decreases the sudden inflow changes that are present in the reconstructed daily discharge time series while maintaining the same yearly sediment transport rates and almost the same yearly discharge distribution. The same approach has been already successfully used for long term morphological simulations by, for instance, Yossef et al. (2008). The downstream hydraulic boundary condition is a flow-dependent water level. In this study, the water level at the end of the model domain is assumed to coincide to the one corresponding to uniform flow.

The upstream boundary condition for the sediment component is a balanced sediment transport (recirculation: input = output). The downstream boundary condition is given by the free sediment transport condition, which allows undisturbed bed level changes up to the end of the model domain. The bed-load transport rate is computed with the Engelund and Hansen (1967) formula, valid for sand-bed rivers. The unvegetated bed roughness is represented by a constant Chézy coefficient, whereas the value of the Chézy coefficient of vegetated areas is computed by applying Eq. (5).

The distribution of floodplain vegetation cover is derived from the reconstructed ecotope maps of 1830 AD provided by Maas et al. (1997) (Fig. 5), and the vegetation parameters used in the model are the ones suggested by van Velzen et al. (2003) (Table 1). The computations are carried out, assuming that floodplain vegetation has not changed much in the 50 years covered by the simulations. The historical Waal River characteristics at an average discharge are listed in Table 2.

The results of the preruns show that a time step of 30 s ensures numerical stability and that after 50 years, the system reaches a state of morphodynamic equilibrium in which the longitudinal slope remains constant and cross-sectional changes of bed elevation are negligible ($T=50\mathrm{years}$). The computational time of a 50-year simulation is 5–6 days.

Fig. 6 shows the results of model calibration. The computed bed topography presents three alternate bars in the main channel [Fig. 6(b)]. These bars are present and have a similar size in the reconstructed bed topography [Fig. 6(a)]. Fig. 6(c) shows both the reconstructed cross-section A-A and the computed one for different values of the calibration parameter. The best result is obtained for $A_{sh}=0.4$. The results confirm what was reported by Le et al. (2018a): the larger $A_{sh}$ (smaller bed slope effects) leads to larger bar amplitudes, shorter bar wavelengths, and smaller longitudinal slopes. The value of $A_{sh}=0.4$ results in the longitudinal slope $i_{\mathit{1800}}=1.03\times {10}^{-4}$, which is equal to the reconstructed longitudinal bed slope.

The calibrated model is then validated on its ability to reproduce the reconstructed bed topography of 1842 AD starting from the reconstructed bed topography of 1800 AD. Fig. 7 shows the results of the model validation. The model reproduces the chute cutoff that occurred in the period 1800–1842, resulting in the formation of a middle bar [Fig. 7(b)], which is similar to the one that can be observed from the reconstructed bed topography of 1842 [Fig. 7(a)]. Computed and reconstructed bars have similar lengths, but the model tends to underestimate the chute cutoff channel depth [Fig. 7(c)]. The computed longitudinal bed slope is $i_{\mathit{1842}}=1.03\times {10}^{-4}$, which equals the reconstructed one. The satisfactory results of validation allow trusting the model in its ability to reproduce realistic bed topographies in the study area.

The starting condition of all scenarios is the reconstructed bed topography of 1842 AD, improved by integrating the measured cross-sections with the computed bed topography. The computations, covering the period of 50 years, produce the hypothetical equilibrium situation at the end of 1892 AD, which is obtained without any morphological acceleration. The discharge hydrograph used as an upstream boundary condition is the schematized yearly hydrograph 1790–1810 (Fig. 4).

The computed 1892 bed topography of the base-case scenario (Fig. 8) shows that the channel to the left of the middle bar tends to silt up, whereas the channel to the right has become the main stream. Because the model does not simulate bank erosion, the main stream narrows to 120 m. This result is supported by the later performed historical works, belonging to the Normalization program (Fig. 1). Due to its decreased size since 1872, the left channel was finally abandoned after the construction of groynes (in 1888 and 1890, see Fig. 1), whereas the right channel was artificially straightened and enlarged to 310 m by removing the middle bar. The averaged longitudinal bed slope at the end of the 50-year evolution is the same as the 1842 slope: $i_{\mathrm{untrain}}=1.03\times {10}^{-4}$.

The results of the two wall scenarios in which the parallel channels have a width of 200 m are shown in Fig. 9. The final system always presents a deeper and a shallower channel. If the training wall starts upstream of point bar I top, the left channel gradually silts up and the right channel deepens [Fig. 9(a)]. Instead, if the wall starts downstream of the same point bar top, the left channel deepens, and the right channel silts up [Fig. 9(b)]. These results confirm the experimental and computational findings of Le et al. (2018a, b). Compared to the base-case scenario, the averaged bed level of the navigation channel lowers by either 3.57 m or 3.66 m for the starting locations at the upstream or downstream of point bar I top, respectively. The bed slopes become milder than in the base-case scenario: $0.91\times {10}^{-4}$ and $0.89\times {10}^{-4}$, respectively, compared to $i_{\mathrm{untrain}}=1.03\times {10}^{-4}$. The bed level of the ecological channel raises by either 3.38 m or 3.76 m, respectively, and its bed slope becomes steeper than the one of the untrained river (base-case scenario): $3.80\times {10}^{-4}$ and $3.77\times {10}^{-4}$, respectively.

In the two groyne scenarios, the channel between groyne heads has either the width of the current Waal River (260 m) or of 200 m. Fig. 10 shows the computed bed topographies. In this figure, the river channel has clearly deepened, as described by Suzuki et al. (1987) and Spannring (1999). Compared to the untrained river scenario (base-case), the averaged bed level lowers by either 2.33 or 4.04 m for the river width of 260 or 200 m, respectively. The bed slope becomes milder: $0.93\times {10}^{-4}$ and $0.87\times {10}^{-4}$, respectively, compared to $i_{\mathrm{untrain}}=1.03\times {10}^{-4}$. The results agree with the general trends that can be expected after river narrowing (Duró et al. 2015).

Jansen et al. (1979) and Duró et al. (2015) show that channel narrowing inevitably leads to river incision in the narrowed reach and upstream. In the narrowed reach, incision is caused by channel deepening and by a longitudinal bed slope reduction. A longitudinal slope reduction lowers the bed and water levels at the upper end of the narrowed reach. This triggers an erosion process propagating upstream in which the final bed level lowering is theoretically equal to the slope reduction multiplied by the length of the narrowed reach. This means that the training setup that results in the smallest longitudinal slope has the highest effect upstream.

Table 3 lists the values of average channel deepening and final equilibrium longitudinal slopes of all considered scenarios. In general, the navigation channel obtained with a longitudinal wall presents less deepening compared to an equally-wide navigation channel obtained with groynes. However, the difference appears negligible for the width of 260 m. Starting the wall just upstream of the point bar top slightly reduces the averaged riverbed lowering compared to starting the wall downstream.

The results show that the groyne scenario with a river width of 200 m presents the largest impact on the longitudinal slope, whereas the wall scenario with a navigation channel width of 260 m presents the smallest one. Even if the differences appear small, these might result in important differences in upstream river incision if the trained reach has a length of the order of 100 km, as it is for the River Waal in the Netherlands.

To check river navigability, the models are used to assess the water depths at the prescribed minimum discharge considering the computed bed topographies of 1892 AD (hydraulic models).

The results, shown in Fig. 11(a), indicate that in the base-case scenario, the river fulfills the criterion for navigation depth only in the right channel around the middle bar where the width is 120 m, which is smaller than the minimum requirement. In all the groyne scenarios [Figs. 11(d and e)], and in the wall scenarios subdividing the river in two 200 m-wide parallel channels [Figs. 11(b and c)], the minimum water depth of 2.8 m is exceeded everywhere. The results of the wall scenario with a navigation channel of 260 m are shown in Fig. 12. In this case, some bars form in local channel expansion areas, where the flow presents larger width-to-depth ratios. As a result, the navigation channel would need a flattening of bar tops or local width adjustments to maintain a depth of 2.8 m over the entire width. The formation of bars leading to local water depths smaller than 2.8 m does not occur in the groyne scenario with the same width [Fig. 11(d)] because of the absence of channel expansion areas. Nevertheless, bar tops flattening is regularly carried out in the current Waal River having the same width of 260 m (M. Schoor, personal communication, 2009–2018).

The hydraulic models are also used to compute the water levels at high flows. The effects of river training on water levels are represented in this study by $\mathrm{\Delta }{Z}_{\mathrm{water}}$. This is the difference between the average water level of the considered training scenario and the one of the base-case scenario at Cross-Section A’-A’, assumed as representative of the study river reach [Fig. 7(a)]. If $\mathrm{\Delta }{Z}_{\mathrm{water}}$ is negative, the training method results in lower water levels, which indicates increased high flow conveyance of the river and vice versa if it is positive. The computed water levels at Cross-Section A’-A’ are listed in Table 4.

The results show that the longitudinal training wall increases the high flow conveyance even compared to the untrained river. This is due to navigation channel deepening. The most appropriate starting location of the wall is upstream of the point bar top. Instead, training by means of a series of groynes decreases the high flow conveyance of the river, potentially creating problems during floods.