Late Holocene flood magnitudes in the Lower Rhine river valley and upper delta resolved by a two‐dimensional hydraulic modelling approach

Palaeoflood hydraulic modelling is essential for quantifying ‘millennial flood’ events not covered in the instrumental record. Palaeoflood modelling research has largely focused on one‐dimensional analysis for geomorphologically stable fluvial settings because two‐dimensional analysis for dynamic alluvial settings is time consuming and requires a detailed representation of the past landscape. In this study, we make the step to spatially continuous palaeoflood modelling for a large and dynamic lowland area. We applied advanced hydraulic model simulations (1D–2D coupled set‐up in HEC‐RAS with 950 channel sections and 108 × 103 floodplain grid cells) to quantify the extent and magnitude of past floods in the Lower Rhine river valley and upper delta. As input, we used a high‐resolution terrain reconstruction (palaeo‐DEM) of the area in early mediaeval times, complemented with hydraulic roughness values. After conducting a series of model runs with increasing discharge magnitudes at the upstream boundary, we compared the simulated flood water levels with an inventory of exceeded and non‐exceeded elevations extracted from various geological, archaeological and historical sources. This comparison demonstrated a Lower Rhine millennial flood magnitude of approximately 14,000 m3/s for the Late Holocene period before late mediaeval times. This value exceeds the largest measured discharges in the instrumental record, but not the design discharges currently accounted for in flood risk management.

