Reassessing Eastern Mediterranean tectonics and earthquake hazard from the AD 365 earthquake

The hallmark of great earthquakes (Mw ≈ 8,3-8.5) in the Mediterranean is the 21 July AD 365 earthquake and tsunami that destroyed cities and killed thousands of people throughout the Eastern Mediterranean. This event is intriguing because most Mediterranean subduction forearcs exhibit pervasive crustal extension and minimal definitive evidence exists for great subduction megathrust earthquakes, consistent with weak seismic coupling. This conundrum has led many to favor rupture of a previously unrecognized upper plate splay fault south of Crete in an Mw 8.38.5 earthquake, uplifting a Holocene paleoshoreline on Crete by up to 9 m. Similar source mechanisms have been adapted and commonly used for seismic and tsunami hazard estimation in the region. We present an alternative model for the uplift of the Cretan paleoshoreline and the AD 365 tsunami that centers on known active normal fault systems offshore of western and southwestern Crete. We use new and published radiocarbon dates, together with historical records, to show that uplift of the Cretan paleoshoreline likely occurred during two or more earthquakes within 2 to 3 centuries. Visco-elastic dislocation modeling demonstrates that the rupture of these normal faults fits observed data equally as well as reverse fault models but requires less slip. Tsunami modeling shows that normal-fault ruptures produce strong tsunamis that better match historical reports than a hypothetical reverse fault. Our findings collectively favor the interpretation that damaging earthquakes and tsunamis in the Eastern Mediterranean originate on normal faults and highlight the potential hazard from tsunamigenic upper plate normal fault earthquakes. Plain Language Summary Most people living and vacationing near the Mediterranean Sea coast are not fully aware of the region's earthquake and tsunami hazard. Here we contribute to understanding the mechanisms for major earthquakes and tsunamis in the Mediterranean by investigating the region’s largest historically documented earthquake. The record of this event is thought to be preserved in part as a fossil beach uplifted by up to 9 m on the island of Crete, Greece. Previous studies assumed that the fossil beach was uplifted during a single earthquake in AD 365. However, our results from the dating of marine fossils that died due to sudden emergence above sea level and an assessment of existing historical and archeological records suggest that there were a series of earthquakes that might have incrementally uplifted the fossil beach. We identify manuscript submitted to AGU Advances and model a previously overlooked source for these earthquakes (normal faults) and tsunamis and find that these sources perform as well as or better than the traditionally assumed earthquake sources when compared to observations. These results highlight the potential importance of considering normal-fault earthquake sources in regions where tectonic plates converge and identify future research directions for more comprehensive hazard characterization.


