Can Modeling the Geologic Record Contribute to Constraining the Tectonic Source of the 1755 CE Great Lisbon Earthquake?

Abstract

The precise location of the seismic source of 1755 CE Great Lisbon earthquake is still uncertain. The aim of this work is to use an onland sedimentary record in southern Portugal to test and validate seismic sources for the earthquake. To achieve this, tsunami deposit thicknesses from over 150 cores collected at Salgados in southern Portugal were compared to the results of a tsunami sediment transport model (Delft3D‐FLOW) that simulates tsunami propagation, inundation, erosion, and deposition. Five different hypothetical seismic sources were modeled with varying bed roughness coefficients to assess how well they reproduced observed patterns of tsunami deposit thicknesses and dune. Modeled and observed historical tsunami arrival times were also used to test different earthquake sources. Based on these comparisons, three modeled earthquake sources were able to reproduce the observed data, suggesting they should be regarded as somewhat more likely sources for the 1755 earthquake in contrast to four other modeled sources. The fault closest to shore (Marquês de Pombal) yielded the best correlations between model and observations.

Introduction

Many different sources have been proposed for the 1755 CE earthquake and tsunami, although to date no single source accounts for the massive energy‐release required to: (1) explain the spatial pattern of earthquake intensity observed along the Cadiz Gulf and both the western and southern mainland Portuguese coast, and (2) agree with tsunami travel times observed around and over the Atlantic Ocean. Some studies suggest that this event was triggered by interconnected faults or landslide movements (e.g., Grácia et al., 2003; Vilanova & Fonseca, 2004). Furthermore, the Cadiz Accretionary Wedge (CAW), Horseshoe Fault (HSF), Gorringe Bank (GB), and Marquês de Pombal Fault (MPF) have all been proposed as primary locations where fault‐rupture might have generated the 1755 CE earthquake (Barkan et al., 2009; Matias et al., 2013; Santos et al., 2009) (Figure 1).

1 Figure. Top image – Modeled earthquake sources locations, area of interest (AOI – Salgados), historical runup and simulated water level source locations. Bathymetry. Levels 0, 1, and 2 are the nested grid limits for numerical modeling. Lower image – Sediment cores (red points size according the sediment deposit thickness and X when no deposit founded), virtual tide station (black cross in white dot), and AB/AC profiles at Salgados beach.

Most of the studies that identify specific sources for the 1755 CE earthquake primarily utilize data derived from reports compiled on Arquivos do Ministério do Reino (1756), which contains information on the locations and times when ground shaking was felt, in addition to reports of damage (Santos & Koshimura, 2015a, 2015b). Other studies were based exclusively on simulations of tsunami travel times, either from proposed earthquake sources to the locations where observed data describes the time of arrival and impacts of tsunami waves (e.g., Baptista, Heitor, et al., 1998; Santos et al., 2009; Wronna et al., 2015), or by using target‐to‐source back‐ray tracing (Baptista, Miranda, Miranda, & Victor, et al., 1996). However, no single triggering mechanism proposed so far has been able to reproduce all of the tsunami travel times inferred from the historical records (Santos et al., 2009).

The generation mechanism of earthquake tsunamis is commonly computed by modeling the initial sea surface elevation through analytical solutions from elastic models that compute crustal deformation according to fault motions (e.g., Mansinha & Smylie, 1971; Okada, 1985). Alternatively, an idealized waveform can be directly introduced into the hydrodynamic model. The tsunami propagation and inundation is therefore simulated by hydrodynamic models mostly based on the principle of continuity and mass conservation. An extensive description of available hydrodynamic and sediment transport numerical models applied to tsunami simulations can be find on the work of Sugawara, Goto, and Jaffe (2014).

Previous works have applied forward modeling for relating modeled tsunami inundation and the spatial distribution of historical tsunami deposits (e.g., Butler et al., 2014; Namegaya & Satake, 2014) or contemporary events (e.g., Grilli et al., 2019). However, this approach is known to be frequently inaccurate when directly comparing these features. Therefore, coupling hydrodynamic and sediment transport modeling recently appears as an optimal approach for estimating paleotsunami hydrodynamic characteristics and constraining earthquake source and/or parameters (Sugawara, Yu, & Yen, 2019).

