The 2018 M_w 6.8 Zakynthos (Ionian Sea, Greece) earthquake: seismic source and local tsunami characterization

SUMMARY

We investigated the kinematic rupture model of the 2018 M_w 6.8 Zakynthos, Ionian Sea (Greece), earthquake by using a non-linear joint inversion of strong motion data, high-rate GPS time-series and static coseismic GPS displacements. We also tested inversion results against tide-gauge recordings of the small tsunami generated in the Ionian Sea. In order to constrain the fault geometry, we performed several preliminary kinematic inversions by assuming the parameter values resulting from different published moment tensor solutions. The lowest cost function values were obtained by using the geometry derived from the United States Geological Survey (USGS) focal solution. Between the two conjugate USGS planes, the rupture model which better fits the data is the one with the N9°E-striking 39°ESE-dipping plane. The rupture history of this model is characterized by a bilateral propagation, featuring two asperities; a main slip patch extending between 14 and 28 km in depth, 9 km northeast from the nucleation and a slightly shallower small patch located 27 km southwest from the nucleation. The maximum energy release occurs between 8 and 12 s, when both patches are breaking simultaneously. The maximum slip is 1.8 m and the total seismic moment is 2.4 × 10¹⁹ Nm, corresponding to a M_w value of 6.8. The slip angle shows a dominant right-lateral strike-slip mechanism, with a minor reverse component that increases on the deeper region of the fault. This result, in addition to the observed possibility of similar mechanisms for previous earthquakes occurred in 1959 and 1997, suggests that the tectonic deformation between the Cephalonia Transform Fault Zone and the northern tip of the Hellenic Arc Subduction zone may be accommodated by prevailing right lateral low-dipping faults, occurring on re-activated structures previously experiencing (until Pliocene) compressional regime. Comparison of predicted and observed tsunami data suggests the need of a better characterization of local harbour response for this type of relatively short-wavelength events, which is important in the context of tsunami early warning. However, the suggested dominantly strike-slip character would in turn imply a reduced tsunami hazard as compared to a dominant thrust faulting regime from this source region.

Tsunamis, Joint inversion, Earthquake ground motions, Earthquake source observations, Seismicity and tectonics

1 INTRODUCTION

A M_w 6.8 earthquake struck the area of the Zakynthos island, western Greece (Fig. 1) on 25 October 2018, at 22:54 UTC. This earthquake occurred ∼36 km to the SW of Zakynthos and ∼40 km to the NW of the Strofades island. Limited structural damages were reported mainly on the dock of the Zakynthos harbour and at the Strofades monastery (Karakostas et al.2018). Tsunami alert messages, based on the earthquake parameters, were issued within ten minutes after the event by both the Italian and the Greek Tsunami Service Providers. The earthquake generated a small tsunami recorded by some tide-gauges, including those of Katakolo and Kyparissia in Greece and Crotone and Le Castella in Italy (Fig. 1).

Figure 1.

Location map of the 2018 Zakynthos earthquake. Black triangles and inverted triangles represent the strong motions and the HRGPS stations, respectively. Cyan circles and blue squares (in the inset) represent the GPS sites and the tide gauges, respectively. The red star indicates the epicentre. Black boxes represent the surface projection of the fault planes corresponding to the NP1 and NP2 USGS focal mechanism solutions. Grey, pink, yellow and green stars, display all the aftershocks, those with magnitudes between 4 and 5, those with magnitude greater than 5 and the foreshock, respectively. Events with M > 4 have been relocated in this study (Table 2). Focal mechanisms of the 1959 and 1997 earthquakes are also displayed. Red dashed line shows the approximate position of the east-dipping Cephalonia Transform fault zone (CTFZ); black line in the inset displays the Hellenic Arc subduction zone [HASZ, Kreemer & Chamot-Rooke (2004)].

Zakynthos is located at the transition between the northwestern tip of the Hellenic Arc subduction zone (HASZ), mainly characterized by low-dipping reverse fault mechanisms (Le Pichon & Angelier 1979; Anderson & Jackson 1987; Papazachos 1990; Papazachos et al.1991; Shaw & Jackson 2010), and the southern tip of the Cephalonia Transform fault zone (CTFZ, Fig. 1), dominated by right-lateral shear on subvertical faults within the upper crust (Scordilis et al.1985; Louvari et al.1999; Sachpazi et al.2000; Kokinou et al.2006). The transition between these two tectonic regimes is highlighted, among other evidences, by a variation in amplitude of the southwest warding velocities, at Global Positioning System (GPS) stations, ranging from ∼30 mm yr^–1 (with respect to stable Europe) close to the HASZ to ∼10 mm yr^–1 or less, nearby the CTFZ (Floyd et al.2010; Pérouse et al.2012).

A detailed description of the earthquake kinematics and of its causative fault is then of great importance to improve the knowledge of this active tectonic region. This, in turn, allows to better constrain the associated seismic and tsunami hazard. Active tectonic region characterization is somehow limited due to the existence of some earthquakes of similar magnitude in the same zone in the last 60 yr (Fig. 1), namely the 1959 (M_w 6.8) and the 1997 (M_w 6.6) events (e.g. Kiratzi & Louvari 2003), for which the ambiguity of the nodal plane has not yet been fully resolved.

