# A numerical study of the February 15, 1941 storm in the Tagus estuary

André Fortunato, Paula Freire, Xavier Bertin, Marta B Rodrigues, Juan Ferreira, Margarida L.R. Liberato

► To cite this version:

André Fortunato, Paula Freire, Xavier Bertin, Marta B Rodrigues, Juan Ferreira, et al. A numerical study of the February 15, 1941 storm in the Tagus estuary. Continental Shelf Research, 2017, 144, pp.50 - 64. 10.1016/j.csr.2017.06.023. hal-01671805

HAL Id: hal-01671805

https://univ-rochelle.hal.science/hal-01671805v1

Submitted on 2 Jul 2025

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1 A numerical study of the February 15, 1941 storm in the Tagus estuary

2

3

4 André B. Fortunato<fnref n="1" />; Paula Freire<fnref n="2" />; Xavier Bertin<fnref n="3" />; Marta Rodrigues<fnref n="4" />; Juan Ferreira<fnref n="5" />;

5 Margarida L.R. Liberato<fnref n="6" />

6

7

8

9 Submitted to Continental Shelf Research: September, 2016

10

Revised: March, 2017

11

2nd revision: June, 2017

<footnote n="1">¹ National Laboratory for Civil Engineering, Av. do Brasil 101, 1700-066 Lisbon, Portugal, afortunato@lnec.pt.</footnote>

<footnote n="2">2 National Laboratory for Civil Engineering, Av. do Brasil 101, 1700-066 Lisbon, Portugal, pfreire@lnec.pt.</footnote>

<footnote n="3">3 UMR 7266 LIENSs, CNRS - Université de La Rochelle, 2 rue Olympe de Gouges, 17000 La Rochelle Cedex, France, xbertin@univ-lr.fr</footnote>

<footnote n="4">4 National Laboratory for Civil Engineering, Av. do Brasil 101, 1700-066 Lisbon, Portugal, mfrodrigues@lnec.pt.</footnote>

<footnote n="5">5 Escola de Ciências e Tecnologia, Universidade de Trás-os-Montes e Alto Douro (UTAD), Vila Real, Portugal, juan.g.ferreira@gmail.com</footnote>

<footnote n="6">6 Escola de Ciências e Tecnologia, Universidade de Trás-os-Montes e Alto Douro (UTAD), Vila Real, Portugal, and Instituto Dom Luiz, Faculdade de Ciências, Universidade de Lisboa, Lisboa, Portugal, mlr@utad.pt</footnote>
12 Abstract

13 On February 15, 1941, a storm caused one of the major natural disasters in the Iberian
14 Peninsula in the past century. The storm made landfall in the north of Portugal, leading to a
15 large surge in the Tagus estuary. Adverse meteorological conditions combined with a high
16 spring tide led to extensive flooding of dry land, causing severe damage and casualties. A suite
17 of regional and local scale models is developed to analyze the event and the relative
18 contributions of the different forcing agents to the extreme water levels. Quantitative and
19 qualitative validations show that the models adequately reproduce this type of events. The
20 models are then used to assess the inundation in the upstream reaches of the estuary where
21 extensive agricultural lands are protected by dikes. Results show that over 25 km² could be
22 inundated today, a value that would increase threefold for a sea level rise of 0.5 m. Then, the
23 relative importance of the different forcing mechanisms on the extreme water levels is
24 investigated through numerical experiments. It is shown that the regional surge and the setup
25 induced by swell are the two main drivers of the inundation. In particular, the modulation of the
26 wave setup by tides induces a semi-diurnal signal which is amplified by resonance inside the
27 estuary.

28 **Keywords:** tide-surge interaction; numerical modeling; SCHISM; flooding; 2010 Xynthia storm;
29 Portugal.
30

# 1. INTRODUCTION

31 Coastal inundations, caused by a combination of tides and storm surges, constitute major 32 natural catastrophes. Recent examples include the devastation of parts of New Orleans by 33 hurricane Katrina in 2005, the inundation of the south-western French coast by the Xynthia 34 storm in 2010 and the ravage of the Philippines by typhoon Haiyan in 2013. Climate change is 35 aggravating this hazard through sea level rise and, in some regions, increased storminess (e.g., 36 Lowe and Gregory, 2005). Simultaneously, the growing settlement of populations in coastal 37 regions increases their vulnerability. There is therefore a strong societal and economic 38 motivation to address this problem.

39 Several avenues have been explored to mitigate this problem. Inundation maps have been 40 developed to determine the hazard and assist in planning (e.g., Dietrich et al., 2010; Fortunato 41 et al., 2013; Bertin et al., 2014; Le Roy et al., 2015; Tomás et al., 2015; Bilskie et al., 2016; 42 Perini et al., 2016). Forecast systems allow the anticipation of extreme events and support first 43 responses (e.g., Zampato et al., 2016; Fortunato et al., 2017). Appropriate solutions are 44 designed to reduce and mitigate the effects of inundations (e.g., Zorndt et al., 2014). On a more 45 fundamental level, a better understanding of the relative importance of the physical drivers 46 responsible for extreme water levels is required (e.g., Melet et al., 2016). This understanding can 47 help us to anticipate worst case scenarios, select the most appropriate tools to study the 48 problem and design better solutions.

49 In this scope, much can be learned from the study of past extreme events. The analysis of their 50 consequences indicates where emergency response is most necessary. The assessment of the 51 effectiveness of the response can help us improve the emergency procedures. Inundation maps 52 for past events can be produced to determine the hazard or to test the effectiveness of different 53 solutions. These events can also help define worst case scenarios, by varying the forcing 54 conditions within reasonable limits (Goutx et al., 2014, Lin and Emmanuel, 2016). Finally they 55 provide severe but realistic settings to analyze the importance of the different physical drivers.

56 In areas where catastrophic events are rare, and thus both the authorities and the general public 57 are less aware of their potential impacts, the search for the most extreme events can lead to a 58 time when little quantitative data were available, and newspapers were the major source of
59 information. In these cases, models are particularly useful, filling in the information gap (e.g. 60 Breilh et al., 2014). Fortunately, the recent availability of long atmospheric reanalyses, such as 61 20CRv2 (Compo et al., 2011), opens the route to simulate and analyze these more ancient 62 events (e.g. Stucki et al. 2015).

63 The Tagus estuary (in the Portuguese nearshore close to Lisbon) is an example of an area that 64 has not suffered any major catastrophic flood in recent years. This estuary is vulnerable to 65 inundations, due to a combination of high urbanization of its margins (Tavares et al., 2015) and 66 extensive low-lying marginal areas. Over the past decades, a few events have inundated the 67 margins and caused moderate damage, but no casualties were reported (Freire et al., 2016). 68 Among them, the 2010 Xynthia storm caused the submersion of some urban areas and the 69 overtopping of dikes. However, these recent storms were relatively mild compared to the one 70 that hit the Iberian Peninsula on February 15, 1941. According to Muir-Wood (2011), this was 71 one of the five most severe wind storms across Europe in the 20th century, from an impact 72 perspective. This storm caused widespread destruction and casualties along the Portuguese 73 coast (Freitas and Dias, 2013), including at least 130 casualties (Muir-Wood, 2011). This author 74 estimated that, if that storm occurred nowadays, the damage would amount to 5 billion euros in 75 Portugal alone. Because this storm occurred at a time when systematic data gathering was still 76 scarce, testimonies from newspapers represent the major source of information.

77 This paper aims at learning from the February 15, 1941 storm. The first goal is to provide new 78 quantitative atmospheric and oceanographic elements about this important event. The second 79 goal is to evaluate the consequences of an extreme storm in the Tagus estuary by determining 80 the inundation extent that a similar storm would cause today, focusing on the extensive low-lying 81 terrains of the left margin of the estuary. Simulations are therefore performed for the 82 atmospheric and ocean conditions of 1941, but for the present bathymetry and mean sea level. 83 The ultimate goal of this paper is to determine the relative importance of the different physical 84 processes on the extreme water levels in the Tagus estuary.

85 The study is conducted by simulating atmospheric, ocean and estuarine dynamics. The 86 atmospheric model WRF-ARW (Skamarock et al., 2008) is used first over the North Atlantic to 87 downscale results from the 20CRv2 reanalysis. The same atmospheric model is used with a 88 finer resolution to provide the forcing for the estuarine wave and circulation model. Results are
then used to force regional wave (WAVEWATCH III® - WW3, Tolman, 2009) and tide/surge
(Semi-implicit Cross-scale Hydroscience Integrated System Model - SCHISM, Zhang et al.,
2016) models to provide the boundary conditions for the estuarine model. Finally, extreme water
levels in the Tagus estuary are obtained with a coupling of SCHISM with a spectral short-wave
model (Wind Wave Model – WWM). The coupled circulation/wave model is denoted SCHISM-
WWM (Roland et al., 2012). All relevant physical drivers are considered in the simulations: tides,
waves, wind, atmospheric pressure and river flow. Numerical experiments are carried out to
determine the inundation extent and the relative importance of the different physical drivers.

This paper is organized in five sections besides this introduction. Section 2 describes the study area and the storms that are studied herein: the 1941 storm and the 2010 Xynthia storm, which is used for model validation. Section 3 presents the data available for model calibration and validation. The models and their application are presented in Section 4. Section 5 presents and discusses the results. The final section summarizes the main conclusions.

## 2. STUDY AREA

### 2.1. The Tagus estuary and its margins

Located in the western Portuguese coast (Figure 1), the Tagus estuary is a mesotidal, coastal plain estuary (Figure 2a). It has a surface area of about 320 km², with a deep, long and narrow tidal inlet connecting the Atlantic Ocean to a shallow, tide-dominated basin. This basin has extensive tidal flats and marshes that cover about 40% of the inner estuary. About 40 km upstream, the estuary markedly narrows at the bay head. In this area, the estuary is bordered by low elevation terrains extensively used for agriculture and generally protected by dikes (Figure 2b). The estuarine bottom is mainly composed of silt and sand, of both fluvial and local origins; marine sands are confined to the mouth and inlet channel (Freire et al., 2007).

Measured tidal ranges vary between 0.55 m and 3.86 m at the coast (Guerreiro et al., 2015) but resonance significantly amplifies the semi-diurnal tidal constituents within the estuary (Fortunato et al., 1999). Simultaneously, the estuary is strongly ebb-dominated due to the large extent of the tidal flats (Fortunato et al., 1999). The estuary is usually partially mixed, although it can be stratified when neap tides are combined with high river discharges (Neves, 2010).
117 The impact of storm surges and waves on the Tagus estuary dynamics is poorly studied. Canas
118 et al. (2009) found atmospheric pressure to be more important than wind in determining the
119 surge within the estuary. Oliveira (2000) applied a wave model to the estuary mouth and found
120 that short waves were dissipated in the ebb-delta sand bars at the entrance of the estuary. Rusu
121 et al. (2011) analyzed the effect of tidal currents on wave propagation at the Tagus estuary
122 mouth, and found the currents to increase the significant wave height and the mean wave period
123 on ebb and decrease them on flood. Inside the estuary, locally-generated waves can reach
124 significant wave heights of 0.8 m. This value was measured between Cacilhas and Seixal
125 (Figure 2a) during a storm event with winds of about 11 m/s (Freire et al., 2009), and simulated
126 by Rusu et al. (2009) for winds of 12 m/s.

127 Included in the Lisbon Metropolitan Area, the estuary comprises along its margins 11
128 municipalities with about 1.6 million of total inhabitants, mostly concentrated in its western and
129 northern sides (Tavares et al., 2015). This densely urbanized water front contrasts with a large
130 extension of productive agricultural areas on the eastern side of the estuary (known as “Lezíria”)
131 and along the three islands formed by alluvial deposition (from South to North: Póvoa, Lombo
132 and Alhandra, Figure 2b). The distinct occupation of the estuary right and left margins is also
133 associated to markedly different topographies. In general, the right margin exhibits a high relief,
134 whereas the left margin is characterized by extensive lowlands. The estuary also hosts a natural
135 reserve (the Tagus Estuary Natural Reserve) that is one of the most important sanctuaries for
136 birds in Europe, covering about 14,000 hectares.

