Late Holocene channel pattern change from laterally stable to meandering – a palaeohydrological reconstruction

River channel patterns may alter due to changes in hydrological regime related to changes in climate and/or land cover. Such changes are well documented for transitions between meandering and braiding rivers, whereas channel pattern changes between laterally stable and meandering rivers are poorly documented and understood. We hypothesize that many low-energy meandering rivers had relatively low peak discharges and were laterally stable during most of the Holocene, when climate was relatively stable and human impact was limited. Our objectives in this work are to identify a Late Holocene channel pattern change for the low-energy Overijsselse Vecht river, to develop and apply a novel methodology to reconstruct discharge as a function of time following a stochastic approach, and to relate this channel pattern change to reconstructed hydrological changes. We established that the Overijsselse Vecht was laterally virtually stable throughout the Holocene until the Late Middle Ages, after which large meanders formed at lateral migration rates of about 2 m yr⁻¹. The lateral stability before the Late Middle Ages was constrained using a combination of coring information, ground-penetrating radar (GPR), radiocarbon (¹⁴C) dating, and optically stimulated luminescence (OSL) dating. We quantified bankfull palaeodischarge as a function of time based on channel dimensions that were reconstructed from the scroll bar sequence and channel cut-offs using coring information and GPR data, combined with chronological constraints from historical maps and OSL dating. We found that the bankfull discharge was significantly greater during the meandering phase compared to the laterally stable phase. Empirical channel and bar pattern models showed that this increase can explain the channel pattern change. The bankfull discharge increase likely reflects climate changes related to the Little Ice Age and/or land use changes in the catchment, in particular as a result of peat reclamation and exploitation.

1 Introduction

Channel patterns describe the planform of a river, which reflects the interaction of the river channel with its floodplain. Channel patterns are classically distinguished: laterally inactive channels consist of straight and sinuous stable planforms, whereas laterally active channels consist of meandering and braiding planforms (Leopold and Wolman, 1957; Nanson and Knighton, 1996). Flume experiments and field data have shown that the channel pattern depends on several variables (Kleinhans, 2010). The first is the potential specific stream power, which is the product of the channel-forming discharge and valley slope (Kleinhans and Van den Berg, 2011; Nanson and Croke, 1992). The second is the bank erodibility (Ferguson, 1987; Friedkin, 1945), which is determined by the presence of bedrock in the valley side (Turowski et al., 2008), the bank cohesiveness (Peakall et al., 2007), and vegetation (Gurnell, 2014; Millar, 2000). The third is the type and amount of sediment supply (Gibling and Davies, 2012; Nanson and Croke, 1992).

Channel patterns can change in response to environmental variations (Ferguson, 1987). Many examples of channel pattern changes from braiding to meandering and vice versa are known to be associated with glacial–interglacial oscillations (Vandenberghe, 1995, 2002). Studies on the last glacial–interglacial transition have especially shown the simultaneous occurrence of channel pattern changes with a changing climate (Kasse et al., 2016; Vandenberghe et al., 1994). Climate change affects the vegetation, sediment availability, and discharge regime and consequently the bank stability, sediment transport, and potential specific stream power, resulting in different channel patterns.

Within the Holocene, several examples of channel pattern changes are documented from braiding to meandering rivers and vice versa (Brewer and Lewin, 1998; Lewin et al., 1977; Passmore et al., 1993; Słowik, 2015). However, channel pattern changes between laterally stable and meandering rivers have rarely been reported (Lewin and Macklin, 2010), except where human intervention transforms meandering rivers into heavily regulated and laterally stable rivers by introducing weirs, dams, groynes, and bank protection measures (Hesselink et al., 2003; Hobo et al., 2014; Słowik, 2013; Surian and Rinaldi, 2003). The partial abandonment of former meandering valleys may also result in underfit, laterally stable rivers like the former Rhine branches in the Niers and Oude IJssel valley (Janssens et al., 2012; Kasse et al., 2005).

Many studies have reported increased fluvial activity (e.g. increased discharge, sediment transport and deposition, and bank erosion rates) in relation to human, environmental, and climatic pressures during the Holocene (e.g. Hoffmann et al., 2008; Lespez et al., 2015; Macklin et al., 2010; Notebaert et al., 2018; Notebaert and Verstraeten, 2010). An example of increased fluvial activity is known from the Pine Creek (Idaho, USA), where mining and deforestation combined with intensive grazing resulted in an increase in discharge and sediment input, followed by river widening and an increase in bank erosion (Kondolf et al., 2002). The reverse change has been observed in settings as a result of afforestation (Kondolf et al., 2002; Liébault and Piégay, 2001) or increases in riparian vegetation fixing the channel banks (Eekhout et al., 2014; Vargas-Luna et al., 2016). The increase in fluvial activity during the Holocene was corroborated by an extensive review of existing studies concerning sediment accumulation in west and central European river floodplains by Notebaert and Verstraeten (2010). They concluded that sedimentation rates increased during the Middle and Late Holocene due to environmental changes. However, it is unknown whether the channel pattern changed simultaneously with the floodplain because no Early Holocene channel deposits representing a stable phase were identified. De Moor et al. (2008) hypothesized that the Geul River in southern Netherlands may have been relatively laterally stable during the Early and Middle Holocene, while it was actively meandering during the past 2000 years. Most of the floodplain deposits from the laterally stable phase have not been preserved, but De Moor et al. (2008) were able to reconstruct the bankfull depth for both periods. They estimated the bankfull depth to be a factor of 2 to 3 higher during the Late Middle Ages compared to the Early and Middle Holocene, and related this change to human and climate impact.

We conjecture that the change from laterally stable to meandering has occurred in some rivers for which increased Holocene fluvial activity was reported. The fact that such changes were not reported in the literature may either mean that critical conditions for channel pattern change were not reached or that evidence of such transitions is poorly preserved or left unnoticed. Both laterally stable and meandering rivers may display sinuous planforms, but the geomorphic processes in both rivers are different. Laterally stable channels are rivers without meandering processes, i.e. helicoidal flows causing bar formation and bank erosion at a significant rate (Kleinhans and Van den Berg, 2011; Nanson and Knighton, 1996; Seminara, 2006). In fact, the bends and channel cut-offs in laterally stable rivers may be the result of random and local perturbations (e.g. falling trees, beavers, bank collapse after heavy rainfall, etc.) leading to very limited and local displacement of the channel. Meandering and laterally stable rivers should therefore be distinguished by their different patterns of bar and floodplain formation, rather than merely by planform (Candel et al., 2017; Kleinhans and Van den Berg, 2011). We suggest that identifying channel pattern changes requires more detailed historic accounts or a much higher resolution of subsurface data than usually gathered because palaeochannels of laterally stable channels poorly preserve in the fluvial archive of meandering channel belts (Van de Lageweg et al., 2016), except when they have been cut off by random and local perturbations prior to the meandering phase. Using numeric (e.g. Oorschot et al., 2016) or scaled (e.g. Van Dijk et al., 2012) river simulation models is problematic for testing these ideas because these have not yet been capable of reproducing channel pattern changes. This reflects the lack of understanding of river processes and patterns (Kleinhans, 2010) and the need to gather such information from field studies.

