Low-Frequency Atmospheric Variability Patterns and Synoptic Types Linked to Large Floods in the Lower Ebro River Basin

: This study analyzes the atmospheric variability that caused the largest ﬂ oods affecting the town of Tortosa, Spain, in the mouth of the Ebro River (northeast Iberian Peninsula). The Tortosa ﬂ ood database and ﬂ ood marks in the nearby town of Xerta are used to de ﬁ ne the more relevant ﬂ ooding episodes (discharges . 2900 m 3 s 2 1 ) of the 1600 – 2005 period. We explore the atmospheric variability based on low-frequency patterns and synoptic types applying a multivari- able analysis to grids at sea level pressure and geopotential at 500 hPa provided by the twentieth-century V3 Reanalysis Project for the instrumental period (since 1836). Output from the Last Millennium Ensemble Project was used to analyze the sea level pressure over the pre-instrumental period (before 1836). Our analysis includes 33 ﬂ ood episodes. Four synoptic types are related to ﬂ oods in Tortosa since 1836, characterized by low pressure systems that interact with the Mediterranean warm air mass and promote atmospheric destabilization. Flooding in Tortosa is related to relative high values of solar activity, positive Northern Hemisphere temperature anomalies, and NAO in positive phase. This result indicates that the major ﬂ oods are related to zonal atmospheric circulations (west-to-east cyclone transfer). During winter, the main impact of the ﬂ oods is located at the western part of the basin, and the Pyrenean subbasins are affected during autumn. The major ﬁ nding is that similar ﬂ ood behavior is detected since 1600, improving our understanding of past climates, enhancing the knowledge base for some aspects and impacts of climate change, and reducing uncertainty about future outcomes. SIGNIFICANCE STATEMENT: A total of 33 large ﬂ oods ( . 2900 m 3 s 2 1 ) were registered since 1600 in Tortosa, Spain, located at the mouth of the Ebro River (northeast Iberian Peninsula). They occur associated with low pressure systems that interact with the Mediterranean warm air mass promoting atmospheric destabilization. The ﬂ oods in Tortosa are also associated with other important processes occurring at signi ﬁ cantly longer time scales: high values of solar activity, positive Northern Hemisphere temperature anomalies, and NAO in positive phase, indicating that the major ﬂ oods are related to zonal atmospheric circulations. The major ﬁ nding is that we detect similar ﬂ ood behaviors since 1600, improving our understanding of past climates, enhancing the knowledge base for some aspects and impacts of climate change, and reducing uncertainty about future outcomes.


Introduction
A centennial-long analysis of floods is fundamental for adequate hydrological hazard assessment. To reduce the uncertainty of future projections, the IPCC (2014) promotes the collection of instrumental and documentary climate and flood data that covers extended periods. Natural and historical evidence on actual severe and catastrophic floods is extremely valuable for expanding the most extreme flooding phenomena to a millennial time scale (Wilhelm et al. 2012;Schulte at al. 2019a). Identifying long-term flood series contributes to the analyses of climatic and meteorological processes that naturally govern heavy precipitation (Glur et al. 2013;Peña et al. 2015).
Based on flood marks and documentary accounts, Ruiz-Bellet et al. (2015 and Balasch et al. (2019) investigated the historical flood discharges in the Ebro River basin during the last 400 years, obtaining the contribution of the main subbasins: upper Ebro, Cinca, and Segre. Hydrological features such as flood magnitude, seasonality, complementarity of the subbasin contributions, and the historical periods with highest flood magnitude explain the hydrometeorological behavior of the whole Ebro River basin, a first-order basin in the western Mediterranean.
To improve meteorological forecasts of future floods, it is crucial to identify meteorological patterns related to past flooding events (Teale et al. 2017;Peña and Schulte 2020). An area affected by flooding can be related to the size of the meteorological system: regional-scale flooding might be associated with synoptic-scale systems, while local flooding is related to mesoscale systems. However, it is important to bear in mind that these mesoscale systems are often embedded within synoptic-scale systems in areas with heterogeneous topography (Pino et al. 2016).
In the case of the lower Ebro catchment, the main meteorological features that cause catastrophic flood events, according to previous works on the studied area and for other European basins (Barriendos and Martín-Vide 1998;Llasat et al. 2005;Wetter et al. 2011;Barriendos et al. 2014), are related to surface, upper air synoptic, and some subsynoptic conditions, mostly related to large-scale atmospheric patterns. According to the literature, two kinds of events are identified: 1) events that require a large forced lift, which occur near the Pyrenees (Llasat and Puigcerver 1994) and are related to a Mediterranean Sea region synoptic pattern, and 2) events with continuous precipitation (3-30 days) affecting a large area of the basin with an Atlantic Ocean origin (Balasch et al. 2019). Local topographical and mesoscale meteorological conditions turn out to play a relevant role in connecting both type of events (Gilabert and Llasat 2018). The role of snowmelt is important only in the western part of the Ebro basin (Balasch et al. 2019).
The aim of this study is to analyze the atmospheric variability connected with the most severe floods affecting the town of Tortosa, Spain (see Fig. 1), located in the lower reaches of the Ebro River, northeast Iberian Peninsula. We use the flood database of Tortosa and the flood marks in the nearby town of Xerta to define the major flooding episodes (peak discharges greater than 2900 m 3 s 21 ) since 1600 (Balasch et al. 2019). The atmospheric variability is explored here by analyzing low-frequency variability modes and synoptic types. For this purpose, we apply a multivariable analysis to grids at sea level pressure (SLP) and 500-hPa geopotential (Z500) obtained from the Twentieth-Century Reanalysis project (20CR, version 3; Compo et al. 2011;Slivinski et al. 2019). Moreover, our study takes a novel approach previously tested in midsize catchments in the Swiss Alps (Peña and Schulte 2020;Schulte et al. 2019b) and in southeastern Spain (Sánchez-García et al. 2019) to analyze the atmospheric variability related to flooding before 1836, when atmospheric instrumental grid data are not available. To do so, we analyzed SLP for the Ebro basin as simulated by the National Center for Atmospheric Research (NCAR) Community Earth System Model-Last Millennium Ensemble Project (CESM-LME; Otto-Bliesner et al. 2016).