137 The detailed land-use cartography of the estuarine fringe (defined by the Portuguese law as the
138 550 m wide fringe above the water line at the highest astronomical tide), which covers a total
139 area of 130 km², shows that the most important occupation types are the agricultural (35%) and
140 the urban (34%) zones, with smaller contributions of industrial, port and airport facilities (24%),
141 green spaces (6%) and natural areas (1%) (Rilo et al., 2012).

## 2.2. The 1941 and the 2010 Xynthia storms' impacts

142 The most severe wind storm to have hit the Iberian Peninsula in the 20th century made landfall in
143 the Portuguese coast on February 15th, 1941 (Freitas and Dias, 2013) coming from the SW
144 (Figure 1). The atmospheric pressure dropped 45 hPa in 18 hours over Portugal (Ganho, 2013),
145 thus it was an explosive cyclogenesis, reaching minimum values of 952.1 hPa at the
146
147 Geophysical Institute in Lisbon and 957.8 hPa in Monte Estoril (Figure 2a). The wind speed
148 reached almost 20 m/s at the same stations, although these values may underestimate the
149 maximum wind speed over the water due to resistance imposed by the urban environment. The
150 wind bursts in Lisbon reached 38 m/s (Nunes et al., 2012). The information about the river flow
151 is scarce. A major river flood occurred in the Tagus River a few weeks earlier, on January 23-25,
152 and flood marks indicate that it led to the second worst inundation on record in the Lezíria. This
153 was the 4th largest flood recorded at Vila Velha de Rodão (located in the Tagus River near the
154 border between Portugal and Spain, about 200 km from the estuary mouth) since the mid-
155 nineteenth century (Salgueiro et al., 2013), so the river flow was likely still high when the
156 February storm hit the coast. The daily-mean river discharge at Vila Velha de Rodão (the
157 hydrological station closest to the estuary at that time) peaked at 3266 m³/s on February 16.
158 This river flow corresponds to a return period of about 3 years, according to the Tagus
159 watershed management plan. Finally, although a tide gauge was already in operation at
160 Cascais, measurements are unavailable between the 13th and the 15th of February due to
161 malfunction. The tidal range obtained from a harmonic synthesis is 3.0 m (Table 1). This value
162 corresponds to 78% of the maximum tidal range, and to the percentile 91 of the tidal ranges.
163 Wave data are unavailable, but a newspaper cited by Freitas and Dias (2013) mentions 20 m
164 waves in the Tagus estuary. While this value cannot be confirmed, it suggests that wave heights
165 were totally uncommon and probably extremely high.

166 The margins of the Tagus estuary were severely hit. Overall, the human losses in the estuary
167 include 28 casualties, 14 wounded, 125 evacuees and 3 displaced (Freire et al., 2016). Freitas
168 and Dias (2013) report flooding in downtown Lisbon, destruction of workshops and warehouses
169 in the harbor, dozens of smashed boats, many sank ships (150, according to Nunes et al., 2012)
170 and the interruption of the railway between Lisbon and Cascais. Similarly, flooding occurred in
171 Cascais causing severe damage. The walls protecting the northern margin of the estuary
172 upstream of Pedrouços were destroyed (Nunes et al., 2012). In the upper estuary, houses in
173 Alhandra were flooded with more than 1 m of water, and 25 of its inhabitants drowned (Muir-
174 Wood, 2011). Since this town is located on the margins of the Tagus estuary across from the
175 Lezíria, it is expected that the Lezíria itself was flooded as well. According to local reports
176 (Companhia das Lezírias, 1941, cited in Madaleno, 2006), this storm destroyed all the channels
177 and dikes in the Lezíria.
A storm following a similar SW-NE path but clearly less intense over the study area, Xynthia, struck the Portuguese coast on the morning of February 27th 2010 (Figure 1; Liberato et al., 2013; Breilh, 2014). Xynthia became known for causing extreme flooding in the central part of the Bay of Biscay on the dawn of February 28 (Bertin et al., 2012). However, in Portugal, Xynthia was less severe than the 1941 storm and the tidal range (2.9 m) was lower. Like in 1941, the storm coincided with a spring tide, exacerbating its consequences. In the days preceding the storm landfall the river flows were high, and floods are reported in riverine towns, such as Santarém (http://otejo.com/2010/02/25/, last accessed on July 11, 2016). During the storm and in the following days, inundations occurred in several places in the Tagus estuary (Freire et al., 2016) and the Aveiro Lagoon (Lopes, 2016) but no casualties were reported. In the Lezíria, there is evidence of flooding during the 2010 February-March storms (Figure 2b). During three consecutive days, witnesses report that the dikes were overflown during the night and dawn. The dikes that were overflown and damaged during the storm were later repaired and raised. Inundations in the Tagus estuary also occurred in the Seixal downtown area. Freire et al. (2016) determined the extent of the inundation in this area through interviews and photographs. The similarities between the two storms (Figure 1 and Table 1) make Xynthia an excellent case to validate the ability of the model to reproduce storm events and, in particular, to reproduce the extent of inundation.

## 3. DATA

### 3.1. Bathymetry and topography

The extreme water levels in an estuary and the extent of an inundation largely depend on the bathymetry and topography. The generation of a digital terrain model was thus a very sensitive task. The lack of a comprehensive, detailed and updated dataset led to the use of multiple sources of data. These sources are briefly described below. The digital terrain model interpolated onto the finite element grid described below is shown in Figure 2.

The only comprehensive bathymetry of the estuary dates back to 1964/67. This dataset was used as a baseline bathymetry and updated wherever possible. In most of the navigated areas (channels and docks) and in the tidal inlet, recent bathymetries were made available by the Lisbon Harbor authority. Topographic data along the margins of the estuary and LIDAR (Light
207 Detection and Ranging) data along the coastal margins, both from 2011, were provided by the
208 Portuguese Environment Agency and the Portuguese institution in charge of the geographical
209 information (Direção-Geral do Território - DGT), respectively. Upstream of Vila Franca, only
210 cross-channel survey data were available. Since the distance between transects (2.5 km)
211 prevents the direct interpolation of the bathymetry, the mean cross-sectional area was evaluated
212 and a linear regression was determined along the axis of the estuary (Figure 3). This regression
213 was then used to interpolate the depths along the riverine part of the domain. Similarly, the
214 spacing between available cross-sectional data along the upper Sorraia River varies between
215 200 and 1500 m. Given the narrow width of this part of the river it was represented with a
216 rectangular cross-section and the bathymetry was interpolated using available cross-sections.
217 The depth along the Risco River was estimated based of the few data points available.

218 The topography of the margins was mostly based on a 2008 dataset provided by DGT and
219 obtained by aerial photogrammetry. Although this dataset was provided on a 2 m resolution grid,
220 there are some doubts on the accuracy of the dikes' crest heights. The height of the crest of the
221 dikes protecting the Lezíria was provided by a farmers' association (Associação dos
222 Beneficiários da Lezíria Grande de Vila Franca de Xira) and verified in a few places by a
223 dedicated topographic survey. The dikes covered by this dataset extend from the margins of the
224 Sorraia River to the Tagus River. These data were measured after the partial reconstruction of
225 the dikes that followed their damage during the Xynthia storm. The dikes represented in the
226 model are therefore higher than the ones that existed during the Xynthia storm.

## 3.2. Dynamic data

228 Atmospheric pressure measurements from 1941 at two stations in the Tagus estuary region
229 were provided by the Portuguese Institute for Sea and Atmosphere. Pressure dropped to 952.1
230 hPa in Lisbon (Table 1) and is used below to validate the atmospheric model.

231 Water levels were measured at 13 tidal gauges installed along the estuary in 1972. These time
232 series were harmonically analyzed to extract the amplitudes and phases of the tidal constituents
233 (Fortunato et al., 1999). Time series for model validation are obtained by harmonically
234 synthetizing the tidal elevations using the same constituents included in the model.

235 Daily averaged river flows measured at Vila Velha de Rodão in 1941 are available from the
236 Portuguese water resources information system (www.snirh.pt). However, this station is located
over 100 km upstream of the estuary. The estimate of the river flow at the model boundary took into consideration that: 1) the peak river flow at Vila Velha de Rodão at the time of the storm, of the order of 3500 m³/s, corresponds to a return period of about 3 years; 2) for that particular return period, the river flow at the estuary's head (station of Almourol) is 4600 m³/s, i.e., about 30% higher. Hence, the river flow at the estuary's head was estimated as 30% higher than the one measured at Vila Velha de Rodão (Figure 4). Cubic splines were used to obtain a smooth time series.

After 1941, several dams were built in this watershed and the reservoir capacity in the Spanish territory now exceeds the average annual water resource volume (Salgueiro et al., 2013). As a result, the same meteorological conditions that occurred in February 1941 would lead to a different river flow nowadays. Hence, an alternative approach to estimate the river flow was followed. Since the trajectories of the February 1941 and the 2010 Xynthia storms were fairly similar, the river flow for the 1941 storm was estimated as the river flow measured at Almourol (about 40 km upstream of the grid boundary) during the Xynthia storm.

Both approaches followed to estimate the fresh water flow have clear limitations: the first approach ignores the changes in the watershed between 1941 and today, as well as the time lag between the river flow at the Spanish border and the model boundary; the second approach assumes that precipitation was similar in the two events. Not surprisingly, the river flows estimated by the two approaches show significant differences, even though the order of magnitude is the same (Figure 4). Hence, a sensitivity analysis of the model results to the choice of the river flow time series is performed below.

Finally, since there are no measurements of river flow in the Sorraia River for the relevant period, this river flow is taken as 5% of the river flow in the Tagus River, based on the ratios between annual averages in the two rivers.

# 4. NUMERICAL MODELS: DESCRIPTION, APPLICATION AND VALIDATION

## 4.1. Atmospheric model

The atmospheric conditions over the North Atlantic region favoring the explosive development of the 15 February 1941 winter storm were simulated using the mesoscale Weather Research &
266 Forecasting Model - Advanced Research model (WRF-ARW) version 3.5 (Skamarock et al., 2008). For this experiment, the initial and boundary conditions for the large-scale atmospheric fields of pressure, wind, humidity and temperature, as well as initial soil parameters (soil water, moisture and temperature) and sea surface temperature, are given by the 2-degree horizontal resolution 20th Century Reanalysis (20CRv2, Compo et al., 2011).
267
268
269
270

271 In order to obtain a good representation of the complete genesis and development of the storm, 272 a two-way interactive nested domains configuration was setup: a parent domain with a 27 km 273 horizontal resolution and a child domain with 9 km horizontal resolution. The outer domain 274 covers an area that includes both the North Atlantic basin extension from the Gulf Stream region 275 to the Iberian Peninsula and most of the subtropical North Atlantic Ocean (Figure 1). Both 276 domains have 32 vertical unequally-spaced levels, with the model top at 50 hPa. Further details 277 on the simulation of explosive developing storms in the North Atlantic are presented in Ferreira 278 et al. (2016).

279 The main physical parameterizations applied herein to all domains follow the work of Ferreira et 280 al. (2014) who performed a sensitivity test to 24 combinations of physical parameterization 281 schemes (short wave radiation, cumulus, microphysics and surface-planetary boundary layer) in 282 36 periods during 2006, comparing model results of temperature, mean sea level pressure, 283 humidity and wind speed with observations along the continental Portuguese area. Thus the 284 parameterization setup used in the present study includes the WRF Single-Moment 6-class 285 microphysics scheme (Hong and Lim, 2006), Duhdia shortwave radiation (Dudhia, 1989), the 286 Rapid Radiative Transfer Model (RRTM) longwave radiation model (Mlawer et al., 1997), the 287 Noah Land Surface Model (Chen and Dudhia, 2001), the Yonsei University (YSU) boundary 288 layer parameterization (Noh et al., 2003) and the Grell-3D scheme (Grell and Devenyi, 2002).

