Evidence of Surface Rupture Associated With Historical Earthquakes in the Lower Tagus Valley, Portugal. Implications for Seismic Hazard in the Greater Lisbon Area

The Lower Tagus Valley Fault, Portugal, has long been associated with the damaging earthquakes that affected the Greater Lisbon Area in historical times. These include a poorly documented earthquake that occurred in 1344, the relatively well-documented 1531 earthquake, and the most recent M6.0 1909 earthquake. In this work, we use a 0.5 m resolution LiDAR-based digital elevation model and a 0.5 cm resolution digital surface model based on UAV photogrammetry to accurately locate the fault scarps in the northernmost portion of the western fault strand and to select sites to perform paleoseimolological investigations. The paleoseismological and geochronological analysis performed in the Alviela trench site document the fault activity in the last 3000 years, including two earthquakes during historical times. We performed ground motion scenarios for 20 km, 40 km, and 60 km ruptures including the trench site. The ground motion fields obtained for the 40 km and 60 km ruptures are in agreement with most macroseismic intensity data available for the 1531 earthquake, implying a magnitude in the range M6.8–7.4. However, the degree of deformation preserved in the trench suggests a value closer to the lower magnitude bound. The intensity level observed in Lisbon in 1531 (IX) is lower than the modeled intensities for all considered scenarios and could be related to a particularly high level of vulnerability of the building stock.


INTRODUCTION
The association of historical events with specific faults is enormously important when defining seismic sources for seismic hazard analysis. The detailed knowledge of the fault's location and kinematics is fundamental to accurately estimate the level of ground motion expected in future earthquakes (e.g., Boore et al., 2014;Chiou and Youngs, 2014).
The Lower Tagus Valley (LTV) displays the highest levels of seismic hazard in Western Iberia (Vilanova and Fonseca, 2007;Giardini et al., 2014;Silva et al., 2014;Woessner et al., 2015). It has been affected by damaging earthquakes in 1344, 1531, 1755, and 1909, with magnitude estimates for local events ranging from 6.0 to 6.9. (e.g., Dineva et al., 2002;Stich et al., 2005;Vilanova and Fonseca, 2007;Fonseca, 2020). The 1531 earthquake was the first known event in Portugal for which the damages were reported in a set of different localities (Moreira, 1984;Justo and Salwa, 1998). The offshore 1755 Lisbon earthquake was the focus of detailed coeval inquiries and has a rich macroseismic database both in Portugal and Spain. Vilanova et al. (2003) proposed that one of the several shocks reported for this event occurred onshore, within the LTV to explain both the large degree of damage far away from the rupture and reports of ground deformation. The 1909 earthquake was subject to detailed geological and macroseismic investigations (Choffat and Bensaúde, 1912). The fault ruptures associated with these events have, however, never been documented.
The Lower Tagus Valley Fault (LTVF) is a left-lateral strikeslip structure that displays two nearly parallel strands along both river margins (Figure 1). The surface evidence of fault-activity in the region is subtle both due to the erosional processes related to the fluvial activity of the Tagus river, and the extensive anthropogenic modifications to the landscape shaped by urban expansion and intensive agricultural activities. Besana-Ostman et al. (2012) and Canora et al. (2015) mapped the fault traces using stereographic analysis of old aerial photographs, and the paleoseismological studies of Canora et al. (2015) at the eastern strand reported the occurrence of 6 earthquakes in the past 10.000 years.
In this work, we use a 0.5 m resolution digital elevation model based on LiDAR and a 20 cm digital surface model based on UAV photogrammetry to identify and document the fault location at the Western Strand, and to select suitable sites to perform paleoseismological studies. We employed paleoseismological and geochronological techniques to identify the earthquake history recorded at the Alviela trench.
To test whether a rupture through the Alviela trench site is consistent with macroseismic intensity data, we performed seismic scenarios for different ruptures on the fault.