Study area
We analyzed the atmospheric conditions during the floods affecting the town of Tortosa (Fig. 1), located in the Ebro River's lower basin, only 40 km upstream from the river's mouth on the Mediterranean Sea. The Ebro River is 930 km long and drains the northeastern part of the Iberian Peninsula (Fig. 1a), which includes most of the southern side of the Pyrenees Mountain range. The catchment area under study represents 99% of the total Ebro basin area.
The river catchment is a Cenozoic foreland basin of the Pyrenees with a northwest-southeast orientation. It is triangular shaped and has an area of 85 534 km 2 , limited by the Pyrenees to the north, the Iberian Mountain range to the south, and the Catalan precoastal ranges to the southeast (Figs. 1b,c). The average rainfall within the whole basin is 622 mm yr 21 (period 1920-2000), with a strong altitudinal gradient: 1000-1500 mm yr 21 at the highlands and less than 400 mm yr 21 in the center of the basin (Vicente-Serrano et al. 2007;Balasch et al. 2019).
The Ebro River basin can be divided into three main subbasins (Fig. 1c). The upper Ebro encompasses the western area of the basin, from the headwaters to the city of Zaragoza. Large peak discharges in this subbasin are linked to Atlantic fronts that typically cause snowmelt and long-lasting rainfalls during winter and spring. The second subbasin encompasses the Segre-Cinca catchments (the two main tributaries of the Ebro) on the southern side of the eastern and central Pyrenees. High peak discharges in this subbasin are caused by either spring snowmelt or intense autumn precipitation lasting few days. The middle and lower Ebro extends from Zaragoza to the Mediterranean Sea, with the meager runoff of streams coming from the eastern side of the Iberian Mountain range and the northern slopes of the Catalan precoastal ranges. Remarkable floods in the Iberian mountains and in this final stretch of the river are related to extended flash floods during summer and autumn, respectively. It is important to note that the lower Ebro River is characterized by meanders and large floodplains where the flow may be laminated (Del Valle et al. 2007;Ollero 2010). Soil use changes and increases in water demand have severely reduced the annual runoff volume in Tortosa from the mid-twentieth century to the present day (Gallart and Llorens 2004;Delgado et al. 2010;Buendia et al. 2016). Furthermore, large reservoirs were built in the Ebro basin between 1920 and 1960, mainly in the Pyrenean tributaries and in the middle and lower reaches of the mainstream. The impoundment runoff index (calculated as reservoir capacity divided by mean annual runoff) is 57%, which causes a reduction of 25% in the expected peak discharges with return periods of between 10 and 25 years (Batalla and Vericat 2011).

a. Peak discharge series
We used the Tortosa flood database to find days with peak discharges greater than 2900 m 3 s 21 . This threshold is the peak discharge that causes overflow in Tortosa, where at least the lower parts of the town, located in the left bank of the river, are affected.
Tortosa has a daily maximum instantaneous flow series (Qci) that covers the period 1952-2017 (Table 1). The instrumental series of Qci for the Ebro River at Tortosa could be extended back from 1951 to 1913 by using the series of daily maximum flows (Qc; see Table 1). The Qci-Qc method for restoring Qci from Ebro data series is used by Ruiz-Bellet et al. (2015) in the following way: (i) a relationship between Qc and Qci must be established with a linear regression (r 2 5 0.99) between 44 pairs of data (Qc, Qci) of the same day for the period 1952-2012; (ii) the equation [Qci (m 3 s 21 ) 5 1.023Qc 1 61] found with the regression is used to estimate the Qci of those years when only Qc is available .
Moreover, we used historical flood information included in the PREDIFLOOD/AMICME database (Barriendos et al. 2014 to extend the flood series back to 1600 (Monserrate 2013;Ruiz-Bellet et al. 2015). Peak discharges were reconstructed from flood marks in the town of Xerta, 14 km upstream from Tortosa (Table 2; Fig. 2); the catchment in Tortosa is only 0.18% larger than that of Xerta. The flood scale located in the church of Xerta gathers information on  Bayliss and Reed (2001): 1 5 unreliable, 2 5 reliable, and 3 5 very reliable. b Precision is the maximum expected difference (cm) between the flood mark and the actual maximum water height.
the height of the floods of the Ebro River between 1617 and the present day, although the last flood noted is from 1961, because the reservoirs have reduced the impact of other important episodes such as that of 1982. The scale was built on Valencian ceramics toward the end of the nineteenth century, integrating the previous marks, and it was destroyed by gunfire during the Spanish Civil War . It had to be rebuilt later and the 1961 episode added to it. The relative and absolute heights of the flood marks at the scale have been corroborated with other flood marks located at other up-or downstream towns such as Miravet, Tortosa, Móra de Ebre, Ginestar, Flix, Benissanet, and Benifallet (Balasch et al. 2019). The scale has been also used to compare water heights in various places, especially Tortosa where there is continuous information since the 14th century (Miravall 1997;Querol 2006;Curto 2007;Boquera 2008).
Using the hydraulic HEC-RAS 4.1.0 model (USACE 2010) we reconstructed the peak flows registered in the Xerta scale. The model works under gradually varied, steady and mixed flow conditions to solve the conservation of mass and energy equations in one-dimensional form. The model needs information about the observed water surface heights, geometry and roughness of the bed-river channel and the floodplain and the hydraulic regime of the flow at the running departure (subcritical flow). The last flood event appearing in Xerta scale (January 1961), with a measured peak discharge of 4580 m 3 s 21 (MAGRAMA 2015), was employed to calibrate the model. Additionally, we used the two-dimensional hydraulic Iber model (Bladé et al. 2014) based on the solution of 2D Saint-Venant equations using a finite-volume method to evaluate the uncertainty of the model results (Ruiz-Bellet et al. 2017).
b. Synoptic types and low-frequency atmospheric variability related to large floods in Tortosa

1) SYNOPTIC PATTERNS
The data used to explore the synoptic patterns were the SLP and Z500. The daily data of SLP and Z500 grids related to the instrumental period  were taken from 20CR (Compo et al. 2011;Slivinski et al. 2019). This project has greatly benefited from international cooperation under the Atmospheric Circulation Reconstructions over the Earth initiative, which undertakes and facilitates the recovery of global historical instrumental marine and terrestrial weather observations. Additional support was provided by the Global Climate Observing System and the World Climate Research Programme. Using a state-of-the-art data assimilation system, this reanalysis assimilates only surface observations of synoptic pressure into NOAA's Global Forecast System. It also prescribes sea surface temperature and sea ice distribution to estimate (from the surface to the top of the atmosphere) temperature, pressure, winds, moisture, solar radiation, and clouds. The 20CR includes atmospheric variables at different vertical levels with 28 of horizontal resolution from 1836 to the present, which allows putting current atmospheric circulation patterns into a historical perspective.
Synoptic patterns related to floods in Tortosa were analyzed by applying multivariable analysis to SLP and Z500 anomalies. The SLP/Z500 daily anomalies from 20CR were calculated by subtracting the arithmetic mean of the 30-yr reference period (1961-90) for each daily SLP/Z500 value. Note that in our study the anomalies are then weighted by the square root of the cosine of the latitude to account for the changing area of the model grid with latitude. The analysis was applied over the eastern North Atlantic and Europe domain: 308-708N, 308W-308E. A principal sequence pattern analysis (PSPA; Compagnucci et al. 2001) in T-mode (temporal pattern: time as the variable) determined the synoptic types related to large floods in Tortosa in the period 1836-2005. The PSPA was applied to the daily SLP and Z500 anomalies, specifically those related to the flood day plus the previous three days. We based the eigendecomposition of the SLP and Z500 anomalies on the correlation matrix and used the scree-test criterion to decide on the number of components. The orthogonal varimax procedure was used to rotate the components.
2) TORTOSA FLOOD VARIABILITY COMPARED WITH LOW-FREQUENCY ATMOSPHERIC VARIABILITY