289 The 20CRv2 dataset provides atmospheric data from 56 ensemble members. For this research 290 a sensitive analysis was made and the atmospheric downscaling was performed to a member of 291 the 20CRv2 reanalysis closer to the pressure evolution observed at the Lisbon meteorological 292 station and corresponded to the member that predicted the most intense storm. The simulation 293 period was between 00 UTC 10/02/1941 and 00 UTC 17/02/1941.

294 The simulations are validated with observational pressure data available at the Lisbon and 295 Estoril meteorological stations (Figure 5). The simulated pressure evolution on both positions
follows the observed values, including the pressure drop on the deepening phase of the cyclone.
However it underestimates the maximum intensification and weakening phase of the storm.
Considering that the reanalysis data are obtained without radiosondes information, the vertical
structure of the conditions during the storm are expected to introduce some errors that for such
nonlinear systems affect the evolution of the storm.

## 4.2. Regional-scale wave model

Time series of directional energy spectra were computed from a regional application of the
WaveWatchIII model (Version 3.14; Tolman, 2009), implemented over the North Atlantic Ocean
with a 0.5° resolution (Figure 1). The spectral space was discretized using 24 equally-spaced
directions and 24 exponentially-spaced frequencies spanning 0.04-0.4 Hz. Sea-ice at high
latitudes was taken into account according to the procedure described in Bertin et al. (2013).
The source terms associated with wind input and dissipation by whitecapping were computed
using the parameterization Cycle4 described in Bidlot et al. (2007). Further details on the model
application and validation are provided in Bertin et al. (2013). The model was forced with time
series of wind fields originating from the dynamic downscaling of the 20CRv2 reanalysis
described above.

Because the validation of the wave model presented in Bertin et al. (2013) focuses on yearly-
mean discrepancies, this hindcast was compared against measured wave heights at Leixões
and Sines, two stations for which have pluri-annual time series are available (Figure 6). This
comparison reveals Root Mean Squared Discrepancies (hereafter RMSD) of 0.32 and 0.38 m,
which yield normalized errors of 16 to 24% respectively. Arbitrarily focusing on wave heights
larger than 4 m, these errors drop to 13 and 19% respectively, with normalized biases smaller
than 5%.

## 4.3. Regional scale tide and surge model

Simulations of tides, surges and waves are performed with the modeling system SCHISM
(Zhang et al., 2016). SCHISM (Semi-implicit Cross-scale Hydroscience Integrated System
Model) is a derivative product built from the original SELFE model (v3.1dc; Zhang and Baptista,
2008), with many enhancements and upgrades including an extension to large-scale eddying
regime and a seamless cross-scale capability from creek to ocean (Zhang et al. 2016). SCHISM
uses semi-implicit finite elements and finite volume methods, combined with Eulerian-
326 Lagrangian methods, to solve the shallow water equations. SCHISM is also fully parallelized. It
327 includes modules for surface waves (Roland et al., 2012), morphodynamics (Pinto et al., 2012;
328 Guérin et al., 2016) and water quality (Rodrigues et al., 2009, 2011; Azevedo et al., 2014).

329 SCHISM was applied in 2D depth-averaged mode to the NE Atlantic (Figure 1). The model is
330 forced by tides from a global tide model (Finite Element Solution - FES2012, Carrère et al.,
331 2012), winds and atmospheric pressure from the WRF-ARW application described above, and
332 tidal potential. Eleven tidal constituents are imposed inside the domain through the tidal
333 potential, and twenty nine tidal constituents are forced at the boundary. The grid has about
334 150,000 nodes and a resolution of about 250 m at the Portuguese coast. The time step is set to
335 5 minutes. A detailed description of the model application and validation is presented in
336 Fortunato et al. (2016).

4.4. Local scale tide, surge and wave model

337 The circulation and wave propagation in the Tagus estuary were simulated with the coupled
338 model SCHISM-WWM (Roland et al., 2012) in depth-averaged mode. Two-way coupling is
339 important, as there is a significant modulation of waves by tidal elevation and currents at the
340 mouth of the Tagus estuary (Rusu et al., 2011). The model was forced by wind and atmospheric
341 pressure from the WRF-ARW simulation, time series of wave spectra from the WW3 simulation,
342 elevations and velocities from the SCHISM simulation and river flows. Radiation stress gradients
343 and wind stresses computed by WWM are used to force the shallow water model, and both
344 elevations and velocities from SCHISM are used in the solution of the wave action equation in
345 WWM. Detailed descriptions on the coupling between SCHISM and WWM are given in Roland
346 et al. (2012) and Schloen et al. (2017).
347

348 The grid was generated with the softwares xmgredit (Turner and Baptista, 1993) and nicegrid
349 (Fortunato et al., 2011), based on grids from previous applications (Guerreiro et al., 2015,
350 Fortunato et al., 2017). The domain is over 110 km long and has a surface area of 1442 km². It
351 extends from the ocean, 27 km away from the estuary mouth, to the river, and includes the
352 agricultural flood-prone areas in the northeastern part of the estuary (Figure 7a). The primary
353 generation criterion was the detailed representation of the geometry. The narrow channels have
354 a finer resolution, to allow the proper propagation of the tide. Nodes and element sides were
355 placed along the crests of the dikes (Figure 7b), to optimize the representation of the elevation
of these sub-grid-scale features (Bilskie et al., 2015). The resulting grid has 140,000 nodes and a resolution between 2 and 800 m (Figure 7a). The time steps were set to 1 minute for SCHISM and 15 minutes for WWM.

Friction was parameterized using a Manning coefficient (Figure 8). This coefficient was selected based on previous experience (Guerreiro et al., 2015), the nature of the estuary bottom and the land cover. In particular, the Manning coefficient was taken as 0.015 m¹³/s in the upper Tagus estuary, where the bottom sediments are muddy. In contrast, values of up to 0.027 m¹³/s were used in the Tagus River, whose bottom sediments are primarily sands. In the dry areas, the Manning coefficient was determined based on the land cover from Chen et al. (2015).

The model's ability to simulate tidal propagation was validated first. A one-year simulation was performed, forced by tides from the regional model of Fortunato et al. (2016) at the ocean boundaries and by average river flows of 517 m³/s and 20 m³/s at the Tagus and Sorraia boundaries, respectively. Biased and unbiased root mean square errors were computed at 13 stations (Table 2). Errors compare favorably with those obtained in previous applications with the same dataset (Guerreiro et al., 2015). Unbiased root mean square errors at high tide were also computed, to assess the model's ability to represent high water levels. Results show that the model reproduces high tides with errors between 10 and 20 cm in the upper estuary (Table 2).

The ability of the model to reproduce flooding during an extreme event was verified by simulating the inundation of the Seixal downtown area during Xynthia, since the inundation limits were available at this location. In this case, the atmospheric model was initialized by results from Interim Reanalyses (ERA-Interim) as obtained from the European Centre for Medium-Range Weather Forecasts Reanalysis (ECMWF). ERA-Interim is the latest ECMWF global atmospheric reanalysis available since 1979 up to the present (Dee et al., 2011). Results show that the model represents inundation accurately, within the limits of grid resolution (Figure 9): ten out of 16 points determined as the limits of inundation fall in the stripe identified by the model as the inundation limit, four other points fall in the last row of wet elements, one point is located in the second row of wet elements and one point is in the dry area. For the same period (February 25 to March 3, 2010), the root mean square elevation error at the Cascais tide gauge is 7.5 cm.
In spite of the good results provided by the validation exercises, some uncertainty remains on the results. Sensitivity analyses were therefore performed on some perceived weaknesses of the model. The errors introduced by the grid resolution were assessed by running the model with a finer grid. In Run 2 (Table 3), all elements were split into four, while preserving the bathymetry. The importance of the uncertainty associated to the river flow was assessed by forcing the model with the river flow measured during Xynthia (Run 3). Finally, the default values of the friction coefficient (0.067 m²/s³) and the wave breaking coefficient (0.78) were used in the wave model. The first value is adequate for fully developed wind seas (Bouws and Komen, 1983), but may be too large for swell dissipation over sandy bottoms (Hasselman et al., 1973), as is the case here. Therefore, a smaller friction coefficient (0.038 m²/s³) was used in Run 4. Similarly, the wave breaking coefficient typically takes values between 0.6 and 0.78. Therefore, a simulation with the lower value was also conducted (Run 5).

The sensitivity of the model was assessed through the maximum water elevations obtained at several tide gauges (Table 4). The station of Cascais was not used. It is located inside a marina, in an area with steep rocky slopes, and the grid resolution is inadequate to properly represent wave breaking. Instead, two stations were considered on both sides of the estuary mouth: Carcavelos and São João (Figure 2). These stations were placed at a depth equal to the mean sea level. Results at these stations are only representative of the region for the same depth, as moving these stations in the cross-shore direction can change the water level by tens of centimeters. Results were also examined at Conchoso, the approximate limit of salinity intrusion under drought conditions (Rodrigues et al., 2016).

In general, the differences between the four simulations are of the order of a few millimeters (Table 4), thereby providing confidence in the results. The maximum differences occur at the edges of the estuary, where they can reach a few centimeters. Refining the grid increases the maximum water level at Conchoso by 5 cm, suggesting that the standard grid may be too coarse in the riverine portion of the domain. Similarly, forcing the model with the river flow from 2010 increases the maximum water level at the same station by 10 cm. These two results suggest that the model may underestimate the extreme water levels and thus the inundation in the upstream part of the domain. However, this problem appears to be localized since the sensitivity of the model at the following downstream station (Vila Franca) is modest. At the coast (Carcavelos and São João) the model exhibits some sensitivity to the friction coefficient of the
416 wave model. However, the sensitivity of the model results inside the estuary to this coefficient is
417 negligible. The sensitivity to the wave breaking coefficient is also negligible.

418 **4.5. Scenarios**

419 Several inundation scenarios were considered. Scenario 1 (Run 1, Table 3) assumes that the
420 1941 storm would occur in 2016, i.e., considering the 2016 mean sea level (2.26 m above the
421 Chart Datum, CD).

422 Scenario 2 (Run 5, Table 3) is a “worst case” scenario, in which the 1941 storm is combined with
423 a high tidal range. The largest predicted tide in a 20 year period (2001-2020) was selected. It
424 occurred on March 1st 2006 and had a tidal range of 3.6 m.

425 Finally, scenario 3 (Run 6, Table 3) takes into account climate changes predicted for the end of
426 the 21st century. Contrary to some areas of the world, the Portuguese coast is not experiencing
427 an increase in storminess. Almeida et al. (2011) did not find any evidence of a change in
428 storminess in the Portuguese south coast between 1952 and 2009. Marcos et al. (2011)
429 suggested that the storm surges in the Iberian coast will decrease slightly in the future.
430 Vousdoukas et al. (2016) reached a similar conclusion for southern Europe and for most climate
431 scenarios. These results are in agreement with the ones obtained for the winter extra-tropical
432 cyclone activity, over this region, where no statistically significant signal is found neither
433 considering all cyclones nor only the most severe systems (Ulbrich et al., 2013). Hence the only
434 change considered herein is the mean sea level rise (SLR).

435 Several values have been proposed for the SLR by the end of the 21st century. Antunes and
436 Taborda (2009) projected a SLR of 0.47 m between 1990 and 2100 based on data from the
437 Cascais tide gauge. Lopes et al. (2011) estimated a local sea level rise of 0.42 m between 1980-
438 1999 and 2091-2100 using results from a climate model for scenario SRES-A2. Also, analyses
439 of the Cascais tide gauge data indicate that the SLR rate is consistent with global estimates: for
440 instance, the mean sea level rose by about 0.15 m in the 20th century (Antunes and Taborda,
441 2009). This consistency may be due to the negligible vertical movement of the area: GPS-based
442 methods indicate that this area has an upward vertical movement which is an order of magnitude
443 smaller than SLR (www.sonel.org/spip.php?page=gps&idStation=648.php). The agreement
between the SLR at Cascais and the global values justifies the use of global estimates of SLR (Guerreiro et al., 2015). The most recent estimates of the mean SLR for the end of the 21st century vary between 0.37 and 0.82 m (IPCC, 2013). Therefore, we adopted the value of 0.5 m as a suitable representative value of SLR between present (2016) and the end of the century.