This research entails a case study focussing on a river in which lateral activity during the past 500 to 600 years caused spectacular meandering: the Overijsselse Vecht in the Netherlands (Fig. 1). Previous work on this system has identified a transition from braiding to meandering during the Late Glacial (Huisink, 2000), while the meandering pattern remained throughout the Holocene until the river was channelized between 1896 and 1914 CE (Huisink, 2000; Neefjes et al., 2011). However, Quik and Wallinga (2018) found that the meanders were relatively young, with the oldest scroll bars dating from ca. 1400 to 1500 CE, by reconstructing meander formation using a combination of optically stimulated luminescence (OSL) dating of scroll bars and planform reconstruction based on historical maps. No fluvial deposits were found dating from before this period, except from a Holocene palaeochannel (here referred to as “Palaeochannel Q”) in a ground-penetrating radar (GPR) profile recorded by Huisink (2000, p. 123) 13 km upstream near Hardenberg (Figs. 1b and 2). Palaeochannel Q is relatively small compared to the meandering channel, seems to lack scroll bars, and was already cut off on the historical map of 1720 CE. Therefore, it is questionable whether the Overijsselse Vecht meandered prior to ca. 1400 CE. Alternatively, the river changed from a laterally stable into a meandering river during the Late Middle Ages. Our aims are (1) to identify whether a channel pattern change has occurred from laterally stable to meandering by collecting and combining detailed subsurface and geochronological data of the river from the pronounced meandering phase and the preceding phase; (2) to develop a methodology to reconstruct bankfull discharge as a function of time using the scroll bar deposits and channel remnants as a geological archive of the former channel dimensions; (3) to test whether palaeohydrological changes may explain the potential channel pattern change; and (4) to elaborate on the potential causes for changes in the discharge and channel pattern.