In this work, we use strong motion data, high-rate Global Positioning System (HRGPS) waveforms and static coseismic GPS displacements in a non-linear joint inversion, to discriminate between the two nodal fault planes and to retrieve the kinematic rupture process of the 2018 Zakynthos earthquake. We also perform a forward simulation of the tsunami generated by our kinematic rupture model and compare it with the observations collected at the four above-mentioned tide-gauges, to preliminarily assess the impact of the model uncertainty on the variability of the predicted tsunami waveforms. Finally, we briefly discuss the causative fault of this earthquake within its tectonic framework.

2 DATA

2.1 Seismicity relocation and moment tensor solutions

The earthquake was recorded by the National Observatory of Athens (NOA) and by the Institute of Engineering Seismology and Earthquake Engineering (ITSAK) seismic stations, among others. The hypocentre (37.38°N, 20.58°E, 11.5 km depth, Fig. 1) was estimated with the NonLinLoc software (Lomax et al.2000).

The moment tensor solution published by USGS indicates either a dominant right-lateral slip on a low dip angle fault plane, with a small thrust component, or a mainly reverse motion on a subvertical fault plane (see Table 1). Similar solutions have been proposed by GCMT (Global Centroid Moment Tensor), QRCMT (Quick Regional Moment Tensor), NOA and GeoForschungsZentrum (GFZ). Some of these solutions show significant non-double couple components, which may be due to the structural Earth's heterogeneity and source complexity. Other earthquakes characterized by a dominant strike-slip mechanism and similar magnitude (Fig. 1) are reported in the same area by Kiratzi & Louvari (2003) for the 1959 (M_w 6.8) and the 1997 (M_w 6.6) earthquakes, whereas lower magnitude dip-slip events are reported slightly to the southwest of the 2018 main shock (Papadimitriou 1993; Kiratzi & Louvari 2003).

Table 1.

Published moment tensor solutions for the 2018 Zakynthos, Ionian Sea, earthquake; last column contains the minimum cost function values of the associated joint inversion (namely ‘joint’) and reached in each inversions, performed in this study, by the single kind of data (‘smH’ stand for dynamic data (strong motion and HRGPS), and ‘gps’ for static data).

		NP2				NP1
Agency	ζ (Strike)	δ (Dip)	λ (Rake)	%DC	ζ (Strike)	δ (Dip)	λ (Rake)	Cost function
USGS	109°	82°	52°	76 per cent	9°	39°	167°	0.133 (joint)
								0.216 (smH)
								0.051 (gps)
NOA	108°	85°	41°	39 per cent	14°	49°	174°	0.405 (joint)
								0.365 (smH)
								0.445 (gps)
GCMT	114°	83°	63°		11°	28°	165°	0.512 (joint)
								0.497 (sm)
								0.527 (gps)
GFZ	107°	85°	68°		5°	23°	167°	0.632 (joint)
								0.505 (smH)
								0.759 (gps)
QRCT	117°	85°	63°		17°	27°	168°	0.714 (joint)
								0.617 (smH)
								0.811 (gps)

The main shock was preceded, 32 minutes earlier, by a M_L 4.9 foreshock, and followed after 15 min by a M_L 5.1 aftershock. The two largest aftershocks of M_L 5.4 and M_L 5.5 occurred on October 30 at 03:00 and 15:12 UTC, respectively. The NOA (NOA 2019) catalogue contains also 1705 aftershocks of magnitudes M > 2 (black points in Fig. 1) occurred in the two months following the main shock. Events with magnitude M > 4 have been relocated (64 events, see Table 2), using the same procedure adopted to locate the main shock. The aftershock locations (Fig. 1) show a diffuse seismicity, mainly concentrated in the uppermost 20 km (Fig. S1, in Supporting Information), and do not univocally depict a fault plane. This might be partly due to the inaccuracy of the locations due to the limited station coverage.

Table 2.

Relocated seismicity, for events with magnitude M > 4, given in terms of origin time, latitude, longitude and depth. ‘No’ is the number of P and S phases; ‘DM’ the minimum station distance; ‘GAP’ is the azimuthal gap of the seismic stations; ‘RMS’ is the root-mean-square of traveltime residuals; ‘ERH’ is the standard error in the epicentre estimate, and ‘ERZ’ is the standard error in the focal depth estimate.