Introduction
On 21 July AD 365, a massive earthquake was widely felt around the Mediterranean. The earthquake generated strong ground motions and produced a tsunami that radiated throughout the Eastern Mediterranean basin, destroying cities and causing numerous casualties (Ambraseys, 2009;Papadopoulos, 2011). In Alexandria, tsunami devastation was so severe that for centuries after the event, the anniversary was commemorated as the "day of horror" (Stiros, 2001). The AD 365 event is generally considered the largest seismic event recorded in the Mediterranean and has been the subject of intense research (e.g., Pirazzoli et al., 1996;Stiros, 2001Stiros, , 2010Shaw et al., 2008;Mouslopoulou et al., 2015b). The event is particularly interesting because subduction zones in the Mediterranean lack the great megathrust earthquakes that typify subduction zones elsewhere in the world (Jackson and McKenzie, 1988). Instead, the main seismic hazards are associated with widespread normal faulting in the forearc crust of the overriding plate. These features-a low degree of seismic coupling on the subduction boundary and widespread extensional deformation in the forearc-are typical for retreating subduction zones and are commonly attributed to the fast rollback of the subducting plate (Royden, 1993).
Historical records and geomorphic evidence of an uplifted Holocene paleoshoreline on Crete, which is found at maximum elevations of ~9 m above sea level at the southwestern tip of the island and spans a ~160 x 80 km area, place the event epicenter offshore of southwestern Crete (Fig. 1). We refer to this prominent paleoshoreline, which is almost continuously mappable around western Crete, as the Krios Paleoshoreline, based on the location of its maximum elevation at Cape Krios in southwestern Crete. Archeological evidence and radiocarbon dating of marine fossils demonstrate that the Krios Paleoshoreline was uplifted in the first centuries AD ( Fig. 1) (Pirazzoli et al., 1996;Dominey-Howes et al., 1998;Shaw et al., 2008). Numerous studies argue for all of the uplift to have occurred in a single event in AD 365, based on the correlation between radiocarbon ages and historical earthquake reports (e.g., Shaw et al., 2008; manuscript submitted to AGU Advances Stiros, 2010;Mouslopoulou et al., 2015b). Hinging on this evidence, and considering that the epicenter lies above the Hellenic subduction zone (HSZ), the largest, fastest, and most seismically active plate boundary in Europe, a logical source for the AD 365 earthquake is the subduction megathrust. However, dislocation modeling demonstrates that Crete is too far landward for megathrust earthquakes to produce significant uplift, as required by the deformed Krios Paleoshoreline (Ganas and Parsons, 2009). manuscript submitted to AGU Advances  (Ganas and Parsons, 2009;Vernant et al., 2014). The red box highlights the extent of (a). (c) Structural map of the Hellenic Subduction Zone as mapped by Chamot-Rooke et al. (2005) Based on the apparent inability of the HSZ megathrust to reproduce the paleoshoreline uplift in a single event, Shaw et al. (2008) and Stiros (2010) argued that the seismic source was a previously unrecognized upper plate reverse fault that potentially splays off the subduction interface and daylights in the Hellenic Trough, a 4-5 km deep bathymetric depression southwest of Crete (Fig. 1). Using dislocation and tsunami modeling Shaw et al. (2008) showed that the reverse fault interpretation would have generated an Mw 8.3 -8.5 earthquake and strong and widely travelled tsunami, suggesting this reverse fault represents an earthquake and tsunami hazard equivalent to hazardous subduction megathrusts globally. Based on this work, subsequent studies have argued that the forearc of the HSZ may have multiple previously unidentified splay thrusts and that all may be capable of generating > Mw 8 earthquakes and strong tsunamis manuscript submitted to AGU Advances (Stiros, 2010;Shaw and Jackson, 2010;England et al., 2015;Mouslopoulou et al., 2015b). Thus, our current understanding of tectonics and seismic hazard of the densely populated Eastern Mediterranean is intimately related to our understanding of the AD 365 event.
While these results are provocative, at present, there is no direct evidence for surfacebreaking reverse faults in the Hellenic forearc. Geological data demonstrate that the upper crust of the Hellenic forearc has been dominated by extension since the Miocene (Angelier et al., 1982;van Hinsbergen and Meulenkamp, 2006), and all known, observed active faults on subaerial forearc highs, such as Crete and Rhodes, are extensional (Angelier et al., 1982;Fassoulas, 2001;Gallen et al., 2014). No active contractional structures have been definitively imaged offshore between Crete and the Mediterranean Ridge Accretionary complex (Lallemant et al., 1994;Chamot-Rooke et al., 2005). Seismic imaging and interpretation of International Ocean Drilling Program IODP drill cores indicate that the Hellenic Trough ( Fig. 1) is a sediment underfilled half-graben bound by a southwest dipping extensional fault (Lallemant et al., 1994;Bohnhoff et al., 2001). Furthermore, wide-angle seismic imaging of the crust down to the plate interface of western and southwestern Crete demonstrate SW-NE and N-S extension, respectively, with no evidence for an embedded thrust fault splaying off the plate boundary and daylighting in the Hellenic Trough (Bohnhoff et al., 2001).
Interpretations in support of a contractional structure southwest of Crete are indirectly derived from several focal mechanisms (Shaw and Jackson, 2010), GPS stations between the forearc and the volcanic arc (Shaw et al., 2008;England et al., 2015), and the assumption that the AD 365 event singularly uplifted the Krios Paleoshoreline. However, GPS data show only minor contraction (1-4 mm/a) (Shaw et al., 2008;Saltogianni et al., 2020) between the forearc and the volcanic arc, yet due to the lack of offshore GPS data, this signal is indistinguishable from elastic deformation associated with minor locking on the subduction interface (Vernant et al., 2014).
Most importantly, there is little definitive support nor requirement for the uplift of the Krios Paleoshoreline in a single event.
Here we contribute to the debate of the AD 365 event by providing an alternative model based on the rupture of known active offshore normal fault systems. Our primary goal is to reconcile late Cenozoic-to-modern forearc deformation with models for the Holocene coseismic uplift as well as historic and geochronologic data. To this end, we present 32 new radiocarbon dates and use them in conjunction with published data to show that the distributions of ages at manuscript submitted to AGU Advances and below the most prominent paleoshoreline on Crete span a period coincident with a series of significant earthquakes. This new evidence implies that the uplift of the Krios Paleoshoreline occurred during at least two but likely several events within two to three centuries. Utilizing this premise we present an alternative hypothesis of Holocene uplift during multiple earthquake events on known active normal faults. These normal fault systems are large (from ~ 50 to ~ 100 km in length), close to the coast (< 5 km), and bound by deep (> 3 km below sea level) bathymetric troughs that represent grabens or half grabens (Fig. 1). We simulate co-and postseismic deformation on these normal fault systems proximal to the Krios Paleoshoreline and model tsunami propagation associated with fault rupture. Based on comparable simulations utilizing a single-event reverse fault, we discuss the feasibility of the two hypotheses with implications for Mediterranean tectonics and regional seismic and tsunami hazard.