PATTERNS AND CLIMATE FORCINGS
The interannual atmospheric variability of the climate for the 1600-2005 was analyzed using SLP anomalies from the climate simulations and validated by comparison with observed and reconstructed data. To this end, the observed atmospheric variability obtained from the SLP grid (Z500 is not used in this period) from the 20CR and reconstructed atmospheric variability from grids provided by Luterbacher et al. (2002) were compared with paleoclimate simulations provided by the CESM-LME (Otto-Bliesner et al. 2016). We examined the 13 full-forcing runs, including greenhouse gases, ozone, tropospheric and stratospheric aerosols, land use, solar radiance, and orbital changes in insolation, each of which was forced with the same observational estimates of historical forcings but initialized by using different atmospheric conditions (however, note that all ocean initial conditions were the same for all ensemble members). The model used was version 1.1 of the CESM Community Atmosphere Model, version 5 (CESM-CAM5), at 28 resolution for the atmosphere and land, and 18 for the ocean and sea ice. The LME spans the years 850-2005, but we analyzed the period 1600-2005 to be consistent with the flood period, and we used the composite of the 13 full-forcing runs.
The low-frequency atmospheric variability patterns were calculated by applying principal component analysis (PCA) in S-mode (spatial pattern: grid points as the variables) to the monthly SLP anomalies from 20CR (Peña et al. 2015). We based the analysis on the covariance matrix and used the scree-test criterion to decide on the number of components. We applied the orthogonal varimax rotation. Furthermore, in line with Peña and Schulte (2020), we examined the accuracy of the PCA representations in the CESM-LME simulations for the pre-instrumental period before 1836 and the Luterbacher et al. (2002) grid for the period 1659-1999.
We based the monthly SLP anomalies for the three analyzed grids on the reference period 1850-1900 [according to Otto-Bliesner et al. (2016)], and anomalies at each grid point were corrected by the square root of the cosine of the latitude. The analysis was applied over the eastern North Atlantic and Europe domain: 308-708N, 308W-308E.
To investigate how flooding in Tortosa compares to radiative forcings, we plotted the peak discharges greater than 2900 m 3 s 21 during the 1600-2005 period and the ordinary (ORD), extraordinary (EXT), and catastrophic (CAT) flood historical series (Barriendos et al. 2014). The flood series were normalized, and the average of them was smoothed using a 22-yr Gaussian filter (INU) to define flood-rich, or flooding, periods (light-brown vertical shading labeled with letters A-D shown in Fig. 5, described in more detail below) and no-flooding periods (white background in Fig. 5, below). We compared these flood temporal distributions with the data series of the paleoclimate records (see Table 3) [i.e., volcanic stratospheric sulfates, simulated Northern Hemisphere (NH) temperature, and atmospheric variability

2355
reconstructed from low-frequency atmospheric variability], described as follows: • For the volcanic response, we adopted the ice-core-based index of Gao et al. (2008). Changes in total solar irradiance (TSI) were prescribed by the composite of the time series (Comp) smoothed by a 22-yr Gaussian filter (f22) from the data series used in the LME: • Wang et al. (2005; the reconstruction is referred to herein as WLS): This is a 1610-2000 CE spectral reconstruction based on a flux transport model of the open and closed flux using the observed sunspot record as the main input. This comes in two versions, (a) a "no background" version that just has TSI variations similar to that seen over a solar cycle today, and (b) a "with background" version with longerterm trends in the solar minimum. First, visual correlations were established to compare and investigate the decadal variability of low-frequency atmospheric patterns, climate proxies, volcanic forcing, and solar activity related to frequency of the major floods in Tortosa.
In a second step, we analyzed the model of regression to evaluate possible correlations (as a synonym of association, not causation). We used INU as the response variable, and the independent variables are total solar irradiance composite (SOL) from seven temporal series (Vieira et al. 2011), volcanic response (VOL), atmospheric variability reconstructed from low-frequency atmospheric variability (NAO), and NH mean temperature composite (TEMP). Furthermore, we used changes in concentration of the principal well-mixed greenhouse gases (GHG) from Schmidt et al. (2011), the KK10 scenario of anthropogenic land cover change from Kaplan and Krumhardt (2011), and the reconstructed agricultural area (pasture and cropland) and associated changes in natural forests, grasslands, and shrublands (referred to as PEA) from Pongratz et al. (2008).
To find small and interpretable models, we would use a selection criterion that explicitly penalize larger models, such as Akaike's information criterion (AIC), Schwartz's information criterion (BIC), Mallows's C p (Cp), and adjusted r 2 (adjr2). We will also look at a selection criterion, cross validated (rss), that implicitly considers the size of the model. We applied the stepwise method to investigate the best models. This method checks going both backward and forwards at every step. It considers the addition of any predictor not currently in the model (forward search), as well as the removal of any predictor currently in the model (backward search). Last, a Monte Carlo simulation was performed to find a set of predictors that best explains the variance in the response variable.
Once the flood pulsations had been visually identified from flood historical time series (light-brown shaded columns in and no-flood periods was calculated to evaluate the variability of these periods with respect to the normal climate.

a. Peak discharge series
In the towns of Tortosa and Xerta 33 large floods (Table 4)  In general, the Atlantic floods of the Ebro take place in winter, whereas in the eastern sector of the basin (Segre-Cinca and Tortosa) floods are more frequent in autumn. The historical documentation consulted makes some reference to snowmelt in most of the Upper Ebro episodes (winter flooding) but is irrelevant in the Segre and Cinca basins (autumn flooding).

b. Synoptic patterns
The PSPA (see section 3) was applied to the SLP and Z500 anomalies on the 112 days related to the 28 flooding episodes during the instrumental period (noting that a flood episode is defined by the flood day and the three previous days). The results showed that four components based on the scree test explain 78% of total variance after the varimax rotation. Figure 3 shows the atmospheric synoptic composites, that is, maps corresponding to typical situations in the atmosphere (synoptic types), which were obtained by averaging the SLP and Z500 grids of dates pertaining to a given cluster. In general, the synoptic types linked to large floods in Tortosa are characterized by low pressure systems from the Atlantic that advect warm and moist air at the low levels of the troposphere with 1) eastern flow from the Mediterranean Sea, which is the principal meteorological cause in flooding in autumn, or 2) southern flow, which is the most predominant cause of flooding in winter. These types of configurations enhance the occurrence of severe mesoscale advective joint to convective systems.  Table 5 show the synoptic types related to major floods in Tortosa for the flood day plus the three previous days.