Date	Origin time	Lat (°)	Lon (°)	Z (km)	M	No	DM	GAP	RMS	ERH (km)	ERZ (km)
18/10/17	22:03:17.09	37.3707	20.6017	9.95	4.2	25	82	248	0.30	2.5	2.5
18/10/19	16:51:56.11	37.4115	20.6455	12.88	4.3	35	37	225	0.24	0.7	0.7
18/10/25	22:22:54.38	37.3763	20.5800	11.29	4.9	42	43	222	0.27	0.6	0.6
18/10/25	22:54:50.09	37.3768	20.5805	11.54	6.6	41	38	226	0.53	2.1	2.0
18/10/25	23:09:21.36	37.1670	20.7065	13.73	5.1	91	61	208	0.58	1.3	1.2
18/10/25	23:17:24.86	37.3538	20.7650	10.12	4.2	45	40	216	0.37	1.2	1.1
18/10/26	00:13:39.03	37.4687	20.6597	7.38	4.5	60	31	179	0.46	1.3	1.0
18/10/26	00:23:13.57	37.3973	20.9408	7.71	4.3	50	36	204	0.54	1.3	0.9
18/10/26	00:32:54.50	37.6877	20.3670	4.35	4.4	48	41	203	0.32	1.0	1.0
18/10/26	01:06:4.080	37.4120	20.8837	8.04	4.5	29	33	206	0.38	2.0	2.1
18/10/26	02:17:34.26	37.5275	20.5832	16.79	4.1	47	30	218	0.49	2.4	2.0
18/10/26	05:48:37.26	37.4035	20.5653	12.82	4.8	73	42	175	0.39	0.8	0.5
18/10/26	06:32:13.84	37.4500	20.7193	9.46	4.2	46	31	208	0.34	1.3	1.1
18/10/26	06:44:8.420	37.4687	20.5265	8.01	4.2	45	38	221	0.23	0.8	0.8
18/10/26	12:11:16.82	37.4663	20.7328	17.04	4.4	34	28	215	0.25	0.9	0.7
18/10/26	12:41:13.77	37.4030	20.5673	13.23	4.6	33	41	226	0.39	2.5	2.0
18/10/26	16:07:9.910	37.4607	20.6240	14.52	4.5	35	33	224	0.25	0.8	0.6
18/10/27	14:33:27.27	37.4843	20.5148	12.20	4.1	44	38	231	0.39	1.3	1.0
18/10/27	10:13:43.21	37.4087	20.6355	12.55	4.1	36	38	221	0.18	0.8	0.5
18/10/27	05:28:46.93	37.5070	20.6735	13.54	4.6	50	27	194	0.33	1.0	0.8
18/10/27	00:05:36.36	37.5768	20.7525	13.97	4.1	43	16	205	0.40	1.0	0.9
18/10/28	20:40:21.47	37.3318	20.7672	11.87	4.1	48	42	222	0.41	1.3	0.9
18/10/29	04:52:12.54	37.5333	20.6302	14.03	4.1	48	27	218	0.38	1.0	0.6
18/10/29	11:29:34.54	37.6547	20.4298	15.70	4.1	25	36	230	0.23	1.2	1.0
18/10/29	15:01:40.53	37.3662	20.6035	14.79	4.3	42	43	228	0.32	1.1	0.8
18/10/29	22:22:14.58	37.4430	20.5580	10.17	4.1	28	38	230	0.29	2.1	2.4
18/10/30	02:59:59.80	37.5993	20.5248	10.04	5.4	49	30	222	0.33	1.1	1.5
18/10/30	06:34:13.51	37.6512	20.5208	12.61	4.2	33	29	221	0.23	0.8	0.5
18/10/30	08:32:25.79	37.4568	20.4275	10.36	4.8	29	46	233	0.24	1.2	1.5
18/10/30	12:49:07.11	37.5592	20.6038	16.35	4.4	66	26	217	0.54	1.5	0.8
18/10/30	14:33:16.72	37.4503	20.4648	12.28	4.1	38	44	229	0.32	1.1	1.0
18/10/30	15:12:03.06	37.5023	20.5127	15.80	5.5	65	37	215	0.60	1.3	1.0
18/10/30	16:39:21.90	37.4480	20.5037	12.71	4.3	40	41	234	0.28	0.8	0.8
18/10/30	18:04:24.40	37.4750	20.5040	14.79	4.3	43	39	226	0.26	0.7	0.6
18/10/31	10:25:12.13	37.4167	20.8413	15.26	4.1	28	32	210	0.44	1.5	2.0
18/11/01	02:44:48.25	37.3818	20.5882	10.34	4.6	61	42	182	0.47	1.3	1.1
18/11/01	05:34:31.16	37.2115	20.6600	8.91	4.3	36	57	238	0.32	1.2	1.3
18/11/02	07:53:14.62	37.6062	20.4360	18.68	4.3	49	37	232	0.35	1.0	0.7
18/11/03	01:10:06.69	37.2322	20.6348	9.43	4.3	47	56	227	0.29	0.8	1.0
18/11/03	08:05:53.83	37.6093	20.3245	15.83	4.3	44	46	208	0.28	0.8	0.7
18/11/04	03:04:31.05	37.1965	20.6205	15.56	4.3	46	60	225	0.30	1.1	1.1
18/11/04	03:12:45.06	37.4052	20.4535	12.63	4.9	55	48	184	0.37	1.2	1.2
18/11/05	02:44:38.22	37.1868	20.6333	9.30	4.4	29	61	241	0.33	2.2	2.4
18/11/05	06:46:13.14	37.6572	20.5322	15.61	4.5	45	27	217	0.24	0.7	0.6
18/11/05	08:31:12.59	37.5712	20.7198	16.38	4.3	47	18	208	0.30	0.8	0.5
18/11/08	22:46:00.54	37.6102	20.4968	10.01	4.2	39	32	226	0.36	1.2	1.1
18/11/10	02:13:38.66	37.6855	20.5332	14.82	4.2	40	27	218	0.37	1.3	1.1
18/11/11	23:38:35.55	37.6610	20.5397	7.79	4.8	50	27	219	0.27	0.6	0.8
18/11/12	06:50:28.87	37.2000	20.6173	5.74	4.7	44	60	229	0.49	1.2	1.0
18/11/15	09:02:05.38	37.5277	20.6877	15.89	4.9	40	24	214	0.26	0.7	0.5
18/11/15	09:09:26.80	37.5207	20.6913	13.04	4.5	31	25	214	0.23	1.2	1.1
18/11/15	11:00:05.14	37.6700	20.5323	14.19	4.4	38	25	219	0.27	1.1	0.9
18/11/18	05:18:03.30	37.4417	20.5348	8.39	4.1	53	40	227	0.40	1.0	0.6
18/11/18	06:06:45.89	37.6030	20.3645	6.04	4.3	53	43	226	0.36	1.0	0.8
18/11/19	05:56:51.26	37.5418	20.6730	13.75	4.1	40	24	215	0.28	0.5	1.0
18/11/19	13:05:55.99	37.2023	20.5800	2.46	5.1	40	61	243	0.37	1.1	1.1
18/11/22	11:11:29.95	37.2013	20.4978	0.05	4.4	34	64	241	0.39	1.5	0.6
18/11/22	16:00:16.36	37.5673	20.3860	2.79	4.1	41	42	237	0.37	1.1	1.2
18/11/29	00:23:00.37	37.6528	20.3007	6.02	4.2	26	41	243	0.21	2.0	3.1
18/12/13	06:26:41.95	37.5502	20.6838	14.60	4.4	54	22	176	0.37	1.2	1.1
18/12/25	01:41:28.21	37.3585	20.8270	11.40	4.6	49	39	216	0.31	1.0	0.8
18/12/26	07:37:53.07	37.3858	20.7900	10.20	4.1	31	36	216	0.23	0.6	0.7
19/01/10	17:30:15.65	37.3215	20.6245	6.95	4.4	33	47	233	0.35	1.1	1.3
19/01/15	01:11:49.20	38.3010	20.4393	0.03	4.2	46	9	121	0.37	1.1	0.4

