Reassessing Eastern Mediterranean Tectonics and Earthquake Hazard From the 365 CE Earthquake

The hallmark of great earthquakes in the Mediterranean is the 21 July 365 CE 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 M w 8.3–8.5 earthquake, uplifting a Cretan Holocene paleoshoreline by up to 9 m. Similar source mechanisms have been adapted for the region, which are commonly used for seismic and tsunami hazard estimation. We present an alternative model for Holocene paleoshoreline uplift and the 365 CE tsunami that centers on known active normal fault systems offshore of western and southwestern Crete. We use new and published radiocarbon dates and historical records to show that uplift of the Cretan paleoshoreline likely occurred during two or more earthquakes within 2–3 centuries. Visco‐elastic dislocation modeling demonstrates that the rupture of these normal faults fits observed data as well as reverse fault models but requires reduced slip and lower cumulative earthquake energy release (∼M w 7.9). 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 can originate on normal faults, highlighting the potential hazard from tsunamigenic upper plate normal fault earthquakes.

event recorded in the Mediterranean and has been the subject of intense research (e.g., Shaw et al., 2008;Stiros, 2001Stiros, , 2010. The event is particularly interesting because subduction zones in the Mediterranean appear to lack the great megathrust earthquakes that typify subduction zones elsewhere in the world (Jackson & 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 × 80 km area, place the 365 CE event epicenter offshore of southwestern Crete (Figure 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 CE (Figure 1) (Dominey-Howes et al., 1998;Shaw et al., 2008). Numerous studies argue that all of the uplift occurred in a single event in 365 CE, based on the correlation between radiocarbon ages and historical earthquake reports (e.g., Shaw et al., 2008;Stiros, 2010). 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 365 CE 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 ).
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 ( Figure 1). Using dislocation and tsunami modeling, Shaw et al. (2008) showed that the reverse fault interpretation would have generated an M w 8.3-8.5 earthquake and a strong and widely traveled 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 reverse faults and that all may be capable of generating >M w 8 earthquakes and strong tsunamis (England et al., 2015;Shaw & Jackson, 2010;Stiros, 2010). Thus, our current understanding of tectonics and seismic hazard of the densely populated Eastern Mediterranean is intimately related to our understanding of the 365 CE event.
While these results are provocative, at present, there is no direct evidence for active, near-surface reverse or thrust faults in the portion of the Hellenic forearc north of the Mediterranean Ridge accretionary complex and south of Crete. 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 & 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;. No active contractional structures have been definitively imaged offshore between Crete and the Mediterranean Ridge Accretionary complex (Chamot-Rooke et al., 2005;Lallemant et al., 1994). Seismic imaging and interpretation of International Ocean Drilling Program IODP drill cores indicate that the Hellenic Trough ( Figure 1) is an underfilled half-graben bounded by a southwest-dipping normal fault Lallemant et al., 1994). 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 .
Interpretations in support of a contractional structure southwest of Crete are indirectly derived from several focal mechanisms (Shaw & Jackson, 2010), GPS stations between the forearc and the volcanic arc (England et al., 2015;Shaw et al., 2008), and the assumption that the 365 CE event singularly uplifted the Krios paleoshoreline. However, GPS data show only minor contraction (1-4 mm/a) (Saltogianni et al., 2020;Shaw et al., 2008) between the forearc and the volcanic arc, yet due to the lack of offshore GPS data, this signal OTT ET AL. 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 365 CE event by providing an alternative model based on the rupture of known active offshore normal fault systems. Our analysis is guided by an attempt to reconcile the record of 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 and below the Krios paleoshoreline on Crete coincide with multiple 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 faults 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 (Figure 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 the fast retreat velocity of this subduction zone (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 (Figure 1). The crust below Crete is between 20 and 35 km thick Snopek et al., 2007) and has been subject to pervasive extension since the Mid-to-Late Miocene (Angelier et al., 1982;Brun et al., 2016;ten Veen & Postma, 1999). Accelerated retreat after 15 Ma is marked by horizontal extension and the formation of 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;. The variable strike of these sedimentary basins indicates horizontal extension both parallel and normal to the subduction zone (Fassoulas, 2001;ten Veen & Postma, 1999). 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 Robertson et al., 2019). Rock uplift is mostly driven by tectonics since erosion rates are only ∼0.1 mm/a, indicating that the component of isostatic adjustment due to erosional unloading is likely negligible (Ott, Gallen, Caves Rugenstein, et al., 2019). The uplift of Pleistocene paleoshorelines is interpreted to be related to a regional-scale uplift signal augmented by local uplift in the footwalls of large normal fault systems Robertson et al., 2019).
OTT ET AL.   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). The subduction thrust is black and bold, reverse faults red, strike-slip faults yellow, normal faults black, and faults with debated or unidentified kinematics blue. Seismic reflection data , focal mechanisms , and active synthetic onshore faults  indicate that two active normal faults lie offshore the western and southwestern coastlines of Crete ( Figure 1, Text S6). We refer to the normal fault offshore of 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 indicate long-term slip rates of 2.9-5.8 mm/a for the Sfakia fault zone . The Sfakia fault zone is also thought to be the offshore continuation of the north bounding fault of the onshore Messara Graben (Figure 1e) Peterek & Schwarze, 2004). Additionally, Late Pleistocene paleoshorelines in western Crete record some of the islands' highest uplift rates with a potential acceleration toward modern , 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 ( Figure 1, Table 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 (SEM) (Quanta 3D FEG) for signs of alteration (Text S1, Figure S8). The samples were calibrated using the Oxcal online calibration (Ramsey, 2009) program with the Marine13 calibration curve (Reimer et al., 2013) and an Eastern Mediterranean marine reservoir effect of 58 ± 85 years as recommended by Reimer and McCormac (2002). Our analysis includes previously published radiocarbon data from Crete and Antikythera (Pirazzoli et al., 1996 and references therein;Shaw et al., 2008;Tiberti et al., 2014;, which were recalibrated by  using the same curves and reservoir offsets applied to our new samples. To test whether radiocarbon samples can distinguish between earthquakes in the centuries before and after the 365 CE event, we generated synthetic radiocarbon posterior distributions of calendar ages for historical earthquakes ( Figure 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 & 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 (Ott, 2020).

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 ( Figure 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 ( Figure S2b), and second, fault parameters including depth, length, starting position, dip, slip, and rake. The coseismic deformation component is calculated following Okada (1992).
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 . 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) above lower crustal material, which mimics the observations of Endrun et al. (2004), suggesting that there might be a thick layer of subducted 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 OTT ET AL.    and wide-angle seismic imaging , 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 ( Figure 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 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-47 km and 30-100 km for western and southwestern Crete, respectively, based on the dimensions of steep bathymetric offshore troughs (Figure 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 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 (Table 1 and Table 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 365 CE 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 & 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 ( Figure 5, Figure S6).
We use EasyWave (Christgau et al., 2014) to generate tsunamis and compute their propagation. EasyWave calculates wave propagation using 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 runup data for comparison. The model takes the fault parameters from our MAP models, calculates the coseismic deformation field following Okada (1985), and propagates the tsunami from this initial displacement.

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 365 CE event, yet there is considerable spread in the data (Figure 2). The probability of all emergence ages gradually increases over 1,000 years before 365 CE, after which it drops off more quickly. The positively skewed distribution suggests uplift in the centuries before 365 CE but could also record organism death before coseismic uplift (Figure 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 ( Figure S1). Several of our samples showed minor growth of secondary calcite during SEM analysis. However, our SEM analysis suggests that minor surficial alteration is effectively removed by sample processing and pre-treatment before radiocarbon analysis ( Figure S8). Nonetheless, secondary growth of carbonate would tend to result in younger ages. Therefore, it is reasonable to assume ages older than 365 CE represent the timing of organism death due to uplift or other means before 365 CE.
Historical records document a series of large earthquakes, including one in 365 CE, that affected Crete in the first centuries CE (Ambraseys, 2009;Papadopoulos, 2011). A large historically documented earthquake and tsunami occurred in 66 CE 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 66 CE, as harbor-basin sediments record marine deposition before and a terrestrial environment after the earthquake (Dominey-Howes et al., 1998;Stiros & Papageorgiou, 2001). Thus, at least locally, there is clear evidence that the Krios paleoshoreline was uplifted in at least two events, rather than the conventional interpretation of a singular uplift event.
Other destructive earthquakes reported on Crete around this time are mentioned for the years BC 368 and 68, and 53 CE, 168, 408, and 620 by Ambraseys (2009) and186, and68 BCE and, 108 CE, 251, 448, 560, 618, and670 by Papadopoulos (2011), some of which are also linked to destruction at archeological sites. However, some of the aforementioned earthquakes are likely the same event with different assigned dates. For example, the 168 CE, 408, and 620 events reported by Ambraseys (2009) are duplicates of the 108 CE, 448, and 618 earthquakes noted by 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 ages collected at and below the Krios paleoshoreline ( Figure 2).
OTT ET AL. Individual probability density functions (PDF) of calibrated radiocarbon samples are added to a summed, normalized PDF (blue-shaded region). "Synthetic" PDF's that are expected for a set of emergence ages of marine organisms that died due to uplift during historically reported earthquakes on Crete are shown in dark red. Dashed lines indicate events with debated timing (Ambraseys, 2009;Papadopoulos, 2011). Several historic events fall into the range of emergence ages (blue). (b and c) Method to generate synthetic PDFs of historical earthquakes. (b) The historically documented calendar age of the earthquake is converted to a radiocarbon age PDF based on the calibration curve and offset by the reservoir correction (dR). (c) Analytical uncertainties from the AMS measurement and marine reservoir effect uncertainty are added to the radiocarbon PDF, subsequently used for the conventional Bayesian radiocarbon inversion to generate a posterior PDF of this calendar age. The synthetic PDFs for the historical earthquakes consider that the sampled data set will include independent variability associated with the uncertainties from 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 observed emergence ages (Figure 2). This analysis suggests that one or more historical earthquakes preceding 365 CE may have contributed to the uplift of the Krios paleoshoreline, consistent with the archeological evidence of uplift and a tsunami impacting the Phalasarna harbor in 66 CE (Dominey-Howes et al., 1998;Papadopoulos, 2011).

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 66 CE and 365 CE, 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 CE. 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 ( Figure 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 that 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 ( Figure S3). The normal faults investigated here do not produce substantial vertical land movement on Gavdos, Antikythera, and Kythera islands. However, from the radiocarbon dating, especially the beta-counting derived ages from Antikythera, it is not possible to conclusively link paleoshoreline emergence ages on the different islands. Therefore, should the 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.
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 ( Figure 3) Shaw et al., 2008), and lead to five critical findings: (1) fault depth, but also dip and slip amount, are dependent on the chosen rheological structure because of changes in postseismic relaxation behavior that scale with the modeled elastic thicknesses ( Figure S2); (2) the fault must rupture the entire elastic layer (≥25 km), and many models rupture deep into the lower crust or mantle ( Figure S2); (3) fault dip is inversely proportional to the elastic layer thickness, but all models require a steeply dipping fault, ≥35°; (4) fault slip must be 30-40 m; and (5) an M w ∼ 8.5 earthquake is required to generate the observed uplift.
The results from the two-event normal fault model fit the paleoshoreline elevation data on Crete as well as the reverse fault model (Figure 3). The results of this analysis indicate: (1) the two normal-sense earthquakes do not need to rupture the entire elastic layer to fit the data, making them less sensitive to the crustal models than the reverse fault scenario; (2) the lengths for both faults are around 45 km and the slip per fault is ∼16 m for our preferred crustal model; (3) the rupture depths range from 15 to 25 km; and (4) the total energy release within 2-3 centuries for both faults combined is M w 7.9 (M w 7.8 Phalasarna fault, M w 7.7 Sfakia fault zone). Fault orthogonal profiles show that the displacement patterns predicted by our models correlate well with the topography and bathymetry (Figure 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 OTT ET AL.    Figure S6). Synthetic mareograms generated by the best-fit reverse fault model have larger peak amplitudes than those generated by the best-fit 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 365 CE 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,378 CE). The tsunami polarity described by Marcellinus is only consistent with a normal fault earthquake offshore Crete (Figure 5b).

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 365 CE event (Ambraseys, 2009;Papadopoulos, 2011). Despite the summed distribution of emergence data showing a peak coinciding with the 365 CE event (Figure 2), the distribution is wide and skewed toward older ages, which may hint at partial uplift in the centuries before 365 CE. 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 and impact is difficult to evaluate (Ambraseys, 2009;Papadopoulos, 2011). Nevertheless, archeological studies have found evidence for a series of destructive earthquake events in the ancient Cretan cities of Kissamos and Gortyn within the first centuries CE (Di Vita, 1985Vita, , 1999Stiros & Papageorgiou, 2001). Altogether, the radiocarbon data are not sufficient to differentiate between uplift in one or several events. However, the consistency of geochronologic data with abundant historical reports, OTT ET AL.

Fault Scalings and Rupture Parameters
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 ( Figure 6). The higher 38 m reverse fault slip for our best-fit model (Figure 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); both assumptions are likely unrealistic on this time scale and may have underestimated the required reverse fault rupture parameters. Meier et al. (2004) showed that earthquakes below western Crete only occur in the upper 20 km; therefore, our 25 km elastic thickness of the crust is a conservative estimate. The amount of slip for the normal-fault scenario is large but substantially lower on individual faults than for the reverse-fault model. Importantly, our model assumes only two events generate uplift of the Krios paleoshoreline. However, uplift likely reflects the cumulative deformation associated with an earthquake cluster, where adjacent fault systems influence each other, perhaps analogous to normal fault earthquake clusters in the central Apennines (Chiarabba et al., 2011). If correct, additional earthquakes would lower the length-to-seismic moment ratios and match empirical relationships of global data while still matching the observed deformation of the Krios paleoshoreline ( Figure 6). Also, the rupture length, especially for the normal fault case, might be underestimated. The Phalasarna fault and Sfakia fault zone extend substantially farther than the 45 km rupture length that the inversion suggests. The shorter rupture length in the inversion is explained by the simplicity of the uniform slip model. A more realistic non-uniform slip model would be able to fit the same paleoshoreline uplift with longer faults and slip that decreases toward the tips.
The rupture parameters proposed here are high relative to those reported in compilations from historical and paleoseismic records, but they are within the range of observed parameters for large normal faulting earthquakes, especially when assuming an earthquake cluster instead of two events. The strongest normal faulting earthquakes are located on the outer trench slopes of subduction zones, with four earthquakes of M w > 8 in the past 100 years and ruptures of >100 km length with >10 m average slip (Lay et al., 2010).
OTT ET AL.
10.1029/2020AV000315 12 of 18  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 M w ∼7.5 (Crone & Machette, 1984;Du et al., 1992;Middleton et al., 2016;Suter, 2008;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 M w 7.5-8 earthquake (Feng et al., 2020;Yuan & Feng, 2010).
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 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 data sets, 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, Figure 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 1303 CE earthquake (estimated ∼M w 8) and tsunami with an assumed hypocenter offshore Crete and Rhodes (El-Sayed et al., 2000;Guidoboni & Comastri, 1997). The initial tsunami polarity described in 1303 CE 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;. However, previous tsunami modeling consistency and our results suggest that perhaps forearc 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 & Postma, 1999) as well as ubiquitous onshore normal faulting, which is in better agreement with our proposed normal faulting origin of the 365 CE 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). We also note that others have made interpretations of extension at the inner graben wall and transpressional features around the outer graben (Chaumillon & Mascle, 1997;. Previous interpretations in support of a contractional structure southwest of Crete are indirectly derived from several focal mechanisms Shaw & Jackson, 2010), GPS stations between the forearc and the volcanic arc (Saltogianni et al., 2020;Shaw et al., 2008), and the assumption that the 365 CE event singularly uplifted the Late Holocene Krios paleoshoreline on Crete. 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 (Saltogianni et al., 2020;Vernant et al., 2014).
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 unnecessary. 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 CE.

The Importance of Postseismic Deformation for the Krios Paleoshoreline Uplift
It has been recognized that the uplift of footwall mountain ranges in extensional settings is in large parts due to postseismic deformation (Thompson & Parsons, 2017). These findings are in agreement with our modeling results from Crete (Figure 4). Our model predicts that, most likely, a significant portion of the 365 CE 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 have 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.

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, this study's findings 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 likely 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 fourth to sixth century CE 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 are 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 data sets are needed to assess interseismic strain accumulation at decadal timescales. With these and existing data sets, 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 365 CE earthquake and uplift of the Krios paleoshoreline centered on the Phalasarna fault and Sfakia fault zone offshore Crete. Radiocarbon data, historical reports, and archeologic 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 fault. The normal fault model predicts a combined energy release from normal-sense earthquakes within the first centuries AD of ∼M w 7.9. 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 365 CE 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.