Reconstructing lateral migration rates in meandering systems – a novel Bayesian approach combining optically stimulated luminescence (OSL) dating and historical maps

Identifying lateral migration rates of meandering rivers is relevant both for fluvial geomorphology and to support river management. Lateral migration rates for contemporary meandering systems are often reconstructed based on sequential remote-sensing images or historical maps; however, the time frame for which these sources are available is limited and hence likely to represent fluvial systems subjected to human influence. Here, we propose to use scroll bar sequences as an archive to look further back in time using optically stimulated luminescence (OSL) dating of sand-sized quartz grains. We develop a modelling procedure for the joint Bayesian analysis of (OSL) dating results and historical map data. The procedure is applied to two meanders from the Overijsselse Vecht, a medium-sized sand-bed river in the Netherlands. We obtained nine samples for OSL dating from scroll bars and combined OSL dating results with historical map data for the period 1720–1901 CE (Common Era). The procedure we propose here incorporates the strengths of both data types for reconstructing fluvial morphodynamics over longer time frames. Using an iterative modelling approach, we translate spatial uncertainty of historical maps into temporal uncertainty of channel position required for Bayesian deposition modelling. Our results indicate that meander formation in the Overijsselse Vecht system started around 1400 CE, and lateral migration rates were on average 2.6 and 0.9 m yr⁻¹ for the two investigated bends, until river channelization around 1900 CE.

1 Introduction

Rivers are one of the world's most important geomorphic agents and represent highly dynamic earth surface systems (Güneralp and Rhoads, 2011; Vandenberghe and Maddy, 2000). Meandering rivers represent the most common river type and can be found all over the world (Hooke, 2013). The mobility of their course has a significant effect on sediment transport and deposition, floodplain development, and landscape change (Kleinhans and Van den Berg, 2011). Related consequences for river management, hazards, and biodiversity pose a practical need to understand the timing and speed of lateral migration (Grabowski et al., 2014; Hooke, 2013; Lespez et al., 2015).

To determine whether changes in meander position and morphology have taken place, it is often necessary to consider relatively long timescales (Hooke, 2013). To establish recent lateral migration rates and reconstruct the planform of contemporary fluvial systems, sequential aerial photographs or topographical maps are often used (e.g. Pišút, 2002; Uribelarrea et al., 2003; Hooke and Yorke, 2010; Eekhout et al., 2013). Unfortunately, the lack of historical sources of this type limits long-term reconstructions of channel position, and reconstructions based on older maps may be strongly affected by spatial uncertainties and mapping inaccuracies. Moreover, as historical maps are often only available for cultivated areas, the fluvial systems displayed are likely to be affected by humans. Such influence can be either through river management (e.g. Hesselink et al., 2003; Frings et al., 2009) or through land-use changes which may change the hydrograph (e.g. Candel et al., 2018) as well as sediment input (e.g. De Moor et al., 2008; Hoffmann et al., 2009). Consequences may be great, as shown in a highly influential paper by Walter and Merritts (2008), who argued that the present-day meandering form of gravel-bedded mid-Atlantic streams in the USA is not a natural phenomenon but rather an artefact caused by mill ponds. Hence, to understand dynamics of fluvial systems and anthropogenic impacts on those systems, there is a need to extend the period of reconstruction to earlier periods, which requires reconstruction methods that do not rely solely on historical sources and/or remote sensing.

To resolve depositional ages and geomorphological process rates, the use of optically stimulated luminescence (OSL) dating to obtain absolute chronologies of fluvial deposits is becoming widespread (Lian and Roberts, 2006; Rittenour, 2008; Wallinga, 2002b). In addition, OSL is increasingly used to date young (<1000 years) sediments as part of late Holocene studies (Madsen and Murray, 2009). Yet, the application of OSL dating to reconstruct lateral migration rates has hardly been attempted (with a few notable exceptions; see, e.g., Rodnight et al., 2005; Rowland et al., 2005; Kemp and Rhodes, 2010). Such limited application of OSL dating for this purpose may well be related to challenges faced when applying luminescence dating to young fluvial sediments. Successful age determination requires the complete resetting of the OSL signal in at least part of the grains, a prerequisite that may not be met due to the limited light exposure during the subaqueous transport of grains in a turbid river.

Bayesian analysis of depositional sequences allows the construction of more precise and robust chronologies by combining ages obtained through a dating method with prior information on the order of events (Bronk Ramsey, 2008). Following pioneering work by Rhodes et al. (2003), Bayesian analysis is now increasingly applied to sequences of OSL ages, often in combination with ages obtained through other methods (e.g. ¹⁴C: Dreibrodt et al., 2010; U-series: Clark-Balzan et al., 2012; dendrochronology: Wallinga et al., 2012) or even multiple methods (e.g. Shanahan et al., 2013; Hobo et al., 2014; Guérin et al., 2017).

Several studies have combined OSL dating with evidence derived from historical maps and documentation as independent age control (e.g. Ballarini et al., 2003; Medialdea et al., 2014), age constraints (e.g. Nielsen et al., 2006; Madsen et al., 2007; Tamura et al., 2011; Kunz et al., 2014), or for historical contextualization of dating results (e.g. Clemmensen et al., 2007). However, the direct integration of OSL dates with legacy data provided by historical evidence into a single geochronology, i.e. where historical evidence is used to help determine the geochronology, is underexplored. An exception to this is the study by Hobo et al. (2014), who reconstructed the floodplain sedimentation of the Dutch river Waal using Bayesian age–depth models based on heavy-metal age estimates, OSL dates and age constraints derived from historical maps.