2.2 Seismic, GPS, and tsunami data

Due to its relatively far offshore location, the main shock was recorded by a limited number of instruments in the near field, that is within ∼80 km from the epicentre three 200-Hz-sampling strong motion stations located in Zakynthos island (LTHK, ZAK2, KRI1), one 10-Hz-sampling HRGPS located in Strofades island (STRF) and nine 30-s-sampling GPS stations (GPS) stations (Fig. 1).

The strong motion stations show maximum peak ground acceleration (PGA) values of ∼0.39 g  on the E–W component at KRI1 and  ∼0.29 g on the N–S component at LTHK (Fig. 2a). For use in the inversion, a bandpass filter in the range 0.02–0.5 Hz was applied and the signals were integrated, to obtain velocity waveforms (Fig. 3).

Figure 2.

Observed data set. (a) Time histories of the observed near field accelerations; (b) and of the high sample rate displacements at the available strong motion and HRGPS sites, respectively (Fig. 1). (c) Horizontal (blue vectors) and vertical (pink vectors) coseismic displacements observed with the associated errors (blue and black dashed ellipsis and bars, respectively) at GPS sites located in the area of interest. (d) Observed tsunami waveforms.

Figure 3.

Comparison between recorded (blue) and predicted waveforms (amplitude in cm/sec) and spectra (amplitude in cm) for each sites’ motion component (panels a) (b), and (c), respectively) and GPS data (panel d) for NP1 (red) and NP2 (green) fault geometry.

The HRGPS waveforms for the north–south, east–west and up–down components were obtained by processing the data in kinematic mode, following the strategy proposed by Avallone et al. (2017). The deformation started ∼13 s after the origin time (Fig. 2b). The ground motion at STRF experienced two main peaks reaching maximum absolute values of ∼15.7 and ∼10 cm in the E–W component, respectively, and two peaks with maximum absolute values of ∼13.5 cm along with the N–S component (Fig. 2b). Two significant peaks on the U–D component (up to ∼9.14 cm) were also observed, before the signal reached its permanent offset. The HRGPS time series in displacement were differentiated to obtain velocity waveforms (Fig. 3).

The static GPS displacements were obtained by comparing 60 d of pre-event coordinates with 60 d of post-event coordinates, calculated by means of Gipsy-Oasis II software (v. 6.4, Bertiger et al.2010, Fig. 2c). STRF is displaced by ∼51.5 mm towards SE, whereas ZAKY and ZAKU are displaced by ∼48 mm towards SW. Significant permanent deformation is observed at AMAL and PYRG (19 and 18 mm, respectively) in the western Peloponnese. GPS sites, which are farther away, do not show appreciable coseismic displacement; nevertheless, they are important to constrain the geodetic seismic moment.

A minor tsunami was recorded by four tide-gauges located at Katakolo and Kyparissia in Greece, and Le Castella and Crotone in Italy (Fig. 1), all with a sampling rate of 1 min. The raw tsunami waveforms (Fig. 2d) were processed by removing the tidal component through a LOWESS local polynomial regression method (e.g. Romano et al.2016). The highest amplitudes are observed at Katakolo (11 cm, ∼70 km to the NE from the epicentre) and Kyparissia (10 cm, ∼100 km to the SE from the epicentre). Lower amplitudes were measured at Le Castella (6 cm) and Crotone (8 cm), both located ∼360 km to the NW from the epicentre.