Background
The HSZ accommodates about 35 mm/a of convergence between Africa and the Aegean, while the Africa-Eurasia convergence is only about 5 mm/a, illustrating fast rollback of the Nubian (African) slab (Reilinger et al., 2006;Vernant et al., 2014). Crete is located about 230 km north of the subduction trench and is a forearc high within the HSZ (Fig. 1). The crust below Crete is between 20 and 35 km thick (Bohnhoff et al., 2001;Makris and Yegorova, 2006;Snopek et al., 2007) and has been subject to pervasive extension since the Mid-to-Late Miocene (Angelier et al., 1982;ten Veen and Postma, 1999b;Brun et al., 2016). Accelerated rollback after 15 Ma induced extension and created numerous Neogene basins in the opening grabens (Angelier et al., 1982;Brun et al., 2016). While the basin bounding normal faults are active and thinning the upper crust, commonly both the footwalls and hanging walls of these structures are uplifted above sea level, indicating a deeper, geodynamic source for long-term regional uplift (Meulenkamp et al., 1994;Gallen et al., 2014;Ott et al., 2019b). Variable strike of these sedimentary basins suggests arc-parallel and arc-normal multidirectional extension (ten Veen and Postma, 1999b;Fassoulas, 2001). Uplift of sedimentary basins also indicates minimum Neogene uplift rates of ~ 0.2 mm/a (Meulenkamp et al., 1994), while Pleistocene paleoshorelines indicate uplift of 0.5 -1 mm/a for the south and west coast and lower, non-uplift, or subsidence in the east and north of the island (Gallen et al., 2014;Ott et al., 2019b;Robertson et al., 2019). Rock uplift is mostly driven by tectonics, since erosion rates are only ~ 0.1 mm/a, indicating that the manuscript submitted to AGU Advances component of isostatic adjustment due to erosional unloading is likely negligible (Ott et al., 2019a). Tectonic uplift of Pleistocene paleoshorelines is interpreted to be related to a regionalscale uplift signal augmented by local uplift in the footwalls of large normal fault systems (Gallen et al., 2014;Ott et al., 2019b;Robertson et al., 2019).
Seismic reflection data Bohnhoff et al., 2001;Alves et al., 2007), focal mechanisms (Bohnhoff et al., 2005), and active synthetic onshore faults (Caputo et al., 2010) indicate that two active normal faults lie offshore the western and southwestern coastlines of Crete (Fig. 1, text S6). We refer to the normal fault offshore western Crete as the Phalasarna Fault, and the fault offshore southwestern Crete as the Sfakia Fault Zone. Sonar data indicate 70-80 m high fault scarps, and the offset of seismic reflectors suggest long-term slip rates of 2.9 -5.8 mm/a for the Sfakia Fault Zone (Alves et al., 2007;Caputo et al., 2010). The Sfakia Fault Zone is also thought to be the offshore continuation of the on-shore Messara Graben ( Fig. 1e) (Peterek and Schwarze, 2004;Gallen et al., 2014). Additionally, Late Pleistocene paleoshorelines in western Crete record some of the islands' highest uplift rates with a potential acceleration towards modern (Ott et al., 2019b), indicating high slip rates on the Phalasarna Fault.