Figure 1 Location of the study area (a), showing the two meander bends Junner Koeland (b) and Prathoek (c). Dotted line in panel (a) indicates the valley side; m a.s.l.: metres above sea level. White numbers along the coring transects in panels (b) and (c) indicate OSL sample locations by abbreviated sample code (all should be preceded by NCL-2415). Samples NCL-2415154–162 were collected in scroll bar deposits; sample NCL-2415163 was collected in a terrace remnant. Lithological cross sections and lithogenetic interpretation are provided in Candel et al. (2018). Digital elevation model (AHN2; horizontal resolution 5 m, vertical resolution 0.2 m): AHN (2018), Van Heerd et al. (2000). Topography (3200 pixels km⁻¹): OpenTopo, Van Aalst (2017).

Here we combine OSL data and historical map evidence using Bayesian chronological modelling to reconstruct planform development and lateral migration rates of the Overijsselse Vecht, a medium-sized sand-bed river in the Netherlands. Our aim is to develop a modelling procedure using the Bayesian depositional sequences tool of OxCal (Bronk Ramsey, 2008) for the joint analysis of OSL dating results and historical map data that (1) allows quantitative analysis of historical maps in geochronological research and (2) supports the development of robust fluvial chronologies with quantified uncertainty that expand into larger time frames. With these methods we determine the age of scroll bar deposits along the Overijsselse Vecht and the average migration rate of two meanders that developed during the past 500 years. This information provides the geochronological basis for a palaeo-hydrological study on the Overijsselse Vecht (Candel et al., 2018).

2 Study area

The Overijsselse Vecht is a medium-sized river originating in Nordrhein-Westfalen in Germany and enters the Netherlands from the east near the town of Coevorden. This rain-fed river has a catchment size of 3785 km² and a rather uniform valley gradient on Dutch territory of $\mathrm{1.4}\times {\mathrm{10}}^{-\mathrm{4}}$ (Wolfert and Maas, 2007). The width of the abandoned channel (i.e. abandoned upon channelization) in the investigated meanders (see below) is ∼30 m, as derived from the present-day digital elevation model (DEM; Fig. 1). According to measurements in the period 1995–2015 of the discharge station in Mariënberg (just upstream of our study area), the average annual discharge and mean annual flood discharge amount to 22.8 and 160 m³ s⁻¹ respectively. Average annual rainfall in the region amounts to 800–875 mm and average monthly temperatures vary between 4.5 and 5.0 ^∘C in January and 22.5 and 23.0 ^∘C in July (KNMI, 2017).

The Overijsselse Vecht was channelized in the late nineteenth and early twentieth century. The original river length of 90 km on Dutch territory was reduced to 60 km by cutting off 69 meanders (Wolfert and Maas, 2007). Weirs were installed to prevent vertical erosion of the river bed, and groynes and revetments were constructed to prevent lateral migration of the river channel (Wolfert et al., 1996; Wolfert and Maas, 2007). In recent decades, water managers have initiated several projects to restore the physical and ecological functioning of the river, which resulted in the need for increased understanding of river morphodynamics (e.g. Wolfert et al., 1996, 2009; Maas et al., 2007; Viveen et al., 2009) and ecological potential (e.g. Duursema, 2004).

The floodplain of the Overijsselse Vecht is flanked on both sides by extensive areas with coversand ridges and drift-sand dunes, overlying fluvio-periglacial deposits of Pleistocene age (Ter Wee, 1966, 1979; Kuijer and Rosing, 1994; geomorphological map of the Netherlands (1 : 50 000); for details, refer to Koomen and Maas, 2004). Its former valley is rather narrow (indicated by the yellow dotted lines in Fig. 1a) and many meanders appear to be laterally confined by relatively stable sides of the former valley. Relief in the area is limited; the floodplain and valley sides are located at ∼5.5 and ∼7 m a.s.l respectively. For this study we focused on two meander bends with large amplitudes that reach outside the former valley side, named Junner Koeland and Prathoek (Fig. 1).

Table 1 Overview of the six historical maps used in this study. All dates are indicated in CE (Common Era). Revision date indicates full or partial revision. The date used for analysis is based on the revision date or, if none is given, on the survey date. For the map sheet of the Hottinger atlas the survey covered 3 years; we therefore took the middle year to use in further analyses. GCPs: ground control points (as used for georeferencing); RMSE: root mean square error (expressed in metres) based on the residuals (or displacement vectors) of the GCPs after using a first-order polynomial transformation in the georectification procedure. JK: Junner Koeland; PH: Prathoek. Map sections are shown in Fig. 2. References list map source and publications on cartographers, survey techniques and depiction methods. References are detailed in the subscript.

1: Algemeen Rijks Archief. Genie-archief, situatiekaart 07. 2: Wolfert et al. (1996). 3: Box (2007). 4: Versfelt (2003). 5: Versfelt and Schroor (2005). 6: CC-BY Kadaster (2018a).
7: CC-BY Kadaster (2018b). 8: CC-BY Kadaster (2018c). 9: Van der Linden (1973). 10: Van der Leest et al. (2005). 11: Stam (2006).

3 Methods

3.1 Lithological survey

