The Seismicity of Ischia Island, Italy: An Integrated Earthquake Catalogue From 8th Century BC to 2019 and Its Statistical Properties

Ischia is a densely inhabited and touristic volcanic island located in the northern sector of the Gulf of Naples (Italy). In 2017, the Mw 3.9 Casamicciola earthquake occurred after more than one century of seismic quiescence characterized only by minor seismicity, which followed a century with three destructive earthquakes (in 1828, 1881, and 1883). These events, despite their moderate magnitude (Mw < 5.5), lead to dreadful effects on buildings and population. However, an integrated catalogue systematically covering historical and instrumental seismicity of Ischia has been still lacking since many years. Here, we review and systematically re-analyse all the available data on the historical and instrumental seismicity, to build an integrated earthquake catalogue for Ischia with a robust characterization of existing uncertainties. Supported by new or updated macroseismic datasets, we significantly enriched existing catalogues, as the Italian Parametric Earthquake Catalogue (CPTI15) that, with this analysis, passed from 12 to 57 earthquakes with macroseismic parametrization. We also extended back by 6 years the coverage of the instrumental catalogue, homogenizing the estimated seismic parameters. The obtained catalogue will not only represent a solid base for future local hazard quantifications, but also it provides the unique opportunity of characterizing the evolution of the Ischia seismicity over centuries. To this end, we analyse the spatial, temporal, and magnitude distributions of Ischia seismicity, revealing for example that, also in the present long-lasting period of volcanic quiescence, is significantly non-stationary and characterized by a b-value larger than 1.