Radiocarbon dating and analysis
To better constrain the emergence of the Krios Paleoshoreline, we collected fossil samples of vermetids (Dendropoma sp.) and corals (Balanophyllia sp.) from at and below the paleoshoreline at eight different sites in western Crete ( Fig. 1, Tab. 1). In Phalasarna, Paleochora, and Chora Sfakion, we sampled vertical transects with six to nine samples per site and ~0.5 m elevation spacing to test for a signal or multi-stage paleoshoreline uplift. Our premise is that organism death and, therefore, our radiocarbon dates are largely related to seismic emergence from the sea. To illustrate this point, we refer to the radiocarbon dates as "emergence ages".
We selected samples from the farthest outgrown bioherm assemblages under the assumption that these are the youngest organisms. We avoided samples that had visual evidence of carbonate recrystallization (e.g., chalky texture and color). Radiocarbon samples were processed using standard techniques and investigated with a secondary electron microscope (Quanta 3D FEG) for signs of alteration (text S1, Fig. S8). The samples were calibrated using the Oxcal online calibration (Ramsey, 2009) program with the Marine13 calibration curve (Reimer manuscript submitted to AGU Advances et al., 2013) and an Eastern Mediterranean marine reservoir effect of 58 ± 85 years as recommended by Reimer & McCormac (2002). Our analysis includes previously published radiocarbon data from Crete and Antikythera (Pirazzoli et al., 1996 and references therein;Shaw et al., 2008;Wegmann, 2008;Tiberti et al., 2014;Mouslopoulou et al., 2015a), which were recalibrated by Mouslopoulou et al., (2015a) using the same curves and reservoir offsets applied to our samples.
To test whether radiocarbon samples can distinguish between earthquakes in the centuries before and after the AD 365 event, we generated synthetic radiocarbon posterior distributions of calendar ages for historical earthquakes (Fig. 2). We use calendar ages, together with the calibration curve Marine13 (Reimer et al., 2013), an atomic mass spectrometer uncertainty of ± 25 years, and a marine reservoir uncertainty of 58 ± 85 years (Reimer and McCormac, 2002) to calculate which distribution of calibrated ages one should expect from these uncertainties, if all organisms died due to uplift during a single earthquake event. This analysis assumes that the variability in the reservoir effect is representative of all collected samples. See acknowledgments for a link to the generated code . manuscript submitted to AGU Advances manuscript submitted to AGU Advances

Inverse fault dislocation modeling
To discriminate between competing hypotheses of forearc deformation, we invert the observed vertical displacement of the paleoshoreline (uplift) for fault rupture parameters (length, location, depth, dip, and slip) using a Bayesian Markov chain Monte Carlo routine (Hastings, 1970). We simulate two scenarios, one model with a single event on a reverse fault and a second with two earthquakes on adjacent normal faults (Fig. 1). The visco-elastic fault dislocation model of Wang et al. (2006) was used to constrain co-and postseismic deformation associated with a hypothetical earthquake. Two types of inputs are required; first, a crustal column with viscosities, seismic velocities, and densities for the respective layers (Fig. S2b), and second, fault parameters including depth, length, starting position, dip, slip and rake.
We used three different crustal models (25mant, 30LC, 40LC) based on different geophysical observations, assumptions, and previous dislocation studies. The crustal model "25mant" is based on the observation that seismicity in western Crete is limited to the upper 20 km of the crust , and the Moho is at ~ 32 km (Bohnhoff et al., 2001).
Therefore, we use a 25 km thick elastic layer on top of a 7 km thick lower crust, sitting on the asthenosphere. The crustal model "30LC" uses both a thicker elastic layer (30 km) and lower crustal material, which mimics the observations of Endrun et al. (2004), suggesting that there might be a thick layer of subduction channel material below western Crete. "40LC" assumes a 40 km thick elastic layer and is mainly used to compare results to Shaw et al. (2008), who employed a 45 km thick elastic layer, which is substantially deeper than observed earthquakes with foci located in the Aegean plate around Crete. Fault slip is assumed to be dip-slip because this is the only component we can resolve within the paleoshoreline uplift dataset. Our models are run for 200 years to capture the full postseismic deformation. Based on geophysical observations of upper crustal earthquake depths  and wide-angle seismic imaging (Bohnhoff et al., 2001), the 25 km thick elastic layer model (25mant) is the most realistic for western Crete and is thus presented in the main text.
The Metropolis-Hastings algorithm (Hastings, 1970) was used to determine the maximum a posteriori (MAP) model and populate the posterior distribution. We calculated 100,000 forward solutions for each of the three different crustal models (Fig. S2b) and manually tuned parameter transition kernels to achieve acceptance rates of 20 -40%. In modeling the normal faults, we assume uniform priors of 50° -80° for the dip, 5 -30 m for slip, and 5 -30 km for depth, while manuscript submitted to AGU Advances for the reverse fault models, we explore a dip range between 25° -55°, a slip of 20 -50 m and depth between 20 -50 km. The length of the two modeled normal faults is 30 to 47 km and 30 to 100 km for western and southwestern Crete, respectively, based on the dimensions of steep bathymetric offshore troughs (Fig. 1). For modeling of the reverse fault, we allow fault lengths between 80 and 200 km. The fault location and strike are predefined by the location and orientation of prominent bathymetric escarpments. The model's likelihood is evaluated based on the fit with interpolated paleoshoreline heights on Crete, Gavdos, Antikythera, and Kythera to remove spatial bias from the measurement locations. For the interpolation, we use the Krios Paleoshoreline heights measured in the field (Tab. 1, S1) and derive a continuous paleoshoreline height raster by kriging. This raster is then sampled in even increments (~ 4.5 km) along the coastline of Crete and neighboring islands.