448

# 5. RESULTS AND DISCUSSION

449

## 5.1. Regional scale atmospheric and oceanographic behavior

Reanalysis data and simulations used for the assessment of the synoptic evolution and dynamical characteristics of the 1941 windstorm allow retrieving the developing cyclone at 00 UTC 15 February on the subtropical North Atlantic near 20°W and 35°N (Figure 1). From this location it experienced a notorious strengthening with a deepening rate over 22 hPa/12h, a value that after geostrophic adjustment to the reference latitude of 60°N increases to 32 hPa/12h, showing how exceptional this event was. In the upper levels of the model grid a meandering, intense jet stream (over 70 m/s at 250 hPa) triggered and propelled the storm to the northeast, similarly to Xynthia (Liberato et al. 2013). The storm made landfall 12 hours later in the north of Portugal, some 350 km to the north of the Tagus estuary. Maximum winds close to 30 m/s hit the south of Portugal coming from the South and Southwest. The cyclone then crossed the Iberian Peninsula and moved northeast, heading to the Bay of Biscay and then to the British Islands, where it weakened.

As previously noted, the 1941 and 2010 Xynthia storms have similar trajectories and lifecycles, though the former storm was even more intense. The reduced number of such extreme cyclones affecting the Iberian Peninsula prevents the establishment of a return period for these events. However, a climatological assessment to characterize the large-scale atmospheric conditions and cyclone tracks during the top-100 potential socioeconomic losses over Iberia associated with wind events, based on 65 years of reanalysis data from the National Centre for Environmental Prediction/National Centre for Atmospheric Research, showed that the lowest minimum mean value for the event day (±1 day) was 966 hPa (Karremann et al., 2016). This analysis indicates that the return period of the 1941 storm is probably higher than 65 years.

462

463

464

465

466

467

468

469

470
471 The cyclonic winds generated energetic waves to the south of the cyclone's track directed to the
472 north-east, moving to the Portuguese coast. The highest waves, with significant wave heights of
473 13-14 m and peak periods around 14 s, hit the coast close to the Tagus estuary mouth. These
474 significant wave heights are exceptional, with a return period of over 100 years (Figure 10a). In
475 contrast, the north of Portugal was hit by smaller waves, with return periods below 10 years.

476 In spite of the severity of the storm and resulting waves, the maximum elevations at the coast
477 are relatively modest. The maximum elevations at the western Iberian coast (neglecting the
478 wave setup) correspond to return periods varying between 2 years in the south and 20 years in
479 the north (Figure 10b). This result is explained by the low tidal amplitude, which corresponds
480 only to the percentile 91 of tidal ranges. The growth of the return period from the South (where
481 onshore winds are maximal) to the North (where the atmospheric pressure is minimal) is
482 consistent with the higher importance of the atmospheric pressure relative to the wind in the
483 generation of storm surges in this region. Indeed, the effect of wind stress on the storm surge
484 decreases with increasing water depth. Since the Portuguese continental shelf is very narrow,
485 the wind setup is modest even for very strong winds.

486 5.2. Inundation in the Tagus estuary

487 The maximum water level obtained at Cascais in Run 1 (hence including wave setup) is 2.08 m
488 relative to present mean sea level. Considering the analysis of extreme water levels based on
489 data from the same station (Guerreiro et al., 2015), this water level corresponds to a return
490 period of just over 20 years. This return period is significantly higher than those obtained in
491 Figure 10b at the two points to the South and North of Cascais (about 2 and 10 years,
492 respectively). This contrast shows that the waves that occurred in 1941 significantly contributed
493 to the storm surge.

494 Results show that a storm similar to the 1941 event could cause today some inundation of the
495 SW area of the Lezíria (Figure 11a), and that the extent of the inundation would be higher than
496 during Xynthia (Table 5). The portion of the dike that was overflown during Xynthia, to the North
497 of the Lombo Island, would be overflown again, even though the crest was raised by about half a
498 meter after the Xynthia event. However, the flooded area would be small (1.6 km², Table 5) due
499 to the presence of a second dike (Figure 2b). The islands of Póvoa and Alhandra would be
flooded as well. On the left margin of the Sorraia River, an area of about 18 km² would be flooded due to the discontinuity and low height of some of the dikes.

The reliability of these results may be limited by the available topographic data, which may be particularly inaccurate around the dikes' crests. To assess the sensitivity of the results, an additional simulation (Run 17, Table 3) was performed in which the height of the dikes' crests was uniformly raised by 20 cm. Results indicate that, in most areas, the inundation extent is not very sensitive to errors of this magnitude, thereby providing further confidence in the estimated inundation areas (Table 5). A related source of uncertainty is the dikes' erodibility, which was not considered in the simulations. The dikes are fragile and require frequent repairs. For instance, in May 2016 a small breach occurred in one of the dikes that protect the Póvoa Island and rapidly grew over the following months. As a result, a significant portion of the island became frequently flooded. In the case of an extreme event, it is likely that the dikes would give in at some places (as was the case during Xynthia), further exacerbating the inundation.

The inundation extent would increase if the tidal amplitude was higher (Figure 11b). In this case, the extreme water level at Cascais increases to 2.39 m relative to present mean sea level, which corresponds to an estimated return period above 1000 years. The second dike in the Lezíria would be overflown increasing the flooded area to 2.5 km², all three islands would be flooded, and the inundated area in the Sorraia left margin would increase by over 50%.

Finally, if the 1941 conditions are combined with a 0.5 m rise in mean sea level, the problem would be significantly aggravated (Figure 11c). Almost 80 km² of land would be flooded, in spite of the presence of the dikes.

## 5.3. Physical drivers of extreme water levels

The effect of different forcing agents on the maximum water levels was determined at several stations throughout the estuary (Figure 2a) in a series of numerical experiments (Run 8 to Run 14, Table 3). The surge was computed by subtracting the time series obtained by forcing the model only with the tide and a climatological river flow (Run 8) from those obtained in the standard simulation (Run 1). The maximum value of the surge was then extracted at each station (black bars in Figure 12). At the upstream stations, analyzing the whole time series could
mask the effect of the different processes because the maximum water level can be higher on February 16 due to the growth of the river flow (Figure 4). Hence, the analysis was only performed for the second high tide of February 15.

In general, the total surge is amplified as it propagates along the estuary, growing from about 1 m at the inlet mouth (P. Arcos) to 1.3 m at Alcochete. This growth may seem due to resonance, which amplifies semi-diurnal tides by 25% between the same two stations (Fortunato et al., 1999). However, the time scale of the passage of the pressure front is about one day. Since the resonance period of the estuary is about 8 to 9 hours, diurnal signals are not affected by resonance (Fortunato et al., 1997). In order to find an alternative explanation for the surge amplification, the individual contribution of each driving process was analyzed.

The importance of each forcing was determined by switching off the forcing agents individually (Table 3), and again evaluating the surge by subtraction of the results of Run 8. Switching off the forcing agents individually, rather than including them individually, maximizes the reproduction of the non-linear interactions. Results are summarized in Figure 12.

The impact of the surge at the ocean boundary on the surge in the estuary is fairly constant in space. This behavior is similar to diurnal tidal constituents, whose amplitude is almost constant throughout the estuary. Given the modest amplitude of this component of the surge (30-40 cm in the estuary) compared to the tide, and its long period, the associated velocities and the corresponding frictional losses should be small. The surge is somewhat amplified in the upstream part of the estuary (Vila Franca and Conchoso), probably due to the convergence of the estuary cross-section.

The swell coming in from the ocean boundary has a significant impact on the estuarine water levels, with a maximum wave setup of 33 cm at the inlet mouth (P. Arcos). The development of a substantial wave-induced setup inside estuaries and lagoons during storms, like the one observed herein, has been identified in previous studies, using both models and observations. During storms, radiation stress gradients can be a leading term in the momentum balance in a tidal inlet (Olabarrieta et al., 2011, Wargula et al., 2014). As a result, the effect of wave-induced setup on extreme water levels can reach tens of centimeters inside estuaries and lagoons (Xie et al., 2008; Bertin et al., 2009; Malhadas et al., 2009; Olabarrieta et al., 2011). Bertin et al. (2015) explain this phenomenon by huge onshore-directed wave forces induced by the breaking
of storm waves over ebb deltas and rocky shoals. Our results show that these large wave-induced setups can also occur in large estuaries, such as the Tagus.

Furthermore, this maximum wave setup of 33 cm at the inlet mouth grows to about 50 cm in the wide part of the estuary. Yet, a simulation forced only by constant 10 m waves at the ocean boundary (i.e., with the tides switched off and the water elevation fixed at mean sea level) shows a homogeneous wave setup of 20 cm within the estuary. Hence, the growth of the wave setup is attributed to the modulation of the wave setup by the tidal elevations and currents. A similar modulation of the wave setup by tides was proposed by Olabarrieta et al. (2011).

This modulation effectively generates a semi-diurnal pattern in the setup. This hypothesis is confirmed by forcing the model only with tides (M2 and harmonics), a constant wave from the West (significant wave height: 10 m; peak period: 15 s) and a constant river flow (Run 15). Results show that the wave setup grows significantly during ebb and decreases on flood (Figure 13). Two physical processes can contribute to this tidal modulation. First, ebb currents increase the wave height and steepness (Dodet et al., 2013), thereby shifting the position of the breaking point seaward. Secondly, a significant fraction of the highest waves break on the terminal lobe of the ebb delta at low tide, generating strong radiation stress gradients and a high setup. At high tide, this fraction decreases. The waves that do not break over the ebb delta propagate into the inlet or to the coast, and contribute less to the setup inside the estuary. Switching off the effect of tidal currents on wave propagation has a relatively small impact on wave setup (Figure 13a), showing that the latter process clearly dominates. This result is confirmed by running simulations forced by waves alone, but with different mean sea levels. The setup increases to 24 cm when the mean sea level is lowered by 1.5 m, and it decreases to 18 cm when the mean sea level is raised by 1.5 m.

While the varying widths of the surf zone during the tidal cycle elucidate the higher setup at low tide, they do not explain the phase lag observed between the tide and surge signals shown in Figure 13. A possible interpretation is the following. The waves generate very strong currents in the inlet mouth. These currents increase the bottom stress, thus reducing the celerity of the tidal wave. Indeed, the high tides in Run 16 occur about 4-5 minutes later than in Run 15. This phase lag is then at least partially responsible for the phase shift between the surge and the tide in Figure 13.
This tidal modulation is amplified by resonance within the estuary (Figure 13a): the maximum wave setup grows from 30 cm at P. Arcos to 38 cm at Alcochete in Run 15. This behavior of the wave setup contributes therefore to the growth of the surge within the estuary. In the upstream part of the estuary, the impact of the swell declines, although it remains non-negligible (above 10 cm) even upstream of the tidal intrusion limit.

Because the higher setup occurs close to low tide, the impact of the wave setup amplification on the maximum water elevations inside the estuary is small. In contrast, the surge coming from the ocean boundary is independent from the tidal phase and plays therefore a more important role in determining the extent of the inundation.

Wind effects can also be significant. In the wide part of the estuary (Alcochete, Ponta da Erva), they can even be the dominant mechanism for the surge. Not only is the estuary particularly wide (exceeding 10 km near Alcochete), but the fetch can exceed 20 km for some particular wind directions. Furthermore, the wind setup peaks at high tide when the fetch reaches its maximum. The wind setup also propagates upstream and is still significant in the riverine portion of the estuary.