INTRODUCTION
The characterization of the seismic activity and of the consequent hazard is largely based on the knowledge that we can gain from past activity, that is, from seismic catalogues. To allow a robust characterization of the seismicity, such catalogues should systematically collect all the known seismic events, as well as characterize at best the completeness of the record in the different periods and the existing uncertainty about seismic parametrization. Especially for historical times and for light to moderate magnitudes (M < 6), this is often very challenging.
Ischia is a volcanic island within the Neapolitan volcanic system. It is located to the southwest of the Campi Flegrei caldera and, along with the volcanic islands of Procida and Vivara, they form the Phlegraean volcanic district (Figure 1). The volcanotectonic framework of Ischia is rather complex: it represents the emerged part of an active volcanic field, which rises more than 1,000 m above the seafloor along the margin of an E-W trending scarp bordering to the south the Phlegraean volcanic district . Seismic hazard is just one of the possible hazards related to the local volcanic system, responsible in the past for landslides, tsunamis, hydrothermal explosions, and volcanic eruptions (see Selva et al., 2019 for a review). Ischia has a long record of eruptions also in historical times, the last one being the 1302 AD Arso lava flow (Vezzoli, 1988;Iacono, 1996), after which it began a long period of volcanic quiescence that nowadays has been lasting for more than 7 centuries. Ischia has been affected by a significant intra-calderic resurgence of approximately 1000 m (Orsi et al., 1991;Acocella and Funiciello, 1999). Seismicity is mainly concentrated on the E-W structures that limit to the north the resurgent block of Mt. Epomeo Cubellis and Luongo, 2018;Trasatti et al., 2019;Cubellis et al., 2020;Carlino et al., 2021;Nappi et al., 2021;and references therein). Local seismicity may be connected to a relevant seismic hazard (Cubellis et al., 2004;Selva et al., 2019) and, in a multi-hazard view, it can be seen as one of the possible events in a multi-hazard chain that may include other phenomena like landslides or tsunami (e.g., Violante et al., 2003;Della Seta et al., 2012;Paparo and Tinti, 2017;Selva et al., 2019).
On the evening of the 21st August 2017, a M w 3.9 earthquake struck the island of Ischia, causing 2 casualties, many injuries, and extensive damage in the village of Casamicciola, on the northern coast of the island (Azzaro et al., 2017). The earthquake was very shallow (depth <1 km), with a source striking E-W in the northern sector of Ischia revealed by ground deformation and surface fracturing (e.g., Gruppo di Lavoro INGV, 2017;Braun et al., 2018;De Novellis et al., 2018;Nappi et al., 2018Nappi et al., , 2019Calderoni et al., 2019;Ricco et al., 2019). Different solutions of the source geometry have been proposed from the joint interpretations of seismic and geodetic data, or field geological data alone, differing in the dip-angle and the fault plane geometry at depth, and its generation mechanism (Braun et al., 2018;De Novellis et al., 2018;Devoti et al., 2018;Calderoni et al., 2019;Cubellis et al., 2020).
The 2017 strong earthquake struck Ischia after tens of years of very low seismicity (D'Auria et al., 2018). However, the same area was strongly hit by a series of heavily damaging earthquakes in the XIX century (Cubellis and Luongo, 2018;Guidoboni et al., 2018;Cubellis et al., 2020), the most devastating of which occurred on July 28 1883, causing more than 2000 casualties and heavy destruction in a large part of the island and following by only two years another destructive event occurred in 1881 (Cubellis, 1985;Alessio et al., 1996;Cubellis andLuongo, 1998, 2018;Cubellis et al., 2004Cubellis et al., , 2020Luongo et al., 2006;Guidoboni et al., 2018). Therefore, the 2017 event renewed the attention on the evaluation of the seismic hazard and related risk at Ischia (e.g., Avvisati et al., 2019;De Natale et al., 2019;Nappi et al., 2021), considering that this small island (less than 50 km 2 ) is densely inhabited with more than 60,000 steady inhabitants, incremented by tourists during the summer season.
The interpretation of the genesis of the seismicity at Ischia and the relationship with its volcanic activity is still debated (e.g., Martinelli and Dadomo, 2017;Cubellis and Luongo, 2018;Calderoni et al., 2019;Selva et al., 2019;Trasatti et al., 2019;Cubellis et al., 2020). The existence of a high-temperature gradient hampering brittle behaviour of the medium limits the seismogenic volume to the first 2-3 km Castaldo et al., 2017;Cubellis et al., 2020). In any case, the existence of long periods of seismic quiescence with very low seismicity rates and of sequences of destructive earthquakes such as the one that occurred in the XIX century shows the complexity of the processes that generate the earthquakes at Ischia.
A unified, extended earthquake catalogue for Ischia does not exist yet. Given the complex character of the seismicity at Ischia, a homogeneous characterization of both historical and instrumental seismicity is fundamental for quantitatively addressing scientific analyses about its spatio-temporal evolution, to better constrain its origin and quantify the seismic hazard and the related risk. Indeed, most of the literature studies have been focused on the main historical events (Cubellis, 1985;Luongo et al., 1987Luongo et al., , 2006Alessio et al., 1996;Cubellis and Luongo, 1998;Cubellis et al., 2004;Carlino et al., 2010), and only some macroseismic datasets (Guidoboni et al., 2007(Guidoboni et al., , 2018 were used to compile the Italian parametric earthquake catalogue (hereinafter CPTI15 1 ; Rovida et al., 2016Rovida et al., , 2019Rovida et al., , 2020. In this frame, the magnitude estimates are controversial, while a systematic collection of the minor seismicity is missing. On the other hand, the instrumental seismicity presented by D' Auria et al. (2018) results limited to the most recent period, with duration magnitude estimations (M d ) incompatible to CPTI15. Therefore the existing catalogues are, in their present form, inappropriate to characterize the seismicity of Ischia due to the very low number of events reported, as well as earthquake parameters not directly comparable to each other and uniform through time.
In this paper, we systematically revise the historical and instrumental seismicity in order to produce the first integrated catalogue of the earthquakes that can be attributed to the local volcano-tectonic activity of Ischia. The analysis of the historical macroseismic records (Section "Macroseismic Catalogue: 8th Century BC -2019") consisted of: (i) scrutinizing the literature data to include all the earthquakes caused by local sources in the island; (ii) homogenizing intensity data; (iii) re-parametrizing the collected data with a standardized procedure that takes into account the relevant uncertainty; (iv) evaluating the completeness of the catalogue from the historical standpoint. As for the recent instrumental seismicity (Section "Instrumental Catalogue: 1993-2019"), we: (i) discuss the evolution of the local instrumental network; (ii) revise the data and the parameter estimations, extending the analysis as far back as the 1990s; (iii) evaluate the completeness of the catalogue. The resulting integrated catalogue (macroseismic plus instrumental records) is then statistically analysed (Section "Statistical Characterization of the Seismicity") through: (i) the characterization of the frequency-size distribution and related uncertainty; (ii) the characterization of the occurrence model of the seismicity, by testing the hypothesis of stationarity in time utilizing a statistical test; and iii) the characterization of the spatial distribution and related uncertainty. The result is an earthquake catalogue spanning over ten centuries, which represents the primary dataset for investigating the long-term behaviour of the volcano seismicity of Ischia and assessing the related seismic hazard at local scale.

MACROSEISMIC CATALOGUE: 8TH CENTURY BC -2019
The long-term seismicity is described in the historical earthquake catalogues by macroseismic intensity data, which are parameterized in terms of epicentral location and magnitude. For Ischia, the CPTI15 catalogue reports the main events above the damage threshold (corresponding to epicentral intensity I 0 > V-VI MCS). However, many studies investigated other minor damaging earthquakes as well as simply felt events. Therefore, there is the possibility to increase the number of records of the historical catalogue taking into account all the results available in the literature. In doing this, we also revised the whole procedure to parametrize earthquakes, that is, the Macroseismic Data Points (MDPs) production, the intensitymagnitude relationship, and the epicentral determinations. This is crucial to produce reliable and homogeneous information through time. Note that, at the present state, there is a large variability of the parameters in the literature even for the most studied events: for example, the magnitude estimations for the largest known event (1883) range in literature from 4.3 to 5.2 (CPTI15, Rovida et al., 2019Rovida et al., , 2020Cubellis and Luongo, 1998;Cubellis et al., 2020;Carlino et al., 2021).
In the following, we describe the procedure adopted to produce the new revised macroseismic catalogue for Ischia through three main steps: (i) producing a comprehensive inventory of the known earthquakes; (ii) homogenising geographically the macroseismic intensity data; and (iii) parametrizing with a standardized procedure the events to determine their magnitude and location, accounting for uncertainties.
The revised catalogue, covering the time-span from the 8th century BC to 2019, is reported as Supplementary Dataset 1.
Step 1: Inventory of the Earthquakes We collected all the information available in the literature into a comprehensive inventory. The starting point is represented by the events included in the Italian Archive of Historical Earthquake Data (hereinafter ASMI; 2 Rovida et al., 2017) the basic tool for collecting, homogenizing, and validating macroseismic data in Italy -updated to 2017 and from which CPTI15 is compiled.
Before the present work, ASMI provided a list of 13 earthquakes for Ischia, known through 18 different studies and catalogues. We integrated this list with a large number of events reported in other studies (for example, Cubellis, 1985;Luongo et al., 1987Luongo et al., , 2006Alessio et al., 1996;Cubellis and Luongo, 2018), and in the CFTI5Med catalogue (Guidoboni et al., 2018), as well as some other unpublished research data. The resulting inventory is reported as Supplementary Dataset 2. It contains a list of 245 studies referred to 102 earthquakes with related parameters and also the primary sources. For each event, we also evaluated the state of knowledge through the critical analysis of data and sources. 2 https://emidius.mi.ingv.it/ASMI/index_en.htm This process mainly focused on the evaluation of the reliability of the primary sources -direct or indirect, coeval or later, local or external to the island etc. -and on the existence of explicit references to the earthquake effects at Ischia. In the inventory, we tracked any specific doubt or issue that emerged, for future investigations. As a result, we distilled the content of the inventory by removing fake and doubtful events potentially generated by other phenomena, and those generated by seismic sources located outside the volcanic system of Ischia. The analysis allowed eliminating 24 events. Among these, it is worth mentioning some events which derived from the misinterpretation of historical sources (e.g., 28-29 July 1762, 1852) and the events of the 1880 seismic sequence located near the Pontine islands (80 km NW of Ischia). The same occurred for more recent earthquakes, such as the 1980 one, now relocated by instrumental data in the Tyrrhenian Sea, and the 1983 event relocated in the Campi Flegrei area.
For several other events, the review of the basic information posed significant doubts about their origin. The most evident situation is when the occurrence of an earthquake is not explicitly reported in the original documents, but it derives from the interpretation of other phenomena to which it could be potentially associated. For example, large landslides may be interpreted as triggered by a strong earthquake (e.g., the one that occurred in 1228), and volcanic eruptions are potentially associated with earthquake swarms. Detailed discussions about these events can be found in Buchner (1986), Civetta et al. (1999), Polara and De Vivo (2011), Guidoboni and Comastri (2005), Gialanella (2013), de Vita et al. (2013, Cubellis and Luongo (2018), and Cubellis et al. (2020). In these cases, we indicated a record in the catalogue reporting "questionable" in the column "Notes".
In a few cases, we also found groups of events not distinguished in terms of time and intensity (e.g., 1883/07/28 or 1932/02/12). In these cases, we indicated a record in the catalogue reporting the information about multiple events in the column "Notes".
In conclusion, the inventory of the historical earthquakes referring to Ischia consists of 102 events, 24 fake ones, and 5 doubtful ones.
Step 2: Homogenization of Intensity Data The rich bibliography available for the historical earthquakes of Ischia includes studies based on different geographic positioning of the localities therein analysed; this makes original intensity data (both values and distribution) not directly comparable to each other. In order to overcome this problem, we applied a homogenization procedure to univocally define a locality and, hence, to produce comparable seismic histories as well as earthquake parameters (see next section).
As a first step, we adopted the same national geographic gazetteer used in DBMI15 , which grants the unequivocal association of a place name with a pair of geographical coordinates (lat, long, ID, place name and other geographical and administrative information). The intensity data available in the literature have been therefore assigned to this geographic reference, which includes municipalities and main inhabited localities. Because of the dense urbanization of the area and the strong decay of intensity in very short distances, we also added new points to the gazetteer in order to better sample the macroseismic effects in the several minor settlements spread in the territory. Details on the implementation of the gazetteer's localities for Ischia are provided in Supplementary Appendices.
The above operation has been accompanied by the homogenization of the MDPs available in the literature. For the best-documented earthquakes (namely 1828, 1881, and 1883), we generally grouped the great number of original MDPs into the reference localities, but we also did the reverse, i.e., we ungrouped information referred to a locality into new sites if adequately supported by the historical sources. When available, alternative studies about the same earthquake were considered. For 49 minor events, the historical information has been instead reviewed and new MDPs were produced. Details of this revision are reported in Supplementary Appendices. The adopted intensity scale is MCS (Cancani, 1904;Sieberg, 1912).
In total, we obtained new (revised or produced ex-novo) MDPs distributions for 54 earthquakes, which significantly improved the ASMI database. The MDPs for the largest known earthquakes in Ischia (1828Ischia ( , 1881Ischia ( , 1883Ischia ( , and 2017) are reported in Figure 2.
Step 3: Parameterization of Historical Earthquakes The seismotectonic features of Ischia required specific procedures for defining the epicentral parameters of the historical events from intensity data, particularly regarding the epicentral location and the magnitude determination. As for the latter, it is widely recognized that the use of regional intensity-magnitude relationships is inappropriate in volcanic areas (for Italy, see for example the Etna case in Azzaro et al., 2011) that are typically characterized by extremely shallow hypocentres like Ischia (depth = 1-2 km; e.g., Cubellis et al., 2020). On the other hand, the low number of events recorded in recent times (D'Auria et al., 2018) prevents the formulation of specific relationships. To better interpret spatio-temporal variations and characterize the seismic hazard and risk, it is critical to homogenize the parameters from different data sources, characterizing as much as possible the existing uncertainty. In the following, we shortly describe how this was achieved for Ischia.

Epicentral Location
The epicentral location from intensity data is routinely calculated in the CPTI15 catalogue through the latest version of the "Boxer" code (4.0; Gasperini et al., 2010). Also for Ischia, we applied the simplest technique for the determination of the epicentrethe so-called "Method 0" -that is calculated as the centre of gravity (truncated average of their coordinates) of the sites with highest intensities, since it proved to be robust in case of poor intensity datasets and provided reliable results at Etna, where the macroseismic features are similarly characterized by high-intensity attenuation and the extreme concentration of damage in small zones.
The results of the location estimation are reported in Figures 3A,B; note that the areas represent only the uncertainty of the epicentres and do not indicate any source geometry.
Given the dense urbanization of Ischia, an estimation of location uncertainty is of primary importance. To this end, we first defined a reliability index for the location of the earthquakes, based on the number of data available: VL (Very Low) for 1 MDP; L (Low) for 2-5 MDPs; H (High) for MDPs >5. For events labelled as VL, we did not provide an estimation of the epicentre while for those labelled as L we calculated only the epicentre, since data are considered not suitable for quantifying the related uncertainty. Finally, epicentre and related uncertainty are reported for all the events classified as H. To do this, we exploited the bootstrap method implemented in Boxer 4.0: the code fits a 2D Gaussian distribution, quantifying its 2 × 2 covariance matrix.
For the earthquakes documented by multiple studies (1828 and 1883), discarding any "unconventional" operation of average among the different intensity values available for a given locality, the solution usually adopted in literature is to reassess intensity starting from the reanalysis of the historical sources and to produce an unified set of MDPs for each earthquake from which an epicentral solution was calculated anew. This approach requires bringing together the authors of all the existing historical studies overcoming the different macroseismic interpretations, to reduce existing epistemic uncertainty. Here, given the purpose of the present paper, we prefer to adopt a different strategy that aims at quantifying the uncertainty arising from alternative interpretations, in order to evaluate its impact on the definition of earthquakes' parameters. Starting from the epicentral solutions and related uncertainty from the geographically homogenized data of the multiple studies, we integrated them with an ensemble model (Marzocchi et al., 2015). Ensemble modelling consists of fitting a parent distribution to the available alternative models to quantify the emerging epistemic uncertainty, weighting each alternative by their credibility. Ensemble modelling represents an evolution of the Logic Tree technique usually adopted in seismic hazard quantification, and provides more robust results when few alternatives are available (for more details, see discussion in Marzocchi et al., 2015 and reference therein). In doing this, we assigned the same weight to the available solutions, merging the best guess location and the related uncertainty. The ensemble model is obtained by fitting a Gaussian distribution through the Maximum Likelihood Estimation using a sample of locations obtained combining two samples with equal size from the individual input uncertainty distributions (Selva et al., 2014(Selva et al., , 2018. The results of the ensemble models are shown in  Ciuccarelli et al., 2018). Thus, the ensemble model extends to almost the larger entire area FIGURE 2 | Intensity maps of the largest earthquake of Ischia. Coloured circles are MDPs; blue triangles, solid and dashed lines represent the location uncertainty distribution (best guess; 1 and 2 σ areas, respectively) for all the earthquakes. When multiple studies are available for the 1828 (A,B) and 1883 (C,D) events, the location uncertainty distribution represents also the ensemble of the available models (in black, see text for details). When only one study is available [1881 in panel (E); 2017 in panel (F)], the ensemble and the study results coincide. EE stands for Environmental Effects.
of uncertainty, with the highest probability in the area in common to the two individual studies. As for the 1883 earthquake, the two studies available (Cubellis and Luongo, 1998;Ciuccarelli et al., 2018; see Supplementary Appendices) produce similar locations, so the ensemble model does not differ much from either individual study, with the zone of highest probability (1 sigma) substantially covering the maximum damage area.

Magnitude
Because of the lack of a specific intensity-magnitude relationship (hereinafter IMR) for the seismicity of Ischia, we built ensemble models taking into account alternative IMRs produced for other Italian volcanic areas, namely Campi Flegrei and Etna, where high values of epicentral intensity are associated to low to moderate magnitude, similarly to Ischia. Given that the reference magnitude for the CPTI15 catalogue is the moment magnitude M w , while for the instrumental magnitude is traditionally M d (see Section "Magnitude"), we produced ensemble models for both the magnitude scales. Considering that CPTI15, in analogy to Etna, estimated the magnitude M w for the Neapolitan volcanoes using first an IMR to estimate M L (Richter, 1935;Gasperini, 2002) from the intensity (Azzaro et al., 2011) and then converted M L to M w , we also considered procedures to quantify M w passing through the estimation of M L .
For Campi Flegrei, Marturano et al. (1988) proposed two IMRs (linear and logarithmic, hereafter MAR88lin and MAR88log, respectively) obtained from few macroseismic data of the bradiseismic crisis of 1982-1984, calibrated on the duration magnitude M d (Branno et al., 1984;Marturano et al., 1988). Since neither model reports uncertainty in the model parameters, we reconstructed the input data to reproduce the best-guess linear and logarithmic IMRs (see Figure 9 in Marturano et al., 1988) with the same methodology adopted in the original paper (standard least-square). We obtained a standard deviation of residuals equal to 0.51 and 0.49, respectively.
A M d -M w relationship has been produced by Petrosino et al. (2008; hereinafter PET08; Figure 4A), calibrated on the instrumental magnitude range 0-4 through a standard leastsquare procedure. Adopting the fit of Petrosino et al. (2008), the uncertainty on M w evaluation from M d results 0.30 (1 standard deviation, evaluated from the original data). With a similar procedure, PET08 produced a M L − M w relationship, with uncertainty on M w evaluation of 0.13.
For Etna, a large set of instrumental and macroseismic data was used by Azzaro et al. (2011;hereinafter AZZ11) to derive the IMRs for M d and M L , adopting a standard least-square procedure. The authors indicate the standard deviation of the residuals as 0.12 for M d and 0.36 for M L . The conversion from M L to M w can be then obtained through the relationship by Saraò et al. (2016;hereinafter SAR16), calibrated on a large set of instrumental data (magnitude range 2.1-4.8) adopting an orthogonal least-square relationship, with a reported uncertainty on M w of 0.2 (Saraò et al., 2016).
Here we considered the available IMRs to define 3 alternative procedures for estimating M d (models D1 to D3) and 5 for M w (models W1 to W5), as reported in Table 1. Note that AZZ11 combined with SAR16 (model W5) and PET08 (model W1) was previously adopted in CPTI15 for estimating M w at Etna and Neapolitan volcanoes, respectively. To account for the uncertainty of each of these procedures, the ensemble takes in input not only best-guess estimates, but also sets sampling the existing uncertainty. Notably, all the models are accompanied by the quantification of uncertainty on the final magnitude, for which we assume a normal distribution with a standard deviation equal to the uncertainty declared in the original study.
To propagate the uncertainty also in the models that consider the sequential application of two relationships (namely the ones to obtain M w : models W1 to W5), we applied a sequential sampling procedure. We first sampled the first uncertainty distribution and then, for each sampled value (M L for W1 and W5, M d for W2, W3 and W5), we applied the second relationship to obtain Mw, sampling also its uncertainty. In this way, we obtained samples of M w that propagate the uncertainty of both models.
These 3 + 5 models represent all the possible procedures, but they cannot be considered equally credible. In particular, MAR88lin and MAR88log models derive from similar data and mainly differ for large magnitudes, which are outside the magnitude range of the original calibration. Given that MAR88lin calculates unrealistically high magnitudes, which for the largest intensities would require too large fault sources for the size of the island and its seismogenic sources, we prefer not to include in the ensemble the models performed with this relationship (Model D2 and W3 for M w ). For M d , the remaining models (D1 and D3) were weighted equally for all intensities. For M w we adopted a more elaborate strategy. Model W2 provides a slope markedly different from the other models, leading to relatively large magnitudes for low intensities and vice versa, and therefore we decided to remove it from the ensemble. The remaining models (W1, W4, and W5) were differently weighted for low (I 0 < VIII) and high (I 0 ≥ VIII) epicentral intensities.
Considering that Campi Flegrei may be considered a volcanotectonic environment more similar to Ischia, but high intensity and magnitude values are available only for Etna, we defined the following weights: for I 0 < VIII, a weight of 2 is given to models including PET08 relationship (Models W1 and W4), and a weight of 1 to the remaining model (Model W5); for I 0 ≥ VIII, equal weight is assigned to all the models (W1, W4, and W5). The ensemble models are then obtained as in Section "Epicentral Location", by fitting a Gaussian distribution through the Maximum Likelihood Estimation using a sample of potential values of M d and M w obtained by combining different samples related to the single models, whose size was proportional to the credibility weight of each model (Selva et al., 2014(Selva et al., , 2018. The results are shown in Figures 4B,C, while the numerical results are reported in Table 2. In general, the ensemble models show a good coherence with the dimensions of seismic sources in such a volcano-tectonic setting, that is, high intensities can be reached with relatively low magnitudes, given that the earthquake hypocentres are very shallow. The seismic cut-off is located at about 2 km depth due to the high geothermal gradient (150-220 • C/km) while causative faults display short lengths because of the extreme fracturing of the shallow crust. These structural conditions limit the maximum possible magnitude estimated between 5.3 and 5.8, considering different assumptions Luongo et al., 2006;Castaldo et al., 2017;Cubellis and Luongo, 2018;Nappi et al., 2018;Cubellis et al., 2020). In Figure 5, we report the magnitude estimations for all the earthquakes in the macroseismic catalogue.
In order to verify the ensemble IMRs, we compared the ensemble models with the instrumental magnitudes of the two earthquakes for which both instrumental and macroseismic intensity estimations are available. This comparison is blind, since these data have not been used for setting up the models. The 2017 earthquake had an epicentral intensity I 0 = VIII (Azzaro et al., 2017) and several magnitude estimations in the literature. In Figure 6A, we compare these data with the probability distributions for an intensity I 0 = VIII; both ensemble models (for M w and M d ) are in good agreement with the instrumental estimates. In Figure 6B, we report the comparison for the 2008 earthquake, with I 0 = V (Cubellis and Marturano, 2009). In this case, even if the ensemble model for M w appears compatible with data, the one for M d appears significantly biased toward higher magnitudes, also considering an uncertainty of 0.3 on the instrumental M d . We note, however, that both the ensembles would be compatible considering a lower epicentral intensity (I 0 = IV).
While for the 1883 earthquake the alternative studies (Cubellis and Luongo, 1998;Ciuccarelli et al., 2018; see Supplementary Appendices) indicate the same epicentral intensity (I 0 = XI), for the 1828 earthquake the two available studies provide two slightly different intensities (I 0 = VIII-IX and IX, respectively). As made in Section "Epicentral Location", we fitted a 1D Gaussian ensemble model with the Maximum Likelihood Estimation from a sample of potential input magnitudes from both models. As expected, the obtained ensemble in Figure 5C shows a distribution centred at an intermediate level, slightly  Completeness" for the historical and the instrumental catalogues, respectively); the dotted red lines report the assumed uncertainty bounds (+/-0.2 M w ) for the completeness magnitude used to check the stability of parameters of the magnitude-frequency distribution (more details in Section "Magnitude-Frequency Distribution"). (C) Ensemble magnitude estimation for the 1828 earthquake compared with ensemble's members (more details in Section "Magnitude").
more dispersed than the original distributions. Noteworthy, the uncertainty of I 0 derived from the alternative MDPs is significantly smaller than the one resulting from the IMR.

Historical Completeness
The long-term seismicity of Ischia appears rather discontinuous and fragmentary, mainly for the most ancient periods. In general, the lack of events in a given historical period does not necessarily indicate a seismicity gap but rather a break in the continuity of information from the contemporary sources (Valensise and Guidoboni, 1995;Stucchi et al., 2004;Cubellis and Luongo, 2018). This situation may depend on several reasons related to both the territory and the social or political context. As for the Neapolitan region and Ischia, a detailed description of the historical situation from Antiquity to recent times is reported in Supplementary Appendices. In the following, we shortly recall the main facts referred to Ischia to evaluate the degree of completeness of its seismic history. During the Greek and Roman periods (VIII BC -V century AD), Ischia was populated mainly along the coast  (Delizia, 1987). The available sources for this period have recorded only earthquakes associated with eruptions, which dramatically impacted the inhabitants. In the early Middle Ages (VI-X century AD), Ischia was scarcely inhabited (Buchner Niola, 1965;Cundari, 1998) and no information is available about local seismic activity. Starting from the XI-XII centuries, the island began to be part of a wider economic and administrative context for the extraction of alum but, unfortunately, the relevant administrative documents were lost during the Second World War bombardments in Naples. In the XIII century, some contemporary sources recorded two natural events of significant impact for the island, a landslide in 1228, and an earthquake in 1275.
The 1302 eruption caused a serious impact on mining, agricultural, and fishing activities, determining a depopulation of Ischia until the second half of the XVI century. As a consequence, in the time-span XIV-mid-XV century, there is a substantial lack of information about facts and events of the island, and probably the apparent lack of seismic events cannot exclude the occurrence of moderate earthquakes (I 0 < VIII MCS), having left no traces in the few written sources of that time. Starting from the late-XV century, the role of thermal baths became increasingly important for Ischia: the settlements grew and the thermal treatments became the main activity, mostly at Casamicciola. Thus in the XVIII-XIX centuries, thanks to the proximity to the Neapolitan area, the island represented an elite destination for Italian and European tourism, so for this period it is very likely that also traces of events of lower intensity should have been reported. This condition of cultural interest for the "geological" events occurring in Ischia was strengthened by the foundation of the Osservatorio Vesuviano in Naples in 1841 and, lastly, of the Osservatorio Geodinamico in Casamicciola in 1885.
In conclusion, the results of the analysis of the historical completeness reported in Table 3 show that the catalogue can be considered complete starting from the mid-XVIII century as regards events of moderate intensity (I 0 ≥ VII MCS), while the completeness for low intensity earthquakes (I 0 ≥ IV MCS) significantly jumps to 1885, when the first local seismic monitoring system was installed at Ischia (Luongo et al., 2012). The first seismic sensor deployed on the island of Ischia dates back to 1885, when G. Grablovitz installed a seismic tank in the newly founded Casamicciola observatory (Grablovitz, 1901;Grablovitz, 1902Grablovitz, -1903Ferrari, 2009;Luongo et al., 2012). The first modern seismic station was installed in the same location by the Osservatorio Vesuviano (hereinafter INGV-OV) (station OC9) in 1993. Since then, the seismic network has been regularly improved up to the present state. The current network is the result of a first set of three analogue stations integrated over the years with digital ones (D'Auria et al., 2018). Since 2015 the network counted 4 sites instrumented with 3 analogue velocimeters (OC9, FO9, and CAI), 3 digital velocimeters (IOCA, IMTC, and IFOR), and 1 accelerometer (IOCA). After the 2017 earthquake, the network was upgraded with 5 sites instrumented with velocimeters and/or accelerometers of the permanent seismic network and up to 6 seismic stations of the temporary network (Galluzzo et al., 2019). The configuration of the seismic network since late July 2018 until present is reported in Figure 1B, while the complete list of stations, time of installation, and sensor type is reported in Table 4.
The most recent instrumental seismic catalogue of Ischia, presented by D'Auria et al. (2018), contains earthquakes located or detected from January 1999 until February 2018. Instrumental earthquake data are routinely produced by INGV-OV and periodically updated, and made available on the web (the Ischia instrumental online catalogue) 3 . The catalogue includes origin time, duration magnitude, and, when possible, location estimation. In some cases, in order to save information on the occurrence area of a seismic event, the catalogue is integrated with the indications of the felt area.
As detailed in the next subsections, to enlarge the instrumental dataset and make it comparable with the macroseismic catalogue,

Revision of the Instrumental Seismic Catalogue
Extension to the Pre-1999 Period The revision of the instrumental seismic catalogue was carried out to reconstruct a robust catalogue of pre-1999 seismicity, since the first modern seismic station was installed in 1993.
The first step was the recovery of the information present in the "Reports on the Surveillance Activity, " periodically prepared by INGV-OV and addressed to the Italian Civil Protection Department.
From a first comparison between the Reports and information reported in the paper archives, we found discrepancies requiring a deeper analysis. As a consequence, the archived seismic traces have been visually verified, when still available. This painstaking analysis allowed us to distinguish between local earthquakes and other transient signals such as anthropic events (explosions by abusive fishermen), atmospheric events (thunders), or regional earthquakes.
The re-examination led to the identification of 42 seismic events dated between 1993 and 1998. The maximum magnitude is M d 1.5 in 1997, recorded by the 3 seismic stations working on the island. Other events of low magnitude were recorded by one or two stations (mainly by FO9; Table 4) and therefore not located. Close to FO9, in the south-western sector of the island, some very low energy shallow earthquakes were identified in a small seismogenic area characterized also by a strong geothermal activity (Selva et al., 2019), in the same area where some events were recently located in 2018, after the improvement of the seismic network.

Epicentral Location
Many events cannot be located, because of their low magnitude and the scarce number of seismic stations deployed on the island.
The first located event was a M d 1.3 earthquake that occurred in 2007. The situation slightly improved in 2015, when the network passed from three to four seismic stations, and further improved in 2018, when the present-day monitoring network became available (Figure 1 and Table 4). The fairly low seismicity rate recorded in recent times in Ischia prevented also the development of detailed tomographic images by using local earthquake recordings. Three velocity models are currently used to locate earthquakes on the island (Figures 7A-C The 1D crustal models (model A in green and B in blue) are very similar (model A resulting slightly faster than B) while model C shows the greatest differences. Models B and C were obtained by using data from the SERAPIS tomographic experiment (Judenherc and Zollo, 2004), aimed at defining the velocity model of Campi Flegrei. Therefore, despite the presence of some seismic stations on the islands of Ischia and Procida, the number of seismic rays crossing the crust below Ischia is relatively low, not allowing a detailed tomography. As a consequence, model C has a resolution of 250 m in the Campi Flegrei area which increases up to 1 km for Ischia (whose size is approximately 9 km × 6 km). Model C is characterized by a strong velocity contrast between Campi Flegrei and Ischia, which, in the discretized velocity model, results in a highly heterogeneous vertical layer.
Model A is the one adopted for the locations reported in the instrumental catalogue (see Figures 3A,C, 7C). Locations are performed using the Hypo71 program (Lee and Lahr, 1972). The epicentral errors, ERH, are estimated through the square root of the sum of the estimated latitude and longitude variance while the vertical error, ERZ, can be interpreted as a 68% confidence interval, assuming a chi-squared value of 1.00 (Boyd and Snoke, 1984;Husen and Hardebeck, 2010). As for the macroseismic catalogue, we classified each event in terms of quality. In particular, we assigned a high-quality class when location error is smaller than the island size (5 km) and a low-quality class when it is larger. When the horizontal and vertical errors were not present in the catalogue we reported the lower limit estimate of 10 km, which is already a relatively large value considering the size of the island and its seismogenic structures (Selva et al., 2019;Trasatti et al., 2019). When the estimation of the location was impossible, we assigned a quality class Null (N).
Notably, this estimation does not account for the uncertainty on the velocity model, which can potentially lead to an important underestimation of the effective uncertainty (Husen and Hardebeck, 2010;Garcia-Aristizabal et al., 2020).
An indication of the impact of the uncertainty of the velocity model on the location is provided in Figures 7C,D, where the largest earthquakes (M d ≥ 0.9) of the 2014-2018 period are located with the A and B velocity models. Note that model B, characterized by higher velocities, tends to spread the epicentres.
To further deepen into the uncertainties introduced by the velocity model and by the use of a small number of stations, we compared the locations obtained with each of the available velocity models for the 3 strongest earthquakes recorded from 2008 to 2018. The analysis, reported in Supplementary Appendices, indicates that whether the uncertainty estimation provided by the location algorithm is often small (<0.5 km), the uncertainty due to the velocity models can be very significant. In addition, the possibility of using more stations is critical to reduce location uncertainty. This is indeed particularly important for Ischia, considering its small size and the complex 3D structure (Selva et al., 2019 and references therein).

Magnitude
The magnitude reported in the catalogue is the duration magnitude M d , based on coda duration, that allows rapid estimates even when the seismic traces are saturated or the signal to noise level is low (Petrosino et al., 2008). M d is still used to compile the seismic catalogue of Ischia because the magnitude of the recorded earthquakes is usually very low (<2.5) and the seismic stations very noisy. In these cases the M L often cannot be estimated.
Duration magnitude is estimated using the relationship derived for the Campi Flegrei caldera (D'Auria et al., 2018) and the duration is estimated through the visual analysis of the seismograms. Before 2017, the seismicity of Ischia was characterized by small and shallow events, most of which were detectable only in Casamicciola Terme; therefore it was not possible to define a magnitude-duration relationship. Thus, the scale created for the Campi Flegrei was adopted also for Ischia on the basis of the similar geological and seismological features of the two volcanoes. The duration-magnitude empirical relation can be found in Orsi et al. (2004).
Based on experience, we consider an uncertainty of 0.3 associated with low M d values. For the events with the highest magnitude (M d 4.0), where the experience cannot be invoked, the uncertainty quantified in the macroseismic catalogue is probably more appropriate (Section "Step 3: Parameterization of Historical Earthquakes").
In order to allow the comparison between the instrumental and macroseismic local and national datasets, we added the estimate of M d and M w to both catalogues. For the instrumental one, following the same approach made for the macroseismic part (Section "Magnitude"), we used the M d − M w relationship defined by PET08 shown in Figure 4A. The uncertainty was estimated by combining the uncertainty of M d (assuming a normal distribution with σ = 0.3) with the one of the M d − M w relationship (normal distribution with σ = 0.3, see Subsection "Magnitude" in Section "Step 3: Parameterization of Historical Earthquakes"), that is, similarly to the procedure adopted for the macroseismic catalogue, we sampled both the uncertainty distribution of M d and, sequentially, the one from the M d − M w relationship for each sampled M d . The M w uncertainty reported in the instrumental catalogue corresponds to 1 σ.

The Sensitivity of the Instrumental Network and Completeness
Considering the low number of events in the instrumental catalogue of Ischia and the potential nonstationarities of local seismicity (Selva et al., 2019; further discussed in Section "Statistical Characterization of the Seismicity"), we evaluated the completeness by analysing the sensibility of the seismic networks operating through time, rather than performing statistical analyses (e.g., Schorlemmer and Woessner, 2008;Tramelli et al., 2013a,b). To this end, we applied the SENSI code (Orazi et al., 2013) 4 to the network composed of 4 seismic stations (as it was on August 21, 2017), to the current permanent seismic network (9 stations), and to the integrated seismic network given by the union of the mobile and permanent seismic networks (14 stations).
The simulation of the detection and location threshold was carried out down to a hypocentral depth of 1,500 m because the high thermal gradient determines the ductile-fragile transition at a depth of about 2 km (e.g., Carlino et al., 2006;Castaldo et al., 2017;Cubellis et al., 2020), inhibiting deeper seismicity.
In Figure 8a-f we show the minimum magnitude to identify an earthquake with a hypocentral depth of 500 m (panels a, b and c) or 1,500 m (panels d,e, and f) b.s.l. at least at one seismic station. Panels a, d display the detection threshold for a seismic network composed of 4 stations, b and e for the permanent seismic network (9 stations), and c and f for the permanent plus mobile networks (14 stations). As expected, the detection threshold decreases when approaching the stations themselves. The detection performance of the 4 stations seismic network was high (M d = 0) in the northern area of the island, leaving the The permanent seismic network, whose deployment terminated at the end of October 2017, is able to detect shallow earthquakes of magnitude M d > 0.5 in the whole island (panels b, e); again, the best coverage remains in the northern area (where the historical main seismicity is concentrated). Note that this network also includes a station on the island of Procida, not existing before. The integrated network (permanent plus mobile stations) has a very high coverage in the epicentral area of the 2017 Casamicciola earthquake in order to identify any aftershocks. Figure 8g-l shows the minimum magnitude to locate an earthquake with a hypocentral depth of 500 m (panels g, h, and i) and 1500 m (panels j, k, and l) b.s.l., picked at 4 seismic stations at least, according to the different network configurations. The localization level of the current permanent seismic network is around M d = 1, reaching M d = 0-0.5 for very shallow earthquakes (depth 500 m) in the central-northern sector of the island. Before the upgrade of the network (concluded in 2018), the location threshold was higher than M d 1.5 on the whole island.

STATISTICAL CHARACTERIZATION OF THE SEISMICITY
Merging the earthquakes of the macroseismic catalogue (up to 1992) with the ones of the instrumental catalogue (from 1993) leads to a unified catalogue of 252 earthquakes, covering the timespan from the 8th century BC to the end of 2019. Among them, 78 events are the ones above the completeness magnitudes. The merged catalogue is reported as Supplementary Dataset 4. This represents the most complete and extended catalogue available for Ischia to date.
With this catalogue, we can attempt to characterize the Ischia seismicity from a statistical point of view, analysing spatial and magnitude-frequency distributions of the earthquakes, and exploring the stationarity process of seismicity.
The completeness of the unified catalogue is defined according to the results of the historical completeness analysis for the macroseismic catalogue (Section "Historical Completeness") and the completeness based on network sensitivity (Section "The Sensitivity of the Instrumental Network and Completeness"), adopting for both the moment magnitude M w ( Table 3). In Figures 5A,B, the variation of completeness through time is compared with the magnitude estimations reported in the catalogue. Given the small number of events and the difficulty in robustly defining completeness, the impact of the uncertainty on the parameter estimations (b-value and annual rates) is investigated.

Magnitude-Frequency Distribution
First, we estimated the magnitude-frequency distribution of the earthquakes assuming a Gutenberg-Richter (GR) law (Gutenberg and Richter, 1944) with a taper for the strongest events. This tapered GR (Vere-Jones et al., 2001;Kagan, 2002) has a soft bound for the highest magnitude, differently from the hard bound of the classical truncated GR, which instead assumes a maximum magnitude that cannot be exceeded. From a statistical standpoint, we prefer the tapered version of the GR law because it is not possible to estimate the maximum magnitude value of the truncated version of the GR in a proper manner solely based on seismic data, even with a catalogue lasting centuries (Holschneider et al., 2011;Geist and Parsons, 2014).
The tapered GR distribution is described by the equation (Kagan, 2002): where F(M) is the cumulative distribution function for the seismic moment M(F = 1-S, survivor function reported by Kagan, 2002); M t is the minimum moment; β is the parameter controlling the slope of the distribution, and M cm is the corner moment that rules the tapering of the right tail of the distribution.
Since we adopt M w in our catalogue, we convert the seismic moment in the previous equation using the Kanamori (1977) formula. Comparing with more commonly used parameters and definitions, M t corresponds to the completeness magnitude, M cm to a "corner magnitude" over which the magnitude-frequency distribution decays quicker than a GR, and β corresponds to 2/3 of the classical b-value. For simplicity, in the following we refer to the classical parameters in terms of magnitude and b-value. To estimate these parameters, i.e., the cumulative annual rate of events (λ), the b-value, and the corner magnitude, in a catalogue with a time-varying magnitude of completeness, we used the method proposed by Taroni and Selva (2021). This approach couples the maximum likelihood estimation (MLE) proposed by Weichert (1980) and adapted for the tapered GR, with a Monte Carlo Markov Chains (MCMC) computation to properly explore and estimate the uncertainty associated with the parameters, allowing for a joint evaluation of potentially correlated parameters (Keller et al., 2014). The smallest explored magnitude (the minimum among completeness levels) is magnitude M w = 1.0. Notably, this method allows accounting for time-variable magnitude completeness and it considers in input the annual rates observed in each magnitude bin (including no observations), evaluated according to the estimated completeness for this magnitude (longer time intervals are available for the larger magnitudes). Consequently, this method is less sensible than other ones (like the classical MLE) to the small magnitudes recorded only in the most recent part of the catalogue. Figure 9A shows the observed and estimated annual rates of the events in the seismic catalogue, while Figures 9B-D display the uncertainty associated with these parameters using 10 4 sampling from the MCMC computation. Figure 9E represents the scatter plot of the joint estimation of the annual rate and the b-value, which are correlated to each other. Overall, the b-value is significantly larger than 1 and it can be constrained in the range 1.0-1.3, with the best guess value of 1.11. The annual rate λ of M w ≥ 1.0 is in the range 4-8 events per year, with the best guess of 5.54 /yr. The corner magnitude uncertainty distribution is almost flat, showing that it cannot be well constrained by the data and demonstrating that the parameter that describes the right tail of any GR distribution is difficult to constrain statistically (Holschneider et al., 2011).
Since several events are close to the magnitude of completeness, we verify the robustness of our results for the uncertainty on its definition. To this end, we simulated 10 4 different sets of completeness, by adding a Gaussian random error with zero mean and a standard deviation of 0.2 (this value can be considered suitable to model the completeness uncertainty) as graphically shown in Figures 5A,B (dashed red lines correspond to +/− 1 standard deviation bounds). The results in Figures 9F-H show the 2D-histogram of the joint estimation of annual rate and b-value, as well as the marginal distribution (i.e., the 1D-histogram) of the annual rate and b-value. Comparing the different panels in Figure 9, we note that the uncertainties resulting from the parameter estimation with an MLE approach are comparable with the uncertainties considering a Gaussian error (with standard deviation of 0.2) on the evaluation of the magnitude of completeness. This means that the obtained results for the magnitude-frequency distribution are robust, being not critically dependent on the selected completeness magnitudes. The strong correlation between the parameters leads to a non-centred joint MLE with respect to the 2D-histogram, whereas the same MLE is compatible with both marginal distributions.
Given the strong correlation between a-and b-values and that the Poisson hypothesis used in Weichert (1980) approach may be challenging in Ischia (see next Section), we further check the stability of results by estimating the b-value also independently from the annual rate. To this end, we use the method described in Taroni (2021) that allows estimating the b-value with a timevarying magnitude of completeness, including the correction for the binning of the magnitudes (in our case 0.1) and the correction for an unbiased MLE (Marzocchi et al., 2020). We obtain a b-value = 1.18 +/−0.13, very similar to the one previously estimated.
Overall, the characterization of the magnitude-frequency distribution appears reliable, and the results, in particular the b-value > 1, appear in line with analogous analysis in other volcanic areas (Vilardo et al., 1991;D'Auria et al., 2013). This value differs from the one estimated by D' Auria et al. (2018). The difference is partially due to the adopted magnitude (M d instead of M w ), whose conversion has a slope smaller than 1 (0.82, Petrosino et al., 2008; see Figure 4A). In addition, the quantification method (Ogata and Katsura, 1993) adopted in D' Auria et al. (2018) estimates, at the same time, the parameters describing both the complete and the incomplete part of the catalogue, fact that, in case of a catalogue with a low number of events (<100), can lead to unstable results. In any case, if we consider the b-value estimation and the associated uncertainty, 0.75+/−0.13, a larger b-value, e.g., 1.1, is still compatible with their findings (i.e., inside the 99% confidence interval). Considering that the largest earthquakes occurred in pre-instrumental times while most of small magnitude events are complete in instrumental times, only the future seismicity detectable by the updated seismic network operating since 2018 will give the possibility to further test our findings (e.g., b-value >1) with a more homogeneous dataset.
Notably, while quite large uncertainties exist in magnitude estimations, in this analysis we did not explore this uncertainty because most likely it has a quite complex structure that could influence the results. For example, the errors on magnitude are probably correlated in the macroseismic part of the catalogue but independent in the instrumental part. This complex uncertainty structure, if not satisfactory described, may bias the results, and may be the topic for future works.

Is the Ischian Seismicity a Stationary Poisson Process?
Since the catalogue shows periods with different seismicity rates, we also estimated the annual rate and the b-value of the tapered GR separately for the two sub-catalogues: the first one spans 1750-1884 and the second one 1885-2019, with a length of 135 years each. The goal is to understand if the higher seismic rate observed in the past (until the 1881-1883 earthquake sequence) is statistically compatible with the present low seismic rate associated with the development of the early seismic monitoring system to date (Luongo et al., 2012). In fact, the observed annual rate of the largest events (from M w 3.6, complete from 1750) in the period 1750-1884 is 7 times higher with respect to the one in the period 1885-2019 (7 and 1 events, respectively).
Note that the 7 events that occurred during the timespan 1750-1884 are not all independent of each other: indeed, applying a classical declustering method to the catalogue (Gardner and Knopoff, 1974), one event is removed.
Assuming a stationary Poisson process, we computed the probability to observe 6 (or more) events in 135 years, i.e., the number of observed events in the first declustered sub-catalogue, using the annual rate and the b-value estimated in the second sub-catalogue (λ = 6.85 /year of M w ≥ 1.0, b-value = 1.34).
The computed probability corresponds to the p-value of a binomial test for the number of observed events (Taroni et al., 2017); we used 10 4 pairs of annual rate and b-value coming from the MCMC estimation in the second sub-catalogue in order to take into account the uncertainty in the parameter estimation. Since 3 over 6 events are near the completeness threshold (see Figure 5A), we also perform the binomial test using only 3 events, in order to check the robustness of our findings. The results are shown in Figure 10A. We obtained a large majority of low probabilities (<0.05), demonstrating that the Poisson hypothesis for the seismic events' distribution can be rejected, independently from the existing uncertainty on the GR parameters or the magnitudes. In practice, the seismicity of Ischia described by our integrated catalogue is a non-stationary process, and significant modulations in the seismogenic process should be invoked to justify the observed long-term oscillations of the seismicity rate.
While the non-stationarity process could not be a surprising feature in volcanic seismicity, it is not so obvious in a volcanic system that did not experience any eruption in the last 700 years (last eruption occurred in 1302 AD) as well as volcanic unrest episodes in recent times (Selva et al., 2019). On the other hand, this significant non-stationarity will challenge the assessment of seismic hazard, as the available data are sufficient to demonstrate that the Poisson hypothesis, commonly adopted in the long-term analyses, does not hold in Ischia.  (B), with σ = 1.5 Km. The colour scale represents the spatial density of events (red higher, blue lower). White dots are reported for all the events in the catalogue, while black circles for events that are also mainshocks, according to the Gardner and Knopoff (1974) declustering method.

Spatial Distribution of the Seismicity
To analyse the spatial distribution of the earthquakes, we built a very high-resolution model of smoothed seismicity by using cells of 0.005 • × 0.005 • . Instead of Frankel's classic smoothed seismicity method (Frankel, 1995), we implemented the innovative method proposed by Hiemer et al. (2014), where the Gaussian smoothing kernel is multiplied by a function that gives more emphasis to the strong past events that occurred when the magnitude completeness of the catalogue was higher, to compensate the lack of low magnitude events in the catalogue. In this way, we can base the smoothed seismicity model catalogue using all the events in the catalogue, with a completeness that varies through time, instead of using the most recent seismicity only. This is particularly important for Ischia, due to the rather limited number of events in the catalogue.
The method uses the following Gaussian kernel: where K ij is the contribution to the j-th spatial cell of the i-th earthquake in the catalogue; r ij is the distance between the i-th event of the catalogue and the centre of the j-th spatial cell; and σ is the so-called correlation distance which regulates the smoothing. To obtain the total spatial rate of the j-th spatial cell, the contributions of all the N earthquakes in the catalogue must be summed:

Figures 10B-D
show the spatial distributions obtained, adopting different σ (0.5, 1, and 1.5 km). We preferred to avoid any optimization procedure, since the total number of data in the catalogue is not large enough to produce a robust inversion of this parameter. The catalogue is declustered adopting the Gardner and Knopoff (1974) method, to avoid a fictitious concentration of the spatial rate where past sequences occurred. The maps, having the same logarithmic scale, represent the spatial probability density function of the events, i.e., the sum of the values in all the cells is 1.
These results show that the most widely adopted method to analyse the spatial distribution of earthquakes (smoothed seismicity with typical settings) recognizes the importance of the well known Casamicciola seismogenic source area, in the northern sector of Mt. Epomeo. Indeed, this small area generated almost all the largest magnitudes in the catalogue (1828, 1881, 1883, and 2017). This area is characterized by E-W faults limiting toward north the most uplifted part of the resurgent block of Mt. Epomeo (Vezzoli, 1988;Tibaldi and Vezzoli, 1998;Funiciello, 1999, 2006;Trasatti et al., 2019), and these events seem to concentrate in the western part of this seismogenic area, possibly due to a more ductile behaviour of rocks in the eastern sector, where most of the volcanic activity took place (Cubellis and Luongo, 1998;Carlino et al., 2010;Selva et al., 2019;Cubellis et al., 2020). However, a non-negligible probability area extends further south, including the majority of the most uplifted part of the resurgent block of Mt. Epomeo (Selva et al., 2019;Trasatti et al., 2019), especially in its western and the central sector. In this area is also located the 1863 event, the only large event showing rather constant intensity values throughout the island and, consequently, with a probable deeper origin. This means that this larger area cannot be completely neglected as a potential source of future seismicity. Noteworthy, adopting the classical Frankel's (1995) method and/or not declustering, we obtain similar results, even if the classical approach -equal weight to all the events in the catalogue -creates a more homogeneous spatial distribution on the island and less emphasis to the stronger past events, mainly clustered in the north-western sector of Ischia.

CONCLUDING REMARKS
This manuscript describes the statistical characterization of seismicity occurring within the volcanic island of Ischia (Italy), based on the analysis of a newly developed and exceptionally long seismic catalogue of the seismicity generated within Ischia. This catalogue is based on an extensive screening of all the information available in the literature, from historical records to instrumental data. This work has been accompanied by a significant effort toward the quantification of existing uncertainty, as well as data completeness. The produced catalogue significantly extends the pre-existing catalogues, producing an unprecedented reach dataset covering several centuries, as summarized in Table 5. This allowed characterizing the statistical properties of the Ischia seismicity, revealing some peculiar properties like, for example, its significant non-stationarity also in a period of no eruptions or unrest episodes.
This study not only provides significant insights into the knowledge of the seismicity of Ischia and its related hazard, but also introduces significant novelties into the quantification and use of the uncertainties in the earthquake catalogues and their statistical characterization. This is particularly important whenever a limited number of earthquakes is available, as in the case of Ischia. The characterization and the management of uncertainty are based on the extensive use of the ensemble modelling, as well as on the developments of tests that quantify the robustness of the statistical characterization accounting for existing uncertainties.
More specifically, the main achievements of this study can be summarized as follows: • Overall, the macroseismic data well describe the seismic history of Ischia and its seismic style, characterized by isolated events, small swarms with few low-energy events concentrated in few months, and sequences with a destructive mainshock accompanied by some minor foreand aftershocks. The events which can be parameterized (location, magnitude) are 57, 16 of which above the damage threshold (I 0 > V-VI MCS). They are being included in the Historical Archive of Historical Earthquake Data (ASMI) and will be used to update the national earthquake catalogue CPTI15 (currently including only 12 events). The standardization of intensity data represents a methodological aspect of interest for future applications. • The integration of results from different studies through ensemble models is also an innovative approach to quantify the epistemic uncertainty of parametric data, providing more realistic uncertainty bounds than any single individual procedure (Taroni et al., 2014;Garcia-Aristizabal et al., 2020). The application to Ischia to the 1828 and 1883 events highlights good compatibility of the results from alternative studies, especially for the 1883 event. Beyond the specific intensity estimates for a given site, we are now aware that the relative difference in the earthquake parameters due to subjective factors (selection of localities, interpretation of historical sources) or objective ones (the most important of which is the building vulnerability) are relatively minor. This is particularly important for the earthquakes having a strong impact on the seismic hazard at the local scale. • The ensemble approach has been also adopted to merge a set of intensity-magnitude relationships, given the lack of data to derive a specific relationship for Ischia. The values of moment magnitude M w obtained for the destructive earthquakes of 1881 and 1883 (epicentral intensity I 0 IX and XI, respectively) are now 4.4 and 5.2, significantly higher than the ones reported in CPTI15 (4.1 and 4.3, respectively). Overall, our results are consistent with the seismogenic volumes proposed for Ischia in the literature (Carlino et al., , 2010Castaldo et al., 2017;Cubellis and Luongo, 2018;Selva et al., 2019;Trasatti et al., 2019;Cubellis et al., 2020;Carlino et al., 2021). • The extension of the instrumental catalogue as far back as 1991 confirms the very low rate seismicity in recent times. The uncertainty on the earthquake parameters, especially the location, is fairly large until October 2017, when the improvement of the seismic network allowed a better definition of the parameters, although a decisive step would derive only from more constrained velocity models, potentially including the entire area Ischia-Procida-Campi Flegrei to take advantage also of the Campi Flegrei seismic network. • By merging the historical and instrumental datasets, we obtained the most complete earthquake catalogue available for Ischia to date, representing a solid base to analyse the statistical features of seismicity and their temporal variations. Regarding the estimation of the parameters of the magnitude-frequency distribution, the best guess b-value obtained from the unified catalogue is 1.11, a value steadily larger than 1. This result is robust also considering the uncertainty in the MLE method or the completeness threshold. This in line with other volcanic zones in Italy and worldwide (Wyss et al., 1997;Murru et al., 1999), appearing also consistent with the high-temperature gradient beneath the seismogenic volume Castaldo et al., 2017;Trasatti et al., 2019;Cubellis et al., 2020), but different from what obtained for Ischia, using only instrumental seismicity, by D' Auria et al. (2018). • The seismicity is not limited only to the well known seismogenic area of Casamicciola, but is extended to the central part of the island, especially along the faults surrounding the most uplifted part of the resurgent block of Mt. Epomeo (e.g., Selva et al., 2019;Trasatti et al., 2019), in particular in its western sector. While a decrease of the uncertainty in the location estimates may better constrain modulations in its spatial distribution, the extension beyond the Casamicciola area is confirmed beyond uncertainty from both the spatial distribution of the instrumental seismicity and the results of standard statistical analyses applied to the entire revised catalogue. Thus, this result is fairly stable and cannot be considered an artefact of the uncertainty on historical or instrumental locations. • The intensity distributions of the largest events (1828, 1881, 1883, and 2017) show an elongated trend in the E-W direction and a rapid decrease with distance indicating a shallow sources. The E-W trend coincides with local fault system, probably indicating the potential contribution of source directivity effects in the spatial distribution of intensities (Vezzoli, 1988;Alessio et al., 1996;Funiciello, 1999, 2006;Carlino et al., 2010;Nappi et al., 2018Nappi et al., , 2021Trasatti et al., 2019;Cubellis et al., 2020). • The observed Ischian seismicity significantly deviates from a stationary process, also taking into account the uncertainty in data. The exceptional higher rates of earthquakes with M w ≥ 3.6 in the years 1750-1884 cannot be explained with the parameters estimated in the years 1885-2019; then a stationary Poisson process is not suitable to describe the Ischian seismicity independently from the application of declustering algorithms. Thus, observations suggest that significant modulations in the seismogenic process have occurred, leading to significant variations of the seismic rate through time.
• In this study, we have considered only events with local origin (within the island of Ischia). External sources may of course contribute to the local seismic hazard. The national seismic hazard map locates Ischia within a single volcanic source area that includes the entire Neapolitan volcanic district and extends to the Apennines, and it reports a relatively low seismic hazard (Selva et al., 2019). The main Apennines sources do not report damages in Ischia, as for example for the 1805 M w 6.6 Matese earthquake or the 1980 M w 6.9 Irpinia earthquakes (Intensities VI and V in Ischia, respectively; Esposito et al., 1987;Gaudiosi et al., 2020). Seismo-volcanic events at Vesuvius and Campi Flegrei are barely felt in Ischia, and mostly are not even felt (Branno et al., 1984;Cubellis et al., 2007).

DATA AVAILABILITY STATEMENT
The macroseismic data produced in this manuscript will be integrated in the same databases (ASMI, CPTI, and DBMI) in their next releases. The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
JS coordinated the overall work, with the help of RA for the macroseismic catalogue, of AT for the instrumental catalogue, and MT for the statistical analyses. JS, RA, GA, CC, EC, SP, and AR developed the macroseismic catalogue of Section "Macroseismic Catalogue: 8th Century BC -2019." JS, AT, MC, DL, and PR participated in the development of the instrumental catalogue of Section "Instrumental Catalogue: 1993-2019." JS, MT, and AT contributed to the development of the statistical analysis of Section "Statistical Characterization of the Seismicity." All the authors reviewed and approved the final manuscript.

FUNDING
This work benefited of the agreement between Istituto Nazionale di Geofisica e Vulcanologia and the Italian Presidenza del Consiglio dei Ministri, Dipartimento della Protezione Civile (DPC). This paper does not necessarily represent DPC official opinion and policies.