Figure 1 Maps of the Overijsselse Vecht. (a) Map showing the location of the Overijsselse Vecht catchment and the location of the study site. (b) Digital elevation map (DEM; Actueel Hoogtebestand Nederland, 0.5×0.5 m) (Van Heerd and Van't Zand, 1999) of the downstream section of the Overijsselse Vecht river, indicating both study sites: Junnerkoeland and Prathoek. DEM of the Junnerkoeland bend (c) and Prathoek bend (d), including locations of cores, OSL samples by Quik and Wallinga (2018), the OSL and ¹⁴C samples from this study, the GPR transects, the grain size samples, and inflection points. The possible historical course of Palaeochannel X according to Maas (1995) is indicated. (e) Zoomed-in figure of Palaeochannel X. (f) Topographical military map (TMK) dating from 1851 CE (CC-BY Kadaster, 2018; Van der Linden, 1973), showing the Overijsselse Vecht during its meandering phase.

Figure 2 Interpretation by Huisink (2000) of subsurface strata from GPR data collected near Hardenberg 13 km upstream of Junnerkoeland (see location in Fig. 1b). Horizontal strata of cover sand deposits (A) on top of the channel deposits of an interpreted braiding system (B). A relatively small, symmetrical palaeochannel is present (C) within the Late Glacial deposits, hereafter referred to as “Palaeochannel Q”. Figure adapted after Huisink (2000).

2 Study area

The Overijsselse Vecht (Fig. 1) is a low-energy, sand-bed river flowing from Germany into the Netherlands, with an average annual discharge (Q_m) of 22.8 m³ s⁻¹ and a mean annual flood discharge (Q_maf) of 160 m³ s⁻¹ derived from the gauging station in Mariënberg for the period 1995 to 2015 (see location in Fig. 1b). The river has a length of 167 km, its catchment covers 3785 km² with the highest point +110 m above sea level (m a.s.l.), and a relatively uniform valley slope of $\mathrm{1.42}\times {\mathrm{10}}^{-\mathrm{4}}$ to $\mathrm{1.7}\times {\mathrm{10}}^{-\mathrm{4}}$ in the Dutch part of its trajectory (TAUW, 1992; Wolfert and Maas, 2007). The Overijsselse Vecht incised its current valley during the Late Glacial within fluvioperiglacial sands, locally covered by aeolian cover sands (Huisink, 2000; Ter Wee, 1966; Wolfert and Maas, 2007). During the Late Holocene, aeolian drift sands formed along the Overijsselse Vecht as a result of agricultural overexploitation (Van Beek and Groenewoudt, 2011). The Overijsselse Vecht was an actively meandering river until 1896, when weirs were constructed and parts of the river were channelized. The river was completely channelized after 1914 CE, with five weirs controlling the water levels. Recently, sinuous side channels bypassing the weirs have been created as part of river restoration aiming to restore past physical and ecological characteristics of the river.

At present the topography of the meandering phase is partly still intact in the floodplain (Maas, 1995). Wolfert and Maas (2007) reconstructed the pre-channelization planform from historical maps of 1720, 1850, and 1890 CE. Large differences in meander development and lateral migration rates were found between different river reaches. In particular in areas where non-cohesive aeolian sands formed the channel banks, large meanders formed and lateral migration reached rates up to 3 m yr⁻¹. In this research we will study two of the large meanders, named Prathoek and Junnerkoeland (Fig. 1), for which Quik and Wallinga (2018) reconstructed the scroll bar development using OSL dating in combination with historical maps.

Here we take advantage of the preservation of a palaeochannel (here referred to as “Palaeochannel X”) with comparable dimensions as Palaeochannel Q (Huisink, 2000, p. 123) preserved in the Junnerkoeland as a sharp bend (Fig. 1c). Maas (1995) interpreted Palaeochannel X to be connected to the oldest swale of the Junnerkoeland scroll bar deposits (Fig. 1c). Palaeochannel X, however, was likely abandoned before the scroll bar formation because large differences in dimensions exist between Palaeochannel X and the meander bend, but the well-preserved nature suggests that Palaeochannel X is relatively young. The small dimensions of both Palaeochannels X and Q would suggest that the river had comparatively less energy and may have been relatively laterally stable prior to the meandering phase.

3 Methods

3.1 Lithological description

Cores were performed in a transect perpendicular to the scroll bars of both meander bends (Fig. 1c, d). An additional transect was cored perpendicular to Palaeochannel X (Fig. 1e). In the case that the deposit consisted of peat we used a gouge auger (Ø: 3 cm); in the case of unsaturated sand we used an Edelman auger, and in the case of saturated sand we used a Van der Staay suction corer (Van de Meene et al., 1979). In total, 68 cores were performed to a maximum depth of 7.3 m. The surface elevation of each coring site was either determined using a GPS combined with a DEM (Van Heerd and Van't Zand, 1999) or with a global navigation satellite system (GNSS) device. A standard method was used to describe the sediment cores in 10 cm thick intervals using the Dutch texture classification scheme, which approximately matches the USDA terminology (Berendsen and Stouthamer, 2001; De Bakker and Schelling, 1966). The median sediment grain size (D₅₀, m) of non-organic, sandy samples was visually checked in the field by comparison with a sand ruler. Grain size analysis was used to estimate a D₅₀ for the entire scroll bar deposit (Sect. 3.3). In addition, the plant macro-remains, any visible bedding, and colour were described. The percentage of gravel (> 2 mm) was estimated in the field using sieves. The lithogenesis was inferred from the lithological properties, facies geometries, and DEM topography, distinguishing fluvial, fluvioperiglacial, cover sand, drift sand, and residual channel-fill deposits (Huisink, 2000; Ter Wee, 1966).

3.2 Ground-penetrating radar

Ground-penetrating radar (GPR) was used to reconstruct the channel dimensions of the scroll bars. GPR measurements were conducted with a pulseEKKO PRO 250 Hz with a SmartTow configuration. The GPR transects were placed along the centreline of the meander bends, perpendicular to the ridge and swale morphology (Fig. 1c, d). The electromagnetic-wave velocity was 0.060 m ns⁻¹, derived by using isolated reflector points (Neal, 2004; Van Heteren et al., 1998) and by comparing depths of recognizable layers with the coring data.

3.3 Grain size analysis

In total 33 samples for grain size analysis were taken from the scroll bar deposits and three samples were taken from Palaeochannel X. The samples of the scroll bar deposits were taken from each 0.5 m interval from the channel lag up to the swale surface at three locations in Junnerkoeland and two locations in Prathoek (Fig. 1c–e). The samples of Palaeochannel X were taken from three locations below the residual channel-fill from the former riverbed. Grain size samples were analysed in a laboratory with a Beckman Coulter LS230 laser particle sizer. This instrument has a measurement range of 0.1 to 2000 µm. Samples were sieved with a 2 mm sieve and prepared with HCl (1 M) and H₂O₂ (30 %). All data were processed using a Fraunhofer.rfd optical model because of the low clay–silt content (Agrawal et al., 1991). Finally, the average and standard deviation were calculated for both the scroll bar deposits and Palaeochannel X and used in the palaeohydrological calculations.

3.4 OSL dating

We used the modelled age–distance relationships determined by Quik and Wallinga (2018) in our calculations. Their obtained OSL ages from the scroll bar deposits were used as priors and combined with historical map data in a Bayesian sequence model using the OxCal software (Bronk Ramsey, 2009). For details on the method see Quik and Wallinga (2018). In this study, we took four additional samples for OSL dating on the inner and outer bank of Palaeochannel X. These samples were collected in an opaque PVC tube (Ø 4.5 cm) mounted on a hand auger, allowing for sampling without light exposure. The analysis in the laboratory followed the same procedure as in Quik and Wallinga (2018). The OSL age was determined at the Netherlands Centre for Luminescence dating, with equivalent doses measured on small aliquots of quartz using the SAR protocol (Murray and Wintle, 2003) and dose rates determined from activity concentrations measured using gamma-ray spectrometry. A bootstrapped version of the minimum age model (Cunningham and Wallinga, 2012) was used to derive the best estimate of the burial dose and deposition age. Given the limited amount of samples associated with Palaeochannel X and the absence of additional age constraints from historical maps, no Bayesian analysis was performed for these samples.

3.5 ¹⁴C dating

A sample was taken in the deepest part of Palaeochannel X, at the sand–peat interface, using a piston corer (Ø: 6 cm). Macro-remains and leaf fragments from terrestrial species were selected from 1 cm intervals in the laboratory using a light microscope. Samples were stored in diluted HCl (4 %) at 5 ^∘C. The sand content was measured for each interval to precisely determine the position of the sand–peat interface. Material with volumetric sand percentages lower than 10 % to 20 % was considered as peat (Bos et al., 2012). The macro-remains from the centimetre above this interface were selected for the ¹⁴C analysis, providing a terminus ante quem date for the abandonment of the channel. The ¹⁴C age was determined by accelerator mass spectrometry (AMS) at the Centre for Isotope Research (Groningen University). For calibration, the IntCal13 curve was used in the OxCal 4.2.4 software (Bronk Ramsey, 2009; Reimer et al., 2013).

3.6 Channel dimensions

The channel dimensions of Palaeochannel X were determined from the lithological cross section. The residual channel-fill was delineated along the sand–peat interface. Bankfull depth (H_bf) was defined from the bottom of the palaeochannel up to the first clear knick point on the bank, which was mapped with a GNSS device such that the width–depth ratio was minimal (Williams, 1986). Relative error of H_bf was assumed to be similar to the relative error of H_bf during the meandering phase (ca. 10 %) and used in the calculations (see details below) because both H_bf values were determined by using coring data. Additional dimensions were measured from the delineated channel, involving the bankfull width (W), cross-sectional area (A), and wetted perimeter (P). These channel dimensions were also measured for Palaeochannel Q from the GPR profile recorded by Huisink (2000, p. 123) (Fig. 2). Here we assumed a similar relative error of W, A, and P as was taken for H_bf.

Figure 3 Sketch of the cross-sectional flow area of a meandering channel used for the bankfull palaeodischarge calculations (Allen, 1965; Hobo, 2015; Leeder, 1973).

The river channel was assumed to have the channel dimensions as shown in Fig. 3 during the meandering phase. This sketch is based on Allen (1965), Leeder (1973), and Hobo (2015). The bankfull depth (H_bf) was estimated from the coring data taken from the bottom of the channel lag up to the surface elevation in the swales (Fig. 4). Small elevation differences were expected to result from local variation rather than real changes in H_bf, and therefore the average H_bf was calculated from the smoothed bottom and surface elevation. The standard deviation of H_bf was calculated from the actual variable bottom elevation over the length of the scroll bar. The transverse bed slope (α) of the inner bend was determined based on the GPR transects (Fig. 5), in which lateral accretion surfaces could be distinguished. The angle was measured on the steepest (middle) parts of the identified lateral accretion surfaces. The average and standard deviation of α were calculated and used in the calculations. The calculations of the channel dimensions follow from Fig. 3. The bankfull width (W, m) and cross-sectional area (A, m²) were determined by Eqs. (1) and (2).

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }W=\mathrm{1.5}{\displaystyle \frac{H_{\mathrm{bf}}}{\tan \left(\mathit{\alpha }\right)}}\text{(2)}& {\displaystyle }& {\displaystyle }A=W{H}_{\mathrm{avg}}\end{array}$$

H_bf is the bankfull depth (m), and $H_{\mathrm{avg}}=\frac{\mathrm{7}{H}_{\mathrm{bf}}}{\mathrm{12}}$ approximates the average water depth (m). The wetted perimeter (P, m) was calculated from the assumed channel geometry (Fig. 3) following Eq. (3).

$$\begin{array}{}\text{(3)}& P={\displaystyle \frac{H_{\mathrm{bf}}}{\sin \left(\mathit{\alpha }\right)}}+{\displaystyle \frac{W}{\mathrm{6}}}+\sqrt{H_{\mathrm{bf}}^{\mathrm{2}}+{\left({\displaystyle \frac{W}{\mathrm{6}}}\right)}^{\mathrm{2}}}\end{array}$$

The hydraulic radius (R, m) was calculated by Eq. (4).

$$\begin{array}{}\text{(4)}& R={\displaystyle \frac{A}{P}}\end{array}$$

For each swale visible on the DEM the sinuosity (s, -), radius of curvature (R_curv, m) and scroll bar surface area (SB_surf, m²) were measured. The former channel sinuosity was estimated by the use of the DEM, measuring the distance along the swales relative to the distance along the valley between the inflection points (Fig. 1c, d). The sinuosity of Palaeochannel X was measured using the same approach (Fig. 1c). The channel slope (S_c, –) was calculated from the sinuosity and valley slope (S_v, –) determined by TAUW (1992) and Wolfert and Maas (2007) following Eq. (5).

$$\begin{array}{}\text{(5)}& S_{\mathrm{c}}=S_{\mathrm{v}}/s\end{array}$$

The volumetric rate of scroll bar growth (SB_vol, m³ yr⁻¹) was determined from scroll bar surface area (SB_surf, m² yr⁻¹) and the thickness between each swale and interpolated time interval following Eq. (6):

$$\begin{array}{}\text{(6)}& {\mathrm{SB}}_{\mathrm{vol}}={\displaystyle \frac{{\mathrm{SB}}_{\mathrm{surf}}{H}_{\mathrm{bf}}(\mathrm{1}-\mathit{\phi })}{\mathrm{\Delta }\mathrm{age}}}\end{array}$$

where φ is the porosity (here 0.3 to 0.35 volume fraction) (Nimmo, 2004), which was included to compare the SB_vol with the sediment transport, and Δage is the age difference between the scroll bars (yr) based on the datings by Quik and Wallinga (2018). Although scroll bar deposits were absent, following Eq. (6) we also calculated the volumetric sediment transport for the fluvial deposits on the inside of Palaeochannel X.

Figure 4 Stratigraphic cross sections of the study sites (for location see Fig. 1). Lithological cross sections of Junnerkoeland (a) and Prathoek (b). Lithogenetic cross sections of Junnerkoeland (c) and Prathoek (d) including the OSL samples by Quik and Wallinga (2018) and OSL and ¹⁴C dating results from this study. The surface and erosive base elevation are indicated with dashed lines, resulting in the inferred bankfull channel depth (H_bf). (e) Zoomed-in lithogenic cross section of Palaeochannel X. The thick dashed line indicates the bankfull level of the palaeochannel.

3.7 Palaeodischarge

The channel dimensions were used to calculate the bankfull discharge (Q_bf, m³ s⁻¹). Bankfull discharge is commonly considered an approximation of the channel-forming discharge with a recurrence interval of 1 to 2 years (Dury, 1973; Wolman and Miller, 1960). We assumed that the bankfull discharge was similar for both Junnerkoeland and Prathoek, regarding the short distance between these river sections (Fig. 1b). Hence the bankfull discharge was presented by combining the bankfull discharges for both meander bends. The bankfull discharge was estimated by applying the Chézy equation, following Eq. (7):

$$\begin{array}{}\text{(7)}& Q_{\mathrm{bf}}=CA\sqrt{R{S}_{\mathrm{c}}},\end{array}$$

where Q_bf is the bankfull discharge (m³ s⁻¹), and C is the Chézy coefficient (m^0.5 s⁻¹). The Chézy coefficient, i.e. flow resistance, was estimated by substituting Eq. (8) in Eq. (7). Equation (8) is an empirical relation (Brownlie, 1983):

$$\begin{array}{}\text{(8)}& {\displaystyle }{Q}_{\mathrm{bf}}={\left({\displaystyle \frac{R}{\left(\mathrm{0.3724}{S}_{\mathrm{c}}^{-\mathrm{0.2542}}{\mathit{\sigma }}_{\mathrm{s}}^{\mathrm{0.105}}{D}_{\mathrm{50}}\right)}}\right)}^{\mathrm{1.529}}W{g}^{\mathrm{0.5}}{D}_{\mathrm{50}}^{\mathrm{1.5}},\end{array}$$

where σ_s is the sorting of the bed material grain size (–) derived from the grain size analysis (Sect. 3.3) and approximated by 0.5($\frac{D_{\mathrm{50}}}{D_{\mathrm{16}}}+\frac{D_{\mathrm{84}}}{D_{\mathrm{50}}})$, D₁₆ and D₈₄ are the 16th and 84th percentile sediment grain size (m), respectively, and g is the gravitational acceleration (m² s⁻¹). As a validation, the calculated Chézy coefficient was compared with average Chézy coefficients of 12 comparable low-energy, sand-bed rivers with scroll bars (S_v < 0.001, 90 < Q_bf < 320 m³ s⁻¹) calculated from a large river dataset (Kleinhans and Van den Berg, 2011; Van den Berg, 1995). The cross-sectional-averaged flow velocity (u_bf, m s⁻¹) was calculated by following Eq. (9).

$$\begin{array}{}\text{(9)}& u_{\mathrm{bf}}={\displaystyle \frac{Q_{\mathrm{bf}}}{A}}\end{array}$$

3.8 Sediment transport

The sediment transport was calculated to compare with the SB_vol (Eq. 6). Sediment transport was calculated in two different ways. The first method was the slightly modified Engelund and Hansen (1967) relation following Eq. (10):

$$\begin{array}{}\text{(10)}& Q_{\mathrm{s},\mathrm{bf}}={\displaystyle \frac{\mathrm{0.05}{u}^{\mathrm{5}}Wti}{(\frac{{\mathit{\rho }}_{\mathrm{s}}}{\mathit{\rho }}-\mathrm{1}{)}^{\mathrm{2}}{g}^{\mathrm{0.5}}{D}_{\mathrm{50}}{C}^{\mathrm{3}}(\mathrm{1}-\mathit{\phi })}},\end{array}$$

where Q_s,bf is the yearly sediment transport derived from the bankfull discharge (m³ yr⁻¹), t is the number of seconds in a year, i is the intermittency assumed to be 0.03 to 0.07 (Parker, 2008), ρ_s is the sediment density (kg m⁻³), ρ is the water density (kg m⁻³), and φ is the porosity assumed to be 0.3 to 0.35 (Nimmo, 2004). The relation of Engelund and Hansen was used because the relation is suitable for sand-bed rivers with relatively low flow velocities (Van den Berg and Van Gelder, 1993), and the input variables required were available.

In the second method the sediment transport was determined for each discharge magnitude and related frequency (Q_s,freq) (Wolman and Miller, 1960) from present-day flow conditions by assuming that the current discharge frequency distribution also applied to the meandering phase. We used the hourly discharge data from 1995 to 2015 of the gauging station in Mariënberg (Fig. 1b). This gauging station is close to the study location and has the lowest amount of data gaps compared to the other stations. The flow duration was calculated for intervals of 10 m³ s⁻¹, and for each discharge interval the sediment transport was calculated using Eq. (10), excluding the intermittency factor. If the discharge was above bankfull, the flow would go across the floodplain. The Chézy coefficient for the floodplain was assumed to be half the Chézy coefficient in the channel because of the higher roughness of the floodplain compared to the channel. We assumed that the floodplain width was 350 m for the start of the meandering phase, which was estimated from the DEM (Fig. 1c), and that the width would increase proportionally with the lateral migration rate for each time step during the meandering phase.

3.9 Potential specific stream power

The potential specific stream power was calculated to plot into a stability diagram. Kleinhans and Van den Berg (2011) distinguished four different stability fields, further building on Van den Berg (1995) and Makaske et al. (2009): rivers with laterally stable channels, meandering rivers with scroll bars, meandering rivers with scroll and chute bars as well as moderately braided rivers, and braided rivers. In this research, only the first two stability fields are relevant. These stability fields are separated by a discriminator that represents the theoretical minimum energy needed for the channel pattern to occur (Kleinhans and Van den Berg, 2011). The potential specific stream power was calculated by applying the relationship presented by Kleinhans and Van den Berg (2011) following Eq. (11):

$$\begin{array}{}\text{(11)}& {\mathit{\omega }}_{\mathrm{pv}}={\displaystyle \frac{\mathit{\rho }g{\sqrt{Q}}_{\mathrm{bf}}{S}_{\mathrm{v}}}{\mathit{\epsilon }}},\end{array}$$

where ω_pv is the potential specific stream power (W m⁻²) and $\mathit{\epsilon }=\mathrm{4.7}\sqrt{\mathrm{s}{\mathrm{m}}^{-\mathrm{1}}}$ for sand-bed rivers (Van den Berg, 1995). The discriminator line was plotted by applying the relationships presented by Makaske et al. (2009) and Kleinhans and Van den Berg (2011) following Eq. (12):

$$\begin{array}{}\text{(12)}& {\mathit{\omega }}_{\mathrm{ia}}=\mathrm{90}{D}_{\mathrm{50}}^{\mathrm{0.42}},\end{array}$$

where the subscript ia refers to the discrimination between laterally stable and meandering channels with scroll bars.

3.10 Bar regime

Bar regime was predicted by applying the relationships of Struiksma et al. (1985) and Kleinhans and Van den Berg (2011). River bends can be seen as an example of a perturbation to both the flow and bed sediment, which have different adaptation lengths over which they return to equilibrium. This difference in response is expressed by the interaction parameter (IP), which is the ratio between the adaptation length of bed perturbation and the adaptation length of flow. The adaptation length of flow was calculated following Eq. (13),

$$\begin{array}{}\text{(13)}& {\mathit{\lambda }}_w={\displaystyle \frac{C^{\mathrm{2}}{H}_{\mathrm{avg}}}{\mathrm{2}g}},\end{array}$$

and the adaptation length of a bed perturbation (m) is calculated following Eq. (14):

$$\begin{array}{}\text{(14)}& {\mathit{\lambda }}_{\mathrm{s}}={\displaystyle \frac{H_{\mathrm{avg}}}{{\mathit{\pi }}^{\mathrm{2}}}}{\left({\displaystyle \frac{W}{H_{\mathrm{avg}}}}\right)}^{\mathrm{2}}f\left(\mathit{\theta }\right),\end{array}$$

where f(θ) is the magnitude of the transverse slope effect (–) calculated following Eq. (15) (Talmon et al., 1995).

$$\begin{array}{}\text{(15)}& f\left(\mathit{\theta }\right)=\mathrm{9}{\left({\displaystyle \frac{D_{\mathrm{50}}}{H_{\mathrm{avg}}}}\right)}^{\mathrm{0.3}}\sqrt{\mathit{\theta }}\end{array}$$

Here, θ is the dimensionless shear stress (–) calculated following Eq. (16):

$$\begin{array}{}\text{(16)}& \mathit{\theta }={\displaystyle \frac{\mathit{\tau }}{\left({\mathit{\rho }}_{\mathrm{s}}-\mathit{\rho }\right)g{D}_{\mathrm{50}}}},\end{array}$$

where τ is the shear stress (Pa) calculated following Eq. (17).

$$\begin{array}{}\text{(17)}& \mathit{\tau }=\mathit{\rho }gR{S}_{\mathrm{c}}\end{array}$$

The interaction parameter (IP, –) was calculated following Eq. (18) to determine the bar regime for the historical and prehistorical Overijsselse Vecht and for comparison with the theoretical thresholds of bar regime (Crosato and Mosselman, 2009; Struiksma et al., 1985) by

$$\begin{array}{}\text{(18)}& \mathrm{IP}={\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{s}}}{{\mathit{\lambda }}_{\mathrm{w}}}}.\end{array}$$