SEISMOTECTONIC SETTING
The Lower Tagus Valley (LTV) is an NNE-SSW trending valley located in the central part of Portugal, which comprises part of the densely populated Great Lisbon Area ( Figure 1). The valley is located within the Lower Tagus Cenozoic Sedimentary Basin (LTB), a NE-SW elongated tectonic depression that includes up to 2000 m of Tertiary sediments (Paleogene to Pliocene), Pleistocene fluvial terraces, and up to 70 m of Upper Pleistocene to Holocene alluvial cover (Cabral et al., 2004;Cabral et al., 2013).
The LTV developed in the Neogene as a compressive foredeep basin linked to the tectonic inversion of the Mesozoic extensional Lusitanian Basin in response to the NW-SE Miocene compression (Ribeiro et al., 1990;Rasmussen et al., 1998). The area is currently affected by the NW-SE trending convergence between the Eurasian and African plates at a rate of 3.8-5.6 mm/ yr (e.g., Fernandes et al., 2003;Nocquet and Calais, 2004;DeMets et al., 2010). In this tectonic context, some of the inherited basement faults in western Iberia are considered active or potentially active under the current stress regime, acting as reverse and strike-slip faults (Arthaud and Matte, 1975). Examples of this type of reactivated faults include the Manteigas-Vilariça-Bragança Fault, the Penacova-Régua-Vérin Fault, and the Lower Tagus Valley Fault (Arthaud and Matte, 1975). Moment tensor inversions for instrumental earthquakes indicate a preponderance of reverse and strike-slip faulting (Stich et al., 2010;Domingues et al., 2013;Custódio, et al., 2016).
Several moderate to large earthquakes affected the LTV during historical times. The 1344 earthquake is poorly documented and has been assigned M6.7 (Vilanova and Fonseca, 2007) based on the maximum Modified Mercalli Intensity of XI-X in Lisbon (Moreira, 1979). The 1531 earthquake caused widespread damage in the LTV with MMI values of IX spanning over 70 km along the valley, from Lisbon to Santarém (Justo and Salwa, 1998). Magnitude estimates for this event range from Mw6.5 (Stucchi et al., 2013) to M6.9 (Vilanova and Fonseca, 2007).
The most recent significant earthquake in the region occurred on April 23, 1909 (Choffat and Bensaúde, 1912;Fonseca and Vilanova, 2010;Teves-Costa and Batlló, 2011). Teves-Costa and Batlló (2011) assigned a maximum MMI intensity of X based on the original report of Choffat and Bensaúde (1912). Greatly damaged localities included Benavente, Samora Correia, and Santo Estevão, settled in the eastern margin of the Tagus river, around 40 km NE of Lisbon. The earthquake caused 47 fatalities, serious property damage, and extensive liquefaction effects within the Tagus alluvial sediments. Based on the scarce number of instrumental records available, Stich et al. (2005) computed a reverse-fault moment tensor solution with a NE-SW trending strike and Mw6.0 for the 1909 event.
An additional earthquake may have occurred within the valley on November 1755 (Vilanova et al., 2003;Fonseca, 2020), triggered by the Great Lisbon Earthquakes (e.g., Reid, 1914;Moreira, 1984;Johnston, 1996). Vilanova et al. (2003) stress that an additional shock would justify both the very high intensities reported for the valley (IX-X MMI according to Moreira, 1984), the reports of ground deformation, and the inconsistencies in the tsunami arrival times.
Despite the occurrence of moderate to large historical earthquakes, the LTV experienced very low levels of seismic activity during the last five decades (Borges et al., 2001;Vilanova and Fonseca, 2004;Custódio et al., 2015). Only a few micro-earthquakes have been recorded within the valley, down to a depth of 20 km (Fonseca and Long, 1991).

The Lower Tagus Valley Fault Zone as a Seismic Source
The LTVFZ is the main tectonic feature within the Lower Tagus Valley (LTV) based on its extension and seismic potential ( Figure 1). This NNE-SSW trending fault zone is considered to be a major structure formed during the Carboniferous Laurasia and Gondwana collision, which was reactivated under the current stress regime (Arthaud and Matte, 1977).
Despite several efforts to characterize the fault zone (Cabral, 1995;Vilanova and Fonseca, 2004;Cunha et al., 2005;Cabral et al., 2013), the location and kinematics for the LTVFZ remained controversial and poorly constrained. Gravimetric and seismic reflection data suggest a complex geometry (Cabral et al., 2003). While Cabral et al. (2003) suggest that this structure is a blind reverse fault, Vilanova and Fonseca (2004) argue in favor of strike-slip kinematics.
The identification of active structures along the LTV is complex due to the high erosion/sedimentation rate of the Tagus River coupled with the moderate tectonic activity of the fault zone. Therefore, the evidence of active tectonics at the surface tends to be concealed by fluvial dynamics. Additionally, the valley is affected by intense anthropogenic modification due to agricultural activities and urban development. Despite the challenges, Besana-Ostman et al. (2012) and Canora et al. (2015) have successfully identified and mapped active tectonic structures along both river margins valley, using stereographic analysis of old aerial photographs and field surveys. Besana-Ostman et al. (2012) produced the first map of active fault traces in the western margin of the Tagus River. They mapped the LTVFZ as an approximately 80 km-long strikeslip structure that transects Miocene to late Quaternary Frontiers in Earth Science | www.frontiersin.org March 2021 | Volume 9 | Article 620778 deposits. Following a detailed quantitative morphological analysis using the mountain-front sinuosity index (Bull and McFadden, 1977), they conclude that the mountain fronts associated with the northern and southern sections of the fault are rated as Class I (active fronts) and Class II (intermediate), respectively. Canora et al.
using geomorphological, paleoseismological, and seismic reflection analysis identified and characterized the Eastern strand of LTVFZ. They mapped the fault for around 70 km and found evidence of six surface fault ruptures in the last 10 ka, with a minimum slip rate of about 0.14-0.24 mm/yr for the fault zone and a mean recurrence interval for surface ruptures around ∼1.7 ka.

