Extreme Floods in Small Mediterranean Catchments: Long-Term Response to Climate Variability and Change

: Climate change implies changes in the frequency and magnitude of ﬂood events. The inﬂuence of climate variability on ﬂooding was evaluated by an analysis of sedimentary (palaeoﬂoods) and documentary archives. A 500-year palaeoﬂood record at Montlle ó River (657 km 2 in catchment area), eastern Spain, revealed up to 31 palaeoﬂoods with a range of discharges of 20–950 m 3 s − 1 , and with at least ﬁve ﬂoods exceeding 740–950 m 3 s − 1 . This information contrasts with the available gauged ﬂood registers (since year 1971) with an annual maximum daily discharge of 129 m 3 s − 1 . Our palaeoﬂood dataset indicates ﬂood cluster episodes at (1) 1570–1620, (2) 1775–1795, (3) 1850–1890, and (4) 1920–1969. Flood rich periods 1 and 3 corresponded to cooler than usual (about 0.3 ◦ C and 0.2 ◦ C) climate oscillations, whereas 2 and 4 were characterised by higher inter-annual climatic variability (ﬂoods and droughts). This high inter-annual rainfall variability increased over the last 150 years, leading to a reduction of annual maximum ﬂow. Flood quantiles ( > 50 years) calculated from palaeoﬂood + gauged data showed 30%–40% higher peak discharges than those using only instrumental records, whereas when increasing the catchment area (1500 km 2 ) the discharge estimation variance decreased to ~15%. The results reﬂect the higher sensitivity of small catchments to changes on ﬂood magnitude and frequency due to climate variability whereas a larger catchment bu ﬀ ers the response due to the limited extent of convective storms. Our ﬁndings show that extended ﬂood records provide robust knowledge about hazardous ﬂooding that can assist in the prioritization of low-regret actions for ﬂood-risk adaptation to climate change. 2 ), including the Montlle ó River basin, shows a di ﬀ erence of 15–25%. These results show the higher di ﬀ erence of ﬂood discharge for a given quantile in smaller catchments, which is likely due to the limited extent of mesoscale convective cells producing large ﬂoods. These ﬁndings have important implications for ﬂood prevention, risk assessment and adaptation in small catchments of the Mediterranean region.


Introduction
Global warming is leading to widespread changes in the hydrological cycle although the chain linking greenhouse gases to outcome impacts is long and complex [1]. A more variable climate and therefore, a more variable water cycle is likely to affect the pattern of rainfall-producing floods [2].
in width) and expands upstream (60-70 m in width) which favours sedimentation during flood events. The analysis of the fluvial morphodynamic changes based on aerial photographs , supported by field check-up, indicates that some channel bed incision occurred mainly since 1976, when in-stream gravel mining started [21,22]. The uncertainty regarding the valley geometry at the time of the deposition of the slackwater flood deposits was considered in the hydraulic modelling as a possible range of channel topography and discharge.

Palaeoflood Record
Quantitative palaeoflood hydrology methods were used to estimate the magnitude and frequency of past flooding on bases of sedimentary evidence [23]. Fine sediments in suspension are accumulated in slackwater depositional environments, such as channel widening, alcoves and caves in bedrock walls, and back-flooded tributary mouths [12,13]. In these sites, sequences of multiple sedimentation layers record multiple flood events, sometimes with intercalated slope deposits, or with contacts marked by disconformities [13]. A detailed stratigraphy with emphasis on contacts between flood-related sediment beds was performed in 8 stratigraphic profiles, along with a The climate is meso-mediterranean at the upper catchment (pine and evergreen oak forests) to thermo-mediterranean (shrublands) in most of the basin, characterised by a mean annual temperature of 9-15 • C, daily winter minima easily below 0 • C and summer maxima above 30 • C. Mean annual precipitation ranges from over 750 mm in Peñagolosa peak to 500 mm in the lower catchment [18]. Rainfall occurs mainly in autumn (September to November), with a second rainfall maximum in late winter to spring (March to May). The most intense rainfall episodes (up to 300 mm in 24 h) occur in autumn, associated to mesoscale convective cells [19], leading to extreme flooding in the region [20].
In the study reach, locally known as La Volta (The Tour), the valley is constrained laterally by bedrock and the channel bed is built on alluvial gravel ( Figure 1). The lower sector is narrow (20-m in width) and expands upstream (60-70 m in width) which favours sedimentation during flood events. The analysis of the fluvial morphodynamic changes based on aerial photographs , supported by field check-up, indicates that some channel bed incision occurred mainly since 1976, when in-stream gravel mining started [21,22]. The uncertainty regarding the valley geometry at the time of the deposition of the slackwater flood deposits was considered in the hydraulic modelling as a possible range of channel topography and discharge.