Tsunami modeling
Historical records of tsunami inundation from Egypt (Alexandria), Sicily, and the Peloponnese for the AD 365 event (Ambraseys, 2009) provide an additional test on fault source kinematics. On a distant shoreline, the first appearance of the tsunami will be either a peak or trough dependent upon the kinematics of the fault rupture and relative position of the observer with respect to the source (Yamashita and Sato, 1974), analogous to the first arrivals of seismic waves, and can aid in the interpretation of vertical rupture kinematics. We simulate the impact of our fault rupture scenarios on tsunami generation in the Eastern Mediterranean (Fig. 4, S6).
We use EasyWave (Christgau et al., 2014) to generate tsunamis and compute their propagation.
EasyWave calculates wave propagation in the linear approximation of the long-wave theory in spherical coordinates (Christgau et al., 2014), which is only valid for water depths > 20 m. The bathymetry utilized here is from EMODNet (http://www.emodnet-bathymetry.eu) and downsampled to 1 km resolution. Using low-resolution bathymetry and a linear approximation for wave propagation is sufficient because we only model open ocean wave height, due to the lack of historic run-up data for comparison. The model takes the coseismic deformation field of our MAP models and propagates the tsunami from this initial displacement. manuscript submitted to AGU Advances

Radiocarbon data and historical records
The summed and normalized probability density function (PDF) of all calibrated radiocarbon-emergence ages (new and previously published) shows a distinct peak coinciding with the AD 365 event, yet there is considerable spread in the data (Fig. 2). The probability of all emergence ages gradually increases over 1,000 years before AD 365, after which it drops off more quickly. The positively skewed distribution suggests uplift in the centuries before AD 365 but could also record organism death before coseismic uplift (Fig. 2). The distribution of emergence ages against notch height at several sites provides evidence for sequential uplift locally but does not conclusively demonstrate evidence for or against a multi-stage uplift among all sites (Fig. S1). Several of our samples showed minor growth of secondary calcite during SEM analysis and ages younger than AD 365. However, our SEM analysis suggests that minor surficial alteration is effectively removed by sample processing and pre-treatment before radiocarbon analysis (Fig. S8). Nonetheless, secondary growth of carbonate would tend to result in younger ages. Therefore, it is reasonable to assume ages older than AD 365 represent the timing of organism death due to uplift or other means before AD 365.
Historical records document a series of large earthquakes, including one in AD 365, that affected Crete in the first centuries AD (Ambraseys, 2009;Papadopoulos, 2011). A large historically documented earthquake and tsunami occurred in AD 66 offshore western Crete (Papadopoulos, 2011). This event is recorded in the archeological record and caused damage along the west coast of Crete, including tsunami inundation. Importantly, archeological evidence indicates the uplift of an antique roman-age harbor in Phalasarna in AD 66, as harbor-basin sediments record marine deposition before and a terrestrial environment after the earthquake (Dominey-Howes et al., 1998;Stiros and Papageorgiou, 2001). Thus, at least locally, there is clear evidence that the Kris Paleoshoreline was uplifted in at least two events, rather than one, as traditionally assumed. manuscript submitted to AGU Advances  Papadopoulos (2011). The precise dates of some events are debated, while the historical texts used to compile these earthquake catalogs might include amalgamations of several events into one (for details, please see the two studies cited and references therein). To avoid some of these complications, we use the dates reported by Papadopoulos (2011) in our analysis.
We generate synthetic radiocarbon-emergence PDFs for the known and debated historical earthquakes for comparison with the radiocarbon-emergence data collected at and below the Krios Paleoshoreline (Fig. 2). The synthetic PDFs for the historical earthquakes consider that the sampled dataset will include variability from uncertainties in the atomic mass spectrometer measurement, the marine reservoir effect, and the radiocarbon calibration curve. The results show that several historical events fall into the range of the distribution of observed emergence ages (Fig. 2). This analysis suggests that one or more historical earthquakes preceding AD 365

Modeling competing earthquake scenarios
Based on the evidence for high slip rates (see Background), along with the radiocarbon data and historical and archeological records discussed above, we hypothesize that the Phalasarna Fault and Sfakia Fault Zone ruptured in at least two events in AD 66 and AD 365, uplifting the Krios Paleoshoreline to its approximate present-day elevation. Consistent with previous studies, we use the modern height of the Krios Paleoshoreline as a deformation marker of coseismic and postseismic uplift in the first centuries AD. However, sea-level change, regional uplift, and interseismic deformation might have affected the elevation of the paleoshoreline. We discuss the amplitude of those effects in the supplement (text S2 and S4), and test scenarios with corrections for sea-level rise and regional uplift but find a limited effect on the parameter distribution of the inversion (Fig. S3); therefore, we use the paleoshoreline height above present-day sea level as our deformation maker. In contrast to previous studies, we include 3.5 m of uplift we observed on the island of Gavdos, 35 km south of Crete (text S3) for the reverse fault scenario; however, simulations that exclude data from Gavdos have little impact on recovered parameters (Fig. S3). The normal faults investigated here do not produce substantial vertical land movement on Gavdos, Antikythera, and Kythera islands. Therefore, should the manuscript submitted to AGU Advances normal faults modeled here be capable of matching the Holocene uplift on Crete, uplift on adjacent islands is likely related to other faults not considered here.

manuscript submitted to AGU Advances
Using the Krios Paleoshoreline as a geodetic marker of deformation, we find that the rupture parameters from our inverse modeling of the reverse fault scenario are consistent with those from previous studies (Fig. 3) (Shaw et al., 2008;Mouslopoulou et al., 2015b), and lead to five critical findings: (1) fault depth, but also dip and slip amount, are dependent on the chosen crustal model because of changes in postseismic relaxation behavior that scale with the modeled elastic thicknesses (Fig. S2); (2) the fault must rupture the entire elastic layer (≥ 25 km), and many models rupture deep into the lower crust or mantle (Fig. S2); (3)   with the topography and bathymetry (Fig. 4). The profiles also highlight that, especially in the normal fault case, postseismic uplift is an important contributor to the observed height of the paleoshoreline (see figure S4 for a map view of co-, post-, and total seismic deformation in all models).

Tsunami modeling
Modeled tsunamis generated by normal faults offshore western and southwestern Crete reach all sites where credible historic reports document sea-level changes and mimic the spatial extent of simulated reverse-fault tsunami models (Fig. 5, S6). Synthetic mareograms generated by the best-fit reverse fault model have higher peak amplitudes than those generated by normal fault rupture. The higher peak amplitudes are due to the greater dip and slip required by a single reverse fault earthquake to match the paleoshoreline observational data. The best historical documentation of the AD 365 tsunami comes from Alexandria, where the Roman scribe Ammianus Marcellinus described the tsunami in detail, writing, "…solid earth was shaken and trembled, the sea with its rolling waves was driven back and withdrew from the land…many men roamed about without fear in the little that remained of the waters, to gather fish and similar things with their hands…For the great mass of waters, returning when it was least expected, killed many thousands of men by drowning." (Marcellinus,378AD). The tsunami polarity described by Marcellinus is only consistent with a normal fault earthquake offshore Crete (Fig   5b).
manuscript submitted to AGU Advances

A single or multiple events?
Was the Krios Paleoshoreline on Crete uplifted in a single event by up to 9 m? There are a number of damaging historical earthquakes that occur within the vicinity of western Crete in the centuries preceding the AD 365 event (Ambraseys, 2009;Papadopoulos, 2011). Despite the summed distribution of emergence data showing a peak coinciding with the AD 365 event (Fig.   2), the distribution is wide and skewed towards older ages, which may hint at partial uplift in the centuries before AD 365. Our synthetic distribution of radiocarbon-emergence ages predicted from historic earthquakes matches well with the observed emergence age distribution. However, the reliability of individual events is, in some cases, debated, and the spatial extent of damage manuscript submitted to AGU Advances and impact is difficult to evaluate (Ambraseys, 2009;Papadopoulos, 2011). Nevertheless, archaeological studies have found evidence for a series of destructive earthquake events in the ancient Cretan cities of Kissamos and Gortyn within the first centuries AD (Di Vita, 1985Vita, , 1999Stiros and Papageorgiou, 2001). Altogether, the radiocarbon data is not sufficient to differentiate between uplift in one or several events. However, the consistency of geochronologic data with abundant historical reports, the uplift of Phalasarna harbor in AD 66, and archeological destruction horizons all support the hypothesis of multi-stage coseismic uplift of the Krios Paleoshoreline. width data for a global compilation (Leonard, 2010) with best-fit and one sigma uncertainty lines.

Fault scalings and rupture parameters
(b) The same data for length vs. seismic moment. In this case, extensional and compressional faults show different scaling as highlighted by the linear fits in grey and black.
To assess the likelihood of the reverse fault and normal fault models, we compare best-fit fault parameters to empirical data for fault ruptures globally (Leonard, 2010). Both scenarios produce faults with length to width ratios that match global observations but generate higher seismic moment to fault length ratios relative to empirical data (Fig. 6). The higher 38 m reverse fault slip for our best-fit model (Fig. 3), compared to recent studies, is explained by a more realistic crustal set-up and incorporation of postseismic deformation. Previous studies have used thicker elastic layers that matched their fault depth (45 km) (Shaw et al., 2008) or only considered elastic deformation (Stiros, 2010;Mouslopoulou et al., 2015b) (Lay et al., 2010). Also, intraplate normal faulting earthquakes have been shown to reach a similar size to events proposed here, with several well-studied normal faulting earthquakes in the USA and China of Mw ~ 7.5 (Crone and Machette, 1984;Du et al., 1992;Suter, 2008;Middleton et al., 2016;Xu et al., 2018). For instance, the 1556 Huaxian, China earthquake produced 10 m of slip over a large portion of its ~ 90 km rupture in an Mw 7.5-8 earthquake (Yuan and Feng, 2010;Feng et al., 2020).
Another possibility to reconcile the high slip for the normal fault events with empirical scalings might be a more complex earthquake rupture. Recent observations from the 2016 Kaikoura and 2010 El Mayor-Cucapah earthquakes suggest complex fault interactions and manuscript submitted to AGU Advances ruptures with several faults at a time; this type of faulting behavior may be more common than previously thought (Fletcher et al., 2014;Hamling et al., 2017). As the results displayed in figure   6 indicate, both the normal and reverse fault models, when compared to global datasets, generate larger seismic moments as a function of modeled fault lengths in order to reproduce the uplift of the Krios Paleoshoreline, signaling that the true earthquake rupture parameters may be oversimplified in the model.

Tsunami Modeling
Our tsunami modeling suggests that only normal fault tsunamis fit the wave polarity described in Alexandria. A north dipping reverse fault offshore Crete also creates an initial negative arrival but with a much lower magnitude (~ 10 cm in our best-fit model, fig. 5).
However, one could argue that despite Ammianus Marcellinus's detail, a single report is not conclusive. Moreover, earthquake-triggered slumping offshore Alexandria might have generated a tsunami with a regressive polarity. Nevertheless, our findings and the Marcellinus report are consistent with reports from the well-documented 8 August AD 1303 earthquake (estimated ~ Mw 8) and tsunami with an assumed hypocenter offshore Crete and Rhodes (Guidoboni and Comastri, 1997;El-Sayed et al., 2000). The initial tsunami polarity described in AD 1303 in Alexandria was regressive, and the polarities reported from several locations around the Eastern Mediterranean basin were only consistent with the rupture of a south-dipping normal fault offshore southeastern Crete (El-Sayed et al., 2000). This event was subsequently often linked to another large reverse fault rupture and tsunami (England et al., 2015;Mouslopoulou et al., 2015b). However, the consistency of previous tsunami modeling and our results suggest that perhaps normal fault rupture is a common tsunami-trigger in the Eastern Mediterranean.

Evaluation of competing hypothesis
The interpretation of structural data and Neogene sediments indicate ongoing forearc extension since the Late Miocene (Angelier et al., 1982;ten Veen and Postma, 1999b;Gallen et al., 2014;Ott et al., 2019b) as well as ubiquitous onshore normal faulting, which is in better agreement with our proposed normal faulting origin of the AD 365 event. Normal faults are better supported by seismic and borehole data indicating the Hellenic Trough is a sediment underfilled half-graben (Lallemant et al., 1994;Bohnhoff et al., 2001). We also note that others have made interpretations of extension at the inner graben wall and transpressional features around the outer graben Chaumillon and Mascle, 1997).