The other forcing agents are generally unimportant. The effect of the river flow decreases downstream. In particular, it is negligible from the wide part of the estuary downstream. This behavior was previously identified by Vargas et al. (2008). Local waves are small, only exceeding 1 m in the deeper channels of the upper estuary. Hence, their effect is modest (up to about 5 cm). Similarly, the effect of atmospheric pressure gradients is negligible at these small scales.

Non-linear effects are seen in Figure 12 as the difference between the total surge and the sum of the individual contributions. These effects are significant and contribute to the surge everywhere inside the estuary. They are particularly important in the riverine part of the domain. Guerreiro et al. (2015) analyzed the effect of sea level rise on the Tagus tidal propagation and concluded that it shifts the resonance period to lower frequencies, thereby increasing the resonance amplification of semi-diurnal tides. The non-linear effects shown in Figure 12 probably have the same origin: each factor increases the mean water level inside the estuary, and thus amplifies the tidal amplitude.
617 **6. SUMMARY AND CONCLUSIONS**

618 A hindcast of the strongest storm that occurred in the 20th century in the Iberian coast was
619 carried out using a suite of atmospheric, oceanic and estuarine models. The models were
620 carefully validated with available data, and the validation was complemented by sensitivity
621 analyses. Overall, the models are shown to provide an accurate simulation of the coastal
622 hydrodynamics, and, in particular, of the water levels and inundation in the Tagus estuary. Since
623 the peak of the atmospheric pressure is overestimated by about 8 hPa, the surge due to the
624 atmospheric pressure alone should be underestimated by at least 8 cm. The ability of the model
625 to reproduce a storm that occurred so many years ago, in 1941, opens interesting perspectives
626 to investigate other old extreme events.

627 The 1941 storm was one of the most extreme cyclones in terms of impacts on southwestern
628 Europe in the 20th century (Muir-Wood, 2011; Karremann et al., 2016). Its genesis occurred on
629 the subtropical North Atlantic and intensified rapidly on 15 February near 20°W and 35°N. The
630 wind storm made landfall during the period of maximum deepening in Portugal, with winds close
631 to 30 m/s in the south of Portugal. The cyclone then crossed the Iberian Peninsula and moved
632 northeast, heading towards the Bay of Biscay and the British Islands.

633 The cyclone generated waves with a significant wave height close to 14 m to the south of its
634 track. These waves were particularly severe in the Tagus estuary area and further south, where
635 their return period was over 100 years.

636 In spite of the severity of the storm and the resulting waves, the maximum sea levels only reach
637 the thresholds associated to return periods between 2 and 20 years (increasing northward).
638 These relatively modest return periods are due to the tidal range during which the storm
639 occurred, which corresponds to the percentile 91 of the tidal ranges.

640 Several physical processes contribute to the extreme water levels in the estuary. The storm
641 surge at the ocean boundary contributes with 30 to 50 cm to the surge in the estuary. The signal
642 propagates almost unchanged along most of the estuary, and is amplified at the bay head,
643 where the estuary narrows. The wave setup associated with the incoming swell is the second
644 most important source of surge in the estuary. This setup is modulated by tides and amplified by
645 resonance through the following processes. At low tide, the highest waves break on the terminal
lobe of the ebb delta. The surf zone is narrow, the radiation stress gradients are high and so is the setup. At high tide, a larger fraction of the waves do not break over the ebb delta, but rather propagate to the coast or into the inlet, and their contribution to the wave setup inside the estuary is smaller. This variation of the wave breaking conditions along the tidal cycle thus generates a periodic pattern in the wave setup at the estuary mouth. As the wave setup propagates along the estuary, it is amplified by resonance, much like the tidal wave, because it shares the same period as the tide. Finally, wind effects can also be important in the upper estuary, where the fetch is longest and depths are small. The effect of the river flow is restricted to the narrow, riverine-like part of the estuary.

The models were used to determine the inundation that a major storm would cause today in this estuary. It is shown that there is a submersion hazard of the fertile lands in the upper estuary, in spite of the extensive dikes that surround them. The total area threatened by inundations grows from 25 km², if the 1941 event were to happen today, to 77 km² if sea level rises by 0.5 m. A present-day worst case scenario, obtained by combining the 1941 storm with an extreme tidal range, would flood over 40 km² of agricultural land. These results will support risk analyses and the definition of adaptation measures to cope with future extreme events.

662 ACKNOWLEDGEMENTS

ABF, PF and MR were partially funded by the European Commission through the H2020 project BINGO (Grant Agreement Number 641739). MR was also partially funded by a post-doctoral grant from FCT – Fundação para a Ciência e a Tecnologia (SFRH/BPD/87512/2012). XB received funding from the European Commission through the FP7 project Risc-Kit (Grant Agreement Number 603458). JF and MLRL were partially supported by FEDER funds through the COMPETE (Programa Operacional Factores de Competitividade) programme and by national funds through the FCT project STORMEx FCOMP-01-0124-FEDER-019524 (PTDC/AAC-CLI/121339/2010). This work took advantage of the Portuguese National Distributed Computing Infrastructure (INCD). The authors would like to thank Drs. Elsa Alves and Rui Rodrigues for useful discussions about the river flow. The constructive comments of Diogo Mendes and two anonymous reviewers are also acknowledged. The authors thank the various institutions that provided data, models and model results mentioned in the text.
675 REFERENCES

676 Almeida, L.P., Ferreira, O., Vousdoukas, M.I., Dodet, G., 2011. Historical variation and trends in storminess along the Portuguese South Coast, *Natural Hazards and Earth System Sciences*, 11/9: 2407-2417. DOI:10.5194/nhess-11-2407-2011.

679 Antunes, C., Taborda, R., 2009. Sea level at Cascais tide gauge: data, analysis and results. *Journal of Coastal Research*, Special Issue 56: 218 – 222.

681 Azevedo, A., Oliveira, A., Fortunato, A.B., Zhang, J., Baptista, A.M., 2014. A cross-scale numerical modeling system for management support of oil spill accidents, *Marine Pollution Bulletin* 80, 1-2: 132 – 147. DOI: 10.1016/j.marpolbul.2014.01.028

684 Bertin, X., Fortunato, A.B., Oliveira, A., 2009. A modeling-based analysis of processes driving wave-dominated inlets, *Continental Shelf Research*, 29/5–6: 819–834, DOI: 10.1016/j.csr.2008.12.019.

687 Bertin, X., Bruneau, N., Breilh, J.-F., Fortunato, A.B., Karpytchev, M., 2012. Importance of wave age and resonance in storm surges: The case Xynthia, Bay of Biscay. *Ocean Modelling* 42/1: 16-30. DOI: 10.1016/j.ocemod.2011.11.001

690 Bertin, X., Prouteau, E., Letetrel, C., 2013. A significant increase in wave height in the North Atlantic Ocean over the 20th century, *Global and Planetary Change*, 106: 77-83. DOI:10.1016/j.gloplacha.2013.03.009.

693 Bertin, X., Li, K., Roland, A., Zhang, Y.J., Breilh, J.-F., Chaumillon, E., 2014. A modeling-based analysis of the flooding associated with Xynthia, central Bay of Biscay, *Coastal Engineering*, 94: 80-89. DOI: 10.1016/j.coastaleng.2014.08.013

696 Bertin, X., Li, K., Roland, A., Bidlot, J.R., 2015. The contribution of short waves in storm surges: two recent examples in the central part of the Bay of Biscay, *Continental Shelf Research*, 96: 1-15. DOI: 10.1016/j.csr.2015.01.005

699 Bidlot, J.R., Janssen, P.A.E.M., Abdalla, S., 2007. A revised formulation for ocean wave dissipation and its model impact, *Tech. Memo*, **509**, 29

701 Bilskie, M.V., Coggin, D., Hagen, S.C., Medeiros, S.C., 2015. Terrain-driven unstructured mesh development through semi-automatic vertical feature extraction, *Advances in Water Resources*, 86/A: 102-118. DOI: 10.1016/j.advwatres.2015.09.020

704 Bilskie, M.V., Hagen, S.C., Medeiros, S.C., Cox, A.T., Salisbury, M., Coggin, D., 2016. Data and numerical analysis of astronomic tides, wind-waves, and hurricane storm surge along the northern Gulf of Mexico, *Journal of Geophysical Research: Oceans*, 121, DOI: 10.1002/2015JC011400.

708 Bouws, E., Komen, G.J., 1983. On the balance between growth and dissipation in an extreme, depth-limited wind-sea in the southern North Sea. *Journal of Physical Oceanography*, 13, 1653-1658.

711 Breilh, J.-F., 2014. *Les surcotes et les submersions marines dans la partie centrale du Golfe de Gascogne : les enseignements de la tempête Xynthia*, Ph.D. Thesis, Université de la Rochelle, 225 pp.
714 Breilh, J.-F., Bertin, X., Chaumillon, E., Giloy, N., Sauzeau, T., 2014. How frequent is storm-
715 induced flooding in the central part of the Bay of Biscay? *Global and Planetary Change*,
716 122: 161-175. DOI: 10.1016/j.gloplacha.2014.08.013

717 Canas, A., Santos, A., Leitão, P. (2009). Effect of large-scale atmospheric pressure changes on
718 water level in the Tagus estuary, *Journal of Coastal Research*, Special Issue 56, 1627-
719 1631.

720 Carrère, L., Lyard, F., Cancet, M., Guillot, A., Roblou, L., 2012. FES2012: A new global tidal
721 model taking advantage of nearly 20 years of altimetry, *Proceedings of meeting "20
722 Years of Altimetry"*.

723 Chen, F., Dudhia, J., 2001. Coupling an advanced land surface–hydrology model with the Penn
724 State–NCAR MM5 modeling system. Part I: model implementation and sensitivity. *Monthly
725 Weather Review*, 129: 569–585. DOI: 10.1175/1520-
726 0493(2001)129<0569:CAALSH>2.0.CO;2

727 Chen, J., Chen, J., Liao, A., Cao, X., Chen, L., Chen, X., He, C., Peng, S., Lu, M., Zhang, W.,
728 Tong, X., Mills, J., 2015. Global land cover mapping at 30 m resolution: A POK-based
729 operational approach, *ISPRS Journal of Photogrammetry and Remote Sensing*: DOI:
730 10.1016/j.isprsjprs.2014.09.002.

731 Companhia das Lezírias, 1941. Relatório do exercício de 1940, presente à assembleia-geral de
732 accionistas de 1 de Abril de 1941. Lisboa, 1941. AHCL. Torre do Tombo (in Portuguese).

733 Compo, G.P., Whitaker, J.S., Sardeshmukh, P.D., Matsui, N., Allan, R.J., Yin, X., Gleason, B.E.,
734 Vose, R.S., Rutledge, G., Bessemoulin, P., Brönnimann, S., Brunet, M., Crouthamel, R.I.,
735 Grant, A.N., Groisman, P.Y., Jones, P.D., Kruk, M.C., Kruger, A.C., Marshall, G.J.,
736 Maugeri, M., Mok, H.Y., Nordli, Ø., Ross, T.F., Trigo, R.M., Wang, X.L., Woodruff, S.D.,
737 Worley, S.J., 2011. The twentieth century reanalysis project, *Quarterly Journal of the
738 Royal Meteorological Society*, 137: 1–28. DOI:10.1002/qj.77

739 Dee, D.P., Uppala, S.M., Simmons, A.J., Berrisford, P., Poli, P., Kobayashi, S., Andrae, U.,
740 Balmaseda, M.A., Balsamo, G., Bauer, P., Bechtold, P., Beljaars, A.C.M., van de Berg, L.,
741 Bidlot, J., Bormann, N., Delsol, C., Dragani, R., Fuentes, M., Geer, A.J., Haimberger, L.,
742 Healy, S.B., Hersbach, H., Hólm, E.V., Isaksen, L., Kållberg, P., Köhler, M., Matricardi, M.,
743 McNally, A.P., Monge-Sanz, B.M., Morcrette, J.-J., Park, B.-K., Peubey, C., de Rosnay, P.,
744 Tavolato, C., Thépaut, J.-N., Vitart, F., 2011. The ERA-Interim reanalysis: configuration
745 and performance of the data assimilation system. *Quarterly Journal of the Royal
746 Meteorological Society*, 137, 553–597. DOI: 10.1002/qj.828.