The IP is strongly related to the width–depth ratio and was therefore separately calculated for the meander bends Junnerkoeland and Prathoek. A low IP means that when a bar forms in response to a local perturbation, such as local curvature, the bar disappears within a short distance of the perturbation (Struiksma et al., 1985). This is called an overdamped regime and occurs in channels with a low width–depth ratio. The threshold between overdamped and underdamped can be calculated following Eq. (19):

$$\begin{array}{}\text{(19)}& \mathrm{IP}\le {\displaystyle \frac{\mathrm{2}}{n+\mathrm{1}+\mathrm{2}\sqrt{\mathrm{2}n-\mathrm{2}}}},\end{array}$$

where n is the degree of non-linearity of sediment transport versus depth-averaged flow velocity (–). Following Crosato and Mosselman (2009) we chose n=4, which corresponds to values for a sand-bed river. A higher IP, and hence a higher width–depth ratio, results in an underdamped regime associated with bars that also form further downstream of the perturbation. The thresholds can be calculated following Eq. (20).

$$\begin{array}{}\text{(20)}& {\displaystyle \frac{\mathrm{2}}{n+\mathrm{1}+\mathrm{2}\sqrt{\mathrm{2}n-\mathrm{2}}}}<\mathrm{IP}<{\displaystyle \frac{\mathrm{2}}{n-\mathrm{3}}}\end{array}$$