manuscript submitted to AGU Advances
Previous interpretations in support of a contractional structure southwest of Crete are indirectly derived from several focal mechanisms (Bohnhoff et al., 2005;Shaw and Jackson, 2010), GPS stations between the forearc and the volcanic arc (Shaw et al., 2008;Saltogianni et al., 2020), and the assumption that the AD 365 event singularly uplifted the Late Holocene Cretan paleoshoreline. It is difficult to distinguish earthquake focal mechanisms generated on faults embedded in the forearc from those occurring on the subduction interface due to hypocenter depth uncertainties in the regions of thin crust south of Crete. However, focal mechanisms probably remain the best argument in support of a reverse fault scenario. GPS data show only minor contraction between the forearc and the volcanic arc (1-4 mm/a) (Saltogianni et al., 2020), yet due to the lack of offshore GPS data, this signal is indistinguishable from minor locking of the subduction interface (Vernant et al., 2014;Saltogianni et al., 2020).
These data and arguments collectively suggest that when the assumption of a single earthquake generating ~9 m of coseismic uplift on Crete is relaxed to allow for two or more events, a large splay reverse fault is not necessary. Due to the consistency of the normal faulting scenario with geologic and tsunami observations, more parsimonious rupture parameters, and the lack of imaging of the hypothesized reverse fault systems, we favor a normal faulting origin for the anomalous uplift and earthquake reports of the first centuries AD.