3 RESULTS

3.1 Kinematic rupture model

The rupture history of the 2018 Zakynthos earthquake was obtained by jointly inverting the geodetic and seismic data through a two-stage non-linear inversion procedure (Piatanesi et al.2007; Romano et al.2010; Cirella et al.2012). In the first stage, a heat-bath simulated annealing algorithm builds up the model ensemble by collecting all models and their cost function values. During the second stage, the algorithm performs a statistical analysis of the ensemble, providing the best-fitting model along with fitting statistics, including the average model and marginal distributions of the parameters. Our final kinematic rupture model is obtained by averaging a subset of the whole explored model ensemble (Ensemble Subset, ES hereinafter), including only the ones that provide a sufficiently good fit to the observations.

We inverted the first 50 s of the velocity waveforms. The Green's functions are computed with a discrete wavenumber/finite element technique (Spudich & Xu 2003), and by adopting a local 1-D crustal velocity model (Haslinger et al.1999; Fig. S2, Table S1).

The fault plane is discretized following the USGS solution (Table S2). Four parameters characterizing the rupture at each node (spaced by 4.5 km in both strike and dip directions): peak slip velocity, rise time, rake angle and rupture time. The slip distribution is obtained by integrating the slip velocity over the rise time.

We set a priori bounds variability for each parameter (details in Table S3). Those range of variability are constrained by looking at both empirical scaling relations (e.g. Leonard 2014), and studies on source dynamics (Bernard et al.1996; Mai et al.2005; Custodio & Archuleta 2007; Cultrera et al.2010).

During the inversion, all parameters are simultaneously inverted at all grid nodes. The rupture time at each grid node is constrained by the arrival time from the hypocentre of a rupture front whose speed values are allowed to vary in the range given in Table S3. To avoid unphysical conditions, all models featuring a-causal local rupture velocity larger than P-wave velocity are discarded by the algorithm. The local velocity is computed from the gradient of rupture times, following the procedure described by Cirella et al. (2012, in particular their eqs 1–3).

We performed several preliminary inversions by using, as input values, the nodal planes of all the moment tensor solutions reported in Table 1, and as nucleation point the hypocentre location determined in the previous section, to find the best fault plane location (strike and dip) that we will use afterward for the inversion of the rupture history.

The best agreement between observed and modelled data is reached by using either of the two conjugate planes of the USGS solution (hereinafter, NP1 and NP2), as shown by the cost function values (Table 1).

The fault plane ambiguity inherent in the moment tensor solution is also clearly resolved by the inversion. The cost function of the NP1 fault plane is 47 per cent lower than that of the conjugated NP2 (Table S4). This is also visually evident by comparing observations with synthetics, both in time and in frequency domains, for each motion component and at each station (Figs 3a–c). NP1 is also strongly supported by the difference between the two corresponding forward predictions of the GPS static offsets at STRF (Fig. 3d). We then adopted this fault plane for the next step.

The average model (Fig. 4) was obtained by averaging the ES of nearly 250k ‘good’ models. ES includes all models with a cost function value exceeding by less than 2.5 per cent the absolute minimum cost function value reached during the inversion. This is our preferred rupture model (Figs 4a and b, Table S2), characterized by two asperities: a main slip patch located 9 km northeast from the nucleation and extending between 14 and 28 km depths; and a smaller and shallower slip patch located 27 km southwest from the nucleation. The main asperity features a maximum slip value of 1.8 m, associated to a rise time value of 3.4 s. This is in agreement with the finite fault model proposed by USGS, estimated through a teleseismic inversion. However, the latter is characterized by one asperity around their epicentre, and located about 20 km westward with respect to our main asperity. Our smaller patch is characterized by lower slip and rise time values (∼1.0 m and 2.6 s, respectively).

Figure 4.

Inverted rupture model (average model from ensemble inference) for NP1 focal mechanism solution. (a) Upper and lower panels show the total slip and the rise time distributions on the fault system, respectively. Rupture time shown by black contour lines (each 1 s); black arrows displayed in upper plot represent the slip vector. Red star displays the hypocentre. (b) Preferred total slip distribution, displayed in panel (a), projected on the Earth surface. (c) Snapshots at different times of slip velocity on the fault plane. White dashed lines show the regions that slipped more than 0.2 m. Grey contours indicate the rupture times (each 1 s). Red star displays the hypocentre.

On the main patch, the rake angles at the nodes show a dominant right-lateral mechanism, with minor reverse component that increases on the deeper part. The overall mechanism is consistent with the moment tensor solutions (Table 1). The rupture has a total duration of ∼14 s, reaching a maximum local velocity of 3.5 km s^–1 . Slip velocity time snapshots (Fig. 4c) highlight clear bilateral rupture propagation. The rupture starts breaking the main NNE slip patch after 4 s from the origin time, involving a total duration of 8 s; while the smaller SSW asperity fails about 4 s later, between 8 and 14 s. Therefore, the maximum energy release occurs between 8 and 12 s, when both patches are breaking simultaneously. The retrieved rupture evolution agrees well with the moment rate function published by USGS. The seismic moment inferred by our model is M₀ = 2.4 × 10¹⁹ Nm, corresponding to a magnitude M_w 6.8. For completeness, we show in Fig. S3 the rupture model obtained by adopting the NP2 fault geometry; the slip is smeared over a large portion of the fault, with a significantly increased value at very shallow depth.