747 Dietrich, J.C., Bunya, S., Westerink, J.J., Ebersole, B.A., Smith, J.M., Atkinson, J.H., Jensen, R.,
748 Resio, D.T., Luettich, R.A., Dawson, C., Cardone, V.J., Cox, A.T., Powell, M.D., Westerink,
749 H.J., Roberts, H.J., 2010. A High-Resolution Coupled Riverine Flow, Tide, Wind, Wind
750 Wave and Storm Surge Model for Southern Louisiana and Mississippi: Part II – Synoptic
751 Description and Analysis of Hurricanes Katrina and Rita. *Monthly Weather Review*, 138:
752 378-404. DOI: 10.1175/2009MWR2907.1

753 Dodet, G., Bertin, X., Bruneau, N., Fortunato, A.B., Nahon, A., Roland, A., 2013. Wave-current
754 interactions in a wave-dominated tidal inlet, *Journal of Geophysical Research:
755 Oceans* 118, 3: 1587 - 1605. DOI: 10.1002/jgrc.20146.
756 Dudhia, J., 1989. Numerical study of convection observed during the winter monsoon experiment using a mesoscale two-dimensional model. *J. Atmos. Sci.*, 46, 3077-3107.
757
758 DOI: 10.1175/1520-0469(1989)046<3077:NSOCOD>2.0.CO;2

759 Ferreira, J.A., Carvalho, A.C., Carvalheiro, L., Rocha, A., Castanheira, J.M. 2014. On the influence of physical parameterizations and domains configuration in the simulation of an extreme precipitation event. *Dynamics of Atmospheres and Oceans*, 68: 35-55, ISSN 0377-0265. doi: 10.1016/j.dynatmoce.2014.08.001

763 Ferreira J.A., Liberato M.L.R., Ramos A.M., 2016. On the relationship between atmospheric water vapour transport and extra-tropical cyclones development. *Physics and Chemistry of the Earth*, Parts A/B/C doi: 10.1016/j.pce.2016.01.001

766 Fortunato, A.B., Baptista, A.M., Luettich, Jr., R.A., 1997. A three-dimensional model of tidal currents in the mouth of the Tagus Estuary, *Continental Shelf Research*, 17/14: 1689-1714. DOI: 10.1016/S0278-4343(97)00047-2

769 Fortunato, A.B., Oliveira, A., Baptista, A.M., 1999. On the effect of tidal flats on the hydrodynamics of the Tagus estuary. *Oceanologica Acta*, 22/1: 31-44.
770
771 DOI:10.1016/s0399-1784(99)80030-9.

772 Fortunato, A.B., Bruneau, N., Azevedo, A., Araújo, M.A.V.C., Oliveira, A., 2011. Automatic improvement of unstructured grids for coastal simulations, *Journal of Coastal Research*, Special Issue 164: 1028-1032.

775 Fortunato, A.B., Rodrigues, M., Dias, J.M., Lopes, C., Oliveira, A., 2013. Generating inundation maps for a coastal lagoon: a case study in the Ria de Aveiro (Portugal), *Ocean Engineering*, 64: 60-71. DOI: 10.1016/j.oceaneng.2013.02.020

778 Fortunato, A.B., Li, K., Bertin, X., Rodrigues, M., Miguez, B.M., 2016. Determination of extreme sea levels along the Iberian Atlantic coast, *Ocean Engineering*, 111/1: 471-482. DOI: 10.1016/j.oceaneng.2015.11.031

781 Fortunato, A.B., Oliveira, A., Rogeiro, J., Tavares da Costa, R., Gomes, J.L., Li, K., Jesus, G., Freire, P., Rilo, A., Mendes, A., Rodrigues, M., Azevedo, A. (2017). Operational forecast framework applied to extreme sea levels at regional and local scales, *Journal of Operational Oceanography*. DOI: 10.1080/1755876X.2016.1255471

785 Freitas, J.G., Dias J.A., 2013. 1941 windstorm effects on the Portuguese Coast. What lessons for the future? *Journal of Coastal Research*, Special Issue 65: 714-719. DOI: 10.2112/S165-121.

788 Freire, P., Taborda, R. Silva, A., 2007. Sedimentary characterization of Tagus estuarine beaches (Portugal). A contribution to the sediment budget assessment. *Journal of Soils and Sediments*, 7/5: 296-302. DOI:10.1065/jss2007.08.243.

791 Freire, P., Ferreira, Ó., Taborda, R., Oliveira, F.S.B.F., Carrasco, A.R., Silva, A., Vargas, C., Capitão, R., Fortes, C.J., Coli, A.B., Santos, J.A., 2009. Morphodynamics of fetch-limited beaches in contrasting environments, *Journal of Coastal Research*, Special Issue 56, 183-187.
795 Freire, P., Tavares, A., Sá, L., Oliveira, A., Fortunato, A.B., Santos, P.P., Rilo, A., Gomes, J.L., 796 Rogeiro, J., Pablo, R., Pinto, P.J., 2016. A local scale approach to estuarine flood risk 797 management, *Natural Hazards*, doi:10.1007/s11069-016-2510-y

798 Ganho, N., 2013. Risco de ventos tempestuosos de escala sinóptica em Portugal continental: 799 análise causal, in *Riscos naturais, antrópicos e mistos. Homenagem ao Professor Doutor 800 Fernando Rebelo*, Universidade de Coimbra, pp. 251-266 (in Portuguese).

801 Grell, G.A., Devenyi, D., 2002. A generalized approach to parameterizing convection combining 802 ensemble and data assimilation techniques. *Geophysical Research Letters*, 29 (14). DOI: 803 10.1029/2002GL015311

804 Guérin, T., Bertin, X., Dodet, G., 2016. A numerical scheme for coastal morphodynamic 805 modelling on unstructured grids, *Ocean Modelling*, 104, 45-53. DOI: 806 10.1016/j.ocemod.2016.04.009

807 Guerreiro, M., Fortunato, A.B., Freire, P., Rilo, A., Taborda, R., Freitas, M.C., Andrade, C., Silva, 808 T., Rodrigues, M., Bertin, X., Azevedo, A. (2015). Evolution of the hydrodynamics of the 809 Tagus estuary (Portugal) in the 21st century. *Revista de Gestão Costeira Integrada*, 15/1: 810 65-80. DOI: 10.5894/rgci515

811 Hasselmann, K., Barnett, T.P., Bouws, E., Carlson, H., Cartwright, D.E., Enke, K., Ewing, J.A., 812 Gienapp, H., Hasselmann, D.E., Kruseman, P., Meerburg, A., Muller, P., Olbers, D.J., 813 Richter, K., Sell, W., Walden, H., 1973. *Measurements of wind-wave growth and swell 814 decay during the Joint North Sea Wave Project (JONSWAP)*, Dtsch. Hydrogr. Z. Suppl., 815 12, A8, 1-95.

816 Hong, S.-Y., Lim, J.-O.J., 2006. The WRF Single-Moment 6-Class Microphysics Scheme 817 (WSM6). *Journal of the Korean Meteorological Society*, 42: 129-151.

818 IPCC, 2013. Climate change 2013: the physical science basis. Contribution of Working Group I 819 to the fifth assessment report of the Intergovernmental Panel on Climate Change. 820 Cambridge University Press, Cambridge, U.K. and New York, N.Y., USA, 1132 pp.

821 Karremann, M.K., Liberato, M.L.R., Ordóñez, P., Pinto, J.G., 2016 Characterization of synoptic 822 conditions and cyclones associated with top ranking potential wind loss events over Iberia. 823 *Atmospheric Science Letters*, 17: 354-361. DOI: 10.1002/asl.665

824 Goutx, D., Baraer, F., Roche, A., Jan, G., 2014. Ces tempêtes extrêmes que l'histoire ne nous a 825 pas encore dévoilées. *La houille Blanche*, 2, 27-33 (in French). DOI 10.1051/lhb/2014013.

826 Le Roy, S., Pedreros, R., Andre, C., Paris, F., Lecacheux, S., Marche, F., Vinchon, C., 2015. 827 Coastal flooding of urban areas by overtopping: dynamic modelling application to the 828 Johanna storm (2008) in Gavres (France), *Natural Hazards and Earth System Sciences*, 829 15/11: 2497-2510. DOI: 10.5194/nhess-15-2497-2015

830 Liberato, M.L.R., Pinto, J.G., Trigo, R.M., Ludwig, P., Ordóñez, P., Yuen, D., Trigo, I.F., 2013. 831 Explosive development of winter storm Xynthia over the subtropical North Atlantic Ocean. 832 *Natural Hazards and Earth System Science*, 13, 2239-2251. DOI: 10.5194/nhess-13- 833 2239-2013.

834 Lin, N., Emanuel, K., 2016. Grey swan tropical cyclones. *Nature Climate Change*, 6: 106-111. 835 DOI: 10.1038/NCLIMATE2777
836 Lopes, C.L.B., 2016. Flood risk assessment in Ria de Aveiro under present and future scenarios,
837 Ph.D. thesis, University of Aveiro, 224pp (in Portuguese).

838 Lopes, C.L., Silva, P.A., Dias, J.M., Rocha, A., Picado, A., Plecha, S., Fortunato, A.B., 2011.
839 Local sea level change scenarios for the end of the 21st century and potential physical
840 impacts in the lower Ria de Aveiro (Portugal), *Continental Shelf Research* 31, 14: 1515 -
841 1526. DOI: 10.1016/j.csr.2011.06.015

842 Lowe, J.A., Gregory, J.M., 2005. The effects of climate change on storm surges around the
843 United Kingdom, *Philosophical Transactions of the Royal Society A*, 363: 1313–1328. DOI:
844 10.1098/rsta.2005.1570.

845 Madaleno, I.M., 2006. Companhia das Lezírias – o passado e o presente. Hispania Nova.
846 Revista de História Contemporânea, 6, http://hispanianova.rediris.es. (in Portuguese).

847 Malhadas, M.S., Leitão, P.C., Silva, A., Neves, R., 2009. Effect of coastal waves on sea level in
848 Óbidos Lagoon, Portugal, *Continental Shelf Research*, 19/9: 1240–1250, DOI:
849 10.1016/j.csr.2009.02.007.

850 Marcos, M., Jorda, G., Gomis, D., Perez, B., 2011. Changes in storm surges in southern Europe
851 from a regional model under climate change scenarios, *Global and Planetary Change*,
852 77/3-4: 116-128. DOI: 10.1016/j.gloplacha.2011.04.002

853 Melet, A., Almar, R., Meyssignac, B., 2016. What dominates sea level at the coast: a case study
854 for the Gulf of Guinea, *Ocean Dynamics*, 66/5: 623–636. DOI 10.1007/s10236-016-0942-2

855 Mlawer, E.J., Taubman, S.J., Brown, P.D., Iacono, M.J., Clough, S.A., 1997. Radiative transfer
856 for inhomogeneous atmosphere: RRTM, a validated correlated-k model for the long-wave.
857 *Journal of Geophysical Research*, 102 (D14), 16663-16682.

858 Muir-Wood, R., 2011. The 1941 February 15th Windstorm in the Iberian Peninsula. *Trébol*, 56: 4-
859 13.

860 Neves, F.S., 2010. *Dynamics and hydrology of the Tagus estuary: results from in situ
861 observations*. 210p, Ph.D. thesis, Universidade de Lisboa, Portugal. Unpublished.
862 Available at: http://repositorio.ul.pt/handle/10451/2003

863 Noh, Y., Cheon, W.G., Hong, S.-Y., Raasch, S., 2003. Improvement of the K-profile model for
864 the planetary boundary layer based on large eddy simulation data. *Bound-Layer
865 Meteorology*, 107, 401-427.