Our approach initially models tsunami propagation from proposed seismic source areas (initial boundary conditions – Figure 2 and Table 1) to selected coastal target locations. Second, travel times are derived and validated with the documentary data (Table 1). Finally, we model patterns of onshore inundation including inland sediment transport and effects on coastal morphology at Salgados lowland (Figure 1). This coastal lowland contains high‐resolution geological and geomorphological datasets that provide objective information on deposition and erosion induced by the 1755 CE Lisbon tsunami (Costa, Andrade, Dawson, et al., 2012; Costa, Andrade, Freitas, et al., 2012). This allows for rigorously testing a number of proposed earthquake and tsunami sources by expanding the number and diversity of metrics used as validation criteria.

2 Figure. Initial condition (sea level displacement) for the five different hypothetical fault (source) tested for the 1755 CE earthquake. Four sources represent uniform slip on faults: Cadiz Accretionary Wedge (CAW), Horseshoe Fault (HSF), Gorringe Bank (GB), and Marquês de Pombal Fault (MPF). Parameters for these sources were derived from previous studies (e.g., Ramalho et al., 2018). Scenario 1 is a rearrangement of the 1969 Lisbon earthquake source and a possible combination with a seismogenic structure combining GB and HSF.

1 TableFault Parameters for the Modeled Hypothetical Tectonic Sources (Values Obtained From Averaging Data Provided in Selected References (Baptista, Miranda, Chierici, & Zitellini, 2003; Baptista, Miranda, Omira, & Antunes, 2011; Barkan et al., 2009; Gjevik et al., 1997; Gutsher et al., 2006; Lima et al., 2010; Omira et al., 2009; Ramalho et al., 2018; Silva 2017; Wronna et al., 2015)

Source	CAW	HSF	GB	MPF	Scenario 1
Depth (km)	5	5	5	5	5
Length (km)	170	170	120	110	114/120
Width (km)	200	100	60	70	100/60
Slip (m)	20	14	10	12	14/10
Strike (°)	349	42	57	23	42/57
Dip (°)	6	35	35	35	35
Rake (°)	90	90	90	90	90
Magnitude	8.8	8.5	8.2	8.2	8.8
Sines historical TTT (30 min)	40	34	34	29	34
Sines historical runup (no information)	–	–	–	–	–
São Vicente historical TTT (16–17 min)	21	18	25	15	20
São Vicente historical runup (11–15 m)	13	19	8	12	17
Lagos historical TTT (23–30 min)	27	27	39	28	28
Lagos historical runup (10 m)	10	9	5	12	12
Huelva historical TTT (45–50 min)	40	50	>60	50	55
Huelva historical runup (no information)	–	–	–	–	–

Note. Comparison of tsunami travel time (TTT), in minutes, and runup, in meters, at coastal locations along the broad Gulf of Cadiz retrieved from the historical record and yielded by modeling different epicentral areas.

Geologic Evidence of the 1755 CE Tsunami in the Salgados Lowland

In this study we use geological data from over 150 cores obtained at the Salgados coastal lowland, Portugal (Figure 1). The lowland is a sediment‐filled lagoon separated from the ocean by a sandy beach backed by a multiple‐ridged dune. Landward of the dune, the 1755 CE tsunami deposit has been characterized as a massive to normally graded, sheet of marine‐facies shell‐rich sand with an erosive base sandwiched in lagoonal mud (Costa, Andrade, Dawson, et al., 2012; Costa, Andrade, Freitas, et al., 2012; Costa, Costas, et al., 2016). The tsunami deposit is roughly 50 cm thick closer to the sea and thins in both landward and alongshore directions (Costa, Costas, et al., 2016). Costa, Andrade, Dawson, et al. (2012), Costa, Andrade, Freitas, et al. (2012), Costa, Costas, et al. (2016), and Moreira et al. (2017) used paleoecological, geochemical, mineralogical, microtextural and grain size data from tsunami and modern surface sediments from Salgados lowlands to show that the primary source of the tsunami sediments were the dunes and secondarily the beach.

Costa, Costas, et al. (2016) present data from a ground‐penetrating radar (GPR) investigation of the dunes at Salgados. Two cross‐shore profiles (AB and AC in Figure 1) extending from the upper beach toward the backbarrier area provided information on the architecture of the dune complex, sediment packages and erosional features. Profile AB is 210 m long and located 300 m westward from Salgados inlet channel, where the dune crest reaches 8.5 m above mean sea level (MSL). Profile AC is 245 m long and is located 120 m westward of AB. Profile AB presents a dune crest that reaches 11 m above MSL. Both profiles contain a clear image of an erosional surface within the dunes at approximately 6 m above MSL. Optically Stimulated Luminescence (OSL) dating of dune sands immediately below and above that surface constrained an episode of erosion to the mid‐seventeenth century (Costa, Costas, et al., 2016). Regional tsunami historical records, however, suggest that wave heights at the coast were higher, up to 12 m above MSL (Costa, Costas, et al., 2016).

Tsunami Modeling Methods

To validate tsunami hydrodynamic and sediment transport models we tested five different hypothetical fault (source) areas for the 1755 CE earthquake (Figure 2). All source areas considered herein have been previously proposed in the literature and include the MPF (Baptista, Miranda, Chierici, & Zitellini, 2003; Lima et al., 2010; Omira et al., 2009; Wronna et al., 2015), GB (Baptista, Miranda, Chierici, & Zitellini, 2003; Omira et al., 2009; Wronna et al., 2015), CAW (Gutscher, 2006; Lima et al., 2010; Omira et al., 2009; Ramalho et al., 2018; Wronna et al., 2015), HSF (Baptista, Miranda, Chierici, & Zitellini, 2003, Baptista, Miranda, Omira, & Antunes, 2011; Barkan et al., 2009; Lima et al., 2010; Omira et al., 2009; Ramalho et al., 2018; Wronna et al., 2015) and one hypothetical composite scenario combining some of these shallow faults (Figures 1 and 2; Table 1).

The initial sea surface perturbation generated by the sources considered herein has been computed using Mansinha and Smiley (1971) elastic deformation approach through Mirone software (Luis, 2007) and on the OpenEarth MatLab script from Deltares (2016). The first four sources represent uniform slip on faults: CAW, HSF, GB, and MPF. Parameters for these sources were derived from previous studies validated against historically observed tsunami arrival times and backward ray tracing, but not against the geological record (e.g., Ramalho et al., 2018; Wronna et al., 2015). The next hypothetical source, Scenario 1, is a rearrangement of the 1969 Lisbon earthquake source (Fukao, 1973) by combining GB and HSF sources occurring simultaneously.

Tsunami propagation, inundation, and sediment transport were modeled using Delft3D‐FLOW, which solves the nonlinear shallow water equations using a finite difference scheme and has been validated against analytical, laboratory, and field measurements of tsunami hydrodynamics (Apotsos et al., 2011). Three nested grids were constructed with spatial resolutions of 232 m (Level 0), 100 m (Level 1), 50, 25, and 5 m (Level 2 – varying spatial resolutions on a single grid) (top image on Figure 1). Also a synthetic tide gauge was added 500 m offshore southward of Salgados near the 10 m isobath to monitor tsunami water levels (lower image on Figure 1).

A combined bathymetric‐topographic DEM was created from three different datasets (from European Marine Observation and Data Network, Copernicus Program, the Planning and Management of Coastal Zones Program of Portugal) with vertical datum adjusted to MSL at the Cascais tide gauge, 25 km west of Lisbon. The DEM was adjusted by using lithostratigraphic data from the 150 sediment cores retrieved from the lowland to reconstruct the approximate surface prior to the 1755 CE event. A final correction of −1.5 m was applied to the DEM to account for the tide level observed at the time of the earthquake.

A depth‐averaged (2DH) model was run using the weakly reflective Riemann boundaries on all grid levels in order to calculate tsunami‐induced hydrodynamics. Runups (i.e., height reached on the observation points) and tsunami travel times were compared to observations at four sites (Sines, Cabo de São Vicente, Lagos, and Huelva) (Figure 1 and Table 1). The tsunami sediment transport model uses the formulation from van Rijn (2007) with 10 vertical layers on the Level 2 grid (3D) in order to include the effects of suspended‐sediment induced density stratification on the vertical turbulent mixing. Furthermore, changes in the bed level caused by erosion and sedimentation processes are updated at each time step of the simulation (Lesser et al., 2004). An unlimited erodible sediment source is represented in the model as a 10–15 m thick sand extending from the offshore to the back of the foredune, with no sand available in the muddy lowland area. In all simulations, the median grain size sediment parameter [D50] used was 250 µm with a density 2,650 kg/m³, based roughly on the D50 observed in the observed tsunami deposits. In order to test the sensitivity of model outputs regarding bed roughness, we adjusted the Manning’s n roughness coefficient between 0.025 and 0.080 in the dune and lowland areas.

Model Results

Studies on tsunami deposit taphonomy (Spiske et al., 2019; Szczucinski, 2012) suggest that the extent and thickness of sandy tsunami deposits generally decrease with time. Therefore, the spatial extents and thicknesses of the observed 1755 tsunami deposit in Salgados are assumed to be minimum representations of the original deposit, making direct comparisons with the modeled tsunami deposit difficult.

Sediment transport simulations for GB source, using a conservative low bed roughness of 0.025, was not able to reproduce a tsunami deposit volume of more than 25% of the measured volume in Salgados (Figure 3). The modeled volumes from the CAW and Scenario 1 sources reach or exceed 100% of the measured volume from the cores over a range of bed roughness values ranging from 0.025 to 0.065. Using a Manning’s roughness coefficient of 0.080 with the CAW source, the simulated volume of sediment deposited was 121% of the volume calculated from the cores samples (Figure 3).

3 Figure. Modeled sediment volume (m3) using the five different hypothetical fault compared to measured sediment volume (m3) from core samples. The black horizontal line represents 100% of the measured sediment volume (m3) from core samples. The variation of Manning’s roughness coefficient shows the value needed to approximate the calculated volumes to the measured volumes.

Results of sediment transport simulations balance (erosion or deposition) for all sources using a Manning’s value of 0.025 for the entire domain are shown in Figure 4. Figure 4 shows profiles (AB and AC black lines on lower image of Figure 1) sediment transport simulations balance. The original surface is the present‐day topography from LiDAR. The GPR surface represent the erosional surface from 1755 CE tsunami. The follow lines represent the erosional surface from the five different hypothetical fault simulations. Erosion on southwest and northeast flanks on the dunes along profiles AB and AC (Figure 4) was reasonably well reproduced using CAW, HSF, MPF, and Scenario 1 sources. Results for simulations related to GB source fail to fully reproduce this erosional pattern.

4 Figure. Modeled sediment transport simulations balance along AB and AC profiles (see Figure 1 for precise location) using 0.025 m−1/3s Manning’s roughness coefficient. Nowadays topography from LiDAR in black line. The GPR surface representing the erosional surface from 1755 CE tsunami in red line. The others lines represent the erosional surface from the five different hypothetical fault simulations.

Comparison between historical data for arrival times and modeled results at Sines, Cabo de São Vicente, Lagos, and Huelva is presented in Table 1. The best overall match between documentary and modeled arrival times were obtained using the MPF and HSF sources, which yielded better correlations with the observed data (<3% and 4% differences respectively) than other sources. The worst correlation corresponds to CAW and GB sources with mean errors of 23% and 28%, respectively.

Gorringe Bank

The simulations results (time travel tsunami and sediment transport) from the sampled cores do not agree well with GB as a probable source, because the modeled runup were far too small to generate significant inundation and consequently unable to produce a tsunami deposit. Compared to historical observations, the modeled runups were smaller and the modeled tsunami travel time from GB was too large.

Cadiz Accretionary Wedge

For the CAW source, the modeled runup elevation agrees well with the historical data. However, the modeled tsunami travel time is shorter than indicated in the historical documents (i.e., modeled tsunami waves traveled ∼20% slower) (Table 1, Figure 4). Furthermore, the modeled sediment deposit volume is 8 times larger than the volume calculated from geological data retrieved from cores (Figure 5). In order to achieve an erosion depth compatible with the GPR (cross‐dune) profiles described by Costa, Costas, et al. (2016) it is necessary to use an unrealistical high roughness coefficient (>0.08) that is not in agreement with the land cover observed in the Salgados lowland (Chow, 1959). Furthermore, the resulting (modeled) deposition/erosion profile does not agree well with field data described in Costa, Costas, et al. (2016).

5 Figure. Erosion and deposition thickness using Manning’s roughness coefficient of 0.025 m−1/3s to all AOI. (a) CAW, (b) HSF, (c) GB, (d) MPF, (e) Scenario 1. Yellow dots indicate the location of observed sediment deposits from the 1755 CE tsunami. Green dots indicate the cores where the 1755 CE tsunami deposit was not observed. Therefore, the region between yellow and green dots roughly approximates a minimum estimate of the inland extent of the 1755 CE tsunami deposit.

Horseshoe Fault

The HSF model results do not agree well with the observed data. The modeled tsunami arrives 4 min later in Sines than in the historical record. Likewise, the modeled runup for Cabo de São Vicente was 5 m higher than reported in historical records. The modeled sediment volumes deposited in the lowland are compatible with objective observations when using Manning’s n coefficient values between 0.025 and 0.040, which are in broad agreement with the land cover. The model predicted larger amounts of erosion than observed on the seaward section of the dune along profile AB, and less erosion on the landward section (Figure 5). On profile AC, the model predicted no erosion.

Marquês de Pombal Fault

Both modeled tsunami travel times and runup magnitude correlate well with the observed historical data. There is also a close correspondence between the observed and modeled volume of the tsunami deposit when using a realistic Manning’s n coefficient of 0.04–0.05. However, the modeled dune erosion only partially matches the GPR data. Actually, the model fails to correctly reproduce erosion along profile AC on the landward region of the dune, possibly because the largest simulated wave could not overtop the dune; this obstructed representation of flow along its landward slope and computation of the corresponding erosion (Figure 5).

Scenario 1

Simulated tsunami travel times obtained from this setting yields good results in all target locations with the exception of Sines and Cabo São Vicente, where simulated travel times are higher. Furthermore, the modeled runup is ∼20% higher than reported in the historical records. However, it is noteworthy that the modeled patterns of dune erosion and the volume of sediment deposited in the lowland predicted by the models correlate well field data when a relatively high Manning roughness coefficient of 0.060–0.070 is used (Figures 4 and 5).

Discussion and Conclusions

The geological record at Salgados of the 1755 tsunami has the potential to provide an independent data set to validate tsunami inundation and sediment transport models and thereby to test hypothetical earthquake sources.

Of the five hypothetical scenarios presented above, the Marquês de Pombal and Scenario 1 provide the best overall match with both source‐to‐target tsunami travel time and runup taken from the documentary record sources. In addition, they also provide the best overall match in terms of predicted erosion/deposition patterns (e.g., total volume) obtained from field (geological evidences).

The source closest to the shore (Marquês de Pombal) yielded the best correlations between modeled and field data. This broadly confirms the region southwest of Cabo São Vicente as the most likely source area of the 1755 CE earthquake. This region has been previously proposed by Baptista, Miranda, and Victor (1992) based on the location of the February 1969 earthquake and also based on back‐ray tracing, travel times and wave heights (Baptista, Heitor, et al., 1998). In contrast, all simulated sources located further south in the Cadiz region (CAW and Scenario 1) over‐predict the volume of the tsunami deposit in Salgados lowland, the magnitude of runup reported in the documentary record, thus suggesting that a CAW source model is an (highly) unlikely source of the 1755 CE.

Although the GB source has been favored by Santos et al. (2009), its use in the context described herein leads to unacceptable mismatch with both sedimentary and hydrodynamic results as well as with the documentary record in terms of travel time and runup. In fact, it presented the poorest overall agreement results among all tested sources. This was mainly due to the large distance traveled by the tsunami waves (>200 km) before impacting the coast. Tsunami travel times and runup inferred from this source were consistently longer and smaller, respectively, when compared with field and historical data.

The numerical modeling approach used in this study that incorporates the geological record was able to partially constrain proposed 1755 CE earthquake sources. It is important to stress that this exercise does not unequivocally resolves the age‐old question about the 1755 CE source, nevertheless it points future directions for other fields of geoscience to pursue and, hopefully, it will contribute to the establishment of more reliable hazard assessments for Iberia and for the mid‐North Atlantic.

Acknowledgments

The authors thank the Fulbright Commission for Junior Faculty Member Award (2018) for the financial support, the Pacific Coastal and Marine Science Center/United State Geological Survey (USGS), the Foundation for Research Support of the Rio de Janeiro State (FAPERJ) for a Young Scientist of Our State Grant (2016). Authors would also like to acknowledge FCT funded project OnOff – Coupling onshore and offshore tsunami record: complementary tools for a broader perspective on tsunami events – PTDC/CTA‐GEO/28941/2017.

Data Availability Statement

The authors confirm that the data supporting the findings of this study are available within the article and on a public open database 4TU. Research Data with the https://doi.org/10.4121/uuid:99cf3b07-43e4-4a3b-986b-835811aad485. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the US Government. The authors would like to thank Profs. Luís Matias and M. A. Baptista and Dr. Rachid Omira for stimulating discussions on the previous version of this manuscript. The authors acknowledge the comments provided by the Editor Dr. Helen Graves, by Prof. David Tappin and one anonymous reviewer that greatly contributed to improve this work.