Many studies on past floods aim to quantify the peak discharge reached during specific events (Herget & Meurs, 2010;Hu et al., 2016;Toonen et al., 2013) or a non-exceeded discharge over a prolonged time period (Enzel et al., 1993;England et al., 2010).
Traditionally, palaeoflood modelling studies have been conducted in geomorphologically stable fluvial settings, such as bedrock canyons, where calculating a representative discharge is relatively straightforward (Kochel & Baker, 1982;overviews in Baker, 2008 andDíez-Herrero, 2015). In settings where the fluvial landscape has changed due to natural fluvial processes and anthropogenic activities, however, the floodplain topography and channel bathymetry must be reconstructed prior to calculations of palaeoflood discharges (Herget & Meurs, 2010;Machado et al., 2017;Toonen et al., 2013).
Although major advantages of deploying 1D hydraulic models are computational speed and little required input, these models cannot resolve spatial aspects of flooding such as floodwater dispersal patterns over the floodplain. Therefore, 1D modelling produces firstorder calculations of palaeoflood peak magnitudes rather than actual palaeoflood simulations. Introducing a two-dimensional (2D) model component lessens the assumptions in the modelling step of palaeoflood research and yields more realistic output than a 1D model (see discussion in Herget et al., 2014), providing stages and discharges at different points in time over the entire model area. Many studies discuss the performance of 2D hydraulic models for present landscape situations (Bomers et al., 2019b;Czuba et al., 2019;Dottori et al., 2013), and the importance of applying these in palaeoflood research is often emphasised, but the actual use is usually outside time and budget constraints (Benito & Díez-Herrero, 2015;England et al., 2010).
Besides computing power limits, the main problem associated with 2D modelling of past floods is a requirement for accurate highresolution input data, most notably a reconstruction of the terrain.
Generating this input involves the use of geological and geomorphological field data, and experimental GIS techniques (van der Meulen et al., 2020). A few studies have applied 2D models to simulate relatively recent floods, occurring in the nineteenth century (Bomers et al., 2019;Calenda et al., 2003;Hesselink et al., 2003;Skublics et al., 2016) and the twentieth century (Ballesteros et al., 2011;Bomers et al., 2019a;Denlinger et al., 2002;England et al., 2007;Masoero et al., 2013;Stamataki & Kjeldsen, 2020). For such recent events, abundant data are available for model schematisation and validation. Changes in morphology can be regarded as either insignificant or well documented and readily incorporated. Simulations of older floods in larger alluvial reaches, such as the Lower Rhine, are lacking due to terrain reconstruction demands. However, these reaches are usually densely populated, warranting detailed assessments of river flood risk. For such regions, it would be of interest to simulate extreme events in earlier historic or pre-historic times, thereby constraining peak discharges over longer time periods that are usually the focus of palaeoflood research.
The Lower Rhine is a large and economically important meandering river which transitions into the Rhine delta near the German-Dutch border (Figure 1). Here, it divides into three main distributaries: the Waal, the Nederrijn and the IJssel. These river branches have been embanked since late mediaeval times (Hesselink, 2002). Across the valley and delta, the river and floodplain have been heavily modified by anthropogenic activities (Hudson et al., 2008;Kalweit et al., 1993;Lammersen et al., 2002;van der Meulen et al., 2020).
Throughout the Holocene, Lower Rhine floods have left scattered geological and geomorphological traces in the landscape (Gouw & Erkens, 2007;Minderhoud et al., 2016;Pierik et al., 2017;Toonen et al., 2013;Willemse, 2019). In historic times, large floods were a recurring threat for the inhabitants of the area (Tol & Langen, 2000), and at present, considerable resources are dedicated to developing and maintaining flood protection structures. A major challenge lies in determining appropriate safety standards for river management works and underlying magnitude-frequency relationships of extreme events (Hegnauer et al., 2014;Hegnauer et al., 2015;van Alphen, 2016).
Previous palaeohydrological research in the delta apex region has provided estimates of Lower Rhine flood frequencies and relative magnitudes over the past centuries to millennia based on sedimentary archives obtained from oxbow-lake fills Toonen, 2013;Toonen et al., 2015;Toonen et al., 2017). These studies have shown that between 4.2 and 0.8 ka (kiloannum = thousands of years ago), that is, in the Late Holocene up to embankment along the downstream reaches, five flood events with recurrence times close to 1,000 years occurred. The largest of these was dated to 785 CE based on a tentative correlation between radiocarbon-dated flood deposits and historical sources Toonen, 2013). Recently, a palaeo-DEM has been constructed for the Lower Rhine river valley and upper delta in the first millennium CE (van der Meulen et al., 2020). The present research builds upon these studies by applying a hydraulic model to the palaeo-DEM, resolving absolute magnitudes of extreme floods in the Late Holocene period preceding the onset of river embankment in late mediaeval times.
The purpose of this study is to constrain past flooding patterns and magnitudes in the Lower Rhine river valley and upper delta ( Figure 1). Specifically, we aim to quantify the peak discharge of the most extreme ('millennial') floods that occurred. Likely, this value is in between the largest measured flood discharge (12,000 m 3 /s) and the largest discharges generated by stochastic weather simulations coupled to hydrological models (GRADE project: 24,000 m 3 /s; Hegnauer et al., 2014;Hegnauer et al., 2015). We set up a 1D-2D coupled hydraulic model with reconstructed terrain and roughness as input to simulate the propagation of extreme discharge waves in the past landscape setting. Subsequently, we compare the output water levels with geological, geomorphological and archaeological observations informing on minimum and maximum palaeoflood levels in the study area. In addition to calculating a Late Holocene millennial flood magnitude for the Lower Rhine river, this study demonstrates for the first time the potential and the difficulties of spatially continuous palaeoflood modelling in large lowland areas. At the least, such research requires (1) an accurate reconstruction of the river and floodplain, (2) a numerical model with a two-dimensional component and (3) a collection of spatially distributed palaeoflood levels.