Figure 5 Example of a ground-penetrating radar (GPR) profile (250 Hz) in the Prathoek bend. (a) Original GPR profile and (b) interpreted GPR profile with lateral accretion surfaces and the channel lag, indicated by yellow lines.

3.11 Errors and uncertainty

The above described calculations (Eqs. 1 to 11 and 13 to 20) were run 10 000 times to take into account the random errors of the input parameters, following a stochastic approach by using Monte Carlo simulations. The uncertainty of these parameters was described above, relating to the transverse bed slope, bankfull depth of the meanders, valley slope, porosity, grain size, intermittency, and the measured channel dimensions of Palaeochannels X and Q. Systematic errors were not taken into account because the palaeohydrological reconstruction was meant to distinguish relative differences between fluvial phases, rather than reconstructing absolute hydrological parameters. Results are plotted with average values from the Monte Carlo simulations when normally distributed, or median values when not-normally distributed, including the 16th and 84th quantile representing the uncertainty margin. All formulas and example data used are made available in the Supplement.

4 Results

4.1 Lithogenetic units

Several lithogenic units were distinguished (Fig. 4), following similar interpretations of the sedimentary units as Huisink (2000). The descriptions of the lithogenic units are summarized in Table 1. The cover sand deposits were sometimes difficult to distinguish in borehole descriptions from the fluvioperiglacial deposits when the latter had a relatively fine grain size. Because our interest is in the delineation of the scroll bar and residual channel-fill deposits, we combined both the fluvioperiglacial and cover sand deposits into one unit. The fining upward sequence within the scroll bar deposits (Table 1) can be recognized in the grain size analysis done for the scroll bar deposits at Junnerkoeland and Prathoek (Fig. 6). The depth-averaged grain size for both scroll bar complexes is 0.28±0.05 mm. Commonly, at the base of the scroll bar deposits, a sharp transition occurs to the brightly coloured substratum of fluvioperiglacial deposits below, which lack organic material (Table 1). Cores that did not reach the fluvioperiglacial deposits below the scroll bar deposits indirectly indicate the boundary between these units because strongly consolidated layers are present in the fluvioperiglacial deposits that were difficult to core into. An example of a consolidated clay layer can be found directly below the southern part of the scroll bar deposits at Prathoek (Fig. 4b).