3.2 Tsunami analysis and modelling

The amplitude and the general features of the small signals observed at the four selected Greek and Italian tide-gauges around the Ionian Sea (Fig. 2d) are qualitatively consistent with those of a tsunami generated by the coseismic seafloor displacement for an earthquake of this magnitude and focal mechanism. This is shown by preliminary simulations performed by CAT-INGV based on the readily available moment tensor solutions from USGS.

To perform an independent test against data not used in the inversion, we numerically modelled 100 synthetic tsunamis, generated by a further uniform sampling of the 250k good models belonging to the ES (see Section 3.1).

Seafloor vertical coseismic displacements from the 100 static slip models were computed with the analytical Okada formulas (Okada 1992), that is considering the final slip distribution at the end of the kinematic rupture process. The contribution of the horizontal coseismic deformation projected onto the oceanic slope is also modelled (Tanioka & Satake 1996). The resulting seafloor vertical displacement was then low-pass filtered through a function of shape 1/cosh(kh) (where k is the wave number and h is the water depth) to model the attenuation through the water column (Kajiura 1963), finally obtaining the vertical static sea-surface displacement. This is used as initial condition for tsunami propagation modelling, performed with the non-linear shallow water multi-GPU code HySEA (de la Asunción et al.2013), benchmarked in the framework of the US NTHMP program (Macías et al.2016; Macías et al.2017). For all the 100 simulations in the ES subensemble, the simulation duration was fixed at 2 hr, saving the tsunami time-series each 30 s.

We compare the resulting ensemble of 100 synthetic tide-gauge signals with the observed tsunami signals at four tide-gauges (Figs 5a–d), with the purpose of qualitatively addressing how the uncertainty of source parameters is mapped onto the tsunami waveforms (e.g. Lorito et al.2008). The general tsunami energy pattern flowing away from the source, generated with the average kinematic rupture model (Fig. 4a), is also illustrated in Fig. 5(e), which shows the maximum sea level elevation during the whole tsunami simulation.

Figure 5.

Tsunami modelling. (a)–(d) Upper panels show the comparison between observed (black line) and predicted (red line) tsunami waveforms for the slip model in Fig. 4(b); yellow lines represent the predicted tsunami waveforms resulting from the 100 simulations in the ES subensemble. Lower panels show the time–frequency analysis of the observed tsunami signals; the vertical white line indicates the separation of the signal before and after the tsunami arrival with the two different standard deviations used to normalize each portion of the wavelet spectrum, which helps emphasizing the pre-event spectral structure. Panel (e) shows the maximum wave elevation during the numerical simulation; the yellow triangles and the red star represent the tide-gauge positions analysed in this work and the epicentre of the 2018 Zakynthos earthquake, respectively.

It is generally expected that a good enough matching between observed and predicted tsunami waveforms can be achieved only within the first 1–2 signal cycles, whereas the agreement should worsen as the time progresses, because of local possibly not well-modelled propagation complexity around the tide-gauges (e.g. Romano et al.2016). The tsunami simulations were indeed performed using a relatively coarse bathymetric computational grid, with spatial resolution of 15 arcsec (SRTM15±), as sufficiently accurate higher resolution bathymetry around the tide-gauges was unavailable. Moreover, the tsunami periods are in this case quite short, on the order of maximum 15 min. Such short period signals are due to the moderate earthquake size (as compared for example to the typically modelled great megathrust earthquakes and tsunamis), and to the steepness of the initial displacement due to strike-slip faulting (Heidarzadeh & Satake 2014). Hence, it is likely that some unmodelled local signal components such as multiple reflections, or even resonances growing in time, may progressively obliterate the tsunami signal, making the coarsely modelled tsunami smaller and smaller in comparison to observations as time advances. This allows only for a limited consistency between the simulated and observed signals, which is higher at the Greek stations, and lower at the Italian ones. For example, as time passes, this seems to be the case at the farthest Italian tide-gauges (Figs 5c and d), characterized by a lower signal-to-noise ratio than the Greek ones, which are closer to the tsunami source. Also, an unmodelled reflection or oscillation likely appears about 50 min after the earthquake origin time at the Kyparissia tide-gauge (Fig. 5a).

A time–frequency analysis of the tsunami signals, performed by means of the Morlet wavelet (Torrence & Compo 1998), corroborates the above considerations, then suggesting to not pushing the conclusions of this comparison too far. This analysis shows in fact that, as anticipated, some energetic features at specific frequencies that exist prior to the event persist or are reinforced after the tsunami arrival, which is particularly evident at Kyparissia, Katakolo and Crotone tide-gauges (Figs 5a, b and d). The signal at Le Castella behaves differently (Fig. 5c), but we refrain to further commenting this issue since this tide-gauge is a non standard, innovative, one (Annunziato et al.2016), whose response function should be probably better investigated. We point out that the wavelet spectra of the pre- and post-event signals were normalized to the respective time-series variance values, in order to make them comparable. In particular, this serves for illuminating the structure of the pre-event signal spectrum, which, if normalized to the variance of the entire series, would have been much smaller. As time passes, the excitation of the harbour frequencies becomes in general more and more evident mixing up with the tsunami ones. Melgar et al. (2017), noted a similar behaviour, through an analogous analysis, at the Crotone tide-gauge in correspondence of the tiny tsunami generated by the M_w 6.5, 2015 Lefkada earthquake.