866 Nunes, A., Pinho, J., Ganho, N., 2012. O Ciclone de Fevereiro de 1941: análise histórico-
867 geográfica dos seus efeitos no município de Coimbra, *Cadernos de Geografia*, 30/31: 53-
868 60 (in Portuguese).

869 Olabarrieta, M., Warner, J. C., Kumar, N., 2011. Wave-current interaction in Willapa Bay,
870 *Journal of Geophysical Research*, 116, C12014, DOI: 10.1029/2011JC007387.

871 Oliveira, F.S.B.F., 2000. Numerical simulation of wave propagation in the entrance of the Tagus
872 estuary, M.L. Spaulding, H.L. Butler (Eds.), Proceedings of the Sixth Conference on
873 Estuarine and Coastal Modeling, ASCE, pp. 510–525.

874 Perini, L., Calabrese, L., Salerno, G., Ciavola, P., Armaroli, C., 2016. Evaluation of coastal
875 vulnerability to flooding: comparison of two different methodologies adopted by the Emilia-
876 Romagna region (Italy), *Natural Hazards and Earth System Science*, 16: 181-194. DOI:
877 10.5194/nhess-16-181-2016.

878 Pinto, L., Fortunato, A.B., Zhang, Y., Oliveira, A., Sancho, F.E.P., 2012. Development and
879 validation of a three-dimensional morphodynamic modelling system for non-cohesive
880 sediments, *Ocean Modelling*, 57-58: 1-14. DOI: 10.1016/j.ocemod.2012.08.005

881 Rilo, A., Freire, P., Ceia, R., Mendes, R.N., Catalão, J., Taborda, R., 2012. Human effects on
882 estuarine shoreline decadal evolution. *Geophysical Research Abstracts*, Vol. 14, EGU2012-
883 10863, 2012, EGU General Assembly 2012.

884 Rilo, A., Freire, P., Mendes, R.N., Ceia, R., Catalão, J., Taborda, R., Melo, R., Caçador, M.I.,
885 Freitas, M.C., Fortunato, A.B., Alves, E., 2014. Methodological framework for the definition
886 and demarcation of the highest astronomical tide line in estuaries: the case of Tagus
887 Estuary (Portugal), *Revista de Gestão Costeira Integrada*, 14/1: 95-107 (in Portuguese).
888 DOI: 10.5894/rgci450.

889 Rodrigues, M., Oliveira, A., Queiroga, H., Fortunato, A.B., Zhang, Y.J., 2009. Three-dimensional
890 modeling of the lower trophic levels in the Ria de Aveiro (Portugal), *Ecological
891 Modelling* 220, 9-10: 1274–1290. DOI: 10.1016/j.ecolmodel.2009.02.002

892 Rodrigues, M., Oliveira, A., Guerreiro, M., Fortunato, A.B; Menaia, J., David, L.M.; Cravo, A.
893 2011. Modeling fecal contamination in the Aljezur coastal stream (Portugal), *Ocean
894 Dynamics* 61/6: 841-856. DOI: 10.1007/s10236-011-0392-9

895 Rodrigues, M., Fortunato, A.B., Freire, P., 2016. Salinity evolution in the Tagus estuary relative
896 to climate change, *Actas das 4as Jornadas de Engenharia Hidrográfica, Instituto
897 Hidrográfica*, Lisbon, 179-182.

898 Roland, A., Zhang, Y.J., Wang, H.V., Meng, Y., Teng, Y.-C., Maderich, V., Brovchenko, I.,
899 Dutour-Sikiric, M., Zanke, U. (2012). A fully coupled 3D wave-current interaction model on
900 unstructured grids. *Journal of Geophysical Research*, 117, C00J33. DOI:
901 10.1029/2012JC007952.

902 Rusu, L., Bernardino, M., Guedes Soares, C., 2009. Influence of wind resolution on the
903 prediction of waves generated in an estuary, *Journal of Coastal Research*, Special Issue
904 56: 1419-1423.

905 Rusu, L., Bernardino, M., Guedes Soares, C., 2011. Modelling the influence of currents on wave
906 propagation at the entrance of the Tagus estuary, *Ocean Engineering*, 38: 1174-1183.
907 DOI: 10.1016/j.oceaneng.2011.05.016

908 Salgueiro, A.R., Machado, M.J., Barriendos, M., Pereira, H.G., Benito, G., 2013. Flood
909 magnitudes in the Tagus River (Iberian Peninsula) and its stochastic relationship with daily
910 North Atlantic Oscillation since mid-19th Century, *Journal of Hydrology*, 502: 191–201. DOI:
911 10.1016/j.jhydrol.2013.08.008

912 Skamarok, W.C., Klemp, J.B., Dudhia, J., Gill, D.O., Barker, D.M., Duda, M.G., Huang, X.Y.,
913 Wang, W., Powers, J.G., 2008. *A description of the advanced research WRF Version 3.*
914 NCAR Technical Note.
915 Schloen, J., Stanev, E.V., Grashorn, S., 2017. Wave-current interactions in the southern North
916 Sea: The impact on salinity, *Ocean Modelling*, 111: 19–37. DOI:
917 10.1016/j.ocemod.2017.01.003

918 Stucki, P., Brönnimann, S., Martius, O., Welker, C., Rickli, R., Dierer, S., Bresch, D.N., Compo,
919 G.P., Sardeshmukh, P.D., 2015. Dynamical downscaling and loss modeling for the
920 reconstruction of historical weather extremes and their impacts: a severe Foehn storm in
921 1925. *Bulletin of the American Meteorological Society*. 96: 1233–1241. DOI:
922 http://dx.doi.org/10.1175/BAMS-D-14-00041.1

923 Tavares, A.O., Santos, P.P., Freire, P., Fortunato, A.B., Rilo, A., Sá, L., 2015. Flooding hazard in
924 the Tagus estuarine area: the challenge of scale in vulnerability assessments.
925 *Environmental Science & Policy*, 51: 238-255. DOI: 10.1016/j.envsci.2015.04.010.

926 Tolman, H.L., 2009. User manual and system documentation of WAVEWATCH III, version 3.14.
927 NOAA/NWS/NCEP/MMAB Technical Note 276, 194 p.

928 Tomás, A., Méndez, F.J., Medina, R., Jaime, F.F., Higuera, P., Lara, J.L., Ortiz, M.D., Álvarez de
929 Eulate, M.F., 2015. A methodology to estimate wave-induced coastal flooding hazard
930 maps in Spain, *Journal of Flood Risk Management*, DOI: 10.1111/jfr3.12198.

931 Turner, P., Baptista, A.M., 1993. ACE/gredit User's Manual. Software for Semi-automatic
932 Generation of Two-Dimensional Finite Element Grids. Center for Coastal and Land-Margin
933 Research, Oregon Graduate Institute of Science & Technology.

934 Ulbrich, U., Leckebusch, G.C., Grieger, J., Schuster, M., Akperov, M., Bardin, M.Y., Feng, Y.,
935 Gulev, S., Inatsu, M., Keay, K., Kew, S.F., Liberato, M.L.R., Lionello, P., Mokhov, I.I., Neu,
936 U., Pinto, J.G., Raible, C.C., Reale, M., Rudeva, I., Simmonds, I., Tilinina, N.D., Trigo, I.F.,
937 Ulbrich, S., Wang, X.L., Wernli, H., and The IMILAST team, 2013. Are greenhouse gas
938 signals of northern hemisphere winter extra-tropical cyclone activity dependent on the
939 identification and tracking algorithm? *Meteorologische Zeitschrift*, 22/1: 61-68. DOI:
940 10.1127/0941-2948/2013/0420.

941 Vargas, C., Oliveira, F.S.B.F., Oliveira, A., Charneca, N., 2008. Flood vulnerability analysis of an
942 estuarine beach: application to Alfeite sand spit (Tagus Estuary). *Revista de Gestão
943 Costeira Integrada / Integrated Coastal Zone Management Journal*, 8/1: 25-43. DOI:
944 10.5894/rgci26

945 Vousdoukas, M.I., Voukouvalas, E., Annunziato, A., Giardino, A., Feyen, L., 2016. Projections of
946 extreme storm surge levels along Europe, *Climate Dynamics*, DOI: 10.1007/s00382-016-019-5.

948 Xie, L., Liu, H., Peng, M., 2008. The effect of wave-current interactions on the storm surge and
949 inundation in Charleston Harbor during Hurricane Hugo 1989, *Ocean Modelling*, 20/3: 252-
950 269, DOI: http://dx.doi.org/10.1016/j.ocemod.2007.10.001.

951 Wargula, A., Raubenheimer, B., Elgar, S., 2014, Wave-driven along-channel subtidal flows in a
952 well-mixed ocean inlet, *Journal of Geophysical Research: Oceans*, 19: 2987–3001, DOI:
953 10.1002/2014JC009839.

954 Zampato, L., Bajo, M., Canestrelli, P., Umgiesser, G. (2016). Storm surge modelling in Venice:
955 two years of operational results, *Journal of Operational Oceanography*, 9: S46-S57. DOI:
956 10.1080/1755876X.2015.1118804
957 Zhang, Y., Baptista, A.M., 2008. SELFE: A semi-implicit Eulerian-Lagrangian finite-element
958 model for cross-scale ocean circulation. *Ocean Modeling*, 21(3-4), 71-96. DOI:
959 10.1016/j.ocemod.2007.11.005.

960 Zhang, Y.J., Ye, F., Stanev, E.V., Grashorn, S., 2016. Seamless cross-scale modeling with
961 SCHISM, *Ocean Modelling*, 102: 64-81. DOI:10.1016/J.OCEMOD.2016.05.002.

962 Zorndt, A., Goseberg, N., Schlurmann, T., 2014. Influence of retention areas on the propagation
963 of storm surges in the Weser estuary. *Coastal Engineering Proceedings*, 1(34),
964 management.51. DOI: 10.9753/icce.v34.management.51

965
## Tables

<table><caption>Table 1. Characteristics of the two storms. The tidal range is computed from the results of the tidal model of Fortunato et al. (2016) at Cascais.</caption><thead><tr><th>Storm</th><th>Tidal range (m)</th><th>Minimum atmospheric pressure (hPa)</th><th>Maximum wind speed (m/s)</th><th>River flow</th></tr></thead><tbody><tr><td>1941</td><td>3.0</td><td>952.1 (Geophysical Institute)<br>957.8 (Monte Estoril)</td><td>18 (Geophysical Institute)<br>19 (Monte Estoril)</td><td>There was a major river flood in January 21-23. River flow in the Tagus river peaked at 3266 m³/s at Vila Velha de Rodão.</td></tr><tr><td>2010</td><td>2.9</td><td>976.2 (Lisbon)<br>986 (Cascais)</td><td>14 (Lisbon)</td><td>River flow in the Tagus river peaked at 3917 m³/s on February 27 at Almourol</td></tr></tbody></table>