Table 1 Description of lithogenic units.

The GPR profiles clearly show the lateral accretion surfaces of the scroll bar deposits (see example in Fig. 5). The GPR results were poor only where the scroll bar deposits are relatively loamy or clayey on top (i.e. northern parts of Prathoek and Junnerkoeland). The bottom of the scroll bar deposits is mostly unrecognizable because of a low GPR reflection at this depth. In Fig. 5 the bottom of the scroll bars is visible because this part is located in the southern part of Prathoek where the above-mentioned clay layer was present (Fig. 4), which caused a strong reflection of the GPR signal. The well-preserved Palaeochannel X is a relatively symmetrical palaeochannel (Fig. 4e) similar to Palaeochannel Q of Huisink (2000) (Fig. 2). The outer bank consists of Weichselian–Early Holocene deposits (Fig. 4c). The average grain size of the Palaeochannel X bed sediments is 0.23±0.12 mm. No lateral accretion surfaces can be observed in the GPR profile that was placed along the centreline of the Palaeochannel X bend (Fig. 1e and Supplement).

Figure 6 Cumulative grain size distributions of the scroll bar deposits in (a) Junnerkoeland and (b) Prathoek. Three series were made for Junnerkoeland and two for Prathoek, each indicated by a different line type. Each sample within a series is indicated by a different grey tone. The averages of D₁₆, D₅₀, and D₈₄ are plotted. Figure 1c and d indicate the locations of the grain size samples.

4.2 Dating results

The channel deposits on the inside of Palaeochannel X date from 850–320 and 1408–918 BCE. Palaeochannel X was cut off at 739–117 BCE (Figs. 1e, 4c, e and Table 2).

Table 2 OSL and ¹⁴C dating results from Palaeochannel X. Locations are indicated in Fig. 1c, d and Fig. 4c.

4.3 Meander and channel geometry

The reconstructed transverse bed slopes do not show a spatial trend (Fig. 7a, b), and hence the mean and standard deviations were used in the palaeohydrological calculations. The transverse bed slope at Prathoek is higher ($\mathrm{4.5}\pm \mathrm{1.0}{}^{\circ }$) than at Junnerkoeland ($\mathrm{3.3}\pm \mathrm{1.3}{}^{\circ }$), but much lower than the transverse bed slope of Palaeochannel X ($\mathrm{16.9}\pm \mathrm{1.9}{}^{\circ }$) and of Palaeochannel Q ($\mathrm{28.8}\pm \mathrm{3.8}{}^{\circ }$). The age as a function of distance of lateral accretion follows from Fig. 7c, d. This relation was used for the meander and channel geometry calculations (Fig. 8). The bankfull depths of Palaeochannels X and Q are comparable to the bankfull depths of the meanders Prathoek and Junnerkoeland ca. 1500 CE (Fig. 8a, b) (3 to 4 m). The bankfull depths at Junnerkoeland decreased relatively fast ca. 1800 CE because the erosive base elevation rises towards the cut-off channel (Fig. 4c). At Prathoek, the bankfull depth decreased more gradually over time. The reconstructed bankfull width of Palaeochannels X and Q is much lower compared to the meandering phase (Fig. 8c, d), resulting in a relatively small cross-sectional area of Palaeochannels X and Q (Fig. 8e, f).

Figure 7 Transverse bed slope derived from GPR cross sections from the inner point bar to the outer bend for Junnerkoeland (a, c) and Prathoek (b, d) as well as lateral migration distance plotted against age for both bends. Panels (a) and (b) show the transverse bed slope of lateral accretion surfaces measured in the GPR profile (example in Fig. 4), including the mean and standard deviation of all measurements. Panels (c) and (d) show the relation between age and migration distance of the bends. Shading indicates standard deviation of the Bayesian deposition model determined by Quik and Wallinga (2018) for the OSL and historical map dates.

4.4 Palaeohydrology

The reconstructed Q_bf is 3 to 9 times higher at the start of the meandering phase (85–194 m³ s⁻¹) compared to the preceding phase represented by Palaeochannels X and Q (19–32 m³ s⁻¹) (Fig. 9). The difference in Q_bf between 400 BCE and 1500 CE is significant despite the relatively large uncertainty. A similar discharge in 400 BCE compared to 1500 CE would require a cross-sectional area 5 times larger than currently estimated (Fig. 8e) or a 50 times higher valley slope, which falls outside the uncertainty ranges of these parameters. The Q_bf eventually declines over time and drops to 32–70 m³ s⁻¹ ca. 1850 CE. The calculated Chézy coefficients for the meandering phase ($\mathrm{47.5}\pm \mathrm{0.9}{\mathrm{m}}^{\mathrm{0.5}}{\mathrm{s}}^{-\mathrm{1}}$; Eqs. 7 and 8) were comparable to average Chézy coefficients derived from 12 low-energy rivers ($\mathrm{44.8}\pm \mathrm{13}{\mathrm{m}}^{\mathrm{0.5}}{\mathrm{s}}^{-\mathrm{1}}$) from the river dataset by Kleinhans and Van den Berg (2011).

Combining the frequency of each discharge interval with the sediment transport rate (Fig. 10a) results in a histogram of the sediment transport contribution as a function of discharge (Q_s,freq, Fig. 10b). The highest measured discharge at the gauging station Mariënberg between 1995 and 2015 is 185.5 m³ s⁻¹. The most frequent discharge occurring in the channelized Overijsselse Vecht is 0 to 10 m³ s⁻¹, with a frequency of 8.2 % (Fig. 10a). When discharge is still below bankfull, sediment transport increases relatively fast with an increasing discharge. Above bankfull, additional discharge largely flows across the more flow-resistant floodplain, and hence the sediment transport rates increase less. The effective discharge (Q_eff) is 29 m³ s⁻¹, represented by the highest sediment transport contribution (Fig. 10a, b).

Figure 8 Reconstructed meander and channel geometry over time, assuming the date–distance relations (see Fig. 7c, d) over the scroll bars. Panels (a) and (b) show the bankfull depth (H_bf) derived from the coring data taken from the bottom of the channel lag to the inferred bankfull water surface (Fig. 4c, d), for both the Junnerkoeland (a, c, e) and Prathoek (b, d, f). Panels (c) and (d) show the bankfull width derived from the bankfull depth and reconstructed transverse bed slope (Eq. 1). The river width data from Wolfert and Maas (2007) observed on historical maps and the bankfull river width data from Staring and Stieltjes (1848) were included for comparison. Panels (e) and (f) show the cross-sectional area derived from the bankfull width and water depth (Eq. 2). Shading indicates the 16th and 84th quantile. X&Q refers to Palaeochannels X and Q.

