Spatial and temporal variations in river terrace formation, preservation, and morphology in the Lower Meuse Valley, The Netherlands

Abstract The Lower Meuse Valley crosses the Roer Valley Rift System and provides an outstanding example of well-preserved late glacial and Holocene river terraces. The formation, preservation, and morphology of these terraces vary due to reach-specific conditions, a phenomenon that has been underappreciated in past studies. A detailed palaeogeographic reconstruction of the terrace series over the full length of the Lower Meuse Valley has been performed. This reconstruction provides improved insight into successive morphological responses to combined climatic and tectonic external forcing, as expressed and preserved in different ways along the river. New field data and data obtained from past studies were integrated using a digital mapping method in GIS. Results show that late glacial river terraces with diverse fluvial styles are best preserved in the Lower Meuse Valley downstream sub-reaches (traversing the Venlo Block and Peel Block), while Holocene terrace remnants are well-developed and preserved in the upstream sub-reaches (traversing the Campine Block and Roer Valley Graben). This reach-to-reach spatial variance in river terrace preservation and morphology can be ascribed to tectonically driven variations in river gradient and subsurface lithology, and to river-driven throughput of sediment supply.


INTRODUCTION
The behavior of fluvial systems is in part controlled by internal processes determining hydraulic gradient, sediment mobility, channel dimensions, and rates of migration (Kleinhans and Van den Berg, 2010) and for the other part by response to external forcing such as climate change, base level fluctuations and tectonic conditions (Vandenberghe, 1995(Vandenberghe, , 2003Holbrook and Schumm, 1999;Blum and Tornqvist, 2000;Wang et al., 2014). For example, erosion by a river occurs whenever the sediment in transport is less than the transport capacity and the remaining stream power exceeds the erosion resistance of the subsurface (e.g., Van Balen et al., 2010). A river terrace can be defined as a patch of higher elevated former floodplain that was abandoned due to incision of the river (Leopold et al., 1964). In tectonically uplifting areas (or areas with another form of base-level lowering) fragments of such produced terraces, with a range of ages, tend to preserve along the sides of valleys. Climatedriven changes (either abrupt or gradual) in precipitation, vegetation cover, sediment supply, or the presence of permafrost can accentuate terrace formation, preservation, and style of dissection (e.g., Antoine et al., 2000;Blum and Tornqvist, 2000;Vandenberghe 2002Vandenberghe , 2003Bridgland and Westaway, 2008;Kleinhans and Van den Berg, 2010).
The Meuse river provides an outstanding example of preservation of different fluvial styles in late glacial and Holocene river terraces in an active rift structure (Roer Valley Rift System, RVRS;Van Balen et al., 2000) and will, therefore, serve as a case study on variations in river terrace formation, preservation, and morphology due to tectonics. Over the last decades, several studies have been performed on the Weichselian and Holocene terraces of the Lower Meuse Valley (LMV) between Maastricht and Nijmegen ( Fig. 1; Pons and Schelling, 1951;Pons, 1954;Van Den Broek and Maarleveld, 1963;Van den Berg, 1989Berendsen et al., 1995;Kasse et al., 1995Kasse et al., , 2005Huisink, 1997Huisink, , 1999Tebbens et al., 1999;Rensink et al., 2015;Isarin et al., 2017). These studies show a general consensus that the youngest terraces in the study area are climate-controlled and hold remnants of cold-stage braided river types, dating from the Weichselian Late Pleniglacial and Younger Dryas, and warm-stage meandering river types of Bølling-Allerød and Holocene age. Most of these existing studies, however, only cover restricted parts of the Meuse valley and/or are restricted to certain time periods (Fig. 2, see section LMV terrace stratigraphic schemes).
An integrated temporal and spatial reconstruction of the geomorphic activity is however, lacking. Such a reconstruction is needed, as a river's response to an external forcing may vary along the river system, owing to both catchment and reach-specific characteristics (Schumm, 1981;Starkel et al., 2007;Erkens et al., 2009). It is likely that these regional and temporal variations in fluvial response are, in turn, represented in the formation, preservation, and morphology of river terraces, which has been insufficiently considered in the reconstructions produced so far.
In this paper we, therefore, compile and synthesize existing information on dated features such as meander-scar fills, point bars, and aeolian dunes acquired from published and unpublished data, into a GIS database. In addition, the database is extended with the wealth of new data, mainly from archaeological exploration, which has become available in the past years. Furthermore, new palynological interpretations and absolute dating of abandoned channel fills are provided for an area of specific interest. The results are presented as a terrace map as well as a series of palaeogeographic maps for the whole of the LMV. This detailed palaeogeographic reconstruction of the Lower Meuse Valley will then provide the basis for improved insight into regional and temporal variations in river terrace formation, preservation, and morphology.

549
Spatial and temporal variations in the Lower Meuse Valley river terraces, the Netherlands SETTING The~900-km-long rainfed river Meuse has a catchment size of 33,000 km 2 . The Meuse has its headwaters in north-eastern France. From here it passes through the north-eastern part of the Paris Basin and the Ardennes Massif, where it is incised into bedrock, before entering the alluvial reach through the RVRS in the border region of Belgium and The Netherlands (Fig. 1). This active rift system is located in the southern Netherlands, north-eastern Belgium and adjacent parts of Germany . In the LMV, the annual mean discharge of the Meuse is approximately 250 m 3 / s. The bankfull discharge is~1500 m 3 /s at Maaseik (Belgium; Paulissen, 1973;Bogaart et al., 2003;. This study focuses on the segment of the Meuse (the LMV) between Maastricht and Nijmegen. This area provides the opportunity to study a well-preserved terrace morphology with light detection and ranging (LIDAR) data. In addition, the Meuse crosses several fault-bounded blocks of the RVRS within the LMV (Fig. 1), of which the tectonic activity is well-known (Geluk et al., 1994;Van den Berg, 1996;Houtgast and Van Balen, 2000;Cohen et al., 2002;Houtgast et al., 2002;Michon et al., 2003;Michon and Van Balen, 2005;Van Balen et al., 2005), making it a suitable area for studying the interaction between neotectonics and fluvial morphology. Several modest tributaries (Jeker, Geul, Geleenbeek, Rur, and Swalm) join the Meuse in the LMV, of which the Rur is the largest. A confluence between the Meuse and the Niers-Rhine course was present in the north of the study area (near Gennep, Fig. 1) during the last glacial maximum (LGM) and into the Weichselian late glacial ( Van de Meene and Zagwijn, 1978;Huisink, 1999), which was presumably abandoned during the Younger Dryas/Early Holocene transition (Kasse et al., 2005;Janssens et al., 2012). The study area extends a little further north to enter regions where continuations of the Meuse terraces lay buried below younger deposits of the Rhine-Meuse delta Cohen et al., 2002Cohen et al., , 2009Cohen et al., , 2012Hijma et al., 2009).

Tectonical and sedimentary setting
The horst-graben structure of the RVRS comprises the Campine Block (CB), Roer Valley Graben (RVG), Peel  Block (PB), and Venlo Block (VB), respectively (Fig. 3). The CB and PB are expressed as topographic highs while the RVG and VB are lows within the topography of the region ( Fig. 1 and 3). The Feldbiss Fault Zone (FFZ), consisting of the Heerlerheide (HH), Geleen (GF), and Feldbiss (FF) faults, separates the subsiding RVG from the CB in the south, whereas the Peelboundary Fault Zone (PBFZ) forms the boundary between the RVG and PB in the north. The Koningsbosch (KB) and Beegden (BE) faults delineate a small horst within the subsiding RVG. The boundary between the PB and the relative subsiding VB is situated along the Tegelen Fault Zone (TFZ; Fig. 3). The present-day extension phase of the RVRS started during the Late Oligocene and experienced a change of the extension direction from a WNW-ESE to a NE-SW direction since the Early Miocene (Geluk et al., 1994;Van Balen et al., 2002a, 2002bMichon et al., 2003). Geluk et al. (1994) indicate that the central RVG subsided~1000-1200 m while the PB experienced no more than 200 m of subsidence. The Quaternary subsidence history of the RVRS was studied in detail by Houtgast and Van Balen (2000), who concluded that the average subsidence rate of the RVG was~60-90 mm/ka during the Quaternary, whereas the PB and CB subsided with a rate of 46 mm/ka and 27 mm/ka, respectively. The displacement rate along the FFZ showed to be around 41-71 mm/ka since the middle and late Pleistocene, while the displacement along the PBFZ ranged between 50-78 mm/ka since the Late Pleniglacial Cohen et al., 2002;Michon and Van Balen, 2005).
The uplift of the Ardennes, upstream of the RVRS, has resulted in a relatively complete fluvial terrace sequence of the last 5-10 glacial cycles. Moreover, within the RVRS, active differential subsidence has caused terraces and buried terraces to be preserved over the same time span (Van den Berg, 1996;Van Balen et al., 2000). In the LMV study area, this resulted in marked differences in the lithology of the subsurface where the river crosses different terrace bodies preserved at different tectonic blocks (Fig. 3). The longitudinal profile of the (Holocene) LMV, which is characterized by an upstream part with a relatively high valley gradient (~0.5 m/km) and a flatter downstream part with a gradient of~0.15 m/km, shows a knick in slope that roughly coincides with a change in the subsurface lithology.
Fine marine and coastal deposits of Miocene age are present in the shallow subsurface of the PB and northern part of the CB (Fig. 3). Within the RVG these deposits are found at a depth of over 200 m. The shallow subsurface of the RVG consist primarily of coarse fluvial deposits of Rhine and Meuse. The shallow subsurface of the VB, and to a lesser extent the northern most part of the CB, is composed of fine fluvial sediments. These fine-grained fluvial sediments are, however, absent on the PB (Fig. 3). Finally, a thin cover (~0-15 m) of coarse fluvial sediments of the Meuse is found over the whole course of the LMV, which partly corresponds to the active layer of the river bed (~1 m).