<table><caption>Table 2. Errors in the Tagus for the 1972 validation run (in cm).</caption><thead><tr><th>Station</th><th>RMSE</th><th>U-RMSE</th><th>HT-RMSE</th></tr></thead><tbody><tr><td>Cascais</td><td>4.6</td><td>4.6</td><td>5.7</td></tr><tr><td>Paço de Arcos</td><td>6.4</td><td>3.8</td><td>4.6</td></tr><tr><td>Trafaria</td><td>9.7</td><td>5.4</td><td>7.7</td></tr><tr><td>Lisbon</td><td>11.0</td><td>6.1</td><td>13.5</td></tr><tr><td>Pedrouços</td><td>6.3</td><td>4.0</td><td>7.4</td></tr><tr><td>Cacilhas</td><td>5.8</td><td>5.3</td><td>6.0</td></tr><tr><td>Seixal</td><td>8.9</td><td>6.4</td><td>8.9</td></tr><tr><td>Montijo</td><td>12.4</td><td>6.6</td><td>13.8</td></tr><tr><td>Cabo Ruivo</td><td>14.5</td><td>8.6</td><td>15.4</td></tr><tr><td>Alcochete</td><td>11.7</td><td>10.4</td><td>9.2</td></tr><tr><td>Santa Iria</td><td>16.0</td><td>12.8</td><td>15.0</td></tr><tr><td>Ponta da Erva</td><td>16.8</td><td>16.7</td><td>10.5</td></tr><tr><td>Vila Franca</td><td>18.5</td><td>9.9</td><td>19.5</td></tr></tbody></table>
<table><caption>Table 3. Model simulations</caption><thead><tr><th rowspan="2">Run</th><th colspan="4">Open boundaries</th><th colspan="2">Surface boundary</th><th>Processes</th><th rowspan="2">Other</th></tr><tr><th>Tide</th><th>Surge</th><th>Waves</th><th>River flow</th><th>Wind</th><th>Atmospheric pressure gradient</th><th>Wave generation</th></tr></thead><tbody><tr><td>1</td><td>1941</td><td>Yes</td><td>Yes</td><td>1941</td><td>Yes</td><td>Yes</td><td>Yes</td><td>—</td></tr><tr><td>2</td><td>1941</td><td>Yes</td><td>Yes</td><td>1941</td><td>Yes</td><td>Yes</td><td>Yes</td><td>Fine grid</td></tr><tr><td>3</td><td>1941</td><td>Yes</td><td>Yes</td><td>2010</td><td>Yes</td><td>Yes</td><td>Yes</td><td>—</td></tr><tr><td>4</td><td>1941</td><td>Yes</td><td>Yes</td><td>1941</td><td>Yes</td><td>Yes</td><td>Yes</td><td>Reduced friction in WWM</td></tr><tr><td>5</td><td>1941</td><td>Yes</td><td>Yes</td><td>1941</td><td>Yes</td><td>Yes</td><td>Yes</td><td>Reduced breaking coefficient in WWM</td></tr><tr><td>6</td><td>2006</td><td>Yes</td><td>Yes</td><td>1941</td><td>Yes</td><td>Yes</td><td>Yes</td><td>—</td></tr><tr><td>7</td><td>1941</td><td>Yes</td><td>Yes</td><td>1941</td><td>Yes</td><td>Yes</td><td>Yes</td><td>SLR</td></tr><tr><td>8</td><td>1941</td><td>No</td><td>No</td><td>Climatology</td><td>No</td><td>No</td><td>No</td><td>—</td></tr><tr><td>9</td><td>1941</td><td>No</td><td>Yes</td><td>1941</td><td>Yes</td><td>Yes</td><td>Yes</td><td>—</td></tr><tr><td>10</td><td>1941</td><td>Yes</td><td>No</td><td>1941</td><td>Yes</td><td>Yes</td><td>Yes</td><td>—</td></tr><tr><td>11</td><td>1941</td><td>Yes</td><td>Yes</td><td>Climatology</td><td>Yes</td><td>Yes</td><td>Yes</td><td>—</td></tr><tr><td>12</td><td>1941</td><td>Yes</td><td>Yes</td><td>1941</td><td>No</td><td>Yes</td><td>Yes</td><td>—</td></tr><tr><td>13</td><td>1941</td><td>Yes</td><td>Yes</td><td>1941</td><td>Yes</td><td>No</td><td>Yes</td><td>—</td></tr><tr><td>14</td><td>1941</td><td>Yes</td><td>Yes</td><td>1941</td><td>Yes</td><td>Yes</td><td>No</td><td>—</td></tr><tr><td>15</td><td>M2</td><td>No</td><td>Constant</td><td>Constant</td><td>No</td><td>No</td><td>No</td><td>—</td></tr><tr><td>16</td><td>M2</td><td>No</td><td>No</td><td>Constant</td><td>No</td><td>No</td><td>No</td><td>—</td></tr><tr><td>17</td><td>1941</td><td>Yes</td><td>Yes</td><td>1941</td><td>Yes</td><td>Yes</td><td>Yes</td><td>Dikes raised by 20 cm</td></tr></tbody></table>
<table><caption>Table 4. Maximum elevations for the sensitivity analysis simulations (m, CD). For each of the runs 2-5 the values representing the maximum discrepancy relative to Run 1 are highlighted in bold.</caption><thead><tr><th>Station</th><th>Run 1<br>(standard)</th><th>Run 2<br>(fine grid)</th><th>Run 3<br>(river flow)</th><th>Run 4<br>(reduced<br>friction)</th><th>Run 5<br>(reduced<br>breaking coef)</th></tr></thead><tbody><tr><td>Carcavelos</td><td>4.787</td><td>4.786</td><td>4.787</td><td><strong>4.798</strong></td><td>4.787</td></tr><tr><td>São João</td><td>4.780</td><td>4.779</td><td>4.781</td><td><strong>4.791</strong></td><td>4.780</td></tr><tr><td>Paço de Arcos</td><td>4.506</td><td>4.502</td><td>4.505</td><td>4.506</td><td><strong>4.513</strong></td></tr><tr><td>Trafaria</td><td>4.597</td><td>4.598</td><td>4.604</td><td>4.599</td><td>4.597</td></tr><tr><td>Lisbon</td><td>4.592</td><td>4.591</td><td>4.599</td><td>4.593</td><td>4.592</td></tr><tr><td>Pedrouços</td><td>4.721</td><td>4.739</td><td>4.726</td><td>4.723</td><td>4.721</td></tr><tr><td>Cacilhas</td><td>4.735</td><td>4.755</td><td>4.742</td><td>4.736</td><td>4.735</td></tr><tr><td>Seixal</td><td>4.777</td><td>4.803</td><td>4.782</td><td>4.779</td><td>4.776</td></tr><tr><td>Montijo</td><td>4.832</td><td>4.854</td><td>4.830</td><td>4.833</td><td>4.831</td></tr><tr><td>Cabo Ruivo</td><td>4.873</td><td>4.900</td><td>4.876</td><td>4.875</td><td>4.874</td></tr><tr><td>Alcochete</td><td>5.013</td><td>5.035</td><td>5.014</td><td>5.015</td><td>5.013</td></tr><tr><td>Santa Iria</td><td>5.100</td><td>5.122</td><td>5.099</td><td>5.103</td><td>5.100</td></tr><tr><td>Ponta da Erva</td><td>5.150</td><td>5.121</td><td>5.158</td><td>5.155</td><td>5.152</td></tr><tr><td>Vila Franca</td><td>5.377</td><td>5.387</td><td>5.368</td><td>5.538</td><td>5.376</td></tr><tr><td>Conchoso</td><td>5.588</td><td><strong>5.640</strong></td><td><strong>5.685</strong></td><td>5.589</td><td>5.590</td></tr></tbody></table>

<table><caption>Table 5. Inundated areas (km²).</caption><thead><tr><th>Zone</th><th>Run 1</th><th>Run 6</th><th>Run 7</th><th>Xynthia (2010)</th><th>Run 17</th></tr></thead><tbody><tr><td>Lezíria</td><td>1.6</td><td>2.5</td><td>25.2</td><td>1.1</td><td>1.3</td></tr><tr><td>Póvoa island</td><td>4.1</td><td>4.2</td><td>4.2</td><td>4.1</td><td>4.1</td></tr><tr><td>Lombo island</td><td>0.8</td><td>4.9</td><td>5.1</td><td>0.1</td><td>0.1</td></tr><tr><td>Alhandra island</td><td>2.4</td><td>2.5</td><td>2.5</td><td>0.1</td><td>1.4</td></tr><tr><td>Sorraia left margin</td><td>17.9</td><td>27.7</td><td>40.1</td><td>10.1</td><td>14.3</td></tr></tbody></table>
979 **Figure captions**

980

981 **Figure 1.** General location: North Atlantic, showing the boundaries of the atmospheric model regional and local grids in black, of the wave model grid in red, of the tide and tide/surge model grid in green. The yellow lines represent the trajectory of the center of the 1941 (solid line) and 2010 (dashed line - Xynthia) storms. The trajectories were determined based on the atmospheric model results described below. The inset shows the boundaries of the Tagus estuary grid.

986 **Figure 2.** Tagus estuary bathymetry relative to chart datum (interpolated onto the model grid) and place names: a) global view; b) detail of the inundation-prone areas of the upper estuary. The maximum high water line of Rilo et al. (2014), where available, is represented in blue, and the dikes in black. The red lines in b) represent stretches of the dikes that were overtopped during the 2010 Xynthia storm. The yellow circles indicate atmospheric stations. Real and virtual tidal gauges are represented as red circles.

991 **Figure 3.** Determination of cross-sectional-averaged depth along the estuary axis upstream of Vila Franca. CD stands for Chart Datum.

993 **Figure 4.** Estimates of the river flow reaching the estuary.

994 **Figure 5.** Validation of the atmospheric model: comparison between the measured (circles) and modeled (lines) atmospheric pressures.

996 **Figure 6.** Validation of the significant wave heights obtained with the regional wave model at Leixões and Sines, respectively to the North and to the South of the Tagus estuary, for the period 2007-2012. RMSD: root mean square deviations; NRMSD: normalized root mean square deviations.

999 **Figure 7.** Model grid: a) model domain and grid resolution; b) detail of inundation-prone areas of the upper estuary.

1001 **Figure 8.** Manning friction coefficient ($m^{1/3}/s$).

1002 **Figure 9.** Validation of the modeled inundation of the Seixal downtown area during Xynthia. The red circles represent the inundation limits determined in the field and the blue shape is an interpretation of the inundated area based on those limits and the existing buildings (from Freire et al., 2016). Dark (light) green triangles represent fully (partially) wet elements as determined by the model. Fully dry elements are omitted for clarity.

1006 **Figure 10.** Return period of the 1941 storm effects along the Portuguese west coast. The vertical dashed lines indicate the value reached in the 1941 simulation at three different latitudes indicated by the color. a) Return periods of the significant wave height based on the 1952-2012 hindcasts of Bertin et al. (2013). The solid lines are a GEV distribution fitted to the empirical distribution (circles). b) Return periods for the maximum elevations at the Portuguese west coast based on the 1980-2010 hindcasts of Fortunato et al. (2016), and computed with the approach of Fortunato et al. (2013). These return periods do not take wave setup into account.

1013 **Figure 11.** Inundation of the upper Tagus estuary for the 1941 storm: a) present MSL; b) worst case scenario, obtained by combining the 1941 storm with an extreme spring tide; c) sea level rise of 0.5 m.

1015 **Figure 12.** Maximum surge at different stations throughout the estuary.

1016 **Figure 13.** Tide-surge interaction in the Tagus estuary: a) evolution of the wave-induced setup along the estuary; b) tidal elevations at Cascais. The dashed line in a) corresponds to a simulation in which the effect of velocities on wave propagation was switched off.
<figure><img src="image_2.png" /></figure>

20°N

39.2°N

39°N

38.8°N

38.6°N

38.4°N

9.2°W 8.8°W
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>
<figure><img src="image_5.png" /></figure>

b)
<figure><img src="image_6.png" /></figure>
<figure><img src="image_7.png" /></figure>
<figure><img src="image_8.png" /></figure>
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>
<figure><img src="image_11.png" /></figure>
<figure><img src="image_12.png" /></figure>
<figure><img src="image_13.png" /></figure>
<figure><img src="image_14.png" /></figure>

The blue line represents the current border of the city. Red dots indicate the locations of existing and proposed border markers.

N
60
0
120
meters
<figure><img src="image_15.png" /></figure>
<figure><img src="image_16.png" /></figure>
<figure><img src="image_17.png" /></figure>

a)
<figure><img src="image_18.png" /></figure>

b)
<figure><img src="image_19.png" /></figure>

water depth (m)

> 0.0 - 0.2

0.3 - 0.4

0.5 - 0.6

0.7 - 0.8

0.9 - 1.0

1.1 - 3.0

3.1 - 5.0

5.1 - 7.0

7.1 - 20.0
<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /></figure>
<figure><img src="image_22.png" /></figure>