Calculated sediment transport rates were higher than the inner bank growth or scroll bar growth, suggesting the channel deposition can be explained entirely by the reconstructed sediment transport (Fig. 10c). The Q_s,bf of the laterally stable phase was much lower than for the meandering channels, explaining the large difference between the growth rate of the channel deposits on the inner bank at Palaeochannel X (7.0 m³ yr⁻¹) and the scroll bars of Junnerkoeland and Prathoek at the start of the meandering phase ($\mathrm{1.8}\times {\mathrm{10}}^{\mathrm{3}}{\mathrm{m}}^{\mathrm{3}}{\mathrm{yr}}^{-\mathrm{1}}$). Both the sediment transport and average scroll bar growth decreased during the meandering phase.

Figure 11a shows that the river theoretically had insufficient stream power for meandering ca. 400 BCE, and the bar regime was overdamped (Fig. 11b). The stream power seemed sufficient for meandering ca. 1500 CE, and the bar regime was underdamped. The potential for meandering gradually decreased during the meandering phase and became again insufficient when the potential specific stream power drops relatively fast ca. 1850 CE. The damping regime also gradually decreased, but remained underdamped ca. 1850 CE.

Figure 9 Bankfull discharge over time, combined for Junnerkoeland and Prathoek. Shading indicates the 16th and 84th quantile. X&Q refers to Palaeochannels X and Q.

5 Discussion

5.1 Laterally stable phase

A relatively laterally stable phase existed prior to the meandering phase, which is corroborated by the geochronological and palaeohydrological reconstruction. Palaeochannel X formed by extremely slow channel displacement of ca. 6 cm yr⁻¹, assuming a constant channel displacement rate, shown by the OSL dates taken from the channel deposits on the inside of Palaeochannel X (Figs. 1e, 4c and Table 2). The lateral migration rate of the Junnerkoeland meander bend was ca. 40 times higher (Quik and Wallinga, 2018; Wolfert and Maas, 2007). The outer bank of Palaeochannel X consists of Weichselian and Early Holocene deposits. No Middle Holocene deposits were found in the corings (Fig. 4c, d), reflecting the stable character of the Overijsselse Vecht during this period.

The preservation potential of deposits associated with the laterally stable phase is likely very small. Deposits and dimensions of active channel reaches are not preserved during the stable to meandering transition because channel belt dimensions increase. Hence, channel reaches are only preserved when they were cut off prior to the stable–meandering transition, e.g. due to local perturbations. A channel cut-off probably caused Palaeochannels X and Q of the laterally stable phase to become disconnected from the main river before the meandering phase started. In this way these reaches escaped from later lateral erosion during the meandering phase. Consequently, the lateral stability of the river is not immediately evident from these preserved channel reaches because the perturbations led to very slow channel displacement, as was found for Palaeochannel X. However, scroll bar deposits did not form and lateral accretion surfaces were lacking (Figs. 2, 4e and the Supplement), showing that the displacement was not related to meandering in which helicoidal flows cause bar formation and bank erosion at a significant rate and all along the channel (Seminara, 2006). The laterally stable phase lacked the potential to meander given its low position in Fig. 11a and is characterized by an overdamped regime (Fig. 11b, c) and low sediment transport (Fig. 10c). Consequently, the formation of bars was suppressed and the inner bank deposition was small (Fig. 10c).

Figure 10 Sediment transport budgets calculated from present-day flow conditions and from meander migration. (a) Discharge and sediment transport characteristics of the Overijsselse Vecht derived from hourly discharge data from 1995 to 2015 of the gauging station Mariënberg, including the frequency of each discharge class over a year on a frequency scale from 0 to 1 and the sediment transport as a function of discharge for a randomly selected year (1546 CE) in the Junnerkoeland meander bend. (b) Histogram of the sediment transport contribution as a function of discharge. (c) The sediment transport and average scroll bar growth over time (JK: Junnerkoeland, PH: Prathoek, X&Q: Palaeochannels X and Q). The abbreviations Q_s,freq and Q_s,bf are explained in Sect. 3.8. The inner bank growth X refers to the growth rate of the channel deposits on the inner bank at Palaeochannel X, assuming a constant lateral migration rate. Shading indicates the 16th and 84th quantile.

The bend curvature is also an indication for the channel stability. Palaeochannel X comprises a very sharp bend ($\frac{R_{\mathrm{curv}}}{W}=\mathrm{1.4}\pm \mathrm{0.2}$) compared to the meandering phase ($\frac{R_{\mathrm{curv}}}{W}=\mathrm{2.1}\pm \mathrm{0.4}$), which is often found in low-energy rivers in which lateral migration is limited (Candel et al., 2018; Hickin and Nanson, 1984). Large similarities exist between the laterally stable phase reported here and the laterally stable channels in highly cohesive sediment on the intertidal mudflat, which are mostly laterally stable except for some sharp bends at which bank failure and flow separation result in very limited and local channel migration (Kleinhans et al., 2009).

5.2 Channel pattern change

The Overijsselse Vecht river changed from a laterally stable into a meandering river. Differences in palaeohydrological conditions between the two phases were large enough to distinguish despite the large uncertainties in the palaeohydrological reconstruction. Bar regime changed from an overdamped regime into an underdamped regime (Fig. 11b, c), leading to overdeepening of the outer bend pool and enhancement of the point bars in the inner bend (Struiksma et al., 1985; Crosato and Mosselman, 2009; Kleinhans and Van den Berg, 2011). The significantly higher bankfull discharge (factor of 3 to 9; Fig. 9) explains the potential to meander (Fig. 11a), the high sediment transport, and the high scroll bar growth (Fig. 10c) at the start of the meandering phase.

The exact moment of the channel pattern change is between the cut-off of Palaeochannel X (400±300 BCE) and the reconstructed initiation of scroll bar formation (1504±52 CE). Most likely, the transition occurred shortly before the latter because both point bars had a relatively similar meander start age (Figs. 2 and 4), the surrounding floodplain is formed by Late Glacial or Early Holocene deposits (Fig. 4), and there is no evidence of older scroll bar deposits in the vicinity of the studied meander bends. Mature meandering river systems would always leave traces of older scroll bar deposits, channel cut-offs, or meander scars because these are never completely being removed by the river (Toonen et al., 2012; Van de Lageweg et al., 2016).

The palaeohydrological reconstruction shows that the increasing bankfull discharge likely explains the channel pattern change. The increasing bankfull discharge may reflect an increase in annual discharge, but could also be related to a more irregular discharge regime because the bankfull discharge largely represents the higher discharges in a river (Dury, 1973; Wolman and Miller, 1960). Consequently, the discharge may have been constant over a year with low peak discharges and a relatively high base flow during the laterally stable phase, changing into a more peaked discharge regime with a relatively low base flow at the start of the meandering phase.