c. Low-frequency atmospheric variability
The modes of the low-frequency atmospheric variability in the monthly SLP anomalies calculated for Europe  from the PCA in S-mode [see section 3b(2)] were summarized in annual terms as well as for the winter, spring, summer, and autumn seasons (see Table 6 along with Table S1 of the online supplemental material).
We used the loading matrix of PCA analysis to map the atmospheric modes shown in Fig. 4. This loading matrix defines the spatial structure, that is, the weight of the factors, represented by the Pearson coefficient of correlation, at each point of the grid. Climatically, this is the degree of relationship between the different pressure centers. Four components were chosen from the scree test that explains nearly the same seasonal variance, approximately 85% of the low-atmospheric variability in the eastern North Atlantic and Europe. The maps obtained here are like those defined in the comprehensive studies of the low-frequency atmospheric variability, and these are remarkably similar to those presented by Barnston and Livezey (1987). For this reason, we use their nomenclature. The main annual modes are described below.
Remarkable seasonal differences can be observed in Fig. 4. The annual modes are very similar to those observed in winter and spring, but they are different from the summer and autumn modes: the principal mode of atmospheric variability for the warm season differs from the NAO mode in the cold season insofar as it presents a northwest-southeast pressure gradient with positive anomalies extending over the Scandinavian Peninsula and the British Islands, as well as slight low pressure anomalies over the Mediterranean area. The warm NAO mode explains nearly 40% of the summer and autumn total variability. Table 6 shows the yearly indices of the principal modes of atmospheric variability prevailing during the flood events in Tortosa. The simulated and reconstructed atmospheric variability indices before 1836 are captured by the ensemble mean of the SLP from CESM full-forcing experiment and indices since 1836 are from the 20CR. Each yearly value is the average of indices taking the value of the flood month plus the previous 11 months.
Moreover, Table 6 can be used to analyze the possible influences of the indices of low-frequency atmospheric modes in the flooding of Tortosa. For the annual period, the positive phase of the NAO was the most frequent for the analyzed floods (19 of 33 floods), followed by the EA in positive phase (8 of 33). The two low-frequency atmospheric modes are coherent with the two principal synoptic types related to western and southern flows. The seasonal analysis (see Table S1 in the online supplemental material) shown for winter and spring floods shows a dominance of the NAO in positive phase, but the second-place position of EA/WR in negative phase indicates that the passage of Atlantic low pressure centers is the main atmospheric trigger for flooding in Tortosa during the cold season. A reverse pattern is observed for the autumn. The principal mode during the warm season shows a northwest-southeast pressure gradient with positive anomalies extending over the Scandinavian Peninsula and the British Islands, which is different from the meridional gradient observed for the cold season [more information in Folland et al. (2009)]. Thus, the low-frequency atmospheric variability related to flooding in Tortosa is characterized by periods with a dominance of low pressure centers related to Atlantic flows, defined by the autumn NAO and the EA/WR modes in, respectively, their negative and positive phases.  Table 7 presents the association between the annual modes of low-frequency atmospheric variability and the synoptic types related to major floods in Tortosa is presented. We used the Pearson product-moment spatial correlation, which shows a strong correlation between the patterns, although seasonal differences are detected. For the annual period, the SEF pattern correlates significantly (p , 0.001) and negatively (20.935) with the NAO mode and positively (10.808; p , 0.001) with SCAND, denoting southeast flows that promote flooding in Tortosa. The winter and spring season also show the SEF pattern correlating significantly (p , 0.001) and negatively with the NAO (20.838) mode, while the SF pattern correlates with EA/WR (20.755; p , 0.001), indicating anticyclonic blocking over Russia and Atlantic low pressure that are associated with western flows that deepen over the Mediterranean Sea and provoke southwest, south, and southeast flows over the Ebro basin. In summer, the SF synoptic pattern correlates significantly (p , 0.001) and negatively (20.658) with NAO, and the EF pattern positively (10.870; p , 0.001) with the EA/WR mode. The autumn season shows the SEF pattern correlates significantly (p , 0.001) and positively (10.913) with NAO, and the EF pattern with the EA/WR mode (10.810), which involves anticyclonic blocking over northern Europe, while the Ebro basin is affected by the warm and humid Mediterranean air mass.