DATA AND METHODS
The Tagus river-dynamics is the main responsible for shaping the landscape in the study area, but so is the enormous amount of anthropogenic activity related to agriculture. Since the tectonic activity of the LTVFZ is moderate, the preservation of surface ruptures along the fault is relatively rare. Additionally, the ground surface is partially covered by vegetation, including evergreen forests and dense crops, further preventing the identification of fault-related features. For this reason, Airborne Light Detection and Ranging (LiDAR) data has been acquired in 2014 for two narrow bands including the northern and central parts of both the Eastern and the Western strands of the LTV. Based on the LiDAR point cloud data a 0.5 m-resolution bare-earth Digital Elevation Model (DEM) was derived for the study area (Foroutan et al., 2016).
In this work, we focused on the northernmost section of the western strand of the LTVFZ (Figure 1). We reviewed the fault traces identified by Besana-Ostman et al. (2012) using the 0.5 mresolution DEM, 1:30,000 scale aerial photos, and 1:50,000 scale geological maps (Serviços Geológicos de Portugal, 1952;Instituto Geológico e Mineiro, 1999). The active faulting geomorphic  features identified in the DEM were further investigated using detailed field surveys. Some of the features identified in the old aerial photos (1946)(1947)(1948)(1949)(1950) were, however, no longer noticeable in the 2014 DEM due to the major landscape modifications that occurred in that time period. Figure 2 shows 21 ± 1 m and 42 ± 1 m left-laterally offset river-channels located within the study area. The age of the channels affected by the fault activity is unknown. However, considering that the drainage network is incised in the deposits of the Q4 terrace (see Figure 2 for location), it is possible to approximate an age for the displacements. Martins et al. (2010) report Luminescence dating ages of 99 ± 6 ka for the Q4 terrace deposits. Therefore, we can estimate a horizontal sliprate ranging from ∼0.21 ± 0.02 mm/yr to ∼0.42 ± 0.03 mm/yr. The uncertainties presented are based on 1) the river offset measurements performed on the DEM and 2) the uncertainties in the terrace deposits dating. However, the actual uncertainties are much larger, as it would be necessary to carry out detailed mapping to ratify the origin of the materials eroded by both rivers and, thus, to have a higher degree of reliability on the assessment of the slip rate.
Other features that can be observed both on the LIDAR-based DEM and in the field are fault scarps, linear valleys, and oriented drainage channels ( Figure 3). Those features while concurrent with left-lateral strike-slip faulting also suggest some degree of vertical displacement. Pliocene and Pleistocene materials display the most abrupt scarps, whilst in the areas filled with Holocene deposits, those structures are very subtle ( Figure 2 and Figure 3). These latter areas are, however, ideal for carrying out paleoseismic studies for characterizing the most recent activity on the fault. In the ensuing section, we detail the sites selected for paleoseismological studies. For these sites, we produced very-high-resolution digital surface models (DSMs) derived from photogrammetry. The DSMs were built from 416 individual images captured with a DJI Phantom 4 PRO, equipped with an RGB camera, using front and lateral overlaps of 80%. The flights were performed using a double-grid design, at heights of 10 m above ground. The images were processed with Agisoft Metashape Professional Software, to build a DSM with a 0.5 cm resolution. Further details on the data processing can be found in Pereira et al. (2020).

PALEOSEISMIC STUDIES ON THE WESTERN LTVFZ
Field campaigns have enabled us to select two sites-Boavista and Alviela-to carry out detailed paleoseismic studies (Figures 4 and Frontiers in Earth Science | www.frontiersin.org March 2021 | Volume 9 | Article 620778 6 5). In both cases, a fault scarp affects the most recent deposits. In the case of the Boavista site, the scarp is very prominent and places Plio-Pleistocene materials in contact with recent alluvial deposits ( Figure 4). In the Alviela area, the fault scarp is located within Holocene alluvial deposits and its height is less than 2 m ( Figure 5). Although the Boavista area was more promising in terms of morphology, the trench exposed a massive unit of clay deposits, 3 m thick, without any sedimentary indicator that could permit the identification of the fault. Due to the proximity of the Tagus River, the sedimentation rate in this area is tremendous, much higher than the deformation rate. The volume of sediments deposited during recent overflows in this area exceeded 2 m in thickness. In contrast, the Alviela trench has shown abundant paleoseismic information, despite the fault surface expression being very subtle, probably because is not located within the Tagus floodplain but in a smaller affluent.