| MATERIALS AND METHODS
The overall 'inverse modelling' approach applied in this study comprises three major steps ( Figure 2). The first step includes the reconstruction of the landscape, which encompasses floodplain topography, river position and bathymetry, and roughness values matching past land cover. Further, an inventory of palaeoflood levels was generated based on earlier geological and archaeological investigations. The second step consists of hydraulic modelling, using the past landscape as input. Next to terrain and roughness, the model set-up requires upstream and downstream boundary conditions, and a rasterization of the model area consisting of a 2D grid on the floodplain and 1D profiles in the channels. Key to the 'inverse modelling' approach is that multiple model runs are conducted with discharge waves of increasing magnitude as upstream boundary condition. In the third step, the output of the successive model runs is compared with the inventory of palaeoflood levels. The match between observed and simulated flood levels (m) is used to evaluate the model output resulting from F I G U R E 2 Overview of the 'inverse modelling' approach applied in this study [Colour figure can be viewed at wileyonlinelibrary.com] F I G U R E 1 Extent of model area and numbered profiles on LiDAR DEM of North Rhine-Westphalia and the Netherlands (modified from van der Meulen et al., 2020). The red crosses indicate the locations of observational data on palaeoflood levels ( Supplementary Table). The grey boxes centred on Cologne and Duisburg-Xanten indicate the extents of Figures 3 and 4 [Colour figure can be viewed at wileyonlinelibrary.com] different peak flows (m 3 /s), thus constraining the maximum discharge that occurred in the studied period.
The upstream boundary of the model area is at river km 638, between Andernach and Bonn, Germany, where the Rhine changes from a bedrock-confined river to an alluvial meandering river in a terraced valley floor ( Figure 1). From here, the model area covers the entire floodplain of the Lower Rhine valley and the delta apex region. The downstream boundaries are placed at river km 908 (Waal), km 912 (Nederrijn) and km 950 (IJssel), where the delta plain widens substantially and water level slopes of the river branches are influenced by tides (Kleinhans et al., 2011;Stouthamer et al., 2011).  (Harnischmacher & Zepp, 2014). In addition, the river position and channel bathymetry were reconstructed by combining available historical, archaeological and geological data, and the natural floodplain elevation along the river was restored using geomorphological interpolations (van der Meulen et al., 2020).
The early mediaeval target age of the palaeo-DEM pre-dates any major direct modifications to the terrain by humans. Therefore, the distal floodplain topography is representative for the entire Late Holocene period up to late mediaeval times (van der Meulen et al., 2020).
The proximal floodplain topography and river position are considerably more variable through time, but within the Late Holocene the fluvial style and channel dimensions of the Rhine have been fairly stable (Klosterman, 1992;Erkens et al., 2011). Furthermore, local morphological change due to fluvial erosion and deposition likely did not alter the overall conveyance capacity for peak discharges at the scale of the entire valley. Similarly, indirect human impacts on the topography, most notably changes in sediment budget (Hoffmann et al., 2007(Hoffmann et al., , 2009Middelkoop et al., 2010;Notebaert & Verstraeten, 2010), affected the area during the Late Holocene, but probably hardly influenced the propagation of extreme discharge waves over the full model area. As precisely these extremes are the focus of our simulations and analyses, we regard the simulation results produced using the early mediaeval palaeo-DEM representative for older Late Holocene times.

| Hydraulic roughness
The next component of the model schematisation is hydraulic roughness, expressed as Manning's n values. Vegetation distribution and land use are broadly known for early mediaeval times in the study F I G U R E 3 LiDAR ground level DEM and palaeo-DEM of model area around Cologne, Germany (modified from van der Meulen et al., 2020). All anthropogenic modifications, such as embankments for roads and railroads, dump sites and quarry pits, have been removed, and the 'original' topography has been restored by various interpolation techniques. The river position has been reconstructed by combining historical and archaeological sources with geologicalgeomorphological insights. The extent of the figure is indicated in Figure 1 [Colour figure can be viewed at wileyonlinelibrary.com] area, but detailed constraints are scarce (Gouw-Bouman, in prep).
Therefore, we developed a simplified land cover reconstruction. We distinguished five landscape classes based on distance to the river and relative floodplain elevation (with respect to a second-order polynomial trend surface), each associated with different natural vegetation and land use characteristics, and hence different average hydraulic roughness values ( Figure 4; Table 1). We used standard values for different vegetation types (Chow, 1959) as starting point to assign n values to landscape classes based on their estimated overall land cover composition.
The areas located more than 5 m above the trend surface (high grounds, class H) were mostly covered by forest, resulting in high n values. Note that class H covers only a small area (Table 1), because we cut off the model area where elevation significantly rises. The river itself (polygon obtained from palaeo-DEM, class R) consisted mostly of bare gravel and sand, resulting in low n values. The proximal floodplain (class P), which we defined as a 1-km-wide zone next to the river (Figure 4), had dispersed riparian vegetation of the type modern floodplain-ecologists aim to restore (van Oorschot, 2017) and was assigned relatively high n values. Realistically, the width of this zone greatly varied along the river, and thus our definition is a major simplification. However, given the large study area, a detailed reconstruction was not feasible. The higher parts of the distal floodplain (class DH) were used most intensively for settlement and agriculture (Pierik & van Lanen, 2019), and thus had little natural vegetation. Agricultural land in early mediaeval times consisted of smaller fields than today, with abundant hedges and local stands of trees. Therefore, we assigned higher n values to this class than used for agricultural land in models of present situations. The lower parts of distal floodplain (class DL) were, and still are, relatively wet and mostly used as meadows for grazing and hay production. In these areas, open grasslands prevailed, but locally swamp forests (carr) occurred as well (Gouw-Bouman et al., in prep). Grasslands have low n values, but the characteristically dense swamp forests have high n values. However, it is not possible to reconstruct the precise spatial distribution of grasslands and swamp forests within this zone. Overall, this class mainly consisted of grass; therefore, we assigned it low n values throughout the model area.
Besides average roughness values (n best ), we assigned lower (n min ) and upper estimates (n max ) to all classes, except for class H with its marginal occurrence (Table 1). The n min and n max were assigned to account for the uncertainties inherent in assigning roughness values (Bomers et al., 2019a;Chow, 1959;Warmink et al., 2013). All landscape classes, expect class R, consist of a mixture of relatively open and forested parts. Because of the cooccurrence of different vegetation types, Manning's n values beyond the lower and upper estimates given in Table 1 are unlikely to be correct for any Late Holocene situation. In the Rhine area, significant catchment-scale deforestation occurred since the Bronze Age circa 3.5 ka. There was a slight increase in natural vegetation cover in the Dark Ages around 1.5 ka, due to population decline and land abandonment following Roman times (Teunissen, 1990;Kaplan et al., 2009;Gouw-Bouman, in prep). For these reasons, lower n values (between n min and n best ) may be appropriate for Middle Holocene to early Late Holocene times, whereas higher values (between n best and n max ) may be typical for Roman times and for late mediaeval times. We used the n best values in the model runs, unless otherwise stated. We used the full momentum equations to solve the system, because these resulted in more accurate model results compared with the simplified diffusive wave equations.

| Model set-up
The upstream boundary condition consists of a discharge time series ( Figure 6). An initial discharge of 1,000 m 3 /s was used in all runs to avoid a dry channel at the start of the simulations. We selected the two largest measured discharge waves as input. These are the floods of December 1925/January 1926 (Bomers et al., 2019a) and January 1995. For the palaeoflood simulations, we used the shape of the 1926 discharge wave as a representative standard hydrograph because it closely resembles the average of modelled extreme flood waves of the Rhine near Andernach (Hegnauer et al., 2014;Hegnauer et al., 2015). By scaling this hydrograph, we obtained a series of input waves with different peak discharges, from 10,000 to 30,000 m 3 /s with intervals of 2,000 m 3 /s ( Figure 6). We simplified the discharge were used as downstream boundary conditions at the locations where water can potentially leave the model domain ( Figure 5). These normal depths were calculated using Manning's equation (Brunner, 2016).
In total, we conducted 15 simulations using the palaeoflood model set-up, one for each different peak discharge and two extra runs for both the 10,000 and 24,000 m 3 /s input waves using the minimum and maximum Manning's n values for all classes (Table 1)

| Evaluation of model output and palaeoflood discharges
As a first step in analysing the model output, we plotted maps showing maximum water depths for each simulation. We used these to extract the maximum inundation extent for each different input discharge. Next, we plotted the simulated water levels and discharges at seven approximately equally spaced profiles across the model area ( Figure 1). From the amounts of water passing these locations (discharge per time step), we determined flood wave propagation. Further, we plotted for all profile locations the relation between input discharge and water level.
To determine the maximum discharge that occurred in the Late Holocene to early mediaeval time period, we compared the simulated flood levels with an inventory of observational data (Supplementary Table) Toonen, 2013). 3 | RESULTS