d. Tortosa flood variability compared with atmospheric variability and climate forcings
Historical flood analysis since 1600 (Fig. 5) illustrates the quasi-cyclic variability of flood intensities and climate proxies. It is important to note that all data constitute a 22-yr Gaussian low-pass filter, including the bars in the reconstructed NAO from CESM-LME and 20CR. Each bar represents the mean of the NAO for the month of the flood plus the previous 11 months. We detected four periods in reference to highfrequency flooding (.2900 m 3 s 21 ) in Tortosa: 1617-43, 1710-87, 1825-84, and 1907-85 (letters A-D respectively in Fig. 5). A visual comparison (see also Table 6) of the Tortosa flood intensities and the NAO provides evidence on the influence of atmospheric circulation dynamics. Furthermore, Fig. 5 suggests that, except for the B event with two minimum temperature peaks, generally high values of TSI, positive  anomalies of the NAO, and positive anomalies of NH summer temperatures correspond to high-frequency flood periods. The lower temperatures in event B responds to volcanic forcing, which causes a cooling effect [see Swingedouw et al. (2017) for details on volcanism effects in climate mechanisms]. This shows the principal divergence between the temperature (negative anomalies) and the TSI (positive anomalies) since 1600. As can be observed in Fig. 5, a 120-yrlong flood disaster gap (.2900 m 3 s 21 ) from 1644 to 1709 can be associated mostly with low solar activity and with negative NAO during the Maunder Minimum (Barriendos 1997;Luterbacher et al. 2001), whereas the gap from 1787 and 1825 coincides at the beginning to low TSI and at the end to relatively high TSI values. This period partially corresponds with the "Maldá" anomaly (1760-1800), characterized by a strong climatic variability period (Barriendos and Llasat 2003). To include additional flooding information, especially during these gaps, Fig. 5 also includes the whole flood chronology of the Ebro River at Tortosa  although for many of the events the peak flow reconstruction is not possible due to the lack of water marks. As can be observed during these two periods ORD and EXT events occurred, but not CAT (Barriendos et al. 2014). Similar conclusions have been obtained when analyzing droughts events in the same area (Tejedor et al. 2019).
We applied a regression analysis to evaluate associations between proxy and floods records [see section 3b (2)]. For this objective, it is necessary to address which models to consider. Model selection involves both a quality criterion and search procedure by considering only models with first-order terms, so no interactions and no polynomials were used. There are seven predictors in the model. So, if we consider all possible models, ranging from using zero predictors to all seven predictors plus the variable response (eight model parameters in total), there are p21 k50 p 2 1 k 5 2 p21 5 2 7 5 128 possible models. For this reason, we often search through possible models in an intelligent way, bypassing some models that are unlikely to be considered good. Table 8 gives us the correlation r and determinant coefficients r 2 between INU and the parameters used for the regression model.
The seven best models are displayed in Table 9 using the stepwise search, indicating that the predictor selected in each step is introduced sequentially, in this case ranging from one to seven predictors.
The predictors selected in each step of the stepwise search are first just SOL, then SOL 1 TEMP in the second step, SOL 1 KK10 1 GHG in the third step, SOL 1 KK10 1 GHG 1 TEMP in the fourth step, SOL 1 KK10 1 GHG 1 TEMP 1 NAO in the fifth step, SOL 1 KK10 1 GHG 1 TEMP 1 NAO 1 VOL in the sixth step, and all predictors in the seventh step. Table 10 shows adjr2, Cp, AIC and BIC selection criteria plus the rsq (coefficient of determination r 2 ) and rss (sum of squares) coefficients and is used to select the best model. Table  10 shows that the best model is model 6, which presents the best adjr2, Cp, and AIC criteria (marked in boldface type in the table). BIC appears as the best criterion in model 4. Model 6 was chosen for the regression analysis, and the predictor PEA was discarded from the analysis.
A summary of regression analysis presenting the parameters from the stepwise regression method is provided in Table 11.
A Monte Carlo simulation was performed to find the parameters of interception model (INT) plus the six predictors that best explain the variance in the response variable. We simulated samples of size n 5 1000 from model 6: where « i ∼ N(0, s 2 5 0.8004); s 2 is the sum of squares error of the original model calculated from stepwise search. The summary of the simulation is shown in Table 12.
From Table 8, the variable that most influences INU is SOL, which accounts for 25% of the variance of the regression model that explains 48% of the variance of INU (see Table 12). These results can be extrapolated to Fig. 5 where visual correlations showed a possible association between flooding in Tortosa and TSI. All variables are positively correlated with INU except KK10. NAO shows no correlation despite being included in the model, and in Fig. 5 flooding seems related to positive phases of the NAO. Table 12 shows that the Monte Carlo simulation improved the stepwise regression model with better results in residual standard error (rse) and adjusted r 2 (cf . Tables 11 and  12). Figure 6a shows a histogram of the simulated values and an overlay of the true distribution. Figure 6b presents the INU (see Fig. 5) composite from historical sources (INUhist) and the fitted values (INUfit) from simulated model 6. The time series show a very similar variability over the time period analyzed. Figure 7 shows the composite of SLP between 1850 and 1900 (reference period; Fig. 7a), flood-rich periods minus the SLP reference period (Fig. 7b), and flood-poor periods minus the SLP reference period (Fig. 7c). The annual composite presents the NAO configuration, with high pressure values over the eastern Atlantic and low pressure in Iceland. The difference of this pattern with high frequency of flooding periods shows a pattern like the NAO with slight differences between the reanalyzed, reconstructed, and simulated SLP grids, while an opposite pattern characterized the poor flooding periods. Note the order-of-magnitude smaller differences for the LME relative to reanalysis and reconstruction in the center and lower panels. If all plots were made with the same scale, the LME would look very different from the other two. This is likely because the LME is an ensemble of simulations; by taking the ensemble mean, a lot of internal climate variability such as NAO variations is averaged out. Only variability that is strongly externally forced will remain in the ensemble mean (e.g., cooling that results from volcanic eruptions).