Figure 11 The potential for meandering with time. (a) The potential specific stream power in a stability diagram (Eq. 11). Two discriminators were plotted for a range of median particle sizes of the bed sediment, which is the range of particle sizes found in the scroll bars and Palaeochannels X and Q (Fig. 6). Panels (b) and (c) show the bar regime for both Junnerkoeland and Prathoek, determined with the interaction parameter (IP) (Eq. 18) and compared to the thresholds (Eqs. 19 and 20). Shading indicates the 16th and 84th quantile. X&Q refers to Palaeochannels X and Q.

A potential cause of the discharge regime and channel pattern change may be the climate change at the start of the Little Ice Age (14th to 19th century) (Grove, 1988) given the overlap in time with the meandering phase (Fig. 2c, d). Although geomorphological responses differ for each river during the Little Ice Age, enhanced lateral migration or incision was generally observed for most rivers in north-western Europe (Rumsby and Macklin, 1996). The increased bankfull discharge in the Overijsselse Vecht may have been caused by higher run-off relative to precipitation due to reduced evapotranspiration rates and frozen soils (Rumsby and Macklin, 1996; Van Engelen et al., 2001) and/or a higher snowfall to rainfall ratio due to lower winter temperatures in the Netherlands and Germany (Behringer, 1999; Lenke, 1968). Higher snowfall rates were also recorded for the United Kingdom (Manley, 1969), where they led to more flooding during the snowmelt period (Archer, 1992). Studies on historical observations of rivers nearby the Overijsselse Vecht (IJssel, Elbe, Lower Rhine, and Meuse) suggested a significantly higher flooding rate during the Little Ice Age compared to more recent flooding rates (Glaser et al., 2010; Glaser and Stangl, 2003; Mudelsee et al., 2003, 2004).

An additional cause for an increasing bankfull discharge may have been land use change in the catchment (Kondolf et al., 2002), which affects the discharge regime due to the direct relation with evapotranspiration (Fohrer et al., 2001). For the Overijsselse Vecht catchment, peat reclamation started in the 12th and 13th century (Gerding, 1995; Van Beek et al., 2015b) and intensified from the 14th century onwards (Borger, 1992; Van Beek et al., 2015a). Reclamation of peatlands partly comprised digging canals to drain the land, and although the reclamation was mainly limited to the margins of peatlands, the hydrological consequences were large. The margins are a natural seal of the peat bog, with a low hydraulic conductivity compared to the remainder of the bog, ensuring peat dome growth. Destruction of these margins will result in drainage of the entire peat bog (Baird et al., 2008; Van der Schaaf, 1999). After several centuries, focus shifted from peat reclamation to exploitation by excavating large peatland areas for fuel during the 17th and 18th century (Gerding, 1995). The largest part of the peat has currently disappeared. Yearly average discharges in peatlands can increase by 40 % in the Dutch climatological setting due to evapotranspiration differences for reclaimed peat areas compared to undisturbed peat areas (Baden and Eggelsmann, 1964; Streefkerk and Casparie, 1987; Uhden, 1967). This increase cannot fully explain the large increase in bankfull discharge in the Overijsselse Vecht (factor of 3 to 9) because peat covered just ca. 27 % of the Overijsselse Vecht catchment area during the 14th century (Casparie and Streefkerk, 1992; Vos et al., 2011), and hence the yearly average discharge of the catchment increased by ca. 11 % due to evapotranspiration differences.

However, several studies have also shown that an increased drainage network in peatlands resulted in higher discharge peaks with a fast discharge response to precipitation (Conway and Millar, 1960; Holden et al., 2004, 2006; Streefkerk and Casparie, 1987). For example, the run-off to rainfall ratio was a factor of 3 higher in a drained Irish peatland compared to an undrained Irish peatland (Burke, 1975), which is comparable to the observed bankfull discharge increase in the Overijsselse Vecht. Finally, canals were not only dug for peat reclamation, but also for shipping and effective generation of water power starting in the 11th and 12th century (Driessen et al., 2000), which may have promoted the higher peak flows even more. New canals resulted in a faster run-off, but also changed the watershed delineation (Driessen et al., 2000). We conclude that both climatic and land use changes were likely responsible for an increase in both total discharge and peak flows, resulting in the transition of a relatively stable river to a highly dynamic meandering system.

5.3 Meandering phase

Interestingly, the bankfull discharge declined during the meandering phase (Fig. 9), leading to decreasing sediment transport relative to the scroll bar growth (Fig. 10c) and insufficient potential specific stream power for meandering after ca. 1850 CE (Fig. 11a). This decline was corroborated by observations of river width from previous studies, which can be compared to the reconstructed widths (Fig. 8c, d). These observations included measurements from historical maps by Wolfert and Maas (2007) and measurements of the bankfull river width over a large river section in 1848 CE by Staring and Stieltjes (1848). The river width data from Wolfert and Maas (2007) largely fall in the range of reconstructed bankfull widths at Junnerkoeland and show a similar decreasing trend (Fig. 8c). However, the historical maps used by them may result in large uncertainties because the water stage that these maps represent is unknown (Quik and Wallinga, 2018). The measured widths by Staring and Stieltjes (1848) are in line with the predicted width at Junnerkoeland, falling within the uncertainty range. The predicted width at Prathoek is underestimated compared to the measured widths by Wolfert and Maas (2007) and Staring and Stieltjes (1848). This underestimation may explain the lower cross-sectional area compared to Junnerkoeland (Fig. 8f) and hence an underestimated bankfull discharge (Fig. 9) and potential specific stream power (Fig. 11a).

The observed decline of bankfull discharge would suggest that the hydrological forcing disappeared or diminished and had a temporary character, which would fit the Little Ice Age that ended in the 19th century as a potential cause. Consequently, it would be expected that the channel pattern reorganized and became laterally stable again. However, the river was still laterally migrating until channelization between 1896 and 1914 CE (Wolfert and Maas, 2007), which may be related to the presence of an underdamped regime enhancing point bar formation in the inner bend (Fig. 11b, c). Additionally, historical bank stability changes may have promoted the river meandering during this period. For example, floodplains were intensively used for cattle grazing, which may have weakened the banks, enhancing meandering after 1850 CE (Beschta and Ripple, 2012; Trimble and Mendel, 1995; Wolfert et al., 1996). Drift sand activity was also initiated by intensive land use since the Late Middle Ages (Fig. 1c, d) (Koster et al., 1993), which may have affected the bank stability.

6 Conclusions

We show that bankfull discharge and the associated river parameters can be reconstructed by following a stochastic approach and through detailed geochronological and lithological analysis of scroll bar deposits and palaeochannels. For the Overijsselse Vecht river we demonstrate that an increase in bankfull discharge ca. 1400 to 1500 CE resulted in a river channel pattern change from laterally stable to meandering. Geochronological data confirmed our hypothesis on the lateral stability of the river prior to the meandering phase, in contrast to previous assumptions that were made of continuous meandering during the Holocene. We show that the reconstructed river parameters are consistent with both the laterally stable and meandering channel pattern by applying empirical channel and bar pattern models. Potential causes for the discharge regime changes include climate change (Little Ice Age) and land use changes (peat reclamation, peat exploitation, digging of canals). We conjecture that the change from laterally stable to meandering has occurred in other rivers for which increased Holocene fluvial activity was reported.