For each meander bend a lithogenetic cross section was constructed based on multiple hand corings from meander base to apex (Fig. 1b, c), following a line perpendicular to the scroll bars, with cores obtained at high and low points in the scroll bar and swale topography. In a separate paper, Candel et al. (2018) provide details of these cross sections (including detailed facies descriptions and an explanation of the lithogenetic interpretation) and use the information to reconstruct palaeo-discharges. Here, we concentrate on those aspects that are directly relevant to the reconstruction of channel position throughout the lifetime of the meander. Most of the cores reached a depth of 4 to 6 m, penetrating the full point-bar deposits (extremely fine to coarse sand) including the channel lag (medium fine to coarse sand with gravel), which occurred at depths varying from 2.5 to 4 m below the surface. Some aeolian reworking occurred locally, as is evident from the DEM (Fig. 1). The lithological and geomorphological information was used to select the suitable depths to collect OSL samples (see below).

3.2 Analysis of historical maps

Maps show the position of the river channel at the time of the map survey; however, the use of this information for quantitative reconstruction of channel position over time is not straightforward. First of all, historical maps are likely less accurate than their modern counterparts. Hence there is an uncertainty in the exact position of the channel. In this section we outline how we quantify this spatial uncertainty. Next, to allow Bayesian analysis of the sequence, this spatial uncertainty must be translated into a temporal uncertainty. The necessity of this step and the approach we take are outlined in Sect. 3.4.

3.2.1 Georeferencing

Six historical maps (dating from 1720 to 1894/1901 CE (Common Era); Table 1, Fig. 2a–f) were selected for analysis. The selected maps vary in extent (e.g. focused on the Overijsselse Vecht or on the Netherlands as a whole), scale (varying from 1 : 14 400 to 1 : 50 000), and original purpose (military or topographical surveys); consequently, they vary in their display of the study area. We included these maps in our analysis based on the degree of topographical quality (i.e. sufficient detail for georeferencing) and a quality check to assess geodetic distortion (see Sect. 3.2.3). Younger maps (example for 1931 CE provided in Fig. 2g) did not provide relevant information as the meanders under investigation were cut-off during river channelization in the early twentieth century and were consequently no longer active. Older maps (>1720 CE) were also available for the area but had to be excluded from the analysis as georeferencing was impossible due to the limited topographical information on the maps. The historical context of the maps and information on the cartographers, survey techniques, and map depiction methods have been discussed in previous publications (see Table 1 for sources and references). For general background information on the use of historical sources and relevant concepts, see, e.g., Hooke and Kain (1982; particularly chap. 3 on accuracy and analysis), Petts et al. (1989), or Gurnell et al. (2003).

Figure 2 (a–e) Sections of historical maps displaying the study area (for explanation refer to Table 1); (g) map showing the same area in 1931 CE (publication date). As section (f) consists of two map sheets two dates are listed (details in Table 1). Image references: (a) Algemeen Rijks Archief, Genie-archief, situatiekaart 07; (b) Versfelt (2003); (c) Versfelt and Schroor (2005); (d) CC-BY Kadaster (2018a); (e) CC-BY Kadaster (2018b); (f) CC-BY Kadaster (2018c); (g) CC-BY Kadaster (2018d).

The relevant map parts were used with the given digital image quality or scanned and saved at a resolution of 600 dpi. The maps were georeferenced using Esri's ArcMap (versions 10.3.1 and 10.5) based on static landscape features (e.g. churches, road crossings) that were displayed on both historical and reference maps. For the latter we used a digital elevation model (AHN2, 5 m horizontal resolution, 0.2 m vertical resolution) and a recent topographical map (OpenTopo) (AHN, 2018; Van Heerd et al., 2000; Van Aalst, 2017). Georeferencing maps of river landscapes can be challenging due to limited presence of stable landscape features in floodplains. As the valley of the Overijsselse Vecht is rather narrow (Fig. 1a), we were able to adequately georeference the map sections based on GCPs (ground control points) located close to, but not inside, the floodplain. We predominantly made use of clearly delineated features such as road crossings but occasionally included features with more irregular edges such as arable fields (ground control points are available in the Supplement). Maps were georectified using a first-order polynomial transformation (following, e.g., Downward et al., 1994; Leys and Werritty, 1999) and nearest-neighbour pixel resampling as this resulted in the sharpest display of the river channel (based on visual assessment). Subsequently all maps were displayed using the Dutch RDnew projection. The ground control points were exported from ArcMap as a *.txt file to use in the freely available software MapAnalyst (see below). Map quality was analysed using (1) a quantification of geospatial error or planimetric accuracy and (2) visualizations of geodetic distortions as generated by MapAnalyst. Both aspects are explained below.

3.2.2 Quantifying geospatial error

The root mean square error (RMSE) of the displacement vectors of the GCPs is often used as a quantification of geospatial error and was automatically calculated by ArcMap (Esri, 2018) according to

$$\begin{array}{}\text{(1)}& {\displaystyle }\mathrm{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{n}}\left(\sum _{i=\mathrm{1}}^n\left({\mathit{\epsilon }}_{x_i}^{\mathrm{2}}+{\mathit{\epsilon }}_{y_i}^{\mathrm{2}}\right)\right)},\end{array}$$