Climate change and vegetation development
Vegetation is a major controlling factor in river activity and morphology, as it influences bank/bar stability, sediment load, and discharge ;

551
Spatial and temporal variations in the Lower Meuse Valley river terraces, the Netherlands Vandenberghe, 2003). Furthermore, vegetation cover changes influence the availability of aeolian sediments on a local scale. Additionally, biostratigraphy can be used as a relative dating tool when the possibility for absolute dating is lacking. For this study, the regional late glacial-Early Holocene vegetation history for the northern LMV (Hoek et al., 2017), together with the biostratigraphy of a representative Holocene site (i.e., Well-Aijen; Bos and Zuidhoff, 2015), were used as a reference for (re)interpretation of palynological data. Figure 4 shows an overview of the late glacial and Holocene vegetation development in the LMV. Since vegetation development is linked to climate change, estimations of the mean annual temperature (˚C) and total annual precipitation (mm/yr) are given in Figure 4 for each "time-slice" of the late glacial and the Holocene based on Bogaart and Van Balen (2000). Figure 4 also gives an indication of the mean annual peak flood discharge (Q af ), as such bankfull discharges are assumed to be responsible for most of the morphodynamic changes in a river system. Modelled discharges based on GCM output and data on land cover change showed a 6.6% increase in mean annual discharge and a decrease in recurrence time of high flow events (Q > 3000 m 3 s -1 ) for the Late Holocene (compared to the Middle Holocene; Ward et al., 2008). Deforestation and the start of large-scale agriculture are indicated to be responsible for the modelled change in discharge regime and a threefold increase in suspended sediment yield (Ward et al., 2008). . Overview of the regional late glacial (after Hoek et al., 2017) and Holocene (after Bos and Zuidhoff, 2015) vegetation history of Lower Meuse Valley together with the reconstructed mean annual temperature, T (˚C); total annual precipitation, P (mm/yr); total annual evapotranspiration, E (mm/yr); and mean annual flood discharge, Q af (m 3 /s) (after Bogaart and Van Balen, 2000). Late glacial biochronostratigrahy after  and Holocene after Zagwijn (1986).