Alviela Trench
At the Alviela paleoseismological site ( Figure 5), a 1-m NE fault scarp offsets the alluvial deposits of the Alviela River, an affluent of the nearby Tagus River. The scarp slightly uplifts the northwestern block creating an area of ponding, an optimum context for developing the fine sedimentation patterns that are critical to paleoseismological studies ( Figure 5). The scarp is clearer in the old aerial photos than in the LiDAR DEM. Both in the field and in the 0.5 cm-resolution UAV photogrammetrybased DSM the scarp is extremely subtle, probably due to persistent plowing of the land ( Figure 6).
A trench of 30 m long and 2 m deep was excavated at the Alviela site ( Figure 6) perpendicular to the fault scarp identified on the LIDAR digital model. The trench exposed a sequence of alluvial layers (units A-E, Figure 7) corresponding mainly to clays and clayey sands typical of alluvial floodplain deposits, and a system of channels, filled by sands, clays, and silts, entrenched on top of it (units F and G; Figure 7). Unit H corresponds to the topsoil developed in different materials and has been greatly altered by anthropic activity (Figure 7).

Age Control
We used the radiocarbon technique on charcoal samples to constrain the age of the units exposed in the trenches. The results are presented in Figure 8. We dated eleven samples, ten of which gave consistent ages for the materials. Only the relationship between samples As3 and As5 is irregular since the stratigraphic order indicates that unit A (2146-1996 cal BP) is older than unit B (2854-2755 cal BP). The reworking of the carbon in sample As3 may have caused that inconsistency. Otherwise, contamination could be affecting the age of one of the samples. At any rate, the ages of deposits A and B do not affect the seismic history recorded in the trench, as we will discuss next.

Evidence of Paleoearthquakes
The Alviela trench exposed a sequence of deformed clay and sand strata, interpreted as deposited in an alluvial floodplain environment, and a sequence of channelized units entrenched into the top of the section (Figure 7). The 2 m-wide fault zone is mainly composed of vertical fault planes affecting all the strata (Figure 7). The cumulated deformation is larger for the older exposed units, suggesting the occurrence of consecutive surface ruptures. A minimum cumulative vertical displacement of 0.25 ± 0.05 m has been estimated for these units, which includes the offset observed on each fault strand. Tracing the fault location was straightforward on the youngest sandy sedimentary units and particularly difficult through the clayey layers. We have detected Frontiers in Earth Science | www.frontiersin.org March 2021 | Volume 9 | Article 620778 slickensides, N155°E/75°E, and 35°N rake, affecting the most recent trench deposits ( Figure 7). Here, the fault offset is attenuated with a minimum vertical displacement of 0.12 ± 0.05 m and a net displacement of 0.21 ± 0.09 m for the most recent units (units F and G). The fault direction provided by the striations is different from the general fault trend (N45°E) observed in the LIDAR digital terrain model and the field probably due to the complexity of the surface rupture produced by the earthquakes in this type of poorly consolidated floodplain deposits. This direction is similar to that observed along the fault trace at other scales ( Figure 6), as can be seen in Figure 2, related with a releasing step-over part  We infer the occurrence of at least two ruptures using structural (fault displacement), sedimentological (changes in the depositional environment), and stratigraphic criteria. The timing of these two events is very well constrained thanks to the good recent stratigraphy exposed in the trench, and a large amount of organic material available to date the deposits.
The oldest deposits in the trench (units A to D) bracket the oldest event E1. Although vertical offsets affecting the deposits are always small, these units are more deformed than the most recent ones (units E to G; Figure 7). The activity of the fault related to this first event (E1) generated a fault scarp that prevented the more recent deposits (units E to G) from sedimenting to the right side of the trench. We consider that event E1 is younger than unit D (664-552 cal yr BP) and older than unit E (514-428 cal yr BP; see Figures 7, 8). Event E2 is the latest that we can identify in the trench. It is affecting all units, including unit H (Figure 7), and shows a fairly small vertical displacement. This event is restraining and offsetting the small sand layers interbedding in unit F. The age of this earthquake is younger than unit H, which is dated around 500-315 cal yr BP (see Figure 8). The two events described above have been correlated with the historical earthquakes recorded in the region (Figure 8). Event E1 could eventually correspond to the 1344 earthquake. For event E2 there are three possibilities in the historical period: the 1531 historical earthquake, the 1909 Benavente earthquake, or the eventually triggered earthquake in 1755. There are no further historical records of damaging earthquakes affecting the LTV between those two events, despite many earthquake damages being reported all over Portugal in particular since 1700AD. The surface-rupture mean-recurrence interval is roughly estimated at 350 years, considering that two events occurred after the deposition of unit D (664-552 yr cal BP).