| Inundation depth and extent
The effect of topography on flooded area is extensive. In the present landscape setting, an extreme flood such as that of January 1995 only inundates the embanked floodplains, amounting to about 450 km 2 , given that no dyke breaches occur ( Figure 7A). In contrast, a similar flood inundates about 2,370 km 2 in the pre-embanked landscape setting ( Figure 7C). Conversely, inundation depths are much higher in the flooded area between the embankments in the present situation ( Figure 7B) than in the past situation when overbank flow was unconfined ( Figure 7D).
The water depths in the pre-embanked landscape are generally larger in the upstream part of the model area than in the downstream part ( Figures 7C and 8). This is because elevation differences decrease in downstream direction, resulting in wider and relatively lower floodplains (Figure 1). This is especially clear in the delta apex where the floodplain divides, circa 50 km upstream of the channel bifurcations ( Figure 8).
Both water depth and inundation extent increase with discharge ( Figure 8; Figure 9). Interestingly, the increase in inundation extent is nearly linear for extreme discharges (Figure 9). Most of the model area consists of a terraced valley, where additional surfaces flood with each step in simulated discharge, although water depths at the newly flooded locations are small. This is different in the present embanked situation, where the inundation extent ceases to increase when all areas between embankments are flooded ( Figure 7A). Despite this terrace effect, the graph does not appear step wise (Figure 9). The expected incremental trend is obscured due to the large size of the area (>4,500 km 2 ) with many small elevation differences and local terrace relief (Figure 3).
When hydraulic roughness values are raised, the inundation extent increases (Figure 9). The total spread in extent resulting from varying the roughness in all landscape classes between n min and n max is approximately 15%. This value is rather stable for different discharges, amounting to 15.5% for a peak discharge of 10,000 m 3 /s, and 14.5% for a peak discharge of 24,000 m 3 /s ( Figure 9).