The importance of postseismic deformation for the Krios Paleoshoreline uplift
It has been recognized that uplift of footwall mountain ranges in extensional settings is in large parts due to postseismic deformation (Thompson and Parsons, 2017). These findings are in agreement with our modeling results from Crete (Fig. 4). Our model predicts that, most likely, a significant portion of the AD 365 uplift was postseismic, confirming studies that highlighted the importance of postseismic far-field uplift on Crete (Shaw et al., 2008). It is essential to note that we run our models for 200 years to capture the full postseismic deformation component. In the normal fault scenario, postseismic uplift is especially high in the southwestern corner of Crete, where the highest uplift was documented, due to the proximity of both normal faults in this area.
These findings illustrate that postseismic deformation may be an important contributor, along with coseismic and regional uplift, to the topography of the up to 2.5 km high coastal mountain ranges of western and southwestern Crete. manuscript submitted to AGU Advances

Implications for geodynamics, seismic hazard, and ways to resolve the debate
Our study proposes an alternative hypothesis to explain the historical observations and paleoshoreline data and reconciles these with long-term upper-crustal extension in the Hellenic forearc. If correct, the findings of this study suggest reduced earthquake and tsunamigenic hazards relative to the single-event splay thrust model, yet they still indicate that upper plate structures represent a significant hazard in the Eastern Mediterranean. If the Krios Paleoshoreline uplifted during multiple events within several centuries, it is likely that normal fault earthquakes in the Hellenic forearc are clustered in time, suggesting temporal variability in earthquake and tsunami hazard. This idea is not new; historical records and radiocarbon dating of uplifted Holocene shorelines indicate an Eastern Mediterranean earthquake cluster during the 4 th to 6 th century AD (Pirazzoli et al., 1996;Stiros, 2001). These findings also highlight the potential role of normal faults for generating strong earthquakes and tsunamis in subduction zone settings.
Nevertheless, future work is needed to conclusively assess the tectonic models and seismic and tsunami hazards of this densely populated region. Current, publicly available offshore seismic data is of poor quality, hence the reliance on onshore data to infer offshore fault kinematics. Improved three-dimensional seismic data would clarify the kinematics and activity of potential offshore seismic sources. Our model attempts to reconcile observed deformation kinematics on Neogene, Pleistocene, and Holocene time scales (see supplement for discussion), but additional InSAR and vertical GPS measurement datasets are needed to assess interseismic strain accumulation at decadal timescales. With these and existing datasets, a more robust picture of the earthquake-tsunami hazards of the Eastern Mediterranean and the geodynamics of the HSZ will emerge.

Summary and Conclusions
We propose an alternative hypothesis for the observations of the AD 365 earthquake and uplift of the Krios Paleoshoreline centered on the Phalasarna Fault and Sfakia Fault Zone offshore Crete. Radiocarbon data, historical reports, and archaeologic findings are insufficient to distinguish between single and multi-stage uplift of the Krios Paleoshoreline. However, the data are consistent with an earthquake cluster on normal faults proposed herein. Inverse modeling of fault parameters suggests that normal faults can generate uplift of the Krios Paleoshoreline with less slip, on faults that do not need to penetrate as deeply, as is required by an offshore reverse manuscript submitted to AGU Advances fault. Our models also highlight the importance of postseismic paleoshoreline uplift. Tsunami modeling of our best fit models shows that normal faults can generate strong tsunamis that impact the entire Eastern Mediterranean basin and are more consistent with the reported tsunami polarity. Based on these findings and the better consistency with the long-term record of crustal extension in the region, we favor a normal faulting origin for the AD 365 and earlier earthquakes. However, we note that more research and especially geophysical imaging, is required to adequately understand the tectonics and seismic hazard of the Hellenic Subduction Zone.