CONSTRAINTS FROM GROUND MOTION SCENARIOS AND MACROSEISMIC INTENSITY DISTRIBUTIONS
Earthquake scenarios are powerful tools for investigating the ground-shaking patterns of both future and past earthquakes. In this section, we develop a set of earthquake scenarios for different ruptures of the Western strand of the LTV fault. We compare the motions with the macroseismic intensity data available from the historical earthquakes in the time frame of the deformation identified in the trench.
We assume that the northernmost position of the rupture is located 7 km northeast of the trench site, at a small a stepover (Mato de Miranda), as it is generally acknowledged that geometric discontinuities may act as barriers to the rupture propagation (e.g., King and Nábělek, 1985;Scholz, 1990). We tentatively consider the southernmost position of the rupture to be located in the vicinity of Santarem, Azambuja, or Alhandra.
The ruptures considered are around 20 km, 40 km, or 60 km long, corresponding to magnitudes in the range Mw6.4, Mw7.0, and Mw7.2 respectively, according to the model of Leonard (2014) for stable continental regions. Leonard's (2014) model uses the most extensive database on earthquake ruptures and relies on investigating empirical relationships between different fault rupture parameters (rupture width, rupture length, and displacement). Substituting those expressions into the definition of seismic moment, he retrieves scaling relationships that are self-consistent with each other and with the definition of seismic moment. Table 1 shows the estimates of Mw based on the length of the fault. The corresponding width values are also presented.
Based on Besana-Ostman et al. (2012) and the field evidence discussed previously, we consider the fault to be vertical and use a dip value of 0-15°for modeling purposes. We assume a seismogenic depth of 24 km, based on the 95% percentile of the depth distribution for earthquakes in the region (Vilanova et al., 2014). When the vertical projection of the rupture width reaches the seismogenic depth the earthquakes become widthlimited and scale differently according to Leonard (2014). For a seismogenic depth of 24 km, this corresponds to M7.5. Therefore all ruptures modeled are not limited by width.
Besides the geometric details of the fault, one needs to select which ground-motion models are most suitable to represent the level of shaking in the region. Vilanova et al. (2012) compared the regional ground motion data ) with a set of ground motion models using the residual analysis proposed by Scherbaum et al. (2004). They conclude that, in general, models developed for stable continental regions better represent the ground shaking levels experienced in regional earthquakes. However, since the database comprehends mostly moderate magnitude earthquakes, (∼M5-6) recorded at large distances (larger than 150 km) they acknowledge that a considerable degree of uncertainty prevails regarding the shaking intensity of future earthquakes in the region.
Earthquakes in stable continental regions display higher levels of stress-drop and lower anelastic attenuation (high quality-factor Q) with respect to earthquakes in regions characterized by more active tectonics such as plate boundary regions. Given the limited amount of data available, it is unclear whether the relatively high motions recorded in SW Iberia are related to high stress-drop, low anelastic attenuation, or both. Chen et al. (2018) use a global dataset of 1 Hz Lg coda quality factor (Q0) for mapping the anelastic attenuation of ground motion attenuation. However, the coverage of SW Iberia in the underlying dataset is scarce (Mitchell et al., 2008). Vales et al. (2020) investigate the Lg coda Q for West Iberia. They find Q0 values varying from 72 to 220. However, as stressed by Havskov et al. (2016) the processing parameter choice significantly influences the results (e.g., the Q0 value computed by Havskov et al., 2016, for Norway is 124 whereas that calculated by Mitchell et al., 2008, is around 750). Therefore, it is difficult to perform global comparisons between studies performed with regional data.
Given these considerations, we evaluate the seismic scenarios using ground-motion models developed for both stable continental regions and for active tectonic regions. We selected ground-motion models for PGA and PGV whose domain of application includes the Mw range of interest, and that consider the influence of site-effects. We analyzed three sets of ground-motion models with regard to the underlying database. The NGA-West2 models are mostly based on California earthquakes (84% of the data). That dataset is complemented with records coming from other global regions characterized by active tectonics. European and the Middle East models are based on the RESORCE dataset (Akkar et al., 2014), which includes mainly data from Turkey, Italy, and Greece. Finally, models developed for Stable Continental Regions are based on either stochastic simulations or hybrid methods. Table 2 summarizes the models considered in the analysis.
We computed the following ground motion fields for the following parameters using the OpenQuake Scenario tool with 1-sigma standard deviation: peak ground velocity (PGV), peak ground acceleration (PGA), and response spectra acceleration (SA) at 0.2 s, 0.5 s, 1 s, and 2 s. We used the Vs30 siteconditions model for Portugal developed by Vilanova et al. (2018), based on correlations between the 1:500,000 geological map and the regional database of shear-wave profile data ( Figure 9).
In general, the models based on a specific dataset produce similar ground motion fields for the fault ruptures analyzed. This suggests that those magnitudes are well represented in the databases. Sensitivity studies demonstrated that slight variations in dip and moment magnitude (68% confidence interval) do not significantly affect the motions. The rupture width has a limited impact on the ground motions given the distance metrics used by the ground motion models (closest distance to rupture or closest distance to the surface projection of the rupture) and the fact that the fault is nearly vertical.
The results show that although Atkinson and Boore's (2006) model predicts higher motions, this is restricted to the near field. Atkinson and Boore's (2006), Atkinson and Boore (2011) model reaches PGA∼2 g for M6.4 in the vicinity of the fault, while for the remainder models the maximum PGA is in the range 0.5-0.7 g at short distances. For PGV, Atkinson and Boore's (2006) model predicts values around 135 cm/s 2 , while the models developed for active tectonic regions estimate values in the range 50-70 cm/s 2 . At distances 30-100 km, the differences attenuate significantly and the model of Boore et al. (2014) for high Q predicts the lowest attenuation out of the models analyzed. Figure 10 shows the PGV ground motion fields predicted by different ground-motion models and for different rupture lengths. We present results for PGV because that parameter is better correlated with macroseismic intensity for higher intensity levels than PGA (e.g., Boatwright et al., 2001;Caprio et al., 2015). We show the influence of a 68% confidence interval in the site conditions model for the Boore et al. (2014) model with high Q. Despite the qualitative nature of macroseismic intensities, these are often the only available data for studying earthquakes in the historical and early instrumental periods. The geochronological analysis of the Alviela trench indicated that both events in the stratigraphic records had a possible correlation with earthquakes that occurred in the historical period.
The oldest earthquake recorded might be the 1344 earthquake, which has partially destroyed the main chapel of Lisbon's cathedral (Moreira, 1984). However, with such limited data concerning a single location, it is difficult to draw further conclusions about this event. Moreira (1985) estimated Modified Mercalli Intensities (MMI) for the 1531 earthquake at 18 locations and outlined the first macroseismic map for a seismic event in Portugal. Justo and Salwa (1998) re-evaluated the data and added few new records, both in Portugal and Spain, and estimated the intensities in the Medvedev-Sponheuer-Karnik scale. Finally, Baptista et al. (2014) performed a reevaluation of the data using the EMS-98 scale. The intensity database produced by FIGURE 10 | PGV Scenario calculations for hypothetical ruptures along the LTVFZ, including the Alviela trench site (located at the northernmost part of the rupture). The macroseismic intensity dataset for the 1531 according to Justo and Salwa (1998) converted to PGV through the Caprio et al. (2015) relationship for Italy is included.
Frontiers in Earth Science | www.frontiersin.org March 2021 | Volume 9 | Article 620778 11 Baptista et al. (2014) is, however, strongly model-driven since the final intensity values are assigned taking into account the source they propose for the earthquake (decreasing the intensities both in Lisbon and in Santarém). For this reason, we prefer to use in this study the dataset by Justo and Salwa (1998).
We converted the 1531 earthquake's macroseismic intensities to PGV through the relationships proposed by Caprio et al. (2015) for Italy. These authors explore the regional dependence of ground motion to intensity conversion using mixed intensity scales and derive invertible relationships for Italy, Greece, and California, besides that for the global dataset. We use the relationships for Italy due to the similarity of construction standards for localities with a significant amount of old buildings. Those relationships predict the highest degree of damage for the lowest value of ground motion (high vulnerability of infrastructures). The results are shown in Figure 10.
The lowest rupture analyzed, from Mato de Miranda to Santarém, which corresponds to M6.4, predicts PGV values significantly lower than those corresponding to the macroseismic intensities for most localities in the dataset with any of the models used in the analysis.
Ruptures of 40 km and 60 km, corresponding to M7.0 and M7.2 respectively, are in general, compatible with most observed intensities. The main exceptions to this are the Alviela and Alverca localities, which stand out in the 60 km scenario, displaying lower intensities that expected had the fault ruptured through those locations. On the other hand, the intensity for Lisbon is higher than estimated by all scenarios analyzed.
In Figure 11 we explore the possible influence of site conditions using the model of Boore et al. (2014) for high Q, which predicts the highest intensities at intermediate distances.
The results suggest that an M7.0 rupture could explain both Alverca and Alhandra intensity values, considering higher than median V30 values for those locations.
The geographic distributions of PGV values predicted by the M7.0 scenario for the localities of Santarém, Castanheira do Ribatejo, and Lisbon (see locations on Figure 11), by the ground motion models analyzed are illustrated in Figure 12.
The dispersion includes both the aleatoric nature of the ground motions, and distinct site conditions within the localities, and the horizontal lines indicate the 1-sigma bounds for intensity IX according to Caprio et al. (2015). All the ground motion models used in the analysis predict PGV values that are compatible with intensities IX in Santarém and Castanheira do Ribatejo. However, the justification of that intensity value in Lisbon would require either favorable site conditions (low velocity/high impedance) or high vulnerability of the infrastructures.
The macroseismic intensities for the 1909 earthquake are presented in Figure 13 together with the ground motion model of Kotha et al. (2016). Both the distribution of high intensities around Benavente, 45 km away from the Alviela trench site, and the moderate magnitude of the event (M6.0-6.2 estimated using early instrumental records, Dineva et al., 2002;Stich et al., 2005) make it most unlikely that the earthquake ruptured through the trench site. It is also unlikely that a triggered local shock associated with the Great 1755 Lisbon Earthquake would reach the trench site since the highest intensities (VIII-IX, IX, and IX-X) were located South of Vila Franca de Xira, 40 km away from the trench site (see https://www. emidius.eu/AHEAD, last visited December 2020), and the ground deformation was reported for the left river margin (Vilanova and Fonseca, 2007).

DISCUSSION
The paleoseismic study of the Alviela trench allowed us to establish that the Western strand of the LTVFZ is tectonically active, and to identify two earthquakes that have occurred during the historical period. The oldest event affected sediments could be associated with the 1344 earthquake, but that event is poorly documented. We propose that the most recent deformation detected in the trench corresponds to the surface rupture of the 1531 historical earthquake.
Based on the offset of the sediments recognized in the Alviela trench and using Leonard's (2014) empirical relationships, we have calculated some parameters for the last event recorded in the stratigraphy. For net displacements of 0.21 ± 0.09 m, we obtain a surface rupture length of approximately 6 km and a magnitude of Mw∼5.5. The 1531 earthquake macroseismic intensity dataset would then strongly disagree with any seismic scenario calculated for ruptures within this magnitude range. Asymmetric offset distributions along fault ruptures have been reported for many seismic events (e.g., Stein et al., 1997;Lin et al., 2002;Manighetti et al., 2005;Wesnousky, 2008;Stirling et al., 2017). Hence, we consider that the displacement identified in the walls of the Alviela trench may not be representative of the overall surface displacement produced by the earthquake. Besides, the average displacement per event produced by an earthquake tends to decreases toward the extremes of the rupture (Choi et al., 2012). The location of the Alviela trench, near the end of the rupture (as modeled in the seismic scenarios and consonant with the local intensity data), could be limiting the extent of deformation found in the trench.
Assuming that no earthquake with magnitude large enough to produce significant damages (M ≥ 6) occurred since 1531, as suggested by the historical records, and taking into account the intensity distribution available for both the 1755 and the 1909 earthquakes, the paleoseismic study performed in this work constrains the rupture of the 1531 earthquake to include the Alviela trench-site. For this reason, we incorporate the fault portion containing the Alviela trench in all modeled seismic scenarios. However, we cannot rule out the possibility of an unrecorded lower-magnitude earthquake affecting the trench area since the M4.9 Teil earthquake ruptured to the surface (Ritz et al., 2020), despite that probably being a rare phenomenon.
The distribution of the intensities calculated for the 1531 earthquake is, in general, in agreement with the seismic scenarios produced for a rupture length around 40 km and Mw ∼7.0. The intensity estimated for Lisbon for that scenario is, however, lower than that reported by Justo and Salwa (1998). Baptista et al. (2014) suggested that the macroseismic intensity in the city of Lisbon could have been overestimated due to particularly strong site effects, since most greatly damaged buildings were located in the lower part of the city, in areas of sedimentary filling. In addition, these authors stress the weakness of the foundations of most of the affected buildings and the type of buildings as a possible cause of the extensive damage. This same argument is used by other authors such as Camelbeeck et al. (2014) to justify the extensive destruction suffered by masonry buildings in central and western Europe affected by moderate historical earthquakes. We, therefore, consider that the M7.0 scenario is the most plausible of the three presented here to justify both the paleoseismic data and the distribution of macroseismic intensities for the 1531 earthquake. The scenario calculations also show that it is unlikely that the 1531 earthquake had Mw∼6.5 as proposed by Stucchi et al. (2013) (see https://www.emidius.eu/SHEEC/sheec_1000_ 1899.html, last visited December 2020). Current seismic activity in the Lower Tagus Valley area is remarkably low compared to the historical seismicity recorded in the area, with at least three moderate-to-large-magnitude earthquakes in the past seven centuries (1344, 1531, and 1909 events). The level of seismic activity from historical records may suggest recurrence periods of a few hundred years for moderate to large earthquakes if a quasi-periodic temporal model is assumed. However, clustering periods intercalated with quiescence periods are usual for intraplate faults that display low long-term slip rates and high-stress-drop earthquakes (Kenner and Simons, 2005). It has also been suggested that in intraplate regions, a moderate-tolarge earthquake may weaken the regional lithosphere facilitating subsequent ruptures even in the absence of strong tectonic loading (Li et al., 2009), thus promoting the concentration of seismicity both in time and space.
The paleoseismological studies of Canora et al. (2015) for the Eastern strand indicate much longer mean recurrence periods for surface ruptures (1.7 ka) than that reported for the Alviela trench (0.6 ka). The active tectonic features are clearer for the Western strand of the fault, which could explain the contrasting recurrence values.
The slip-rate calculated assuming that the displaced drainage dates from Pliocene Q4 terrace deposits allowed us to estimate a sliprate in the range 0.2-0.4 mm/yr. The morphometric calculations of Besana-Ostman et al. (2012) for the northern western segment of the LTVFZ indicate a Class-1 type fault, similar to that of the wellestablished Vilariça fault, for which a long-term slip rate of 0.3-0.5 mm/yr was reported (Rockwell et al., 2009).
The complexity of spatial-temporal patterns of the faults in intraplate regions significantly affects the assessment of both the seismic hazard and the seismic risk. The earthquake clustering phenomena may lead to overestimating the seismic hazard in areas where large events occurred in the recent past and underestimating it in regions experiencing long periods of seismic quiescence. It is critical to clarify whether the historical earthquakes affecting the area are part of a period of activity, followed by a long period of quiescence, or the seismic cycle of these faults is quasiperiodic, in which case the earthquake recurrence could be much higher. We consider the first option as the most plausible since, as suggested by Stirling et al. (2012), quasi-periodic seismic cycles are very unlikely for complex intraplate regions.

CONCLUSION
The paleoseismic study of the Alviela trench shed light on the recent seismic activity of the Western strand of the LTVFZ (past 3000 years) and allowed us to identify the occurrence of two earthquakes rupturing to the surface in historical times. The oldest event occurred between 664 and 552 cal yr BP and 514-428 cal yr BP and could eventually correspond to the poorly documented earthquake known to affect the city of Lisbon in 1344. The most recent deformed unit in the trench dates from 500 to 315 cal yr BP and was most probably affected by the surface rupture of the 1531 historical earthquake associated with the LTVFZ. It is unlikely it corresponds to the ruptures of the eventual triggered local rupture in 1755 or in the 1909 earthquake given the spatial distribution of macroseismic intensities. However, we cannot rule out the possibility that an event unreported in the historical catalog, is responsible for deforming the trench deposits. Seismic scenario calculations for ruptures including the trench site and for Mw6.8-7.4 are in agreement with the macroseismic field available for the 1531 earthquake. The degree of deformation found in the trench favors the scenario performed for 40 km-long rupture and Mw∼7.0. In any case, the intensity IX in Lisbon is underestimated, which could result from the fact that our model could not capture the particularly strong site-effects, or that the building stock had a higher degree of vulnerability than that assumed by the intensity scales.
Paleoseismological studies are critical to investigate the longterm behavior of active faults, including those located in difficult environments such as the LTV. We successfully performed a paleoseismological investigation in a challenging area by using a methodological approach that combines stereographic analysis of old aerial photographs, high-resolution DEMs, and field surveys. In the LTV low tectonic-rates are coupled with high degrees of erosion and sedimentation associated with the dynamics of a major river and large anthropogenic disturbances of the natural ground surfaces due to agriculture and urban development. We overcame these difficulties by 1) focusing on a site with condensed sedimentation, e.g., in an affluent of the Tagus River rather than on the Tagus floodplain, and 2) identifying and accurately locating in the field the subtle scarps recognized in the old aerial photographs.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
CC and SV conceived the study. CC, CV, and YD performed the trench studies. CC wrote the first draft of the manuscript. SV performed the scenario calculations and wrote the corresponding sections of the manuscript. PP and SH developed the DSMs for the trench sites. All authors discussed the results and contributed to final version of the manuscript.

ACKNOWLEDGMENTS
The Portuguese Foundation for Science and Technology (FCT) funded this work through the research project Seismic Hazard: probabilistic Assessment at long Return Periods (SHARPE) (IF-EXPLOR). S. P. V. acknowledges FCT for her contract IF/01561/ 2014/CP1214/CT0006 under the IF2014 Program. The UAV flights and processing were done in the scope of Project ALTITUD3 (PTDC/EAM-REM/30475/2017) funded by the FCT. Centro de Recursos Naturais e Ambiente (CERENA) research unit is funded by FCT through strategic project UID/ ECI/04028/2013. We are grateful to Madalena Abecassis and José Infante da Câmara for permitting us to excavate trenches in their properties and for providing logistic support to the fieldwork. We are grateful for the supportive comments received from editor Ramon Arrowsmith and two anonymous reviewers. Their remarks and suggestions increased the clarity and overall quality of the paper.