| Flood wave propagation
Simulated peak discharges in the pre-embanked landscape remain remarkably constant along the Lower Rhine ( Figure 10). Retention in upstream parts of the area is apparently small and balanced by the discharge contribution of tributaries, indicating that nearly the complete floodplain conveys water downstream. The rise in water levels with discharge is largest in the upstream part of the study area (Figure 11), which is related to the floodplain widening in downstream direction.
An increase in peak discharge from 10,000 to 30,000 m 3 /s leads to a 1.5 to 3 m rise in maximum water levels in the river valley (profiles 1 to 3 in Figure 11), but less than 1 m in the delta (profile 6 in Figure 11). Both extent and depth of inundation exhibit slightly declining (concave down) increasing trends with discharge ( Figures 9 and   11), because flow velocities increase as water depths across the floodplain rise.
In all model runs, the westward route in the delta (Waal and Nederrijn Rivers and associated floodplain) carries more water than the northward route (IJssel valley). The division ratio is approximately 2:1 for a 10,000 m 3 /s input discharge wave but reduces to approximately 3:2 for 18,000 m 3 /s ( Figure 10). For discharges lower than 10,000 m 3 /s (not the focus of this study), an increasingly larger share of the flood water takes the westward route. After approximately 6,000 m 3 /s enters the model area at the start of the modelled time period in the 18,000 m 3 /s run, only one-fifth of total discharges flows northward ( Figure 10). The Q-h curve in the upstream part of the westward route (profile 4 in Figure 11; 'Gelderse Poort') is similar to those in the upstream (profile 5 in Figure 11; 'Oude IJssel') and downstream (profile 7 in Figure 11; 'IJssel') parts of the northward route.
However, in the downstream part of the westward route (profile 6 in Figure 11; 'Betuwe'), the Q-h curve deviates and the rise in water level with discharge is significantly lower, owing to widening of the delta plain.  Table). Somewhat in line with the observational data, the simulated water levels for large (10,000 m 3 /s) and extreme (30,000 m 3 /s) discharges differ in general by about 2 to 3 m in the upstream part of the study area, and by about 1 to 2 m in the downstream part ( Figure 12). This indicates that water levels are less sensitive to discharge in the downstream and delta parts of the study area, which is attributable to the wider floodplain in these areas.
The large spread in the observational data implies that, for a simu-  In the Lower Rhine delta, the floodplain splits into two distributive systems (Figure 1), resulting in a larger area for floodwater dispersal than in the valley. In the present situation ( Figure 7A), the water diverges where the river bifurcates, close to Lobith. However, in the past situation ( Figure 7C), floodwater diversion occurred considerably further upstream, close to Wesel. Current river management in the Netherlands mainly accounts for water entering the country near Lobith, but potential dyke breaches in the German-Dutch border region and consequent floodwater diversion upstream of the river bifurcation can greatly affect the discharge distribution of the Rhine river branches (Bomers et al., 2019b;Parmet et al., 2001;te Linde et al., 2011). Our palaeoflood simulations underline this often overlooked importance of the delta apex area in flood risk analyses.