where ${\mathit{\epsilon }}_{x_i}$ and ${\mathit{\epsilon }}_{y_i}$ respectively represent the x and y component of the displacement vector that indicates the difference between the location of GCP i on the reference map and the location of the same point on the historical map after first-order polynomial transformation. The total number of GCPs is indicated by n. The RMSE was used to assess geospatial error or planimetric accuracy (for the range of spatial errors see RMSEs in Table 1), defined as “the extent to which distances and bearings between identifiable objects coincide with their true values” (Jenny and Hurni, 2011).

After georeferencing, the channel centre line was deduced from each map to reconstruct the channel position for each time slice. The RMSE was used to draw an error band around each centre line. The centre line and error band represent the positional mean and related standard deviation respectively. The channel centre lines and error bands were subsequently used in the chronological modelling (see Sect. 3.4). The final position of both investigated meanders was based on the abandoned channel as displayed by the DEM of the present-day situation. Historical evidence indicates that Prathoek was channelized between 1884 and 1901 CE (Fig. 2ef). Junner Koeland was channelized between 1894 (Fig. 2f) and 1906 CE, as indicated by the topographical map of the Netherlands 1 : 50 000 (map sheet Coevorden, surveyed between 1900 and 1906; not shown).

3.2.3 Visualization of geodetic distortion