Vegetation changes during the Weichselian late glacial and Holocene
Pollen analyses show that a sparse vegetation cover was present during the Late Pleniglacial (~19.0-14.7 ka), although many places were barren due to the low mean annual temperatures and total precipitation. This period is characterized by deposition of aeolian cover sands interbedded with deflation horizons (Kasse et al., 2007). Due to a rise in temperature, a dwarf shrub and herbaceous plant cover could develop during the earliest phase of the late glacial, which can be considered as the transition from tundra towards shrub-tundra (Hoek and Bohncke, 2002). As a result of continued warming and increased precipitation during the Bølling (~14.7-14.0 ka), a predominantly open landscape with birch shrubs and birch trees developed (Hoek and Bohncke, 2002). The subsequent short Older Dryas period (~14.0-13.9 ka) is characterized by a relative open landscape with sparse vegetation, which locally led to renewed aeolian deflation and deposition. Although temperatures and precipitation were supposedly not as high as during the Bølling, a rather open Betula forest developed during the Betula phase of the Allerød (~13.9-12.85 ka) after which Pinus started to dominate and the landscape changed into a pine forest during the second phase of the Allerød. As temperature decreased at the start of the Younger Dryas (~12.85-11.7 ka) the vegetation cover opened, leading to a steppe-like landscape (Hoek et al., 2017). During this period, the characteristic aeolian river dunes were formed along the eastern bank of the Meuse. Due to the colder climate and decreased vegetation cover, peak discharges increased substantially during the Younger Dryas. The transition from the late glacial into the Holocene is characterized by a strong increase in temperature and precipitation. This climatic "improvement" resulted in the gradual development of a dense birch forest during the first part of the Preboreal after which it was replaced by a pine forest during the last phase of the Preboreal (~11.7-10.8 ka). The succeeding Boreal (~10.8-8.9 ka) is characterized by the immigration of thermophilous deciduous trees and shrubs (Fig. 4). During the course of the Boreal, the pine forest developed into a relatively open mixed deciduous forest (Bos and Zuidhoff, 2015). The expansion of this mixed deciduous forest continued during the Atlantic (~8.9-5.8 ka). During the latter part of the Atlantic, the landscape gradually opened due to the first occurrence of farmers in the Meuse valley (Bos and Zuidhoff, 2015). The occurrence of Fagus is characteristic for the Subboreal (~5.8-3.1 ka). The start of the Subatlantic (~3.1-0 ka) still experienced a relatively closed vegetation cover, although expending agriculture led to a more open landscape with an increase of grasses, heathers, and cereals later on during the Subatlantic (Bos and Zuidhoff, 2015).

Existing Lower Meuse Valley terrace stratigraphic schemes
Several terrace stratigraphic schemes have been proposed during the last decades for the Weichselian and Holocene terraces of the LMV (Fig. 2). The northern LMV (between Venlo and Nijmegen) is the most intensively studied reach due to the wellpreserved terrace fragments, with clearly recognizable fluvial styles. Terrace maps based on soil development were made by Schelling (1951), Pons and Schelling (1951), Pons (1957), and Miedema (1987). Detailed geomorphological maps were produced by Buitenhuis and Wolfert (1988) and Wolfert and De Lange (1990). Kasse et al. (1995) and Huisink (1997) discriminated terrace levels using a combination of elevation, morphology, sedimentology, and palynological data. Tebbens et al. (1999) based their subdivision on morphology and a large number of radiocarbon dates (40) of basal residual channel infillings. Additionally, Berendsen et al. (1995) introduced a division of the buried terrace surfaces in the most northern part of the study area, based on a set of 14,000 hand borings and dated basal organic channel fills, which was expanded upon by Cohen et al. (2002) and Cohen (2003). These studies discriminated late glacial and Early Holocene terraces and channel belts based on elevation, buried-surface morphology and pedogenesis, architecture, lithology, pollen, and 14 C-dating.
The middle LMV (between Grevenbicht and Venlo), which comprises the RVG and PB, was studied by Van den Broek and Maarleveld (1963) and Van den Berg and Schwan (1996). They based their terrace scheme on differences in elevation, fluvial style, and soil development.
The Dutch part of the southern LMV (between Maastricht and Grevenbicht) has been studied by Van den Berg (1989). He established his results on relief, morphology, and terrace correlation based on gradient lines. Later on, Houtgast et al. (2002) partly updated the terrace stratigraphy for the Meuse terraces around the FFZ. The terraces in the Belgian part of the southern LMV were investigated by Paulissen (1973). In this work, the terraces are discriminated based on a combination of morphological and geological mapping, radiocarbon dating, and sediment analysis.
Thus, although several studies are available for the LMV, an integrated temporal and spatial reconstruction of the Lower Meuse geomorphological activity since the Weichselian last glacial maximum is lacking.
Recently, Rensink et al. (2015) and Isarin et al. (2017) produced a map of terraces with ages starting from the Younger Dryas. Their map covers the entire LMV and is based on geoarcheological investigations, morphology (digital elevation model), and absolute and relative dating. However, due to the archaeological focus, their map is partly inconsistent with fluvial morphologic principles. Moreover, Late Pleniglacial and early late glacial terraces are missing from their palaeogeographical reconstructions.
This study presents an internally consistent and detailed palaeogeographic reconstruction of the whole LMV, based on stratigraphic principles and genetic points of view. Our terrace map starts at the Weichselian last glacial maximum and thus includes several important climatic transitions. Thus, in comparison to Rensink et al. (2015), our terrace map focuses on dynamic fluvial evolution and covers a longer time span. In addition, it is also based on more and partly new data. We use the reconstructions in order to address regional 553 Spatial and temporal variations in the Lower Meuse Valley river terraces, the Netherlands differences in responses to external forcings within a river system, which has been insufficiently considered up to now despite the well-known tectonic activity.

Lithological cross sections and cores
Transects were hand cored over abandoned channels at selected locations, identified in LIDAR imagery (AHN2, 2017). The transects served to find the most suitable sampling locations for radiocarbon dating and/or pollen analysis on abandoned-channel fills. Reconnaissance coring used a soil auger, a 20-mm diameter gouge, and occasionally a suction corer to retrieve water logged sandy sediments from below the groundwater table. Sediments were described in the field for lithology, carbonate content, color, grainsize (NEN 5104), and presence of plant remains. At sites of thickest infill, sediment cores were taken in meter-long sections with a modified Livingstone piston sampler and, in case of sandy or compacted peat, with a 60-mm diameter gouge. Cores were wrapped in plastic foil, transported, and stored in a cool room. Sediment cores Maaseik and Oude Hoeve Meeswijk were subsampled in the field because of the fragile state of the cores.