4 DISCUSSION AND CONCLUSIONS

We inferred the fault plane and estimated the rupture history associated to the 2018, M_w 6.8 Zakynthos earthquake, by jointly inverting geodetic and seismic observations. Several kinematic inversions were performed, starting from the parameters of the published focal mechanism solutions and from a relocated hypocentre. In particular, we discriminated the main fault plane from the auxiliary one, and finally proposed a preferred kinematic rupture model. The consistency of this rupture model with the observed tsunami signals was finally investigated.

Our preferred model features a non-uniform slip distribution and heterogeneous rupture propagation, characterized by two distinct slip patches, with peak slip of ∼1.8 m on the main NNE asperity. Despite the azimuthal station gaps, the analysis of the slip velocity at different instants during the rupture propagation suggest that: (i) the two slip patches are robust features of the inverted models; (ii) the minor SSW slip area is not an artefact from inversion; (iii) the rupture is characterized by a bilateral propagation. These claims are supported by a synthetic (spike-like) resolution test (Figs S4 and S5), which demonstrates that the stations distribution is in principle suitable to constrain the inversion for this specific fault plane, despite the azimuthal gap. In particular, we consider the minor slip area a robust feature of the rupture model since the spike test indicates enough resolution at that location (Fig. S5) and because this specific slip patch is the only one capable of producing a correctly oriented static displacement at the STRF (Fig. S6) when the static and kinematic GPS data are combined in the inversion.

At least six strong events occurred in the same region in the last 60 yr: on 1958 August 27 (M_w 6.4, Kiratzi & Louvari 2003), on 1959 November 15 (M_w 6.6, Baker et al.1997), on 1962 April 10 (M_w 6.2), on 1968 March 23 (M_w 6.0), on 1976 May 11 (M_w 6.3) and on 1997 November 18 (M_w 6.5, Kiratzi & Louvari 2003). For the 1958 and 1962 events, limited information is available, only concerning the isoseismal trends of the earthquakes (Papazachos et al.1997). For the 1976 event, slightly to the southwest of the 2018 main shock, approaching to the trench, a reverse fault mechanism seems to be reported (Papadimitriou 1993; Baker et al.1997). For the 1959 (McKenzie 1972; Baker et al.1997), the 1968 (Anderson & Jackson 1987) and 1997 (Kiratzi & Louvari 2003) earthquakes the moment tensor solutions are roughly similar to the one of 2018. Furthermore, we noticed that the azimuths of the 1997 coseismic displacement at Zakynthos and Strofades (Hollenstein et al.2006; Hollenstein et al.2008) are comparable with the azimuths of the 2018 coseismic displacements at the same sites. This comparison reinforces the possibility of similar mechanisms for these three events.

Regarding the 1997 earthquake, Kiratzi & Louvari (2003) suggest a source with two different pulses, the second one showing larger scalar moment and stronger strike-slip component, consistently with the available Harvard CMT solution at that time and with the solutions for the 2018 earthquake.

The availability of more near source stations that were in operation in 2018 (especially the one in Strofades island) allows us not only to discriminate the causative N9°E-striking and 39°-ESE-dipping focal plane, but also to image in detail the rupture history. Low-angle (δ < 50°) E-dipping faults in the upper crust (depth < 10 km), above the subduction zone, were also observed in several E–W seismic profiles acquired between the Zakynthos and Strofades islands during the SEAHELLARC project (Papoulia et al.2014). These faults are interpreted as the results, in the Pliocene, of the emplacement of the Ionian thrusts over the Apulian Ridge margin due to evaporite diapirism (Wardell et al.2014). However, the focal solutions of the 1959, 1997 and 2018 earthquakes, as well as the fault geometry and rupture model proposed in this study would suggest that these faults could have been more recently (Pleistocene) re-activated with a dominant right-lateral strike-slip mechanism. In particular, in the transition zone between the HASZ and the CTFZ, we suggest that the tectonic regime related to the CFTZ appears to be still dominant in the upper crust faults between Zakynthos and Strofades.

Further investigation of the dynamics of faults and earthquakes in this area will require denser geodetic and seismological arrays, especially with offshore stations. Nevertheless, although the occurrence of dip-slip events cannot be excluded, we note that the suggested change in the tectonic regime from thrust in almost pure strike-slip would correspond to less vertical sea floor displacement and thus lower tsunamigenic potential caused by an earthquake of a given magnitude, which would have in turn a great importance for tsunami hazard assessment. The influence of this hypothesis could be addressed by testing the sensitivity of existing tsunami hazard models (Basili et al.2018; Basili et al.2019; see acknowledgements for related scientific projects) to different assumptions on the local sources.