| Palaeoflood magnitudes in the Lower Rhine valley and upper delta
Our 'inverse modelling' approach has enabled us to translate geological, archaeological and historical data on past flood levels into discharge values (Figures 2 and 12-14). Comparison of the palaeoflood simulations with observational data suggests a Late Holocene to early mediaeval millennial flood magnitude lower than 18,000 m 3 /s, most probably around 14,000 m 3 /s. This value is larger than any measured flood in the instrumental record (approximately 12,000 m 3 /s in 1925-1926. Accordingly, a peak discharge of 14,000 m 3 /s is our best estimate for the geologically recognised event in 785 CE Toonen, 2013) as well as for the flood cluster episode (cf. Toonen et al., 2017) that is dendrologically identified (Jansma, 2020;Sass-Klaassen & Hanraets, 2006) and that is thought to have initiated the IJssel northward deltaic branch in the sixth to seventh century CE (Cohen et al., 2009;Groothedde, 2010;Makaske et al., 2008). Note that the IJssel river channel is absent in the palaeo-DEM and model schematisation, as it had not matured before late mediaeval times (circa 1,100 CE).
In late mediaeval times, the sedimentary record reveals one flood event of greater apparent magnitude than any in the millennia before Toonen, 2013). This event is linked to the year 1374 CE, which is recognised in historical sources as a year of extreme water levels and major flooding damage along the Lower Rhine (Buisman, 1996;Gottschalk, 1971;Weikinn, 1958). A crosssectional calculation of the 1374 CE event near Cologne resulted in an estimated peak discharge of almost 24,000 m 3 /s (18,800 m 3 / s < Q < 29,000 m 3 /s; Herget & Meurs, 2010). Our findings suggest that only the very lower end of that estimate may be realistic, even though the event is outside the period covered by our simulations as it post-dates embankment in downstream parts of the study area.
Future reconstruction and modelling efforts are underway to specifically address the magnitude of the 1374 CE flood.
In the study area, embankments and other river management works are currently designed to accommodate peak discharges with recurrence times varying from hundreds to thousands of years, with differences between German and Dutch safety standards (Hegnauer et al., 2015;te Linde et al., 2011). The associated peak discharge has been set to 16,000 m 3 /s at Lobith after the Lower Rhine flood of 1995 (Chbab, 1996). Although this value is a subject of discussion and the newly adopted risk approach no longer identifies a single peak flow as standard (Deltacommissie, 2008;van Alphen, 2016), the design discharges in the German-Dutch border region exceed our F I G U R E 1 4 Number of exceeded minimum and non-exceeded maximum elevation markers for each input discharge.