Radiocarbon dating
Materials for radiocarbon dating were picked from split core halves on inspection of the sediments in the laboratory, targeting botanic macrofossils. Selected subsamples (slices of core material, 1-2 cm thick) were treated with diluted potassium hydroxide (peat samples) or disodium pyrophosphate (clayey samples), to disperse organic compounds or clay bonds, respectively. Subsequently the samples were wet sieved over a 200 µm mesh. Terrestrial macro remains were extracted from the remaining sample material by use of a binocular macroscope and identified to species level. Table 1 shows the sample depth, sampled material, and the macro remains selected for radiocarbon dating. The macro remains were (AMS) radiocarbon dated at the Beta Analytic laboratory (Miami, Florida).

Digital Mapping
Digital mapping was performed making use of a GIS system originally developed for the Rhine-Meuse delta. The system splits mapping of geomorphological elements (individual meander scars, terrace fragments) from cataloguing naming, dating, and literature referencing results for individual elements (Berendsen et al., , 2007Cohen et al., 2012). The coverage of the GIS database was extended to include the full study area. Existing information on dated geomorphological features, such as fine-grained meander-scar fills, point bars, and aeolian dunes, was acquired from published sources and institutional databases (i.e., from journal papers, archaeological reports, and geological survey mapping activities) and made browsable in GIS (Fig. 5A). Furthermore, geomorphological datasets of continuous cover (e.g., AHN2 [2017] LIDAR data, Fig. 5B), lithological databases (cores and outcrop descriptions; from new fieldwork in this study, but also abundant pre-existing data from institutional databases), and earlier mappings (geological, geomorphological, pedological, archaeological, and historical) were collected. Using GIS visualization and geological analysis, this dataset was manually (re)interpreted to the geomorphological and geological maps featured in this paper. These maps honor the dating information as well as the geomorphic crosscutting relations (Fig. 5C). The embedded geological reconstruction is described in the further sections of this paper. The morphostratigraphic correlation decisions at the level of individual meander/terrace fragment elements (polygons in the GIS layer, newly digitized), are described in the catalogue part of the database (Fig. 5D). Linkage of GIS polygons and dating information stored in the catalogue (using channel belt and terrace fragment identification numbers [#IDs]), allows quick generation of a terrace or age map (Fig. 5E) as well as time-slice map series of the palaeogeographic evolution ( Fig. 5F). In this way, the GIS becomes a tool for iterative partial improving of the reconstruction (Berendsen et al., 2007;Cohen et al., 2012; and not just a storage and visualization system. In order to produce internally consistent geologicalgeomorphological maps and palaeogeographic reconstructions, it is essential to assign an age to each terrace remnant in an architecturally correct manner. In this study, we map and date terrace morphology from successive stages of river abandonment. The final stage of river bed sedimentary activity is therefore considered as the moment that produces the "end date" of a channel. In other words, the last stage in which flow through the channel is capable of bedload transport, is considered to be the last moment at which a channel was morphologically active and featured in the terrace age maps. Ages obtained from organics at the base of a finegrained channel infill, for our study, provide a minimum age of channel activity (a terminus ante quem date for this activity) and are therefore, in most cases, assumed to be the end of river bed activity. In the same philosophy (dated) overbank deposits provide a minimum age for the underlying riverbed sediments. Terrace remnants that are not directly dated are nevertheless attributed with an "end date" based on correlation (i.e., height, morphology, lithology, and historical maps) to locations where dates are obtained. Due to incision of the Meuse since the Late Pleniglacial, relative elevation at local scale provides a means for such correlation. It should, however, be emphasized that correlation of intra-Holocene terraces is, in some cases, difficult based on height differentiation alone. In these cases, other correlation methods (e.g., lithology, spatial continuity of geomorphological features, and historical maps) become more important in the correlation of individual terrace fragments. The mapping uses a customized color legend consisting of 12 age classes for the staircase of Meuse terraces. One unit covers pre-LGM time, eight units LGM to Early Holocene, and three units the rest of the Holocene period. The legend uses age boundaries that for the late glacial follow Hoek et al. (2017) and for the Holocene subdivision follow Zagwijn (1986; IntCal13 recalibrated). In parts of the mapping the dating resolution is still suboptimal and we therefore used broader Holocene age-classes, grouping them into an early-(Preboreal and Boreal; ca. 11.7-8.9 ka), middle-(Atlantic and Subboreal, ca. 8.9-3.1 ka), and Late Holocene (Subatlantic, youngest 3.1 ka) age class. It should be stressed that whenever the adjectives "early," "middle," and "late" are used for the Holocene in this paper they refer to the subdivisions indicated above. The capitalized form "Early," "Middle," and "Late" Holocene is used whenever we refer to the formal subdivision of the Holocene Series/Epoch as defined by Walker et al. (2012).
The map time series (Fig. 5F) distinguishes active channel belts and abandoned terrace fragments, and uses shading to indicate "later reworked" from "preserved up to present" (as in  and Cohen et al. [2012]). This palaeogeographical reconstruction output allows us to calculate (and conceptually optimize through iterative editing) the spatial extent of the active channel belt as it varied through space and time. It is emphasized here that the map series displays the eventually reached width at the end of activity as catalogued for each discriminated unit (the above #IDs). Owing to lateral migration activity in phases prior to abandonment, this displays a somewhat larger width than the river channel would have had at any instantaneous moment (i.e., channel belt width > active channel width). Both the map of the

555
Spatial and temporal variations in the Lower Meuse Valley river terraces, the Netherlands LMV Weichselian and Holocene terraces, as well as the map time series are projected in the Dutch national grid (Amersfoort/RD New). The GIS database and associated age and palaeogeographic maps are freely available for download at: https://doi.org/10.17026/dans-xkk-f29b.

Database of sites with age control
A first aim of this study is to synthesize previous studies and create an overview of the status of data availability regarding the fluvial geomorphology along the Lower Meuse. At present, our database consists of 381 dates from 255 unique site locations, which include the four new sites presented in this study. The distribution of dating methods, as well as the dated time periods, are given in Figure 6. The majority of these dates (71%) are derived from channel fills and had direct dating of channel abandonment and accumulation of infill as their original purpose. Other dates were performed on, amongst others, aeolian dunes, point-bar deposits, bedload deposits, levees, and floodplain deposits.
Of all the channel-fill datings, 44% targeted the basal channel infill, whereas the rest was taken from higher up in the fill sequence (often in addition to sampling the very base). The reason for this, in some cases, has been the difficulty to pick suitable datable (organic) materials from channel fills that trapped clastic deposits in their initial filling phases. Figure 6 shows that only a small percentage of the dates is of (Late) Pleniglacial age. This might be due to a lack of datable organic material, or due to the fact that a large part of the (Late) Pleniglacial channels was reworked later on. The relatively high amount of Atlantic dates can be ascribed to the extensive peat formation, and therefore available climate archives, in many of the (upper) residual channel fills in the study area during this time (Cohen, 2003;Zuidhoff and Huizer, 2015). As over one-third of the dates were performed in an archaeological context, a peak in Subatlantic dates is to be expected as well. Although already a large amount of data is available within the current database, it is designed so that new data/sites can be easily added and is, therefore, a starting point rather than a finished product.
Newly dated channel fill sites Cross sections were made over four newly investigated paleochannels and their (basal) infill was AMS 14 C dated ( Fig. 7A; Table 1 and Supplementary Information), as these channels are positioned in an area of specific interest, i.e., around the Feldbiss Fault zone, which is suspected to have a pronounced effect on the fluvial dynamics of the Meuse. Moreover, only few absolute dates were available for this region. The Vijverbroek meander was selected because it is a very distinct meander at the outer rim of the Holocene plain. It has laterally eroded the Pleniglacial terrace and could therefore provide an indication of the start of lateral erosion by Holocene channels in this region. Two of the four channels, the Wurfeld and Maaseik meanders (Fig. 7B), provided radiocarbon ages of 9720 ± 40 14 C yr BP (Beta-454392) and 9160 ± 30 14 C yr BP (Beta-454391), indicating that these channels were abandoned during the Preboreal. The dates for the Wurfeld and Maaseik meanders also indicate the moment of the abandonment by avulsion of the meander belt of which they are part. The avulsive abandonment of the channels is supported by the lithology of the infillings. The lithological cross sections show that the Wurfeld meander, which is located upstream (Fig. 7), has a stage of clastic infill prior to the peat formation in the abandoned channel. This clastic infill indicates narrowing of the channel before a full cutoff, most likely representing the transitional stage of active abandonment as described by Toonen et al. (2012) for avulsion-abandoned channels. The fact that the Maaseik residual channel does not show such a distinct clastic channel infill may be ascribed to the more downstream position of this meander to the avulsion

Terrace map
The map of the LMV Weichselian and Holocene terraces (Fig. 8) shows the preservation of terrace remnants by age of abandonment in 14 C yr BP. The boundaries of the age classes were also calibrated to cal yr BP (IntCal13) to enhance correlation with e.g., archaeological periods. Four AHN2 (2017) LIDAR-based topographic profiles were drawn over the course of the LMV. These are situated on the four different tectonic blocks to illustrate height differences between the successive terraces ( Fig. 9). To illustrate the contrasts in fluvial style of the terrace fragments over different tectonic blocks along the course of the LMV, we highlight representative subareas in Figure 10. The classification of channel planform is based on Leopold and Wolman (1957), Nanson and Knighton (1996), and Kleinhans and Van den Berg (2010).

Topographic profiles
Overall, the topographic profiles ( Fig. 9) indicate that the height differentiation between the Weichselian Pleniglacial and late glacial terraces (Δ~2-4 m) is greater than that between the different Holocene paleofloodplains (Δ~1-1.5m). The profiles also show that the height differentiation between the successive terrace remnants is more pronounced on the VB (profile A-A') and especially the PB (profile B-B') compared to the RVG (profile C-C') and CB (profile D-D'), where differentiation between successive terrace remnants in the latter two may be as small as~0.5 m.

General overview
The terrace map shows that late glacial terrace remnants are well-preserved on the VB and PB. Holocene terrace remnants are, on the other hand, best preserved in the south, on the CB and in the RVG. The overall trend on the VB and PB is that of an incising and narrowing river plain, preserving a set of successive river terraces ( Fig. 8 and 9). An exception to this is the ca. 1000-yr-long Allerød period, during which individual meanders seem to have widened (by lateral migration) compared to Bølling precursors (Fig. 8). The CB and RVG show 557 Spatial and temporal variations in the Lower Meuse Valley river terraces, the Netherlands a less complete preservation of successive river plains and late glacial terrace remnants are scarce. A broad amalgamation of middle and late Holocene meanders is present (Fig. 8), especially where the river enters the RVG from the CB. The average width of the complete Holocene river plain is approximately 1400 m on the VB ( Fig. 9; profile A-A') and 550 m on the PB (profile B-B') compared to~3000 m in the RVG (profile C-C') and~2000 m on the CB (profile D-D').

Description of terrace remnants
Terrace remnants of LGM and Late Pleniglacial age are preserved throughout the LMV (Fig. 8). The LGM and Late Pleniglacial terrace fragments show a morphology of many low-sinuous and shallow channels separated by multiple bars, indicating a braided river system. The maximum width of the LGM and Late Pleniglacial braidplain varies over the course of the LMV, from~6000 m across the CB to 18,000, 8000 and 12,500 m in the RVG and on the PB and VB, respectively (Fig.  11A). Thus, starting at the CB, the braidplain first widens in the RVG, then narrows over the PB, before diverging again on the VB (Fig. 8). This suggests that the uplifted PB block acted as a bottleneck and choked the LGM braidplain width. Preservation of the LGM and Late Pleniglacial system is the same for the different tectonic blocks, except for the CB in the south. This is likely because here the LGM and Late Pleniglacial braidplain is bound by bedrock and by relatively coarse-grained Pre-Weichselian fluvial deposits. This limited valley width caused  efficient reworking by the Meuse in the subsequent late glacial and Holocene (Fig. 8).
Preservation of the Bølling system is variable over the LMV. Bølling terrace remnants are present on the VB, PB, and in the RVG (Fig. 8). The terrace remnants show a lowsinuosity, multi-channel planform (Fig. 10O-Q). The maximum width of the Bølling plain varies from~7000 m in the RVG to~3750 and~5000 m on the PB and VB, again showing the bottleneck of the uplifting PB. Furthermore, the river system is narrower compared to the LGM and Late Pleniglacial braided plain (Fig. 11B). It is striking that no terrace remnants of Bølling age could be distinguished at the CB.
Allerød terrace remnants are primarily located on the VB and PB (Fig. 8) and contain a high-sinuosity meandering planform with point-bar ridges and swales ( Fig. 10L and M). The lateral extent of the Allerød plain varies betweeñ 3000 m and~5000 m on the PB and VB, respectively. Lateral displacement of the Allerød meanders resulted in reworking of parts of the LGM, Late Pleniglacial, and Bølling terrace remnants (Fig. 11C). The CB holds only one Allerød terrace fragment of a large, low-sinuosity meander, preserved in the very south (Fig. 10N). In the RVG Allerød-age meanders are only present in the north, bordering the PBFZ (Fig. 12A). This hampers reconstruction of the width of the Allerød meander belt (as presumed for this time) in a critical part of the study area that is the CB-RVG transition (Fig. 12C). The Allerød meander complex bordering the PBFZ (in part overbuild by the city of Roermond) has an anomalous (extremely high) sinuosity and scroll morphology. It, furthermore, shows upstream migration while all other Allerød meander forms show downstream migration.
Remnants of the Meuse river system in the Younger Dryas are present throughout the LMV. However, preservation of the Younger Dryas system is marginal in the northern part of the CB and in the central RVG compared to the rest of the LMV. Terrace fragments of the Younger Dryas (Fig. 10H-K) show a morphology with multiple channels and bars, indicating a multi-channel to a braided river plain. At the transition between the CB and the RVG (i.e., around the FFZ), however, the palynologically dated Younger Dryas terrace fragments of Kingbeekdal, Dukkelaar, and Korbusch exhibit a meandering rather than a braided planform ( Fig. 12; CD). The Younger Dryas braid belt width maximizes to~5000 m on the CB, widens a bit to ~5700 m in the RVG, then narrows down to~2000 m on the PB, after which it modestly widens again to~3400 m on the VB. The Younger Dryas braidplain thus was significantly wider on the CB and in the RVG than across the PB and VB. On the PB and VB, the Younger Dryas river incised and narrowed its river plain compared to Allerød meandering and earlier braided stages ( Fig. 9 and 11D).
Early Holocene terrace remnants are primarily observed on the VB. Their preservation is minimal on the CB and in the RVG. Remnants of this age are absent on the PB. Planform fluvial style of the early Holocene terrace patches differs between the VB and the CB-RVG transition zone (Fig. 10E-G). While in the RVG high-sinuosity meandering channels are present on the terrace fragments ( Fig. 10F and 12C), the terrace remnants in the VB and CB show a multi-channel low-sinuosity planform (Fig. 10E, 10G, 11E). It should be noted that in the RVG Early Holocene channels of true meandering style have only been found just downstream of the FFZ. The width of the Early Holocene channel belts is at least~1500 m on the CB (based on preserved remnants). A maximum for the estimated width is~5000 m (similar to the Younger Dryas) if channel belt width of the Middle Holocene meander belt is assumed to be representative for the Early Holocene. The width of the Early Holocene river plain increases up to~6000 m in the RVG. Further downstream it reduces to a maximum width of~500 m on the PB (i.e., the space in-between preserved Younger Dryas fragments). The width of the Early Holocene plain increases to a maximum of~2200 m on the VB. As in previous time periods, this indicates that the PB is a bottleneck in the width of the channel belt in the Early Holocene.
Middle Holocene terrace remnants have primarily been preserved on the CB and in the RVG (Fig. 8). On the PB and VB, they are present as part of a point-bar series. Middle Holocene terrace remnants have a low-sinuosity singlechannel planform at the VB and PB (Fig. 11F), while the terrace fragments in the RVG contain remnants of a highly sinuous meandering channel with lateral channel migration and crosscut channels (Fig. 10A-D). At the CB, the Middle (and Late) Holocene terraces comprise remnants indicating both a low-sinuosity planform (southern CB), with a narrow plain of~700 m, as well as a high-sinuosity planform   561 Spatial and temporal variations in the Lower Meuse Valley river terraces, the Netherlands (northern CB) that has a maximum lateral extent of~4500 m. Several small Pleniglacial and Younger Dryas terrace fragments are present within the Middle (and Late) Holocene river system of the CB. The maximum lateral extent of the Middle Holocene system is~4000 m in the RVG. On the PB it is~600 m and on the VB~700 m.
Late Holocene terrace remnants are well-preserved over the entire LMV. On the VB, PB, and southern CB, a low-sinuosity to straight river pattern is observed for the late Holocene Meuse. The RVG and northern CB, on the other hand, show a high-sinuosity river pattern during the late Holocene. The maximum width of the Late Holocene river plain varies from~1800 m on the CB and~2200 m in the RVG, to~600 m on the PB and~1100 m on the VB. So, while the maximum lateral extent of the river plain is decreasing on the CB and in the RVG (compared to the  preceding Middle Holocene), the width on the PB remains almost the same and is somewhat increasing on the VB.

DISCUSSION
The general consensus on the LMV Huisink, 1997;Tebbens et al., 1999) is that the (abrupt) climate and associated vegetation changes in northwestern Europe since the LGM were the steering forces behind the differentiation in fluvial style of successive fluvial terraces (reconfirmed and further specified in the current study). This was also identified in studies on the rivers Warta and Vistula in Poland and the Seine in France (Vandenberghe et al., 1994;Roblin-Jouve, 1995;Starkel, 1995). As noted earlier in such studies, however, the river morphology may vary for each valley reach (Vistula [Starkel et al., 2007] and the Upper and Lower Rhine [Erkens et al., , 2011), which is also the case in the LMV. Although climatic forcing was uniform for the entire LMV, spatial variation exists in the morphological planform and preservation of the various stages of late glacial and Holocene valley development. Figure 13 summarizes how various reaches of the Meuse river adjusted their morphological planform to the same climate forcing, depending on their tectonic setting (i.e., relative uplift or subsidence), subsurface composition, and upstream-downstream propagation of morphological change. Below, we first review the observed climate-driven changes in fluvial morphology through time.
Hereafter, we will discuss the spatial variations observed within the LMV terrace remnants.

Temporal variations in river terrace morphology
The wide braided river system of the LGM and Late Pleniglacial (Fig. 13) can be explained by a scarce vegetation cover combined with a high sediment load and peak (spring) discharges (due to a nival discharge regime originating from a hinterland with permafrost-affected subsoil). Peak discharges flooded the braided river plain, depositing sediment as channel bars. During falling and low stages, water concentrated in numerous small and shallow, low-sinuosity erosive channels separated by bars/islands . The main channel(s) of this kind of fluvial system are not preserved because they were used as the starting channel for the next lower terrace. Their positions, however, can be Warming and increased precipitation during the Bølling period (Vandenberghe and Bohncke, 1985;Renssen and Isarin, 2001) led to increasing vegetation cover and a more mixed discharge regime that includes rain (storm) discharge peaks. The developing vegetation cover (e.g., Betula, Juniperus, and dwarf scrubs) on the LGM and Late Pleniglacial braidplain stabilized the LMV river banks, while upstream in the Ardennes, establishing vegetation along the valley edges increasingly trapped hillslope delivered sediments. In the LMV, the regime and sediment-supply changes resulted in a concentration of discharge into a decreasing number of channels. At times of peak discharge, these channels inundated the floodplains. Along the edge of the channels, modestly elevated natural levees were formed (with sandy loamy textures), presumably due to the trapping of sediment by establishing vegetation on the banks. Away from the channels, flood sediments were trapped, blanketing abandoned braidplain terraces with floodplain loam. Contraction of discharge in a fewer number of channels raised the stream power in these channels, causing them to incise. The loam cover stabilized the floodplain and river banks, helping to deepen the thalweg of the remaining channels.
Over time, the more regular discharge regime, increasing vegetation cover, and soil development drove further fluvial change, causing the abandonment of the low-sinuosity multichannel system during the Bølling, in favor of a lesser number of larger and higher sinuosity meandering system in the Allerød ( Fig. 13; Kasse et al., 1995;Hoek et al., 2017). This transition is probably a gradual development. Over time, the beginning meandering channels from the Bølling period became increasingly more sinuous. As a result, meanders that became abandoned in the Allerød have a higher sinuosity than those already abandoned in the Bølling period. In this way, the LMV river system developed from a braided to a multi-channel low-sinuosity planform, to a relatively highsinuosity single-channel planform over the Late Pleniglacial to Bølling to Allerød periods (Fig. 13), as concluded previously by Kasse et al. (1995), Makaske and Nap (1995), Huisink (1997), and from sites to the north of the study area by Cohen (2003). The system transition from Late Pleniglacial braided to highly-sinuous meandering planform during the Allerød is observed in lowland rivers throughout Europe (Antoine, 1993;Vandenberghe et al., 1994;Brown, 1995;Gabris, 1995;Kalicki, 1995;Roblin-Jouve, 1995;Rose, 1995;Starkel, 1995;Vandenberghe, 2002;Erkens et al., 2009Erkens et al., , 2011Janssens et al., 2012). As concluded in general terms by Schumm (1973Schumm ( , 1981, such gradual morphological transitions in lowland alluvial valleys might well have occurred faster in some reaches than in others, i.e., the main climatic change at the onset of the Bølling can be considered as the external forcing trigger, and the resulting geomorphological changes as a complex response. Climatic cooling at the onset of the Younger Dryas resulted in a decreased vegetation cover as well as the reestablishment of discontinuous permafrost, both in the study area (Bohncke et al., 1993;Kasse et al., 1995Kasse et al., , 2005Hoek et al., 2017) and in the hinterland (Vandenberghe and Pissart, 1993). As a consequence, spring discharges and sediment transport capacity increased, leading to the observed shift into a braided fluvial system during the Younger Dryas (Fig. 13). However, it is known from the Warta (Poland) and Tisza (Hungary) rivers, for example, that such a shift into a braided planform is not self-evident, as these rivers retained their meandering planform during the Younger Dryas stade (Vandenberghe et al., 1994;Kasse et al., 2010). The observed width of the braidplain during the~1200-yr-long Younger Dryas period is considerably less (~3 times) than that of the LGM and (Late) Pleniglacial (Fig. 8). This might be due to increased floodplain/bank strength by deposition of overbank loams during the Bølling/Allerød or by the incised position of the river during the start of the Younger Dryas, increasing the volume of sediment that needs to be eroded for the river to extent its plain laterally. In general, only the outer edges of the Younger Dryas river plain, which predominantly contain the smaller secondary channels, are preserved over the course of the LMV (Fig. 8). As no remnants of the main Younger Dryas channel are likely to be preserved (because of reworking), it might be proposed that this main channel had a less braiding planform, in which case the Younger Dryas river within the LMV might be better characterized as a wandering river (i.e., transitional between braided and meandering) rather than a fully braided one. The morphology of the Younger Dryas terrace fragments on the VB, containing multiple parallel low-sinuosity braided channels, show that this was most likely the case in the LMV.
The rapid rise in temperature at the start of the Holocene caused the discharge regime to change (again) from nival to rainstorm-dominated. Furthermore, the local permafrost within the LMV degraded quickly and the depth of seasonal frost decreased, both contributing to an enhanced soilinfiltration capacity. Together with an increasing vegetation cover, this resulted in the concentration of discharge into a decreasing number of channels until only one channel remained. In parallel, renewed trapping of overbank loams occurred and bank stability increased. In other words, very similar "stepwise channel abandonment" occurred in the LMV during the Early Holocene, as reconstructed, over much greater valley width in the Bølling period. The avulsion of the Wurfeld and Maaseik meander belt during the Preboreal indicates that multiple (meandering) channels were active during the Early Holocene in the RVG and, most likely, also during the late glacial (Fig. 13). Multi-channel meandering systems, simultaneously active during the Early Holocene, were also recognized by Erkens et al. (2009) for the Rhine river in the Upper Rhine Graben and for the Meuse in the downstream Rhine-Meuse delta (Cohen et al., 2002Hijma et al., 2009).
After climatic amelioration in the Preboreal, the Holocene climate stabilized and gave way to continued vegetation succession. This led to the development of a mixed deciduous forest by the end of the Boreal and further stabilization of river banks and (floodplain) soils. Consequently, concentration of discharge continued, leading to a single-channel meandering system during the Middle Holocene in the LMV (Fig. 13).
The final phase of terrace formation within the LMV, during the Late Holocene, can be ascribed to human impact. Deforestation and wide-spread agricultural practices in the study area and the hinterland in the last 4000 yr, have led to increased discharges (6.6%), higher sediment load (~300% increase), and increased bank stability by overbank deposits (Dambeck and Thiemeyer, 2002;Bos et al., 2008;Erkens, 2009;Erkens et al., 2009). At some locations, lateral channel activity could increase by, for instance, deforestation of levees and floodplain in combination with an increased discharge (e.g., the Stokkem site on the CB [Paulissen, 1973]).

Spatial variations in river terraces
The fact that most preceding studies on LMV terraces have focused on the record in the northern half of the study area ascribes this region to hold a well-preserved terrace record of distinct fluvial styles. Sufficient lateral movement was possible at the VB for the river to adjust its planform to changing climatic conditions, while, at the same time, the river was able to incise and thereby preserve its former plain (Fig. 8, 9, and 13). Additionally, the relatively low gradient (~0.15 m/km) and subsidence, together with the downstream position, i.e., in the Niers-Rhine confluence as it was functioning in Late Pleniglacial, Bølling-Allerød, and Younger Dryas times (Kasse et al., 2005), provided favorable conditions for terrace formation and preservation.
The PB forms a bottleneck in the course of the LMV where the width of the river plain has reduced significantly over the PB compared to the upstream RVG, especially since the Younger Dryas (Fig. 13). This narrowing can be explained by a combination of relative uplift of the PB with respect to the adjacent VB and RVG, and to restricted lateral mobility of the incised river channel, due to the cohesive subsurface of the PB. The intensified narrowing since the Younger Dryas is probably caused by continued lowering of the thalweg in the resistive subsurface, limiting the available energy for lateral displacement (i.e., vertical incision exceeds lateral movement; Fig. 13).
Lateral movement dominates over vertical incision in the subsiding RVG, where the lithology of the subsurface consists of coarse fluvial deposits, which is represented by a Holocene high-sinuosity meandering river with a relatively wide plain (Fig. 13). The amount of lateral displacement during the Holocene is higher in the RVG compared to the other tectonic blocks and is most likely caused by a combination of a relatively high gradient (~0.5 m/km) and subsidence rate together with a high sediment load in this sub-reach. This relatively high sediment load is probably the result of the narrow river plain and erosional character of the Meuse river on the upstream CB and in the hinterland (Ardennes). This combination of factors makes avulsion of river channels a likely process in this area.
On the southern CB, where the Meuse has to cut through Cretaceous and Paleogene limestone, vertical incision clearly prevailed over lateral movement (Fig. 13). The domination of vertical incision over lateral displacement reduces in the northern CB, where the subsurface consist of similar cohesive sediments as at the PB. Here, several Pleniglacial and late glacial terrace fragments lie within the Holocene river plain, indicating that a multi-channel/avulsive system was present during the Middle Holocene and most likely already during the Younger Dryas or Early Holocene. This resulted in the relatively poor preservation of older terraces due to reworking by the laterally active younger channels. The fact that the planform of the Holocene Meuse on the northern CB is sinuous rather than straight, i.e., like on the PB, is probably caused by the more upstream position of the CB along the course of the LMV. This upstream position resulted in a higher gradient, and therefore higher stream power, and greater sediment load, making the river more likely to move laterally by either meandering or avulsion.
The meandering Younger Dryas channels at the Feldbiss Fault Zone (FFZ; Fig. 12C and D) indicate that the Meuse, at this location, either remained in a meandering mode since the Allerød, or that the river started to meander within the Younger Dryas. Meandering Younger Dryas river planforms were also observed in the Warta and Tisza (Vandenberghe et al., 1994;Kasse et al., 2010). In case of the Meuse river, possible causes which may be responsible for this anomalous meandering planform are a (sudden) change in gradient as the river crosses the active FFZ into the RVG or a change in lithology or grain size of the subsurface over the fault zone (passive tectonic control; Fig. 3). The same holds true for the anomalous planform of the Allerød meander at Roermond, where vertical displacement along the active PBFZ (~66 mm/ ka; Michon and Van Balen, 2005) hampered downstream migration of the meander, forcing it to move laterally and in upstream direction. The different morphological development around the FFZ during the Younger Dryas and PBFZ during the Allerød suggests that the normal faults of the FFZ and PBFZ influenced the river dynamics of the Meuse. It strongly suggests that fault zones within the RVRS executed direct controls on local fluvial morphology in the Weichselian late glacial and Holocene, which will be the subject of a separate paper.

CONCLUSIONS
The palaeogeographic evolution and associated changes in fluvial style in the Lower Meuse Valley at the transition of the Weichselian to the Holocene are controlled by changes in climate and associated hydrological, vegetational, and sedimentary responses. The Lower Meuse Valley was braided during the last glacial maximum and had a multi-channel low-sinuosity planform during the Bølling that evolved into a meandering planform during the Allerød interstade. Subsequently, it was braided again during the Younger Dryas cold stage. Finally, it became meandering again at the onset of the 565 Spatial and temporal variations in the Lower Meuse Valley river terraces, the Netherlands