Palaeoflood Record
Quantitative palaeoflood hydrology methods were used to estimate the magnitude and frequency of past flooding on bases of sedimentary evidence [23]. Fine sediments in suspension are accumulated in slackwater depositional environments, such as channel widening, alcoves and caves in bedrock walls, and back-flooded tributary mouths [12,13]. In these sites, sequences of multiple sedimentation layers record multiple flood events, sometimes with intercalated slope deposits, or with contacts marked by disconformities [13]. A detailed stratigraphy with emphasis on contacts between flood-related sediment beds was performed in 8 stratigraphic profiles, along with a stratigraphic panel of the main outcrop with lateral changes of flood layer contacts and unconformities. In these stratigraphic profiles (B1, B2, etc.; Figure 1) particular attention was given to the identification of individual flood layers (unit number; Unit #), texture, colour, organic remnants for age dating (Table 1) and sedimentary flow structures [24]. The stratigraphy was completed with a trench, dug on the inner alcove infill (close to its end), as well as on remnants of high SWD sediments attached to the alcove walls. At selected sites, sediment peels of the stratigraphic profiles, measuring approximately 80 cm × 50 cm in size, were made in the field [25].

Palaeoflood Dating
The age of the palaeoflood events was determined using radiocarbon and optically stimulated luminescence (OSL) dating of samples collected from individual flood units. Radiocarbon AMS analyses were carried out at three laboratories: (1) Spanish Accelerator Centre in Seville (CNA), (2) Poznan Radiocarbon Laboratory (Poz) and (3) Radiocarbon service of Beta Analytic in Florida (Beta). Radiocarbon dates are quoted in the text as the two-sigma calibrated age range (95%). Radiocarbon ages were calibrated to calendar ages by Oxcal 4.3 software [26] based on IntCal 13 [27] calibration data set (Table 1). Luminescence measurements were carried out in two laboratories: (1) the Nordic Laboratory for Luminescence Dating in Denmark, and (2) the Radioisotopes Unit at the University of Seville.
Luminescence dating determines the time elapsed since quartz and feldspar grains were last exposed to daylight [28,29]. For purposes of dating flood deposits, the general assumption is that the sediment was last exposed to daylight (bleached) during transport prior deposition. The OSL dating was completed for a total of 8 sand samples (Table 2) collected in the field using PVC cylinders to avoid exposure of the sediment to white light. In these sand samples, quartz grains of size 180-250 µm were extracted under light-controlled conditions for luminescence measurements. This grain size has been selected as it has shown to be better bleached than smaller grain sizes [30]. In addition, when measured in small aliquots (<30 grains per aliquot), it provides enough resolution on dose to identify incomplete bleaching (insufficient exposure to daylight to reset the luminescence signal during transport) of the sediment [29,30]. A preheat temperature of 200 • C for 10 s and a cut-heat of 180 • C at a heating rate of 5 • C/s was used for all measurements. A dose recovery experiment has shown that this temperature combination does not derive in thermal transfer. In total, 40-60 aliquots of 2 mm (~20 grains each aliquot) have been measured for each sample to generate the corresponding dose distributions. Individual equivalent doses (De) were determined by interpolating the natural luminescence signal into the dose response curved defined using the single aliquot regenerative dose (SAR) protocol [31]. The measurements included a step with infrared stimulation to detect the presence of feldspar [32], which could compromise the reliability of the detected signal assumed to generate from quartz. Burial De for each sample was estimated using the Internal-External Consistency Criteria (IEU) approach [33,34] proposed to identify the population of doses most likely to correspond to well bleached grains. The parameters needed for the IEU approach are based on an ideal dose distribution only affected by intrinsic factors, therefore excluding the effect of incomplete bleaching. Such "ideal" dose distribution was generated by bleaching a subsample of sample MLL-17 and irradiating it with a known dose. The derived dose distribution was characterized by an over-dispersion of 15%.
The annual dose rates are based on the radionuclide activity concentration measured on~200 g of bulk material using high resolution gamma spectrometry. Contribution of cosmic radiation has been calculated based on the burial depth [35]. Total dose rates ( Table 2) were calculated using DRAC v1.2 [36] taking into account attenuation due to moisture and grain size. An uncertainty of 2% has been added to the assigned water contents to account for variability during the burial period.

Documentary Flood Record
Documentary flood data were compiled to complement the palaeoflood data. Archival flood evidence provides direct (e.g., flood date, duration, stage or spatial extent) or indirect (e.g., post-flood damage repair) information of individual events, allowing flood chronologies and socio-economic impacts to be analyzed [37,38]. Compilation of historical floods that occurred in the studied region was obtained mainly from scientific and technical reports, local history books and non-systematic compilations by historians Balbás Cruz [39], Fogues [40], Fontana [41], and from historical flood compilations by Beltrán Manrique [42], Sánchez-Adell et al. [43], and Camarasa and Segura [19]. Flood intensity was classified according to its basic hydrological behavior and the impacts on five categories following Barriendos et al. [37] classification: (1) Ordinary flood (non-overbank flood + disturbance); (2) Ordinary/extraordinary flood (Non-overbank flood + disturbance + damage); (3) Extraordinary flood (overbank flood + disturbance); (4) Extraordinary/catastrophic flood (Overbank flood + disturbance + damage), and (5) Catastrophic flood (Overbank flood + damage + destruction). A recent update of Mediterranean flood cases from documentary archives including our study region was published by Barriendos et al. [44].

Hydraulic Modelling
Palaeoflood discharge estimates follow from the assumption that the elevation of paleostage evidence (sediment layer) relates closely to the maximum stage attained by an identified flood [23]. Therefore, flood stages associated to the elevation of each flood unit in the stratigraphic profile were converted into discharge values. This conversion is an inverse problem, where the discharge is obtained by trial and error using hydraulic modelling by relating the surveyed elevation of the flood deposits to a simulated water surface profile [45].
The discharge estimation associated with the slackwater flood deposits was calculated using the US Army Corps of Engineers River Analysis System computer program-HEC-RAS [46]. The computation procedure is based on the solution of the one-dimensional energy equation, derived from the Bernoulli equation, for steady gradually varied flow. Therefore, the palaeoflood discharge is based on the calculation of a step-backwater profile producing the best correlation with the geological indicator of flow stage. A field survey was carried out in 2009 using a kinematic differential GPS Trimble 4700 with a precision of ±1 cm in horizontal and ±2 cm in vertical on real-time survey performance.
A total of 8 cross-sections were surveyed along a 1 km reach ( Figure 1). Subcritical flow conditions were assumed along the reach, with critical flow selected as the boundary condition at the most downstream cross-section. Assigned Manning's n values were 0.03 for the valley floor and 0.04 for the valley margins. These n values were chosen according to our previous works in Rambla de la Viuda that used the Limerinos equation [47] for D 50 and D 84 grain sizes [20,48]. The hydraulic model was calibrated using field survey of high water marks of the December 2007 flood (20 m 3 s −1 ). A sensitivity test performed on the model shows that for a 25% variation in roughness values, an error of 10-15% was introduced into the discharge results. Uncertainties in Manning's n values and energy loss coefficients had much less impact than uncertainties in cross-sectional data on discharge values estimated from hydraulic modelling [49].
Along the study reach, the potential channel incision, due to recent gravel mining in the channel, was estimated at one-meter average, based on field survey of gravel bars and channel remnants identified in the 1956 aerial photographs. As this channel bed degradation is due to gravel mining that took place from the 1980s, the discharge estimates of earlier floods are likely to overestimate the real peak flow in ca. 10%. However, as discharge estimates based on paleostage indicators represent a minimum discharge and assuming that slackwater deposition occurs within a water depth of 1-1.5 m, the estimated discharge may be close to real peak values. Individual flood layers overtopping the flood bench commonly increase in elevation into the valley side. Two discharge values were assigned to each flood set derived from modelled water surface profiles matching (1) the base of the alcove; and (2) the highest end-point of the flood sediments.

Frequency Analysis
Flood record stationarity for censored samples (systematic and palaeoflood) was confirmed using Lang's test [50]. This test assumes that stationary flood series can be described by a homogenous or Water 2020, 12, 1008 7 of 23 stationary Poisson process. The 95% tolerance interval of the accumulative number of floods above a threshold, or censored, level, is computed. Stationary flood series are those remaining within the 95% tolerance interval [51]. A flood frequency analysis was performed with the software PeakFQ [52] where systematic and non-systematic data are fit to a log-Pearson Type III distribution. The program applies a generalised method-of-moments estimator denoted the Expected Moments Algorithm (EMA [53]) and a generalised version of the Grubbs-Beck test for low outliers [54]. Visual matching to the calculated palaeoflood and systematic peaks plotting positions and the statistical parameters were used to test for goodness of fit.

Uncertainties Affecting Flood Dating
Chronological control of the timing of the main stratigraphic packages was acheived by combined radiocarbon and OSL ages, and supported by documentary flood records. A most likely age was assigned to sedimentary units based on the combined intersection of luminescence and radiocarbon dating age range (e.g., highest probability of calibrated age ranges; Table 1) within a correct stratigraphic order. In some cases, it proved to be difficult to assign an unambiguous calendar age to individual flood events from 14 C and OSL dating. Post-Sixteenth century radiocarbon dating is affected by human influence on radiocarbon production (industrial and nuclear effects) which can make it difficult, in some cases, to assign a contemporary calendar age [55]. Regarding OSL dating, the dose distributions are characterized by large scatter with over-dispersion values >80%, much higher than the 15% over-dispersion observed on the artificially bleached and irradiated subsample of MLL-17. Several authors applied descriptive and robust statistics to increase precision of the true burial dose of a population with incomplete bleaching grains [29,33,56,57]. Here, we applied the IEU minimum dose approach [33,34] which has successfully estimated ages for sequences of flood deposits in this region [29]. IEU identifies samples B1#3, B1 #4 and MLL #2 to be the ones with large (40-50%) population of well-bleached grains. In the remaining samples, only an average of~15% of the De were identified to belong to well-bleached grains indicating a deficient exposition to light caused by reasons such as a short transport, high turbidity water or if the flood occurred at night. Estimating the equivalent doses from the population most likely to be well-bleached resulted in luminescence ages which are stratigraphically consistent (except sample MLL#17), although comparison with the radiocarbon ages still suggest a minor bias toward an overestimation of the burial dose.

Historical and Systematic Flood Records
The compiled documentary flood register from local historical chronicles (Castelló region) report a total of 400 large flood records since CE 1378. In the case of the Rambla de la Viuda and Montlleó rivers, the level of flood perception by its inhabitants was relatively low, mainly due to the location of major villages far from the river channel and flood areas. In addition, the river channel is entrenched in old alluvium, and the overflowing into urban settlements only occurs at the most downstream reaches, near the coast. Historical flood registers in this area report mainly about the damages caused on linear infrastructures crossing the channel such as irrigation canals, roads and bridges. According to the flood magnitude ranking (1 to 5) [37], most documented floods were classified within the extraordinary (level 3) to catastrophic (level 5) flood categories.
The analysis of regional documentary data (i.e., Millars River, Rambla de la Viuda, and other rivers at Castelló region) indicates an alternation of flood-rich and flood-poor periods over the last 500 years. The highest flood frequency occurred between 1580-1620, 1770-1810, 1880-1925, and 1945-1969. In contrast, periods between 1620-1770 and 1810-1880 were characterised by a high irregularity with few reported events or even absence of floods. In the Rambla de la Viuda river, seven large flood events were reported with category 5 over the period 1580-1900, namely in 1580, 1597, 1617, 1776, 1783, Water 2020, 12, 1008 8 of 23 1787, and 1883 ( Figure 2). During the 20th century, large floods occurred in 1900, 1920, 1922, 1949, 1957, 1962, 1969 and 2000. The analysis of regional documentary data (i.e., Millars River, Rambla de la Viuda, and other rivers at Castelló region) indicates an alternation of flood-rich and flood-poor periods over the last 500 years. The highest flood frequency occurred between 1580-1620, 1770-1810, 1880-1925, and 1945-1969. In contrast, periods between 1620-1770 and 1810-1880 were characterised by a high irregularity with few reported events or even absence of floods. In the Rambla de la Viuda river, seven large flood events were reported with category 5 over the period 1580-1900, namely in 1580, 1597, 1617, 1776, 1783, 1787, and 1883 ( Figure 2). During the 20th century, large floods occurred in 1900,1920,1922,1949,1957,1962,1969   In the Montlleó River, continuous gauging began in 1971 at a station located 5 km upstream of the study area, and it is included in the Spanish Automatic Hydrological Information Systems In the Montlleó River, continuous gauging began in 1971 at a station located 5 km upstream of the study area, and it is included in the Spanish Automatic Hydrological Information Systems (SAIH) network. The water authority provides, for this station, the daily mean discharge data. These daily discharge values were converted to peak discharges using the CEDEX approach [58] applied at regional level in the Júcar River basin, and modified in one parameter calibrated for Rambla de la Viuda catchment by Beneyto et al.  [61]. Epigraphic marks at the dam wall show that flood stage during the 1969 flood was halfway between 1962 and 2000 levels.
In the Montlleó River catchment, the largest reported flood occurred on the night of 9 October 1883 reaching ca. 6 m water stage in the small village of La Estrella. In this location, this heavy rainfall produced a debris flow that destroyed 17 houses, damaged another 28 buildings, and caused 27 casualties (Newspaper El Liberal, 17 October 1883). Another extreme documented flood occurred in 1787 (October 8) in the nearby Alcora River (Rambla de la Viuda tributary) that destroyed the old Alcora dam built in 1580 ( [62] cited in [63]).

Palaeoflood Hydrology
The stratigraphic profiles along the study reach are located in Figure 1d. Specific flood units discussed within the present work are referred to by the profile character identification (e.g., B1, MLL) and the flood unit number (unit #) indicated in the stratigraphic columns. Figure 3 presents information regarding flood sequences with characteristic lithofacies, sedimentary structures, and interpretation.

Flood Geomorphology and Sedimentology
Along the 6.5 km study reach, a lateral gravel bar was deposited on the river's eastern margin, attached to a terrace 6 m above the river thalweg (Figure 1d). This alluvial terrace is overlaid by flood sands and it surface shows traces of flooding, namely sand crevasse deposits. The stratigraphy in profiles B1 and B2 showed a similar number of flood events, 4 and 5 flood units respectively ( Figure 4). The most complete stratigraphy (B2) accumulated a 1.1 m thick sequence with evidence of at least five flood units, with clear contacts marked by stone lines and bioturbation. The second lower unit consists of a medium sand trough cross-bedding structure indicative of high energy conditions (sequence PLT; Figure 3). Radiocarbon dating of the unit provided an age of cal. CE 1674-1894. The upper three units are composed of medium to fine sand within sequences of parallel lamination and climbing ripples in-drift, with indication of reverse flow direction (CR and PLR sequences). Profile B1 contains four units, the upper two are made of SWDs correlated to the B2 upper set, whereas the lower one comprises fluvial sands with OSL dates~4000 years ago ( Figure 4).
The best SWD outcrop is found within an alcove located on a 20-m cliff exposing Pleistocene gravels on the right valley side (circled A in Figure 1d). The alcove is filled with ca 5-6 m-thick slackwater flood deposits intercalated with layers of alluvial gravels and colluvium. The length of this remarkable outcrop (~14 m), together with additional transverse exposures dug across the alcove, provide a good three-dimensional composition of the SWD architecture ( Figure 5). The frontal outcrop shows a complex flood stratigraphy due to its position on a valley side expansion, which creates a high energy eddy environment, and the recurrent input of gravelly hillslope colluvium often associated to periods of decreasing flood occurrence.        A composite stratigraphic sequence combining profile MLL and MLR shows at least 19 individual flood units within three SWD packages bound by sharp unconformities (Figure 6). Package I overlies fluvial gravels and comprises eight flood units (5-20-cm thick) made of very fine to fine sand (CR sequence; Figure 3) with clear contacts frequently marked by silt and/or clay laminar and/or bioturbation marks (burrows and/or root casts). Unit 1 yielded a radiocarbon age of cal CE 1672-1953 (Table 1), and units 2 and 8 were OSL dated to 0.39 ± 0.03 ka (CE 1595-1655) and 0.45 ± 0.05 ka (CE 1515-1615), respectively ( Table 2).
Water 2020, 12, x FOR PEER REVIEW 12 of 23 A composite stratigraphic sequence combining profile MLL and MLR shows at least 19 individual flood units within three SWD packages bound by sharp unconformities (Figure 6). Package I overlies fluvial gravels and comprises eight flood units (5-20-cm thick) made of very fine to fine sand (CR sequence; Figure 3) with clear contacts frequently marked by silt and/or clay laminar and/or bioturbation marks (burrows and/or root casts). Unit 1 yielded a radiocarbon age of cal CE 1672-1953 (Table 1), and units 2 and 8 were OSL dated to 0.39 ± 0.03 ka (CE 1595-1655) and 0.45 ± 0.05 ka (CE 1515-1615), respectively ( Table 2). Package II truncates and overtops package I. It is composed of seven flood units (units 9 to 15; Figure 6) that stack on a depression space at the downstream side of the alcove. In this set, unit 10 is the most prominent reaching ~1 m in thickness and consisting of fine sand with well-developed trains of climbing ripples with upstream flow direction indicative of a large eddy (sequence CR). OSL dating of unit 10 provided an age of 0.49 ± 0.08 ka (CE 1443-1603). Moreover, a sample at the base of unit 12 was radiocarbon dated to cal CE 1660-1953 (Table 1) and unit 13 was OSL dated to 0.6 ± 0.08 ka (CE 1333-1493; Table 2). Package II truncates and overtops package I. It is composed of seven flood units (units 9 to 15; Figure 6) that stack on a depression space at the downstream side of the alcove. In this set, unit 10 is the most prominent reaching~1 m in thickness and consisting of fine sand with well-developed trains of climbing ripples with upstream flow direction indicative of a large eddy (sequence CR). OSL dating of unit 10 provided an age of 0.49 ± 0.08 ka (CE 1443-1603). Moreover, a sample at the base of unit 12 was radiocarbon dated to cal CE 1660-1953 (Table 1) and unit 13 was OSL dated to 0.6 ± 0.08 ka (CE 1333-1493; Table 2).
Package III (units 16 to 19) overlies a wedge of gravelly hillslope colluvium along the unconformity that truncates Package II ( Figure 6). Its basal unit 17 is composed of a 1-m-thick sand layer with parallel lamination with occasional development of current ripples with downstream flow direction (sequence PLR). Radiocarbon dating of unit 15 provided an age of cal 1677-1940 whereas the OSL date resulted in 1.6 ± 0.1 ka (CE 313-513). This OSL date denotes a poor bleaching of the luminescence signal of the quartz grains. Package III was partially cut by gravel units and sands (Figure 5), most likely reworked from former flood deposits and coarse material that fell from the alcove walls.
A fourth SWD package (IV) overlies a 30-cm thick hillslope gravel unit that truncates most of package III (Figures 5 and 6). Package IV fills an erosion depression formed at the alcove entrance. Package IV flood units are described in profile MLC, showing two main flood depositional sets. The lower six units (MLC units 1 to 6;~10 cm-thick each) consist of fine to medium sand with parallel lamination and reworked stone lines, mainly at the bottom part of the sequence. This lower set is overlain by hillslope colluvium of massive gravel that thin out towards the depression sides ( Figure 5). The upper set consists of five units (MLC units 8 to 12;~5-10 cm in thickness) composed of fine to very fine sand fining upward with parallel lamination, and silty clay laminae ( Figure 6). The lowest unit contains a piece of green plastic indicating a recent age (likely post 1970s).
A cut trench inside AML2 exposes eight flood units overlying slope wash gravel (Figure 7), most likely equivalent to the gravel unit at the base of MLC. The basal units 6 and 7 are thicker and consist of fine sand with PLR type sequences culminated by silty clay laminae. Unit 7 was radiocarbon dated to cal CE 1690-1924 (Table 1). In the overlying units, no radiocarbon dates were obtained. However, these deposits are likely to correspond to 20th century floods like those in MLC.
Water 2020, 12, x FOR PEER REVIEW 13 of 23 Package III (units 16 to 19) overlies a wedge of gravelly hillslope colluvium along the unconformity that truncates Package II ( Figure 6). Its basal unit 17 is composed of a 1-m-thick sand layer with parallel lamination with occasional development of current ripples with downstream flow direction (sequence PLR). Radiocarbon dating of unit 15 provided an age of cal 1677-1940 whereas the OSL date resulted in 1.6 ± 0.1 ka (CE 313-513). This OSL date denotes a poor bleaching of the luminescence signal of the quartz grains. Package III was partially cut by gravel units and sands (Figure 5), most likely reworked from former flood deposits and coarse material that fell from the alcove walls.
A fourth SWD package (IV) overlies a 30-cm thick hillslope gravel unit that truncates most of package III (Figures 5, 6). Package IV fills an erosion depression formed at the alcove entrance. Package IV flood units are described in profile MLC, showing two main flood depositional sets. The lower six units (MLC units 1 to 6; ~10 cm-thick each) consist of fine to medium sand with parallel lamination and reworked stone lines, mainly at the bottom part of the sequence. This lower set is overlain by hillslope colluvium of massive gravel that thin out towards the depression sides ( Figure  5). The upper set consists of five units (MLC units 8 to 12; ~5-10 cm in thickness) composed of fine to very fine sand fining upward with parallel lamination, and silty clay laminae ( Figure 6). The lowest unit contains a piece of green plastic indicating a recent age (likely post 1970s).
A cut trench inside AML2 exposes eight flood units overlying slope wash gravel (Figure 7), most likely equivalent to the gravel unit at the base of MLC. The basal units 6 and 7 are thicker and consist of fine sand with PLR type sequences culminated by silty clay laminae. Unit 7 was radiocarbon dated to cal CE 1690-1924 (Table 1). In the overlying units, no radiocarbon dates were obtained. However, these deposits are likely to correspond to 20th century floods like those in MLC.  On the right side of the alcove, there is evidence of a former SWDs fill, with pockets of flood stratigraphy preserved at the shelter wall at a higher elevation position ( Figure 5). This stratigraphy (AML1; Figure 7) shows at least five flood units consisting of very fine sand with parallel lamination, ripples and climbing ripples in-phase with flow directions both towards the inner and outer alcove sides. Unit 1 was OSL dated to 0.47 ± 0.07 ka (CE 1475-1615; Table 2) and units 2 and 3 were radiocarbon dated cal CE 1646-1810 and cal CE 1690-1926, respectively (Table 1). In the upper section, two flood units are overlying a gravel bed. Although no dating material was found, the deposition of the upper flood unit could be attributed, considering their stratigraphical position and the documented flood records, to the 1962 flood, the largest in the 20th century flood record.

Discharge Estimation
Water surface profiles were calculated for multiple discharges related to the elevation of the slackwater flood deposits along the study reach. In the central part of the alcove, profiles MLC and AML1, containing flood evidence at least beginning in the 20th century provided minimum discharge ranges of 220 and 300 m 3 s −1 , respectively. The flood unit 8 at profile MLC, that were found to contain plastic material, indicates that since the 1970s at least five floods exceeded a minimum discharge of 250 m 3 s −1 . The uppermost flood units at both profiles MLC and AML1, indicate that the largest recent flood, the 2000 event, reached a minimum discharge of 500 m 3 s −1 in the study reach. This discharge is about three times the mean daily discharge recorded in Montlleó gauging station data. Figure 8 illustrates the combined systematic and palaeoflood data used in the flood frequency analysis, where best-estimate ages are provided for the individual palaeofloods. Palaeoflood data provides information related to floods exceeding a given perception threshold in a period of known duration (Figure 8). The perception threshold is defined by the elevation of upper flood bed, therefore after deposition of each flood layer there is a self-rising of the minimum discharge threshold required to accumulate the next flood unit. Consequently, the exceedance probabilities of floods present on the stacked palaeoflood stratigraphy decrease as beds are built up. The sedimentary record is assumed to be complete for each flood exceeding the threshold of discharge at the palaeoflood site. Therefore, palaeoflood data can be treated as censored information, which can be handled efficiently by appropriate statistical methods [64,65]. Palaeoflood data from profiles AML1, B1 and B2, provided critical information regarding the largest palaeofloods (discharge and age) as lower bound data [66]. Moreover, their relative simple stratigraphy reduces the uncertainty of the missing of flood events exceeding the elevation of the depositional site. Additional palaeoflood information from other profiles (MLL, MLC and AML2) were incorporated as censured data or known flood discharges below threshold values.

Flood Frequency Analysis
Lang's stationary test was passed for the palaeoflood and systematic records over the period CE 1590-2000, considering a threshold discharge above 450 m 3 s −1 , (Figure 9). This confirms the assumption of stationarity in time (i.e., all floods were randomly generated from a single probability distribution with stable moments) for the use of statistical parametric models [50]. The null hypothesis is to assume that the distribution of the number of exceedences can be described by a homogeneous Poisson process (i.e., independent and identically distributed random variables; [50]). The test is passed when the accumulated flood line (red line with black dots) stays inside the 90% tolerance interval (95 and 5% quantiles; dashed lines), implying the compliance of the Poisson process hypothesis. Palaeoflood data from profiles AML1, B1 and B2, provided critical information regarding the largest palaeofloods (discharge and age) as lower bound data [66]. Moreover, their relative simple stratigraphy reduces the uncertainty of the missing of flood events exceeding the elevation of the depositional site. Additional palaeoflood information from other profiles (MLL, MLC and AML2) were incorporated as censured data or known flood discharges below threshold values.
Lang's stationary test was passed for the palaeoflood and systematic records over the period CE 1590-2000, considering a threshold discharge above 450 m 3 s −1 , (Figure 9). This confirms the assumption of stationarity in time (i.e., all floods were randomly generated from a single probability distribution with stable moments) for the use of statistical parametric models [50]. Palaeoflood data from profiles AML1, B1 and B2, provided critical information regarding the largest palaeofloods (discharge and age) as lower bound data [66]. Moreover, their relative simple stratigraphy reduces the uncertainty of the missing of flood events exceeding the elevation of the depositional site. Additional palaeoflood information from other profiles (MLL, MLC and AML2) were incorporated as censured data or known flood discharges below threshold values.
Lang's stationary test was passed for the palaeoflood and systematic records over the period CE 1590-2000, considering a threshold discharge above 450 m 3 s −1 , (Figure 9). This confirms the assumption of stationarity in time (i.e., all floods were randomly generated from a single probability distribution with stable moments) for the use of statistical parametric models [50]. The null hypothesis is to assume that the distribution of the number of exceedences can be described by a homogeneous Poisson process (i.e., independent and identically distributed random variables; [50]). The test is passed when the accumulated flood line (red line with black dots) stays inside the 90% tolerance interval (95 and 5% quantiles; dashed lines), implying the compliance of the Poisson process hypothesis. The null hypothesis is to assume that the distribution of the number of exceedences can be described by a homogeneous Poisson process (i.e., independent and identically distributed random variables; [50]). The test is passed when the accumulated flood line (red line with black dots) stays inside the 90% tolerance interval (95 and 5% quantiles; dashed lines), implying the compliance of the Poisson process hypothesis. Table 3 shows the peak discharge associated to different annual exceedance probabilities (return period) both with palaeoflood with systematic data and using only systematic record. Flood frequency analysis using palaeoflood data showed higher peak discharges than those obtained from the gauge record (Table 3). For instance, the calculated 1% annual probability flood (100-years average return interval) from combined palaeoflood and instrumental datasets is 710 m 3 s −1 whereas using only the systematic record is 460 m 3 s −1 . In the case of the 1000-year flood, conventionally used for hydraulic design of infrastructures, frequency analysis with palaeoflood-gauge records show a discharge of 1525 m 3 s −1 , whereas with only gauged data is 1020 m 3 s −1 . Table 3. Flood quantiles for different exceedance annual probabilities in Montlleó River (657 km 2 ) and Rambla de la Viuda (1500 km 2 ) fitted to, firstly, the annual maximum discharge of the systematic record only, and secondly to the combined systematic and palaeoflood discharges. The difference in % represents changes on flood probability introduced by climate variability related to "flood-rich" periods. Rambla de la Viuda frequency analysis for a two-component extreme value (TCEV) [20]. In Figure 10, the upper tail includes palaeofloods with discharge thresholds of 590-620 m 3 s −1 , from which it is known were exceed 7 times in a period of 423 years, the last time being in year 2000. The plotting positions of the palaeoflood and gauged data information are within the confidence limits of the fitted distribution, an indication of a good fit. The low mean square error (MSE = 0.021) value also suggests a good consistency between the flood plotting positions (green dots in Figure 10) and the expected probability values. The oldest flood may correspond to the known catastrophic events of 1580, 1597 or 1616, the latest known as the year of diluvium [67].

Annual
Water 2020, 12, x FOR PEER REVIEW 16 of 23 Table 3 shows the peak discharge associated to different annual exceedance probabilities (return period) both with palaeoflood with systematic data and using only systematic record. Flood frequency analysis using palaeoflood data showed higher peak discharges than those obtained from the gauge record (Table 3). For instance, the calculated 1% annual probability flood (100-years average return interval) from combined palaeoflood and instrumental datasets is 710 m 3 s −1 whereas using only the systematic record is 460 m 3 s −1 . In the case of the 1000-year flood, conventionally used for hydraulic design of infrastructures, frequency analysis with palaeoflood-gauge records show a discharge of 1525 m 3 s −1 , whereas with only gauged data is 1020 m 3 s −1 . Table 3. Flood quantiles for different exceedance annual probabilities in Montlleó River (657 km 2 ) and Rambla de la Viuda (1500 km 2 ) fitted to, firstly, the annual maximum discharge of the systematic record only, and secondly to the combined systematic and palaeoflood discharges. The difference in % represents changes on flood probability introduced by climate variability related to "flood-rich" periods. Rambla de la Viuda frequency analysis for a two-component extreme value (TCEV) [20]. In Figure 10, the upper tail includes palaeofloods with discharge thresholds of 590-620 m 3 s −1 , from which it is known were exceed 7 times in a period of 423 years, the last time being in year 2000. The plotting positions of the palaeoflood and gauged data information are within the confidence limits of the fitted distribution, an indication of a good fit. The low mean square error (MSE = 0.021) value also suggests a good consistency between the flood plotting positions (green dots in Figure 10) and the expected probability values. The oldest flood may correspond to the known catastrophic events of 1580, 1597 or 1616, the latest known as the year of diluvium [67].

Information Content of Palaeoflood Data
The stratigraphic packages described in the main alcove, each comprising multiple SWD sequences, represent similar flood periods to those recorded in documentary sources. These flood deposit packages are bound by colluvial deposits corresponding to periods of lower flood frequency. The multiple insets of sedimentary bodies inside the alcove preserve a continuous record of the medium and large floods that occurred in the Montlleó River over the past 500 years.
The stratigraphy of the slackwater flood deposits discloses two main flood magnitude series. The first series contains the most complete SWD record corresponding to events of low to medium magnitude. It records up to 19 flood events representing floods with estimated minimum discharges between 20-750 m 3  Here the discharge of the year 2000 event, the largest in the gauge record, was recorded with a mean daily discharge of 129 m 3 s −1 . Field evidence from high-elevation slackwater flood deposits of the 2000 event (AML2#8), identified by plastic and other human materials within the sediments, provide a minimum discharge of 520 m 3 s −1 . The palaeoflood data analysis reveals the often insufficient quality of the gauged data in registering large flood events, and the added value of hydrological information obtained from the slackwater flood stratigraphy [11].

Palaeohydrology in the Climate and Environmental Context
The palaeoflood and documentary flood data presented show the alternation, during the past 500 years, of flood-rich and flood-poor periods involving both flood frequency and magnitude. A combined analysis of the information gathered from the documentary archives and from the palaeoflood sedimentary records made it possible to identify four flood cluster episodes: (1) 1570-1620 (Package I), (2) 1775-1795 (Package II), (3) 1850-1890 (Package III), and (4) 1920-1969 (lower set in Package IV). These flood-rich periods, lasting a few decades (30-50 years), are followed by a longer period characterised by a lack of flood deposits and dominated by slope washout sediment deposits within the stratigraphic profiles. These four flood periods correlate in time with others described elsewhere in the Mediterranean based on historical documentary evidence [44,68,69], lake records [70,71], and alluvial records [72,73].
Climatically, periods 1 and 3 correspond to cooler than usual (about 0.3 • C and 0.2 • C) climatic oscillations within the Little Ice Age (1300-1870) (cf. [44,74]). An accepted explanation is that heavy torrential rainfalls developed when anomalous cold air depressions come into contact with humid warm surficial air in the Mediterranean Sea. This type of Mediterranean cyclogenesis causes 75%-80% of the severe rainfall events in the area (mainly in autumn), the remainder being produced by winter-spring Atlantic fronts associated with zonal circulation [75]. Periods 1 and 3 also agree with decades of positive rainfall anomaly during both autumn (SON) and spring (MAM) and floods occurring during all seasons [76].
Period 2 comprises of a 20-years flood cluster at the end of the 18th century (large floods: 1776, 1783, 1787, 1795). This flood cluster coincides climatically with the "Malda Anomaly" (1760-1800) described for the Western Mediterranean as a period of high inter-annual hydrological variability (floods and droughts, [77]. In this high variability rainfall pattern, large floods occurred predominantly (>70%) in the autumn due to Mediterranean associated mesoscale convective systems [76]. A high frequency of extreme rainfall indicates a predominant meridian circulation pattern (low index circulation) which it may lead to cut-off-low situations from the main westerly trough systems of cold air that often move along the Mediterranean Spanish coast [78].
Our palaeoflood record shows that the largest floods (e.g., 1581, 1597, 1617, 1783, 1787, 1883, 1962) occurred within periods of higher flood frequency, i.e., wetter than normal years. Machado et al. [76] describes that wet years produce higher spring and winter rainfalls, meaning that floods have secondary highs in winter or even spring or summer. In our Montlleó record, catastrophic and extraordinary floods occurred during autumn months (SON), whereas ordinary and a few extreme floods take place in winter and spring.
Large magnitude flooding over the late 18th century correlates in time with palaeovegetation and geochemical evidences important changes to land use (deforestation and grazing; [20]). The thick colluvium deposits, usually culminating in sequences of short-lived continuous slackwater flood units reflects drier climate conditions and intense land use changes, particularly the thick hillslope deposits overlying package III ( Figure 5). Indeed, late 19th century coinciding with historical evidence of major economic and land-use changes (deforestation and increase in cultivated land in hillslopes; [79]).

Impacts of Global Change on Flood Frequency
The understanding of the causes of natural variability of flooding and their link to weather patterns and climate is essential for the assessment of the future changes of extreme events. Previous flood rich periods occurred during cooler and wetter climate conditions and only since the late 18th century flood clusters are often linked to periods with high inter-annual climatic variability (floods and droughts). Indeed, the frequency of continuous wetter phases has decreased since early 17th century, whereas the ones with marked annual variability have increased over the last 150 years [76]. Analysis of documentary records suggest a stationary long-term behaviour of high magnitude floods (>300 m 3 s −1 ) over the last 500 years. Conversely, moderate-to-low magnitude events show a trend towards a reduction in frequency over the 20th and beginning of 21st century, in parallel with a reduction of annual maximum flow [8]. The reduction of mid-low discharge floods does not correlate with an equal negative trend on flood losses that rather show a general upward trend during the period 1971-2010, suggesting an increase on socio-economic vulnerability [80].
An important issue in hydrology is to define the risks of large floods. Conventional statistical analysis is usually restricted to a few dozen observations, which are used to estimate the probability of floods that exceed a chance of at least one in a hundred years. In this respect, palaeoflood data contains a wider range of climate variability impacts on flooding including situations of high hydrological variability exacerbated over the last 200 years. Hence, the impacts of climate variability on flood frequency can be evaluated comparing the flood quantiles resulting from analysis of gauged records with those extended by palaeoflood records. In the case of the Montlleó River, calculated flood quantiles (>50 years) from palaeoflood-gauged data shows 45%-50% higher discharges than those based only on instrumental data, the latter reflecting minor inter-decadal climate variability. Similar analysis for the Rambla de la Viuda catchment (1500 km 2 ), including the Montlleó River basin, shows a difference of 15-25% [20]. These results show a higher difference of flood discharge for a given quantile in smaller catchments, which is likely due to the limited extent of convective cells producing large floods.
The analysis of the palaeoflood record has revealed what flood magnitude is really possible over a long time frame that includes the effects of low-frequency climate variability in flooding. Flood quantile changes are, therefore, helpful indicators for understanding the impacts of climate variability on floods and for supporting adaptation to climate change. Our results, supported by scientific analysis, provide quantitative evaluation of flood change to assist in the decision of flood-risk management strategies and to identify low regret actions to cope with flooding [4]. As future flood projections are highly uncertain [1,3,6,7], the study of past extreme floods occurring during a full range of hydroclimatic variability allows a reliable assessment of flood hazards for a flexible adaptation to climate change.

Conclusions
The slackwater flood deposits (palaeofloods) and documentary flood archives of the Montlleó River (eastern Spain) provide evidence for extreme floods and their response to natural climate variability. The palaeoflood record identifies flood-rich and flood-poor episodes over the last 500 years. Four multi-decadal flood clusters were identified in the stratigraphic record, namely at periods 1: 1570-1620 (Package I), 2: 1775-1795 (Package II), 3: 1850-1890 (Package III), and 4: 1920-1969 (lower set in Package IV). Climatically, periods 1 and 3 were characterised by cooler than usual temperature (about 0.3 • C and 0.2 • C) during climatic oscillations of the Little Ice Age (1300-1870). However, periods 2 and 4 correspond to episodes of high inter-annual hydrological variability (floods and droughts) which became more frequent over the last hundred years. A direct consequence is the reduction of mid to low discharge floods over the 20th-early 21st century, which is also observed in the gauge record [8].
Conversely, our palaeoflood record shows a stationary long-term behaviour of high-magnitude floods (>450 m 3 s −1 ) over the last 400 years. In the Montlleó River, at least five floods were recorded as high indicating minimum discharges of 740-950 m 3 s −1 . This contrasts with the largest reported flood in year 2000, with a mean daily maximum discharge of 129 m 3 s −1 . The impacts of flood probability changes induced by climate variability and change were weighed in relation to the catchment area size: the Montlleó River (657 km 2 ) and Rambla de la Viuda River (1500 km 2 ). In the case of the Montlleó River, calculated flood quantiles (>50 years) from palaeoflood-gauged data shows 45-50% higher discharges than those based only on instrumental data, the latter reflecting minor inter-decadal climate variability. Similar analysis for the Rambla de la Viuda catchment (1500 km 2 ), including the Montlleó River basin, shows a difference of 15-25%. These results show the higher difference of flood discharge for a given quantile in smaller catchments, which is likely due to the limited extent of mesoscale convective cells producing large floods. These findings have important implications for flood prevention, risk assessment and adaptation in small catchments of the Mediterranean region. Funding: This research was funded by the Ministry of Science, Innovation and Universities through the project "Assessment and modelling eco-hydrological and sedimentary responses in Mediterranean catchments for climate and environmental change adaptation" EPHIMED (CGL2017-86839-C3-1-R) co-financed with European FEDER funds.