Distortion grids were generated for each map section based on the GCPs using the software MapAnalyst (Jenny and Hurni, 2011; http://mapanalyst.org/, last access: 12 September 2017). A world file (*.txt) was created for the reference image and, together with the historical map image, loaded in MapAnalyst. In the *.txt file with GCPs that was generated with ArcMap, the first two columns show the map points indicated in pixel rows/columns, using the upper left corner as the origin. The third and fourth column show the points of the reference image, using the units of its coordinate system (in our case the Dutch RDnew projection with units in metres). Using a text editor, the points for the map and reference image were split into two separate files. In addition, each point pair was assigned an ID number in the first column to allow MapAnalyst to link the point pairs from the two files. For the map images MapAnalyst uses the lower left corner as origin instead of the upper left corner as is done in ArcMap; therefore, the pixel values in the y direction for the map points were converted prior to loading them in MapAnalyst. Map distortions were calculated using a 1 km raster and exported as *.png images.

3.3 OSL dating

OSL dating makes use of a low-intensity light signal emitted by mineral grains (e.g. quartz) upon optical stimulation. The intensity of the signal can be linked to the age of the sediment. The (latent) OSL signal builds up as the grains are exposed to natural ionizing radiation in the environment (from uranium, thorium, potassium, and cosmic rays; known as the environmental dose), which causes the accumulation of trapped charge in the crystal structure. When the grains are exposed to daylight during transport, the trapped charge is released and the OSL signal is reduced (bleached) to a low level, often close to zero. When the grains are deposited and buried, they are no longer exposed to light and the signal can build up again.

Through a comparison of the natural luminescence signal with that induced by laboratory irradiation, the environmental dose received by the grains since burial can be estimated (the palaeo-dose). Combined with information on the environmental dose rate experienced by the grains, the age can be calculated by age (ka) = palaeo-dose (Gy) ∕ dose rate (Gy ka⁻¹) (see, e.g., Rhodes, 2011 for more information on OSL dating).

With dating of fluvial sediments, the assumption in luminescence dating that the OSL signal was completely zeroed prior to sediment burial is questionable (Murray et al., 1995) due to the limited exposure of grains to light when transported in turbid water (e.g. Wallinga, 2002a). This means that at the time of sediment burial, some trapped charge that was built up since the previous deposition may remain, which may consequently lead to a remnant OSL signal upon deposition and an overestimation of the age. Poorly bleached sediments mostly consist of a mixture of poorly bleached and well-bleached grains (Olley et al., 1998), which can be detected by broad and skewed equivalent dose distributions (e.g. Wallinga, 2002b). Heterogeneous bleaching in channel site samples may result in apparent doses ranging over several orders of magnitude and may limit the accuracy of dates (Murray et al., 1995; Olley et al., 1998), although methods have been developed to largely overcome these challenges (e.g. Galbraith et al., 1999; Cunningham and Wallinga, 2012).

3.3.1 OSL sampling

In total 10 samples were collected for OSL dating to determine the time of deposition and burial of the sediments. Six samples were collected in the Junner Koeland scroll bar deposits, three samples in the scroll bar deposits of Prathoek, and one sample in an early Holocene fluvial terrace remnant near Prathoek. To reduce the chances of sampling sediments that were reworked after initial scroll bar formation, we avoided sampling in swales. Samples were roughly evenly distributed over the scroll bar sequence (sample spacing ∼200 m; see Fig. 1b and c). The scroll bar samples were collected just above the channel lag deposit because (1) the channel lag can generally be recognized easily, thereby preventing accidentally sampling older deposits below the scroll bar deposit, (2) at this depth recent aeolian reworking of sediment can be excluded, and (3) the channel lag is located below the oxidation–reduction zone, supporting the assumption that the sediments have been water-saturated since deposition (and thus reducing uncertainty in calculating the environmental dose rate).

Samples for OSL dating were obtained through a modified Van der Staay suction corer (Wallinga and Van der Staay, 1999). First a hole was augered down to the desired depth, using an Edelman auger above the groundwater level and the suction corer below groundwater level. Then, the suction corer, with a loosely fitted PVC sample tube of 30 cm length at its end, was inserted in the hole and pushed down to the lower sample depth. After retrieving the suction corer, the sample tube was removed, and both ends were immediately covered with plastic caps and taped with light-impermeable black tape. Sampling locations were recorded with a GPS device with a horizontal precision of around 5 m.

3.3.2 OSL measurements

The samples were analysed in the laboratory of the Netherlands Centre for Luminescence dating in Wageningen, the Netherlands. Under amber safelight conditions, the sample tubes were opened and two subsamples were obtained from each core; material from the light-exposed outer ends of the tube was prepared for dose rate estimation (Sect. 3.3.3), while the material from the inside of the tube was suitable for palaeo-dose estimation. Luminescence measurements were performed on quartz grains of the fraction 212–250 µm. This fraction was obtained by a combination of sieving and subsequent treatment with HCl (to remove carbonates), H₂O₂ (to remove organic material), and HF (45 min at 40 %, to dissolve feldspars and etch the quartz grains). Samples were measured on a Risø TL/OSL DA-20 reader, using blue LEDs for stimulation and a 7.5 mm U340 filter for detection. The single-aliquot regenerative-dose protocol (SAR) (Murray and Roberts, 1998; Murray and Wintle, 2000, 2003) was used, with details as listed in Table 2. Net OSL signals were obtained using an early background subtraction approach (Cunningham and Wallinga, 2010). The preheat temperature was selected based on preheat plateau tests. With the applied procedure, a laboratory dose (i.e. exposure to a known amount of radiation) could accurately be recovered, as indicated by the average dose recovery ratio of 1.03±0.01. At least 22 aliquots, each consisting of ∼75 grains (2 mm mask size), were measured per sample (31 aliquots on average). The equivalent dose distributions were displayed in radial plots, constructed with the Luminescence package (Kreutzer et al., 2015) employed in R (version 3.2.4).

3.3.3 Dose rate

Light-exposed material from the outer ends of the sample tubes was prepared for dose rate estimation. Samples were dried and combusted, allowing the measurement of water and organic content, and then mixed with wax and moulded into pucks. These were measured on a high-resolution broad-range gamma-ray spectrometer, to determine radionuclide concentrations of K-40 and several radionuclides in the U and Th decay chains. From these, dose rates experienced by the quartz grains were calculated, taking into account attenuation due to water, organics, and grain size, and adding a contribution from cosmic rays (assuming instant burial to present depth). Methods are detailed in Wallinga and Bos (2010).

Table 2 The SAR protocol used for equivalent dose measurement, including measurement parameters.

3.3.4 Palaeo-dose estimation

We analysed the equivalent dose distributions of all samples from scroll bar deposits using the unlogged version of the bootstrapped Minimum Age Model (bsMAM; Galbraith et al., 1999; Cunningham and Wallinga, 2012). This approach is specifically suited to our project as it takes into account the uncertainty of σ_b (see below) and allows the construction of likelihood distributions of the palaeo-dose, which can be used as priors for Bayesian analysis (Cunningham and Wallinga, 2012). The unlogged version of the MAM (Minimum Age Model) was used, as the logged version was not applicable to the youngest samples. Being aware of reservations with regard to using the unlogged model for older samples (Arnold et al., 2009), we checked that, for older samples, the results of logged and unlogged models were in agreement.

The (bootstrapped) MAM requires an estimate of overdispersion in the distribution that is not caused by heterogeneous bleaching (σ_b). Ideally, this estimate is based on the observed overdispersion for well-bleached samples from the same lithology, depositional environment, and age (e.g. Chamberlain et al., 2018). As none of the samples obtained from the Holocene scroll bars appeared to be well-bleached, we had to take another approach. This alternative is provided by the overdispersion obtained for the oldest and best-bleached sample (NCL-2415163; D_e distribution available in Fig. 6b), which was collected from an early Holocene terrace remnant. Although the depositional environment and age are different, the lithology of the material is nearly identical as the scroll bar deposits consist mostly of reworked fluvial deposits of Late Pleistocene/early Holocene age. The overdispersion of sample NCL-2415163 was 0.15±0.03, and this value was used as σ_b input for the bsMAM model. The bootstrapped MAM analysis of the equivalent dose distribution is combined with the sample dose rate to provide a likelihood distribution of estimated burial ages for each sample (see Cunningham and Wallinga, 2012). These age likelihood distributions are implemented as OSL priors in the Bayesian modelling (see Sect. 3.4). Sample NCL-2415163 did not suffer from poor bleaching and based on the normally distributed and relatively narrow equivalent dose distribution for this sample, the Central Age Model (CAM; Galbraith et al., 1999) was adopted for palaeo-dose estimation.

3.4 Bayesian reconstruction of river channel position

Bayesian statistical approaches were used to combine the numerical age data obtained from the historical maps and OSL dating, with auxiliary data on the stratigraphic order. This analysis was performed using the freely available OxCal software (version 3.4; main introduction to the programme described by Bronk Ramsey, 1995; see also Bronk Ramsey, 2009). Bayesian analyses of deposition sequences are normally based on age versus depth profiles (Bronk Ramsey, 2008). Here, we analyse age as a function of lateral migration distance, where the location of the oldest swale (meander base) is taken as the origin and distance is measured along the central axis of the scroll bar deposits. For each of the age constraints (OSL sample age or survey age of a historical map), we have both an age and a location along the central axis. In addition, we know that the meander grew from its base to the final position, providing an order for all information. In the Bayesian approach, all this information is combined to obtain a robust reconstruction of the channel position through time.

To jointly analyse the OSL priors and the historical map data, both were added as entries in a Bayesian deposition model using the P_Sequence (Fig. 3). This command fixes the chronology of the sampling points using the known order of events (based on the known meander expansion direction) and makes use of a Poisson distribution. With this depositional model, channel migration is assumed to be random (i.e. varying in rate but always positive) giving approximate proportionality to distance (factor z) (Bronk Ramsey, 2008). The P_Sequence model needs information on the size of the units (deposition events) within the sequence. Here, we chose to use the average dimension of a scroll bar as unit, i.e. the mean distance between two swales. Based on meander amplitude and the visual counting of scroll bars from the DEM, this distance was estimated to be 77 and 42 m for Junner Koeland and Prathoek respectively. This provides an event spacing (k₀ value) of 0.013 events m⁻¹ for Junner Koeland and 0.024 events m⁻¹ for Prathoek. By applying the regular boundary function, the entire distributions were allowed to be used in the calculation.

Likelihood distributions of estimated burial age (see Sect. 3.3.4) were uploaded to the OxCal directory as priors (*.prior file), with the unit in years CE. For this OSL age data, the position is relatively well known (∼5 m), while there is considerable uncertainty in the time domain (age estimate).

Figure 3 OxCal script as used for the fourth iteration for Junner Koeland. A similar script was used for Prathoek.

For the historical map data we assume the age is known exactly (survey date; see Table 1), while reconstructed locations are uncertain (up to 107 m; see RMSEs in Table 1). Incorporating this data into the Bayesian model requires uncertainty to be expressed in the time domain rather than spatial domain. Moreover, a slight correction is needed as the OSL data provide age information on the scroll bar formation (i.e. inner bank), while the channel centre line was reconstructed from the maps. Channel width varied significantly between the historical maps and may not represent the true low-water channel width of past times. The historical river channel was therefore deduced from each map as a low-water channel centre line. To make sure that the age of the map is assigned to the inner bank instead of the channel centre in the chronological analysis, we corrected the distance (z value) for each centre line using half the channel width. For this calculation, we assumed that channel width has been relatively stable throughout meander formation and is represented by the width of the abandoned channel (∼30 m) as derived from the DEM (Fig. 1).

Uncertainty in the spatial domain was calculated for each map based on σ (see Sect. 3.2.2) and assuming a normal distribution. To transpose this spatial uncertainty to an age uncertainty, the migration rate of the channel is required. This provides a problem, as the procedure aims at a reconstruction of channel position through time, which is needed to determine the migration rate. Nevertheless, a crude estimation of migration rate is possible based on the available data, and an iterative approach can be used to improve this estimation with the model output. For the first iteration error bands were converted to years using an estimated migration rate based on the bsMAM age of the OSL sample closest to the meander base and historically known maximum channelization date (1894 and 1906 CE for Prathoek and Junner Koeland respectively), when the channel reached its final position. For subsequent iterations, migration speed was calculated for each map individually by dividing the difference between the modelled mean ages for its error boundaries (in years) by the error width (in m). This migration rate was then used to convert RMSE (available in Table 1) to years. The converted RMSE (i.e. σ) was used together with the map age (i.e. the mean) as input for the next iteration. This process was repeated until the model outcome converged and remained stable. This required four iterations of the Bayesian model for both investigated meander bends.

Queries were built into the model to generate output for the requested distances, including the meander base (calculated by OxCal using extrapolation of the model; consequently, the uncertainty for this point will be larger) and the spatial boundaries of each error band to allow the calculation of a migration speed for each historical map as described above.

Upon convergence and stabilization of the model outcome, the deposition model of the final (i.e. fourth) iteration was used to interpolate ages throughout the scroll bar deposits. Using the graphical options in OxCal, interpolation results were plotted versus distance from the meander base as the origin toward the meander apex using linear interpolation in between the data points. Plots of the input data were created with the ggplot2 package (Wickham and Chang, 2016) employed in R (version 3.4.0).

We performed sensitivity tests to check model sensitivity to (1) the initial migration rate that was used to convert error bands to years, (2) the correction for channel width, and (3) a rigid versus variable k₀ value. Test outcomes indicated that (1) either increasing or decreasing the initial migration rate by a factor of 2 yielded similar results after several iterations, (2) model outcome with and without using the channel correction differed only by a couple of years, and (3) a rigid or variable k₀ value resulted in age differences of no more than 1 year. Based on these tests we concluded that model sensitivity to these factors is minor.

4 Results

4.1 Historical maps

Figure 4 shows the distortion grid for the map of 1851 CE. As the grid appears virtually straight, geodetic distortions in the map section are minimal. Distortion grids for the other maps were similar (provided in the Supplement).

Table 3 Results of the OSL dating. JK: Junner Koeland; PH: Prathoek; RD: Dutch coordinate system (Rijksdriehoekstelsel). Palaeo-doses and ages are based on the bootstrapped Minimum Age Model for all samples from scroll bar deposits (JK and PH) and the Central Age Model for the terrace sample; n/a: not applicable.

The channel centre lines that were deduced from the maps are shown in Fig. 5. On the map of 1894, Prathoek is already channelized (Fig. 2); therefore, this line is only drawn for Junner Koeland. Based on the current DEM and coring evidence, the georeferenced map dating from 1720 had to be rejected from the Bayesian modelling based on obviously erroneous representation of the river course. The 1720 map depicted the Junner Koeland meander where corings proved that river deposits were lacking, i.e. located too far west by up to 80 m. We tried to include the georeferenced map of 1786 in the Bayesian modelling, but this map resulted in poor model agreement (further discussed below). The correctness of this map was doubted because it displays the meander apex at Prathoek directed southwards instead of southwest like the other channel centre lines (Fig. 5). This would indicate that sedimentation towards the inner bend has taken place in order to move the apex westward, which seems unlikely. Based on this improbable meander planform and the statistical issue created by this map in the Bayesian modelling, we excluded it from further analyses.

Figure 4 Section of the georeferenced historical map of 1851 CE (details listed in Table 1), displaying ground control points (white crosses) and a distortion grid (in white; maze size of 1 km). The grid approaches a perfect grid, indicating that map distortion is minimal. Grid generated with MapAnalyst software by Jenny and Hurni (2011). Image reference historical map: CC-BY Kadaster (2018a).

4.2 OSL ages

The results of the OSL dating and bsMAM calculations are presented in Table 3. The spread in the equivalent doses of single aliquots for the scroll bar samples (example in Fig. 6a) shows that light exposure of the grains prior to deposition and burial was not sufficient to completely reset the quartz OSL signal of all grains. Strikingly, the equivalent dose distribution of many samples shows two distinct peaks: a younger population representing the age of recent transportation and deposition by the river and an older population with an age of about 11 ka, likely representing the depositional age of the grains prior to erosion and transport by the late Holocene Vecht river, as indicated by the OSL age obtained on a terrace remnant at Prathoek (circa 11 ka, Fig. 6b).

Figure 5 Channel centre lines as derived from the historical maps listed in Table 1 and shown in Fig. 2. The DEM (digital elevation model) line represents the course of the meanders preserved upon abandonment at the time of channelization (between 1884 and 1901 CE for Prathoek and between 1894 and 1906 CE for Junner Koeland).

Figure 6 Two radial plot examples showing (a) a scroll bar sample (NCL-2415156) from Junner Koeland and (b) the sample collected in a terrace remnant near Prathoek (NCL-2415163). The scroll bar sample in panel (a) is poorly bleached, showing a large population of grains with a D_e of about 10 Gy. Such equivalent doses are similar to those of the terrace sample in panel (b) (CAM: D_e10.43±0.43 Gy).

Figure 7 Comparison of results before (a, c) and after (b, d) the Bayesian modelling for Junner Koeland (a, b) and Prathoek (c, d). Panels (a), (b) and (c), (d) have equal y axes, and panels (a), (c) and (b), (d) have equal x axes. The legend in panel (a) also applies to panel (c). In panels (a) and (c) the historical map data are indicated with error bars displaying RMSE of the georeferencing procedure (in m); the OSL results of the bootstrapped Minimum Age Model indicate the modelled mean and 1σ error bar (in years). In panels (b) and (d) the blue areas in the Bayesian models indicate the 68.2 % confidence interval. The data point indicating the meander base in panels (b) and (d) is calculated through extrapolation by the Bayesian model. The slope of the graphs in panels (b) and (d) indicates lateral migration rate. Red numbers indicate OSL samples (all should be preceded by NCL-2415). Blue numbers indicate the age of the historical maps (years CE) as used in the Bayesian analysis.

4.3 Reconstructed channel position and chronology

The Bayesian model was iteratively repeated until model results converged to similar output, which occurred after four iterations for both Junner Koeland and Prathoek; therefore, the deposition models that resulted from these iterations were used to interpolate ages throughout the scroll bar deposits. Table 4 lists both the un-modelled and modelled ages. The model agreement index (A value) should be above 60 % to indicate proper model functioning (Bronk Ramsey, 1995). The high A values in Table 4 show that the modelled ages fit the data well. The graphical output of the model as generated by OxCal is displayed on the right in Fig. 7b and d and demonstrates the improvement of the geochronology after the Bayesian modelling compared to the un-modelled data points shown on the left in Fig. 7a and c.

Table 4 Results of the Bayesian model for (a) Junner Koeland, iteration 4, and (b) Prathoek, iteration 4.

OSL sample NCL-2415156 suffered from very poor bleaching (Fig. 6a) and therefore had to be rejected from the Junner Koeland model; otherwise the model A value did not rise above 60 %. In addition, the historical map of 1894 had to excluded, as it resulted in poor agreement that could not be resolved by adjusting the model (e.g. k₀ value). In the fourth (and final) iteration, the model had an A value of 143.8 % (Table 4a). For Prathoek, the poorly bleached sample NCL-2415160 resulted in poor agreement; after excluding this entry, the A value of the fourth iteration amounted to 112.2 %.

The Bayesian model also allows extrapolation of the results towards the meander base (Fig. 7). In our study area the positions of the two meander bases appear to be located close to the former valley side (Fig. 1); therefore, the modelled ages of the meander bases indicate the moment when the meanders sufficiently eroded the valley side to break free from lateral confinement. Model results indicate that this took place during the fourteenth century CE for Prathoek and about 65 years later during the fifteenth century for Junner Koeland. Subsequently, the meanders migrated outside the former river valley, reworking unconsolidated Pleistocene fluvial deposits and coversands, locally with Holocene drift sand on top. This occurred at an average rate of 2.6 and 0.9 m yr⁻¹ for Junner Koeland and Prathoek respectively (calculated by dividing the distance from meander base to apex by the difference between their modelled mean ages). For both meanders, there seems to be an increase in lateral migration rate in the final stages prior to abandonment, roughly between 1825 and 1900 CE.

5 Discussion

Our study demonstrates that the presented Bayesian chronological modelling approach combining OSL and historical map data improves the robustness of fluvial chronologies and can elongate them into longer time frames, compared to what could be achieved by using either data type separately (Fig. 7).

The use of historical maps to determine fluvial migration rates is often challenging due to limited options for georeferencing and uncertainty regarding planimetric accuracy. Historical maps can offer a diachronic perspective on the development of meander shape and channel position (Figs. 2, 5). For our case study this was, however, limited to recent history (after 1720 CE). When going further back in time, the availability of historical maps was limited and map quality was generally too coarse for a detailed reconstruction of meander migration. Doubts regarding the correctness of the maps dating from 1720 and 1786 resulted in excluding these maps from the Bayesian analysis in our study.

With a limited number of GCPs RMSE might underestimate the geospatial error (Esri, 2018; Hughes et al., 2006). When used as a rigid buffer around channel centre lines to determine whether channel changes larger than the assumed error buffer occurred (i.e. deviations that are so large that they cannot emerge from error) (e.g. Urban and Rhoads, 2003; Rhoades et al., 2009), this may result in not detecting small channel changes and could lead to false interpretations of channel stability. Other approaches have been proposed to account for this problem by analysing independent test points or using cross-validation techniques to develop spatially variable error buffers (Hughes et al., 2006; Lea and Legleiter, 2016). However, these methods were developed for aerial imagery where temporal uncertainty is absent and do not convert uncertainty in the spatial domain to uncertainty in the temporal domain. Additionally these methods generally require higher numbers of GCPs than could be registered for georeferencing the historical maps (e.g. 35 per image as used by Lea and Legleiter, 2016). In contrast, the Bayesian modelling methodology that we have presented permits a flexible interpretation of the RMSE buffer, which can easily be calculated even with lower numbers of GCPs, and the iterative procedure with the Bayesian model allows the conversion of the geospatial error into temporal uncertainty. This temporal uncertainty is included in the model as a probability density function, i.e. allowing all channel positions but each with a specific probability of occurrence.

In our analysis there are two potential sources of spatial uncertainty. The first is error in the GPS-measured location of OSL samples (∼5 m). The second is the distance between channel centre line (derived from the historical maps) and scroll bar deposit formed in the inner bend. As the historical maps are likely to display the low-water channel, we made use of the width of the low-water abandoned channel (as depicted on the DEM of the present situation) for deriving the location of scroll bar deposition from the channel centre lines (assuming a constant channel width of ∼30 m). There are indications that bankfull width varied over time (Candel et al., 2018), which would result in an additional spatial error in the location of scroll bar deposition. However, our sensitivity analyses indicate that minor changes in the location of model entries (i.e changing their z value by 15 m) has very limited effect on model outcome (i.e. only a couple of years), indicating model robustness against channel width changes and possible errors due to GPS accuracy of OSL sample locations.

The use of OSL to date young fluvial sediments is challenging, as limited light exposure during the subaqueous transport of grains in a turbid river may lead to an incomplete resetting of the OSL signal in at least part of the grains (Wallinga, 2002b). As demonstrated by Fig. 6 the sediments were indeed affected by poor bleaching, in some cases so extreme that the equivalent dose determined on most aliquots reflected the burial dose of the source material rather than the deposit of interest. The bootstrapped MAM approach yields very large uncertainties in such cases, and the Bayesian analysis indicates that these uncertainty estimates are reasonable. Only for two samples (NCL-2415156 and NCL-2415160, i.e. one from each meander) was the OSL age rejected by the Bayesian model. The Bayesian modelling significantly enhanced the interpretation of the OSL dating results and the robustness of the geochronology.

Through extrapolation, the Bayesian model could estimate the date of meander expansion beyond the former valley sides. From this pivotal moment in time, the meanders were no longer confined and large meanders formed. Our results indicate average lateral migration rates of 2.6 and 0.9 m yr⁻¹ for Junner Koeland and Prathoek respectively. It seems likely that migration rates may have been much higher during brief periods, alternating with periods of relative standstill (Gurnell et al., 1994). However, our chronological data lack the precision and resolution to identify such fluctuations, possibly with the exception of a period with rapid migration in the period 1825–1900 CE.

The methods outlined in this paper provide the chronological constraints that are needed to relate meander expansion beyond the former valley side and changing river dynamics to climate and land-use changes. Such analysis is beyond the scope of this contribution and is reported separately by Candel et al. (2018).

6 Conclusion

Our study demonstrates that combining historical map data and OSL dates in the presented Bayesian modelling approach yields an integrated geochronology with a higher robustness than the sum of its parts. The analysis requires conversion of geospatial error from the historical maps into a temporal uncertainty, and we show that this is possible through an iterative approach. Our method incorporates the strengths of both data types (historical maps and OSL dating) for reconstructing fluvial morphodynamics and supports elongating fluvial sequences into larger time frames. Our study indicates that the meanders of the Overijsselse Vecht expanded beyond the former valley sides from about 1400 CE onwards and that the channels migrated at an average rate of up to 2.6 m year⁻¹ during their lifetime (1400–1900 CE).