Discussion
The flood database of Tortosa and the reconstructed floods using the flood marks in the nearby town of Xerta from 1600 to 2005 allow us to analyze the atmospheric variability during well-defined flood periods. We have shown that both the high-and low-frequency atmospheric variability, among other factors, influence flood events (peak flow higher than 2900 m 3 s 21 ) in the lower Ebro basin. This fact has also been observed in neighboring catchments: Llobregat, Cinca, and Segre, among others ( Fig. 1; Pino et al. 2016;Merino et al. 2017;Gilabert and Llasat 2018).
The prevailing synoptic conditions (high-frequency atmospheric variability related to day-to-day changes) show a good correspondence with the hydrological reconstruction. The synoptic types are characterized by Atlantic low pressure systems that interact with the Mediterranean warm air mass, promoting the destabilization of the atmosphere. The atmospheric configuration at the midlevels of the troposphere is dominated by negative geopotential anomalies, showing the presence of a low pressure center or an Atlantic cold front. These configurations enhance the occurrence of convective events due to the temperature difference between the surface (warm and moist air from the Mediterranean Sea) and the middle levels of the troposphere. Furthermore, the stagnation of the synoptic configuration due to the presence of a high pressure system in Europe can provoke long-lasting rainfall over the Ebro basin.
Therefore, we can differentiate between two different types of heavy rainfall that affected the lower Ebro basin since 1836: 2361 1) Advective precipitations during winter and spring with low or medium convectivity linked to a zonal disposition of the synoptic configuration (SF and WF types), which favored the passage of frontal systems from the Atlantic region and produced heavy accumulated precipitation over several days (see also Merino et al. 2017).
2) Highly convective autumn storms associated with a blocking anticyclone in northern Europe (EF and SEF types), which promoted the passage of low pressure centers over the northeastern Iberian Peninsula. Depressions are usually linked with the Atlantic cyclones that develop or become more intense over the Mediterranean Sea, and  (Barriendos et al. 2014, NAO values from Table 6 (this paper; from CESM-LME simulated data for the pre-instrumental period and 20CR for the instrumental period), TSI (Steinhilber et al. 2009), Northern Hemisphere temperature (Peña and Schulte 2020; from CESM-LME simulated data), and volcanic eruptions (Gao et al. 2008) for the period AD 1600-2005. The high flood frequency is marked on the chart in light-brown shading. These periods are defined by floods with discharges greater than 2900 m 3 s 21 and by the 22-yr Gaussian filter of the composite INU (for details see the text). ORD (ordinary): non-overbank flood 1 disturbance; EXT (extreme): overbank flood 1 disturbance 1 damage; CAT (catastrophic): overbank flood 1 damage 1 destruction (see Barriendos et al. 2014). Letters A-D mark the different high flooding periods. J O U R N A L O F C L I M A T E VOLUME 35 they produce long-lasting and intense rainfall due to 1) the high-water vapor content of the Mediterranean air mass, 2) the orographic uplift of these air masses, and 3) the reinforcement suffered by the negative anomalies of temperature and geopotential height that occur at the low-and midlevels of the troposphere (Peña et al. 2015).
Similar results were found by Gilabert and Llasat (2018) in the analysis of 56 catastrophic floods in Catalonia over the 1900-2010 period. Their results showed that flooding was caused mainly by cyclonic (78.6%) and advective weather systems (19.6%).
In the Ebro River basin, the peak discharges have two distinct kinds of behavior within the same series (Ruiz-Bellet et al. 2015;Pino et al. 2016;Gilabert and Llasat 2018). These differences are linked to the two types of precipitation. A more frequent behavior related to flat-peak discharges might be caused by winter/spring advective rainfall of low or medium convectivity, while a less frequent one associated with high-peak discharges might be caused by highly convective summer and autumn mesoscale rainstorms. The flooding in Tortosa might be caused by the synergistic effect of two different types of rainstorm occurring in two different types of catchment within the Ebro basin (Balasch et al. 2019): the floods coming from the western part of the catchment (the ones that affect the upper Ebro basin up to Zaragoza) may be caused by a persistent rainfall precipitating gently for several days on catchments with a slow hydrological response, which then disperse over a great area and produce flat hydrographs but important runoff volumes. Conversely, the floods coming down from the Pyrenees (the Cinca and Segre basins) into the Ebro River downstream from Zaragoza may be caused by heavier rainstorms falling from few days up to a month on very mountainous, soil-saturated, quick-response catchments that are located in a relatively small area and, thus, are more efficient in terms of runoff generation (Balasch et al. 2019).
The modes of the low-frequency atmospheric variability of the 1-month mean SLP anomalies calculated from the PCA in S-mode were summarized in annual terms (month of flood plus previous 11 months) as well as for the winter, spring, summer, autumn seasons. In general, the flood events are related to those years with positive phases of the NAO (19 events) and EA (8 events). Despite the location of Tortosa, at the northeastern part of the Iberian Peninsula a few kilometers from the mouth of the Ebro River on the Mediterranean Sea, the atmospheric variability that triggers flooding in the town is not directly connected to the atmospheric causes related to flooding on the eastern side of the Iberian Mountain range and both slopes of the Catalan precoastal ranges, which are modulated essentially by Mediterranean flows. The contributions of these Mediterranean flows at the study site can act only to strengthen the normal autumn atmospheric situation originating in the Atlantic Ocean region. Thus, this Atlantic regime explains the fact that flooding in Tortosa is more probable during the cold season of the year (December-May) because of the southern migration of the zonal flow in the Northern Hemisphere's atmospheric general circulation. The Ebro basin encompasses 85 000 km 2 , of which more than 40 000 km 2 correspond to the western area upstream of Zaragoza, and it clearly demonstrates Atlantic (not Mediterranean) behavior, which explains why an important part of the floods (around 50%) occur in winter and the rainfall over the western half of the basin is from the Atlantic. Most of the Mediterranean floods are in autumn and affect the eastern half of the basin (Balasch et al. 2019).
It is important to note that the signal of the SLP pattern related to flood-rich periods is very weak (Fig. 7), with very distorted centers of action, especially the low pressure center over Iceland. This may be related to the poor positive correlations (no statistical significance) that have been observed in the north of the Iberian Peninsula between the NAO and rainfall (Martín-Vide and Fernández 2001;Benito et al. 2021). Instead, these correlations in the central and southern part of the Iberian Peninsula are significant but negative. This fact is in concordance with the results of the regression model that show a poor influence of the NAO in the regression model although it is one of the predictors chosen by the stepwise method.
Positive phases of NAO related to floods in Tortosa in our results are sensible if we relate these findings with other  TRUE  TRUE  TRUE  FALSE  FALSE  5  TRUE  TRUE  FALSE  TRUE  TRUE  TRUE  FALSE  TRUE  6  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  FALSE  TRUE  7  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE radiative and climate forcings. Based on visually analyzing the time series of climate proxies, the historical and contemporary flood analysis suggests that flooding periods with atmospheric conditions marked by a positive NAO phase are, in general, related to positive temperature anomalies in the NH and high values of TSI (see Fig. 5). Similar findings were reported by Corella et al. (2014). These authors used an annual reconstruction of extreme rainfall since the fourteenth century in Montcortés Lake (northeastern Spain, 1027 m above sea level), and they found that extreme rainfall variability positively correlated with solar activity and autumn reconstructions of the NAO. Research based on the paleoclimate simulations performed by Haigh (2003) and Gray et al. (2010), among others, provides evidence that the general atmospheric circulation pattern and particularly El Niño-Southern Oscillation (ENSO) in a warm phase (El Niño), can be influenced by high solar activity. The major implications of the solar activity in the atmospheric variability modes were already observed at the end of the twentieth century (Rogers 1984;Polonsky and Sizov 1991): an intensification of the atmospheric circulation in the North Atlantic (positive NAO) is produced just before and during the maximum development of the ENSO warm phase (high solar activity), and a relaxation (negative NAO) is observed during the ENSO cold phase (low solar activity). The configuration of atmospheric general circulation is like that described by various authors, who refer to it as a paleo-NAO pattern (Wanner et al. 2008(Wanner et al. , 2011. This is defined as the decade state of the atmospheric variability in the North Atlantic, and it seems to be the dominant atmospheric variability mode during warm climate pulsations and high-solar activity (Wirth et al. 2013).
Nevertheless, opposite trends are observed in other regions in central Europe and Iberian Peninsula. Paleoflood analysis sited in the Swiss Alps (Peña et al. 2015;Schulte et al. 2015Schulte et al. , 2019bPeña and Schulte 2020) and Duero basin (Benito et al. 2021) have shown that high summer flood frequencies in the Alps and central part of the Iberian Peninsula are mostly related to the negative SNAO [see Folland et al. (2009) for more details about this atmospheric mode] and low solar activity. This discrepancy is in accordance with the findings presented by Corella et al. (2014Corella et al. ( , 2016 on Montcortés Lake, as well as other findings on the lake, fluvial, and historical records in the Iberian Peninsula (Benito et al. 2003;Vaquero 2004;Moreno et al. 2008) and the western Alps (Schulte et al. 2009(Schulte et al. , 2015(Schulte et al. , 2019bWilhelm et al. 2012;Peña et al. 2015). We suggest that complex solar-induced changes in atmospheric circulation patterns reflect a strong regional contrast in the occurrence of flooding events and that the mechanisms related to the correlations between solar activity and NAO dynamics are not sufficiently understood. By developing a multi-archive integration approach in the Bernese Alps, Schulte et al. (2019b) credit the phenomenon of inhomogeneous response of local neighboring catchments to 1) their physiographic parameters, including size, altitude, storage capacity, and connectivity of basins, and 2) their climate parameters, including type, spatial distribution, duration, and intensity of precipitation.
Additionally, our analysis seems to support the hypothesis that the increased frequencies and magnitudes of floods occurred in periods of rapid climate change when, according to Knox (1993), changes in atmospheric circulation patterns can result in changes in magnitude and the recurrence rates of  TABLE 11. Summary of regression analysis, giving results from the stepwise regression method. Significance is indicated as follows: p , 0.001 is indicated by three asterisks, p , 0.01 is indicated by two asterisks, p , 0.05 is indicated by one asterisk, and p , 0.1 is indicated by a Maltese cross. The rse is 0.7984 on 394 degrees of freedom (DF), the multiple r-squared is 0.3721, the adjusted r-squared is 0.3625, and the F statistic is 38.9 on 6 and 394 DF, with p value , 2.2 3 10 216 . However, the construction of reservoirs during the second half of the twentieth century in the Ebro River reduces the runoff peak discharge, which helps to explain the reduction of peak flow during the 1982 flood in Tortosa: discharge changes from 5200 m 3 s 21 just after the Segre-Cinca contribution to only 3780 m 3 s 21 in Tortosa (Batalla and Vericat 2011). These issues play an unquestionably important role in understanding flood processes and managing flood risk but in other catchments this anthropogenic impact does not mask the climatic imprint on floods during the twentieth and twenty-first centuries (Peña et al. 2015;Schulte et al. 2019a;Blöschl et al. 2017;Benito et al. 2021). The dynamics of the North Atlantic contribute with valuable information or risk assessment and management as also indicate recent studies (Schulte et al. 2020;Peña and Schulte 2020;Benito et al. 2021). The systematic evaluation of paleoclimates simulated by models helps to improve our understanding of past and future climates, as well as their impacts. Whereas the modeling of paleoclimates has improved during the past two decades, the simulation of synoptic pattern of past hydrological extreme events is novel and under development (Peña et al. 2015;Schulte et al. 2019a,b;Peña and Schulte 2020). In this respect, our study contributes to understanding the atmospheric variability related to flood frequency over the past 400 years in a large western Mediterranean basin. Moreover, the method used here has been useful for obtaining the low-frequency atmospheric variability modes in the eastern North Atlantic and Europe, as well as for obtaining the synoptic types related to the largest floods in the lower Ebro basin. This has improved both our understanding of floods and their forecasting in the future.

Conclusions
We study the atmospheric variability associated with Tortosa major flood episodes (peak discharges greater than 2900 m 3 s 21 ) since 1600 using a historical flood reconstruction database from the Ebro River in Tortosa and Xerta, data from the Twentieth Century Reanalysis Project , sea level pressure reconstructions (Luterbacher et al. 2002;1659-1999, and the Last Millennium Ensemble Project simulations . The prevailing synoptic conditions (high-frequency atmospheric variability related to day-to-day changes) that favor heavy rainfall explain well the hydrological reconstruction of the flood events. The synoptic types are characterized by Atlantic low pressure systems that interact with the Mediterranean warm air mass that destabilizes the atmosphere due to temperature differences between the surface (warm and moist air from the Mediterranean Sea) and the middle levels of the troposphere. Furthermore, due to high pressure systems located in central Europe causing stagnation of the synoptic configuration, long-lasting rainfall can occur over the Ebro basin.
We detected four periods (or clusters) of high-frequency flooding (.2900 m 3 s 21 ) in Tortosa: 1617-43, 1710-87, 1825-84, and 1907-85. Most of these are related to high solar variability phases (high to low values of solar activity or vice versa), which altogether highlight atmospheric and hydrological instability during periods of rapid climate change.
In the lower Ebro River, the low-frequency atmospheric variability connected to these flood periods is related to the positive phase of the NAO, relative high values of solar activity, and positive Northern Hemisphere temperature anomalies. This provides evidence that complex solar processes might provoke changes in temperature and variability in the atmospheric circulation. The NAO mode shows that the major floods in the region are related to the zonal atmospheric circulation (west to east cyclone transfer). These atmospheric disturbances have a winter effect in the western part of the basin (upstream of Zaragoza), while the Pyrenean subbasins (the Cinca and Segre Rivers) are affected during autumn.
The results from our study help to improve our understanding of the past, present, and future climates, as well as their impacts, thereby enhancing the knowledge base for addressing some aspects and impacts of climate change so as to reduce uncertainty about future outcomes. Furthermore, future investigations seeking to detect and prevent extreme events will find it particularly useful to establish relationships TABLE 12. Summary of regression analysis, giving results of the Monte Carlo simulation results to improve the stepwise regression model parameters. Significance is indicated as follows: p , 0.001 is indicated by three asterisks, p , 0.01 is indicated by two asterisks, p , 0.05 is indicated by one asterisk, and p , 0.1 is indicated by a Maltese cross (but none of the entries fall into that category). The rse is 0.7753 on 394 DF, the multiple r-squared is 0.4767, the adjusted r-squared is 0.4414, and the F statistic is 228.4 on 6 and 394 DF, with p value , 2.2 3 10 216 .  2367 entities: MAGRAMA, who provided the peak discharge data for Tortosa (Ebro River); NSF/CISL/Yellowstone; the DOE INCITE program; the Office of Biological and Environmental Research (BER); NOAA; and the Last Millennium Ensemble Community Project of CESM1 (CAM5), who provided the sea level pressure, geopotential height at 500 hPa, and Northern Hemisphere temperature datasets from the Twentieth Century Reanalysis Project. We also thank the editor Matt Elmore, who checked the spelling and grammar of the submitted paper. In addition, we thank the reviewers, and we thank the editor of the journal Dr. Andrew Hoell for his comments that have helped to substantially improve the paper.