| Additional applications of model output
The model output supports our understanding of geological and geomorphological characteristics of the study area. For example, the coarsening of the river sediment since the onset of embankment (Frings et al., 2009)  Future modelling efforts may incorporate sediment transport and morphological processes to explain overbank deposition patterns in the natural landscape setting, such as oxbow-lake infilling (Ishii & Hori, 2016;Toonen et al., 2012) and natural levee formation (Johnston et al., 2019;Pierik et al., 2017). However, modelling floodplain sedimentation in detail is largely limited to local studies for present landscape situations (Middelkoop & Van Der Perk, 1998;Nicholas et al., 2006;Thonon et al., 2007).  Table), which is similar in other parts of the Roman Empire (Obrocki et al., 2020;O'Shea & Lewin, 2020 Flooding may have been both a taphonomic factor affecting preservation and an archaeological factor affecting post-Roman land use and settlement continuity. The upper delta and especially the Betuwe area (between the Waal and Nederrijn rivers) flooded rapidly in all palaeoflood simulations before arrival of the discharge peak. This explains the habitation patterns that occurred in the first millennium CE, which were mostly limited to former natural levees (Pierik & van Lanen, 2019). Further, it illustrates the importance of 'zijdewendes' (local dykes perpendicular to the river), which were constructed already prior to dykes along the river and have continued to play a role in flood mitigation (Alkema & Middelkoop, 2005;van de Ven, 1993).

| Modelling advancements in palaeoflood hydrology
Our study shows that upscaling palaeoflood reconstructions to multiple dimensions provides important information that cross-sectional or longitudinal approaches cannot resolve (Baker, 2008;Benito & Díez-Herrero, 2015;Webb & Jarrett, 2002). Models with a 2D component can account for spatial flow patterns such as velocity and direction (Horrit & Bates, 2002;Liu et al., 2015;Tayefi et al., 2007). Especially in a delta setting, the incorporation of multiple dimensions is crucial, as discharge diversions can otherwise not be properly resolved.
Another aspect not incorporated by 1D models is that the maximum water level is rarely horizontal in a direction perpendicular to the main flow. The two-dimensional nature of flood flows may cause some deposits to be placed above the cross-sectional average water surface, resulting in an often-ignored inaccuracy factor in palaeoflood studies (Benito & O'Connor, 2013). This offset is significant (up to 1 m; Figure 7D) and thus relevant when comparing observational levels with simulated output (Figure 12).
A major benefit of 1D-2D modelling compared with crosssectional 1D analyses (Herget et al., 2014;Herget & Meurs, 2010;Toonen et al., 2013) is that not only the peak discharge but also the shape of the discharge wave may be varied ( Figure 6). This is important in palaeoflood research since hydrograph characteristics such as volume affect the flooding patterns and water levels (Bomers et al., 2019;Vorogushyn et al., 2010). Choosing a discharge wave requires identifying a shape that is representative for an extreme flood in the river system under study, as we did in section 2.3. Due to its large size, our study area covers different geomorphological settings, which respond differently to extreme floods. An increase in discharge results in a larger water level rise in the valley than in the delta (Figure 11). This is related both to the lowering of the gradient and to differences in cross-sectional morphology (visualised in Figure 12). Because of the low sensitivity of water level to discharge, determining palaeoflood magnitudes from geological or his- mediaeval times (before the first embankments) had a peak discharge lower than 18,000 m 3 /s, with a best estimate of 14,000 m 3 /s. This millennial flood magnitude is larger than the most extreme measured discharge in the instrumental record, but lower than the magnitudes accounted for in flood risk management of the modern, engineered Rhine.

ACKNOWLEDGEMENTS
This work is part of the research programme 'Floods of the Past,