Concerning the observed tsunami, we noted that its features are generally consistent with the rupture model obtained in this study. However, this preliminary tsunami analysis deserves further consideration. Higher resolution simulations should be used for the characterization of ports and harbours, especially for relatively small, steep/short wavelength earthquakes, which may resonate with local site periods. Harbour features for example played a role in enhancing the tsunami impact in the Balearic Islands during the 2003 M_w 6.9 Boumerdes-Zemmouri earthquake (e.g. Alasset et al.2006). Modelling tide-gauge observations was challenging also for that tsunami, possibly because of local propagation complexity.

In the framework of tsunami early warning this implies further challenges that should be addressed with high-accuracy bathymetric grids and time-consuming high-resolution simulations.

SUPPORTING INFORMATION

Figure S1. Seismicity associated to the 2018 Zakynthos earthquake. Pink star shows the events with magnitude M > 4, relocated in this study. Black points are all the events with magnitude M > 2.

Figure S2. 1-D P- and S-wave velocity profiles adopted in this study (Table S1).

Figure S3. Inverted rupture model (average model from ensemble inference) for NP2 fault plane geometry of the focal mechanism solution. Upper and lower panels show the total slip and the rise time distributions on the fault, respectively. Rupture time shown by black contour lines (in seconds); black arrows displayed in upper plot represent the slip vector. Red star displays the hypocentre.

Figure S4. Spike-like test target models. Each panel display the target slip distribution assumed for each of the 35 kinematic inversions performed. Black contour lines show the region of the fault plane that slipped more than 0.2 m during the earthquake (Fig. 4a). Red star is the location of the hypocentre. The number on each panel identifies the Spike-like Test’ number.

Figure S5. Spike-like test inverted models. Each panel display the retrieved slip distribution for each of the 35 kinematic inversions performed. Black contour lines show the region of the fault plane that slipped more than 0.2 m during the earthquake (Fig. 4a). Red star is the location of the hypocentre. The number on each panel identifies the Spike-like Test’ number.

Figure S6. Comparison between the observed (blue) and forward modelled (red, whole slip distribution as in Fig. 4(a); green, only the main slip patch) velocity time histories (a) and GPS (b) data.

Table S1. P- and S-wave velocity profiles adopted in this study (extracted from table 1 in Haslinger et al.1999).

Table S2. Fault parametrization adopted in this study and kinematic source parameter values retrieved by the inversion. Each grid node on the fault plane is defined by geographical coordinates (Lon, Lat), depth (Z) and the along strike and along downdip fault system coordinates (U and V, respectively). The subfault length and width are both equal to 4.5 km. Strike and dip angles are kept fix and equal to 9° and 39°, respectively. For each node, corresponding inverted rake, slip, rise time and rupture time values, are given.

Table S3. A priori bounds variability for each kinematic parameter we invert for; given in terms of min, max and step values (Min, Max, Step).

Table S4. Cost function values resulting from the inversions performed by adopting the NP1 and NP2 USGS fault plane solutions. costGPS, costSM and costJoint stand for the minimum cost function value associated to the GPS, Strong Motion and HRGPS, and the combination of the two data sets during the joint inversion, respectively.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

ACKNOWLEDGEMENTS

The data for this paper can be accessed at the doi:10.6084/m9.figshare.8313785.

We thank Roberto Basili and Nicolas Chamot-Rooke for helpful discussion. We thank University of Patras (Greece) to make the LTHK strong motion data available and ITSAK-EPPO (Greece) staff for the maintenance of the KRI1 and ZAK2 strong motion stations (http://www.itsak.gr/en/page/infrastructures/networks/acc_network). We also thank scientific (NOANET, http://www.gein.noa.gr, and RING, http://ring.gm.ingv.it) and commercial (METRICA, https://www.metrica.gr, and URANUS, http://www.uranus.gr/) GNSS networks for the dissemination of their data to the scientific community. We are particularly grateful to the RING (INGV, Italy) and NOANET (NOA, Greece) technical staff for the installation and maintenance of these networks and, in particular for the recent STRF GNSS station installation, so crucial for this work. We thank IOC/UNESCO and VLIZ (Nederland) (http://www.ioc-sealevelmonitoring.org/) for distributing the tsunami data. We also acknowledge the Institutions who own and maintain the tide-gauge stations we used which are: NOA (Greece), ISPRA (Italy), INGV (Italy), JRC (EU). We finally acknowledge the Italian (CAT-INGV) and Greek (NOA) Tsunami Service Providers for providing us the information about the tsunami alert messages issued for the Zakynthos earthquake. Global CMT solutions available at http://www.globalcmt.org.

USGS focal mechanism and finite fault slip model from https://earthquake.usgs.gov/earthquakes/eventpage/us1000hhb1/executive.

Preliminary tsunami simulations by CAT-INGV can be found at https://ingvterremoti.wordpress.com/2018/10/26/il-maremoto-del-26-ottobre-2018-nel-mediterraneo-aggiornamenti-delle-ore-14–00/. Bathymetric model used for the tsunami simulation can be downloaded at https://topex.ucsd.edu/WWW_html/srtm15_plus.html. For the European seismic and tsunami hazard map refer to SHARE (http://www.share-eu.org) and TSUMAPS-NEAM (http://www.tsumaps-neam.eu) projects. AA thanks his daughter Beatrice for allowing him using her laptop for drafting the paper.