A common source for the destructive earthquakes in the volcanic island of Ischia (Southern Italy): insights from historical and recent seismicity

The island of Ischia, located in the Gulf of Naples, represents an unusual case of resurgent caldera where small-to-moderate magnitude volcano-tectonic earthquakes generate large damage and catastrophic effects, as in the case of 4 March 1881 (Imax-VIII-IXMCS) and 28 July 1883 (Imax X-XI MCS) historical earthquakes, and of the recent 21 August 2017 MW = 3.9, event. All these earthquakes struck the northern area of the island. With about 65,000 inhabitants, Ischia is a popular touristic destination for thermals baths, hosting more than 3,000,000 visitors per year, thus representing a high seismic risk area. Assessing its seismic potential appears a fundamental goal and, to this end, the estimate of the magnitude of significant historical events and the characterization of their source are crucial. We report here a reassessment of historical data of damage of 1881 and 1883 earthquakes to evaluate the main source parameters of these events (obtained with the BOXER and EXISM software) and quantitatively compare, for the first time, the results with the source characteristics, obtained from instrumental data, of the recent 2017 earthquake. The results allowed us to assess the location, as well as the possible dimension and the related maximum magnitude, of the seismogenic structure responsible for such damaging earthquakes. Our results also provide an additional framework to define the mechanisms leading to earthquakes associated with the dynamics of calderas.


Introduction
The seismicity of active volcanic areas is generally characterized by low magnitude earthquakes (McNutt and Roman 2015). Nevertheless, the occurrence of moderate (4 < M < 5) volcano-tectonic events is not rare and can generate damage and fatalities (Zobin 2012), mainly due to their shallow hypocentral depth (e.g. Convertito and Zollo 2011). Volcanotectonic earthquakes are mostly associated with the magma dynamic but, in some cases, could be not straightforwardly correlated to primarily magmatic processes. A well-known case is the active volcanic island of Ischia (located in the Gulf of Naples), where recent and historical earthquakes have caused heavy damage and thousands of fatalities (Cubellis and Luongo 1998;Carlino et al. 2010) (Table 4, Appendix) (Fig. 1). These events appear to be associated with a phase of subsidence of the central part of the island (Trasatti et al. 2020).
Ischia is a 46 km 2 island located a few km west of Naples. It is an active resurgent caldera, whose central part underwent large uplift, about 900 m since ~ 55 ka (Sbrana et al. 2009), and subsidence in historical time (Buchner et al. 1996;Manzo et al. 2006). The resurgence generated the Mt. Epomeo block, an approximately ~ 2 × 2 km 2 squared structure, bordered by a system of faults mainly oriented NW-SE, E-W, and N-S (Fig. 1). The resurgence was accompanied by volcanic activity outside the block (last eruption in 1302) and produced the exhumation of a high-temperature hydrothermal system, with geothermal gradients > 150 °C km −1 (Vezzoli 1988;Sbrana et al. 2009;Sbrana and Toccaceli 2011;Carlino 2018).
During the nineteenth century, the island of Ischia was a very important site for the Earth Science scholars and one of the most famous places for spas in Europe (Carlino, 2019). For these reasons, many accounts and reports were produced after the earthquakes of 4 March 1881 and 28 July 1883. In particular, many coeval papers, technical reports, and various accounts described the damage distribution of these two events (see Cubellis and Luongo 1998, and references therein), leading to the assessment of macroseismic intensity of I 0 VIII-IX MCS and I 0 X-XI MCS, respectively, for the 1881 and the 1883 earthquakes (Cubellis and Luongo 1998;Cubellis et al. 2004;Carlino et al. 2010). Intensity data (EMS98) for both earthquakes are included in the Parametric Catalogue of Italian Earthquakes (CPTI15, Rovida et al. 2019). In particular, the CPTI15 lists 17 (11 on the island) and 27 (18 on the island) intensity-points (localities) for the 1881 and the 1883 earthquakes, respectively, for which the epicentral intensity (I 0 ) IX EMS98 (1881) and X EMS98 (1883) are reported.
With more than 2300 victims and the whole destruction of the town of Casamicciola Terme, the 1883 event was the most devastating (Carlino et al. 2010). Considering the small-to-moderate magnitude inferred from previous studies (e.g. 4.6 ≤ M ≤ 5.2; Cubellis and Luongo 1998), the high damaging level of this earthquake can be primarily related to the shallowness of the seismogenic source, which is located between 1 and 2 km of depth (Cubellis and Luongo 1998;Carlino et al. 2010).
As for the whole of Italy (Crescimbene et al. 2015), seismic risk has not been adequately perceived by the inhabitants of the island, until 21 August 2017 when a M d = 4.0 (M W = 3.9) earthquake stroke the island. The occurrence of this recent earthquake, after 134 years of almost complete seismic silence, has brought out again the problem of the relatively high seismic risk for the island (Briseghella et al. 2019;De Natale et al. 2019;Marotta et al. 2019). Based on the current national seismic hazard map provided by the Istituto Nazionale di Geofisica e Vulcanologia (http:// esse1. mi. ingv. it/), for the whole island the acceleration with probability of exceeding equal to 10% in 50 years ranges between 0.125 and 0.175 g (being g the gravity acceleration and the values corresponding to the 50th percentile) with a range of variability (0.1, 0.2) g (corresponding to the 16th and 84th percentile, respectively). In fact, given the grid at which the national hazard map is computed there are only two points for the entire island. In this regard, the 2017 earthquake was the first significant event at Ischia recorded instrumentally, for which a PGA of 0.25 g was recorded at about 800 m from the epicentre (De Novellis et al. 2018). Although the comparison between the recorded value and the expected values indicated in the national seismic hazard should take into account the probabilistic nature of the latter and its reference to bedrock soil condition, we note that the recorded value is larger than the expected value at the 84th percentile.
The source mechanism of the 2017 earthquake has been differently interpreted by various authors (De Novellis et al. 2018;Braun et al. 2018;Calderoni et al. 2019;Nazeri et al. 2020). As for the 1881 and 1883 events, albeit to a limited extent, this last earthquake produced serious destruction in spite of its modest magnitude, causing two fatalities, more than 40 injuries and serious damage in the upper part of the town (Nappi et al. 2018). Before its occurrence, only about 50 very-low magnitude earthquakes (M < 2.3) had been recorded since 1999 by the seismic network installed on the island. Despite the seriousness of the situation, the 2017 event gave us the opportunity to deepen our understanding of seismic processes occurring in the island, comparing the latter event to the historical ones, which seem to have similar characteristics.
In this study, we compare for the first time the macroseismic data of two main historical earthquakes of Ischia (1881, 1883) to those of the first instrumental earthquake occurred in the island in 2017. This represents a first important step to quantitatively compare macroseismic and instrumental events for the study of the seismic source of this volcanic area. We use uniform and standard criteria to represent the macroseismic field and to retrieve the fault parameters of the three events, in order to analyse the location and the kinematic of the seismogenic structures generating highly damaging earthquakes in the island, verifying their relation and/or their possible coincidence. In particular, we elaborated the intensity data of both the 4 March 1881 and the 28 July 1883 earthquakes, which are mainly inferred from the analysis of historical archive reports. We reassessed the data for the 1881 and 1883 earthquakes following macroseismic criteria and using the EMS98 scale, to obtain macroseismic fields. We inferred new magnitude estimations and retrieved the geometry of the faults using different computational methods: BOXER (Gasperini and Ferrari 2000;Gasperini et al. 1999) and EXSIM (Motazedian and Atkinson 2005;Boore 2009;Assatourians and Atkinson 2012). Then, we compared the results to the fault plane solution and the damage of the 21 August 2017 earthquake. The results suggest the presence of a unique seismogenic structure, which can be periodically reactivated during the present quiescent phase of the volcano, providing new elements to assess the processes that could generate earthquakes in subsiding calderas. Noteworthy, the fault geometries obtained in the present study can also help to refine the seismic hazard estimates for the island, in particular for those approaches based on single faults rather than areal seismic source zones (e.g. Akinci et al. 2009;Convertito et al. 2006;Pace et al. 2006).

Historical Data
The description of the damage of the 1881 and 1883 earthquakes and the felt report data had been previously reported in various papers and databases (Cubellis and Luongo 1998;Cubellis et al. 2004;Carlino et al. 2010;Rovida et al. 2019). Those data had been 1 3 elaborated by various authors to provide an estimation of epicentral intensity, epicentral location, and intensity-related magnitude (see for instance Cubellis and Luongo 1998;Carlino et al. 2010, andRovida et al. 2019).
The macroseismic intensities of the 1883 event published by Cubellis and Luongo (1998), Cubellis et al. (2004) and Carlino et al. (2010) have been obtained by assigning the intensity value to very small localities, represented by a group of edifices located in a single street or district. This procedure does not match with the standard macroseismic method (see, for instance, Grünthal 1998, Musson et al. 2008, because in many cases the intensity was assigned according to damage suffered by a too low number of edifices that is not representative for an intensity estimate. On the other hand, in the CPTI15 , the number of intensity points represented on the Ischia Island for each one of the three considered earthquakes differs significantly (11 points for the 1881, 17 for the 1883, 24 for the 2017). In particular, the CPTI15 reports a lower number of points in the highest damage (and densely inhabited) areas for the 1881 and 1883 earthquakes with respect to 2017. Given the very shallow hypocentre depth of these events (generally associated with very fast spatial variation of the intensity) and pointing at retrieving more information on their source characteristics, our first intent was to homogenize the areal coverage of the intensity points of 1881 and 1883 events, to get a distribution comparable with that of 2017. In particular, the macroseismic intensities for the 1881 and 1883 earthquakes can be reviewed on the basis of the large amount of available data.
In general, the analyses of seismic intensity data represent the only way to retrieve information about the location and the dimension of the seismogenic source of seismic events occurred in the pre-instrumental era (e.g. Panza and Cuscito 1982;Gasperini and Ferrari 2000;Gasperini et al. 1999;Tertulliani et al. 2012). However, these methods are affected by a number of uncertainties, also due to the local geological conditions, affecting the seismic wave propagation and amplification, and to the different seismic response of the buildings. However, the above uncertainties can be reduced in the case of the 1881 and 1883, in consideration of the quality and quantity of available historical data.
Considering the moderate magnitude of the events, the observed very fast decay of the damage level (that is, the fast spatial variation of intensities) would require a higher level of detail of the intensity pattern with respect to larger earthquakes, to retrieve information on the source characteristics. Thus, we collected the original descriptions of the damage, for both the earthquakes, and increased the number of data points in the northern zone of the island, which was affected by larger damage.
For the 1883 event, we used the original description of the damage, derived from the consultation of the report of the Comitato Centrale di soccorso pei dannegiati dell'isola d'Ischia (1885) (Rescue Committee for the damage of the island of Ischia), which represents the most detailed data source for the damage caused by this event. Immediately after the events, the Committee provided technical forms with the description of damage suffered by each building (with indication of the street and, when available, of the building number) (see also Cubellis et al. 2004 andCarlino et al. 2010).
We grouped the original damage data by location, prior to assessing intensity, as recommended by Musson et al. (2008). We used the criteria for which a location is represented by multiple settlements (MS) , which is a significant number of buildings grouped in a well-known site. The locations are distributed within the six different municipalities of the island (Casamicciola, Lacco Ameno, Forio, Serrara Fontana, Barano, Ischia). Apart from a few exceptions regarding non-damaged zones, we considered only the MS having a number of buildings > 50 units. The report of the Comitato Centrale pei Dannegiati dell'isola d'Ischia (1885) also includes damage data of the island for the 1881 earthquake. However, due to the limited information, we also consulted further literature data (see Cubellis and Luongo 1998, and references therein). We adopted the same criteria used for the 1883 event in the identification of the localities, finally obtaining intensity data for the 1881 earthquake.
In the end, damage data for about 3000 edifices (total) located in the island were analysed. Furthermore, surface effects, such as major fractures, open-cracks and landslides have been also considered in the evaluation of the intensity, taking into account the contemporary chronicles and technical reports (see Cubellis and Luongo 1998, and references therein). On the base of the damage descriptions, we evaluated the EMS98 intensity for the 1881 and 1883 earthquakes, obtaining intensity data for 20 (Table 5, Appendix) and 26 ( Table 6, Appendix) locations on the island, respectively, assigning a name to each MS on the basis of the present toponymy (Grünthal 1998). In Fig. 8 of Appendix, we also report the different intensity maps (1881 and 1883) in which the distribution of our intensity points is compared with those reported in the CPTI15 .
As for the vulnerability classification (Grünthal 1998) for Ischia Island, apart for a few noble buildings, at the end of nineteenth century, the building heritage of the island was mainly made by the same local tuff masonry, with local rough stone and/or squared stony blocks (Polverino 1996(Polverino , 1998Cubellis and Luongo 1998), all belonging typologies to vulnerability class B. Also, recent studies report that (at least in the most damaged area) the building heritage erected before the 1919 mainly falls in the vulnerability class B and C and only in minor part to class A (~ 12% and ~ 25%, respectively, for the Lacco Ameno and the Casamicciola Terme municipalities; Verderame et al. 2017;Del Gaudio et al. 2018). The poor building typology might have increased the overall damage. In fact, many locations of the northern sector of the island (between Casamicciola Terme and Lacco Ameno municipalities) fall in the damage classification 3 to 5 (damage from "significant" to "very heavy"), while the damage were heavy at the epicentre (Casamicciola Terme) for both the earthquakes (1881, 1883). In the municipality of Casamicciola Terme, 250 edifices out of about 800 suffered heavy damage during the 1881 event, and about 670 were collapsed during the 1883 earthquake (Comitato Centrale pei dannegiati dell'isola d'Ischia, 1885) Furthermore, the similarity of the poor architectural technique in the whole island, possibly contributed to a general increase of damage, and not in its difference in the spatial distribution (Polverino 1998). Finally, it should be stressed that the damage produced by the 1881 earthquake had probably weakened the buildings stock, contributing to an increase of the effects during the following 1883 event. Anyway, due to the lack of punctual information of the 1881, this cannot be quantified in terms of increase of vulnerability. Fig. 2 depicts the macroseismic field of 1881 and 1883 earthquakes, obtained from the interpolation of the intensities data listed in Appendix Tables 5 and 6, and the 2017 macroseismic field from CPTI15 ) data (we used the kriging method for the interpolation of data; Kerry and Hawick, 1998).

Direct modelling of the 2017 earthquake
To infer the geometry and the orientation of the faults responsible for the two historical 1881 and 1883 earthquakes, we implemented a procedure similar to that proposed by Convertito and Pino (2014) to study the 1908 Messina Straits earthquake. We first selected a set of parameters for model faults in terms of orientation (strike and dip), Mach number values (i.e. the ratio between rupture velocity and shear-wave velocity), and static stress drop. Then, for each model, we computed peak-ground acceleration (PGA) and peakground velocity (PGV) by using the EXSIM code (Assatourians and Atkinson 2012; Atkinson and Assatourians 2014) at the geographic locations where the observed macroseismic intensities are available. The computed PGAs and PGVs are finally converted to macroseismic intensities by using the relationship proposed by Faenza and Michelini (2010) and compared with the observed ones. The best model corresponds to the one that minimizes the misfit function: In Eq. (1), m is the model parameters' vector, whose components are the length L and the width W of the fault, its strike and dip, the stress drop (Δσ), and the Mach number (α). N is the number of intensity data points, I obs i are the observed intensities and I cal i are the intensities obtained from PGAs or PGVs.
Before analysing the two historical earthquakes, in order the check the performance of the whole procedure, we tested the method on the recent 21 August 2017 (M w = 3.9) earthquake that stroke the same area and for which a reliable model exists for its source geometry and slip distribution (De Novellis et al. 2018). We examined both the direct and the inverse modelling. As for the direct problem, we first computed the predicted intensities from PGAs and PGVs obtained for a realistic source model and checked the results against the observations. Then, as for the inverse problem, we computed intensities for distinct simplified source models and used the observed intensities to determine the source characteristics better simulating the data and compared the best solution with the actual source.
In all the EXSIM computations, we assumed a crustal model characterized by average S-wave velocity value of 1700 m/s, density of 2700 kg/m 3 , and a frequency-dependent anelastic attenuation model Q(f) = 21f 0.6 (being Q the quality factor and f the frequency) (Nunziata and Rapolla 1987;Petrosino et al. 2008;Del Pezzo and Bianco 2013;Capuano et al. 2015). Given the relatively short source-to-site distances considered in this study, a 1/R geometrical spreading is considered (where R is the minimum site distance from the fault). The fault plane was assumed to be rectangular and was subdivided into an appropriate number of 0.1 × 0.1 km 2 sub-faults, which are modelled as point sources characterized by a ω −2 spectrum. The upper left corner of the fault is used as reference fault point (Assatourians and Atkinson 2012; Atkinson and Assatourians 2014). Moreover, in the simulations, we always investigated three stress-drop values (0.1, 1, and 3 MPa), as typical ranges for volcanic areas (e.g. De Natale et al. 1987Natale et al. , 1988 and three Mach number values (0.6, 0.7, and 0.8), as reasonable range for light-to-moderate earthquakes (e.g. Seekins and Boatwright 2010;Convertito et al. 2012).
As for local effects, previous studies suggest that in the northern sector of the island the effect of seismic waves propagation, as well as site amplification, can be considered less important with respect to source effects, because of both the very shallowness of the hypocentre and the relatively short hypocentral distances (Cubellis and Luongo 1998;Gasperini and Ferrari 2000). Also, recent seismic measurements performed in the most damaged area of the 2017 earthquake conclude that H/V spectral analyses performed on seismic noise do not show important peaks related to site amplification (Vassallo et al. 2018). On the other hand, microzonation investigations in three municipalities, in the northern and western areas of the island, determined amplification factors in the range 1 to 3 for the frequency band of interest (1.25-10 Hz) (http:// www. commi ssari orico struz ionei schia. it/ Esiti-Micro zonaz ione. html). Unfortunately, the three districts only cover a limited extent of the island and less than 50% of the data points. On these grounds, we cannot account quantitatively for the site effect, unless assigning arbitrary values to all the sites with unknown site amplification factor. Thus, we simulated PGAs and PGVs at the bedrock.
Firstly, for the direct modelling of the 2017 earthquake, we adopted the fault geometry (L × W = 3.0 × 1.8 km 2 ; fault strike 86° ± 5°; dip angle 70° ± 7°) inferred from the analysis of DinSAR data by De Novellis et al. (2018) and their slip distribution, characterized by a single patch with most slip confined approximately in a 0.5 km-radius circle. We used the observed macroseismic field available in the CPTI15 parametric catalogue , consisting of 24 intensity data points ranging between III and VIII. For both PGA and PGV, we obtained the best solutions for Δσ = 1.0 MPa and α = 0.7 (Fig. 3 and Model M1-A-B in Table 1), and maximum PGA and PGV value of about 0.5 m/s 2 and 0.16 m/s, respectively. Notably, PGA-derived intensities definitely provide a better fit with the hypocentral distance with respect to PGV. This result can be ascribed to the very short sourceto-site distances and to the fact that the considered event is a low magnitude one, with a major relative contribution of high frequencies to the peak-ground motion. Based on this test, we considered the PGAs better able to reproduce the data and, thus, to provide a more robust link between the source fault and the observed intensities.
To test the effects of our assumptions on the local amplification, we modified the PGAs and PGVs obtained for the best model and the "true" slip distribution and fault geometry, considering the relevant amplification factors. In particular, we applied amplification factors in the frequency bands 2-10 Hz and 1.25-2.5 Hz, respectively, to PGA and PGV. As obvious, this could be done only for the sites for which the factor is available (squares in Fig. 3). Remarkably, only for three of the tested sites the resulting intensity is 0.66 larger than the value obtained from PGA and PGV at the bedrock and, most important, the modified values provide comparable (0.120 vs. 0.117, for PGA) or worse (1.2259 vs. 0.268, for PGV) fit to the observed intensities.

Inverse modelling of the 2017 earthquake
In order to investigate how much of a realistic source can be retrieved from the modelling of macroseismic data points, we applied an inverse approach by assuming for the 2017 earthquake the observed intensities. We initially used the CPTI15 epicentre ) and the fault azimuth and extent resulting by using BOXER code (Gasperini and Ferrar 2000;Gasperini et al. 1999Gasperini et al. , 2010, obtaining 80.50° (± 31°) and 0.8 × 1.9 km 2 , respectively. Then, we assumed uniform slip distribution and, by slightly varying the parameters, searched for the best fault by converting the EXSIM predictions to intensity and comparing these latter with the observed ones. In particular, we used the EXSIM results to test the sensitivity of the intensity data points to fault length and also checked the possibility of the whole procedure to retrieve the fault dip. We investigated three fault surfaces extents (1 × 1, 2 × 2, 3 × 3 km 2 ), with the top in correspondence of the Earth surface, and two opposite fault strike directions (80.50° and 260.50°). In the lack of any constraints about the hypocentre, we fixed it at half-length along the strike and 0.2 km above the fault bottom (e.g. Mai et al. 2005). As for the dip angle, the major faults on the island are all vertical or subvertical (Acocella and Funiciello 1999;Vezzoli 1998; Sbrana and Toccaceli 2011) thus, for 1 3 each combination of fault strike and fault extent we computed synthetic intensities for two fault dips, 70° and 90°, with the two strikes giving coincident faults for the vertical planes.
The best model is illustrated in Fig. 4 and listed in Table 1 (MB17). As expected from the previous test, the PGA-derived intensities simulate the data significantly better then PGVs. Both PGAs and PGVs indicate 80.50° as preferred strike and stress drop Δσ = 0.1 MPa. However, while PGAs prefer a small, vertical fault (1 × 1 km 2 ) with Mach number α = 0.7, PGVs are better simulated by using a larger fault (2 × 2 km 2 ) dipping 70°  and with slightly lower rupture velocity (α = 0.6). Considering the assumption of uniform dislocation and the slip distribution derived from the DInSAR analysis (De Novellis et al. 2018), the smaller fault extent appears to be in better agreement with the geodetic solution. Finally, the difference in the fault dip angle between the two preferred solutions can be considered as an estimate of the uncertainty for this parameter. The maximum PGA and PGV values are 0.5 and 0.32 m/s, respectively. Based on the above observations, we deemed the modelling of the PGA-derived intensities more reliable in determining the source geometry and used the PGVs anyway to evaluate the robustness of the preferred solution.

Magnitude estimation of 1881 and 1883 events
In order to proceed with the simulations for the 1881 and 1883 earthquakes, we first assessed the magnitude of these events, which have been used together with the other parameters, as input data of BOXER. Different intensity-magnitude relations have been already used for the 1883 earthquake and provide a magnitude range between 4.6 and 5.2 (Cubellis and Luongo 1998 and references therein). After the 2017 earthquake, we could have a first instrumental assessment of the magnitude for a significant earthquake at Ischia (M d = 4.0, M W = 3.9; http:// cnt. rm. ingv. it), which can be compared with its maximum intensity (I 0 EMS98 VII-VIII) . Thus, among the most used intensity-magnitude relationships for volcanic areas (Patanè et al. 1986;Marturano et al. 1988;Azzaro et al. 2011) (Table 2), we used the ones that provide the most consistent result for the 2017 event. Considering I 0 EMS98 VII-VIII, the intensity-magnitude relationships that provides the closer values to M W = 3.9 (or M d = 4.0) are that derived by Patanè et al. (1986), giving M = 3.8, and the one proposed by Azzaro et al. (2011), giving M d(average) = 3.8. On the other Table 2 Magnitude estimation for the 1881, 1883 and 2017 Casamicciola earthquakes, using different I 0 -M (maximum intensity-magnitude) relationships (Patanè et al. 1986;Marturano et al. 1988;Azzaro et al. 2011) See text for details *According to Gasperini and Ferrari (1990), Gasperini et al. (1999) and Gasperini et al. (2010)  1 3 side, the result of the magnitude-intensity relation of Marturano et al. (1988) (Table 2) provides a difference with respect to the 2017 M d and M W of 0.4 and 0.3, respectively. This difference is larger or at the limit of the uncertainty in the evaluation of the magnitude of 2017 earthquake (± 0.3; http:// cnt. rm. ingv. it), thus we excluded this relation from the following computations. By applying the relations proposed by Patanè et al. (1986) and Azzaro et al. (2011), we obtained M = 4.2 and M d(average) = 4.1, for the 1881 event, and M = 5.0 and M d(average) = 4.8, for the 1883 event, respectively (Table 2). However, we notice that, compared to that inferred by Patanè et al. (1986), the relation proposed by Azzaro et al. (2011) is better constrained, being inferred from an about fivefold number of earthquakes. Thus, in the next steps, we utilized the relation derived by Azzaro et al. (2011).

Retrieved faults parameters for the 1881 and 1881 events
Concerning the 4 March 1881 earthquake, we used the macroseismic intensities obtained in the present study consisting of 20 intensity points ranging between IV and IX. We first used BOXER to infer the best macroseismic epicentre and surface fault projection extent. We obtained 40.744°N (± 0.2 km) 13.903°E (± 0.4 km) for the epicentral location, very close to that obtained for the 2017 event, a best azimuth of 88° (± 30°), and fault extent 1 × 2 km 2 . Next, following the same approach used above for the 2017 event, we computed PGAs and PGVs and converted them to macroseismic intensities. We tested two opposite strikes (88° and 268°), two dips (70° and 90°), three fault geometries (1 × 1, 2 × 2, 3 × 3 km 2 ), three stress-drop values (0.1, 1, 3 MPa), and three Mach number values (0.6, 0.7, 0.8). For the historical events, considering the uncertainty in the BOXER epicentral location, we also tested different reference points, by moving it ± 500 m along the N-S direction and, for each one of these two locations, also shifting it ± 500 m along the direction identified by the azimuth obtained from BOXER. We used uniform slip distribution for all the investigated models. The best models for PGA and PGV are represented in Fig. 5 and listed in Table 3 (MB81A and MB81V). Like for the 2017 event, a lower misfit results for PGA. Both PGA-and PGV-derived intensities indicate a 2 × 2 km 2 fault plane, with Mach number α = 0.6 and stress drop Δσ = 1 MPa. The stress drop is larger than what obtained for the 2017 event, possibly reflecting the higher maximum intensity of the 1881 earthquake with respect to what estimated for 2017. In contrast, different fault dips result for PGA and PGV, respectively, with the former dipping 70° and the latter being a purely vertical plane. The two planes have the same fault centre, which is slightly shifted (500 m) to the south and to the west relative to the BOXER solution. This disparity may result from the different assumption of the two inversion schemes: BOXER uses only the  points corresponding with the major intensities, while the synthetic intensities obtained through EXSIM are compared to all the data points on the island. In fact, the observed macroseismic field of the 1881 earthquake is characterized by a clear asymmetry, with rapidly decreasing intensity moving eastward from I 0 , while higher values are located in the opposite direction, at the same distance. This makes EXSIM procedure to prefer a slightly larger plane, with the centre shifted south-westward. As for the 28 July 1883 earthquake, we used the observed macroseismic intensities obtained in the present study consisting in 26 data points with intensities ranging between VI and XI. The macroseismic epicentre obtained by using BOXER is 40.746°N (± 0.5 km) 13.893°N (± 1.2 km), the best azimuth 86° (± 30°), very similar to those resulting for the 2017 and 1881 events. The surface fault projection extent is an almost squared plane with dimensions 2.5 × 3.4 km 2 , whose larger extension is the N-S direction. Following the same modelling strategy described above, we used EXSIM code and obtained the best models listed in Table 3 (MB83) and represented in Fig. 6. For this earthquake as well, PGA-and PGV-derived intensities indicate distinct values for all the parameters and, also in this case, the PGA simulation provides better fit to the data. The PGA best solution is associated with a 1 × 1 km 2 vertical plane, centred 500 m to the south and to the east with respect to the BOXER epicentre, while a 2 × 2 km 2 fault, dipping north 70°, results for PGV, with fault centre coincident with the BOXER one.
It should be remarked that, although associated with a higher magnitude, the 1883 fault dimension resulted to be smaller than 1881 but with the same stress drop. This result derives from the apparently different intensity distribution for the two events (Fig. 8,  Appendix), displaying a definitely stronger gradient in the epicentral area for the 1883 (Fig. 2). We speculate that this might be related to differences in the slip distribution, possibly smoother for the 1881 earthquakes.
As for the maximum PGA and PGV values for the 1881 event, we obtain 6.7 m/s 2 (0.7 g) and 0.30 m/s, while for the 1883 event 14.3 m/s 2 (1.5 g) and 0.80 m/s, respectively. For both the events, the simulated PGA values, while reproducing the observed macroseismic intensities, significantly exceed the range of values (0.1, 0.2) g (corresponding to the 16 th and 80 th Fig. 4 Best model solution for the August 27, 2017, M W = 3.9, earthquake for uniform slip distribution, using distinct fault reference points. The lower left panel depicts the observed macroseismic data and the surface fault projection obtained by using BOXER (grey box). The black star in all the panels identifies the macroseismic epicentre. The upper left panel shows the intensity values as obtained by converting PGAs in intensities and the inferred fault geometry (black box) (see Table 1), the upper right panel those obtained by using the PGVs. The lower right panels show the observed intensities as function of the hypocentral distance (grey circles) and the computed intensities (crosses) (upper for PGVs and lower for PGAs) (See also Table 1 model MB17-A-V) percentile, respectively) prescribed in the official hazard map as reference for the design of civil structures. It should be noted that the g values obtained from the PGA and PGV values are theoretical and also represent the maximum at the epicentre for the two historical events. In particular for the 1883 event, the g theoretical value is six times the value for the 2017 (0.25 g), recorded 800 m away from the epicentre, while the inferred magnitude is about 1 degree larger than the 2017 event.  also Table 3 model MB81-A-V) 1 3

Discussion
The above elaborations indicate that the seismogenic sources associated with the 1881, 1883, and 2017 earthquakes are located along the northern rim of the Mt. Epomeo resurgent structure with strike roughly in the E-W direction. The fault dimensions resulting for the three events ranges between 1 × 1 km 2 and 2 × 2 km 2 , associated with vertical or almost vertical planes, possibly dipping south. According to the island volcano-tectonics, the system of faults along the northern sector was moved with a thrust mechanism during the resurgence, while the latter was inverted as a normal one, during the subsidence phase (Acocella and Funiciello 1999;Acocella et al. 2001).
Considering the uncertainties in the intensity estimate-either for the observed data points or for the PGA and PGV converted values (± 0.35 and ± 0.26, respectively;Faenza and Michelini 2010)-and in the modelling procedure (propagation and site effects), we consider that, based on the available data, stress drop and Mach number estimates are hardly constrainable with a lower degree of uncertainty. On the other hand, also taking into account the local geology and independent results from analysis of the 2017 earthquake (De Novellis et al. 2018), we conclude that the results about the fault geometry and extent can be considered realistic and reliable within the sampling step of the parameters used in our analysis. The seismogenic source resulting from our study is compatible with the faults system bordering the northern rim of the resurgent structure. The lateral extension of this fault system is structurally limited by the extension of the active part of the resurgent block, about 2 km in length (Acocella and Funiciello 1999;Carlino et al. 2006;Di Giuseppe et al. 2017;Sbrana and Toccaceli 2011). This limit encloses the fault plane (De Novellis et al. 2018) and the observed coseismic fractures (Nappi et al. 2018) of the 2017 events and seems to represent the maximum possible extent of seismogenic sources in the northern sector of the island. Moreover, the maximum depth at which the fragile shear failure can occur is limited by the brittle-ductile transition of rocks, which is located at a depth of about 2 km beneath the Mt. Epomeo (Carlino et al. 2012;Castaldo et al. 2017;Carlino 2018). This depth also corresponds to the cut-off depth of seismicity recorded in the island (D'Auria et al. 2018).
Furthermore, the macroseismic field for the 1881, 1883, and 2017 events shows that the most damaged areas for the three events have similar shape and are strikingly coincident (Fig. 7). The obtained 1881 and 1883 epicentre locations and that of 2017 are very close to each other, with maximum separation of ~ 600 m. These observations suggest a single seismic source as responsible for the destructive earthquakes of Ischia Island, whose reactivation is likely associated with a local stress field variation (Trasatti et al. 2020). The latter hypothesis is also supported by other observations: (i) the absence of significant (M > 1.5) seismic events outside the northern sector of the island; (ii) the absence of regional seismogenic structures crossing the island; (iii) the occurrence of very few aftershocks associated with the 2017 event (less than 30) and localized along the slipped zone (De Novellis et al. 2018); (iv) the presence of faults system around the Mt. Epomeo structure, which does not show a preferential regional pattern orientation, and whose formation is associated with the resurgence process and to the progressive failure of the brittle curst (Sbrana et al. 2009;Carlino 2012;Di Giuseppe et al. 2017).
On the basis of the above observations and of the result of the simulations, the seismogenic fault of 1881 and 1883 does not exceed 2 × 2 km 2 . Given that this source area is actually confined by the structural evolution of Mt. Epomeo and the brittle-ductile transition below the island, we suggest that this specific seismogenic fault can provide a magnitude that is unlikely to exceed the range 5.0 ± 0.3, assuming empirical relations between rupture area and magnitude (Wells and Coppersmith 1994;Somerville et al. 1999).
Furthermore, the results obtained from the BOXER and EXSIM elaborations, which agree with a roughly E-W strike of the seismogenic fault, locate this latter along the most strained zone of the Mt. Epomeo (Acocella and Funiciello 1999), where the maximum load of the resurgent block is exerted during the present subsidence phase. The subsidence possibly took place since ancient Roman Age at least (Buchner et al. 1996), therefore the historical earthquakes likely occurred in a geodynamic context similar to the present one. This makes the vertical loading of the Mt. Epomeo the predominant stress (σ 1 ) of the area (Manzo et al. 2006;De Martino et al. 2011;Castaldo et al. 2017;De Novellis et al. 2018;Trasatti et al. 2020). Moreover, the minimum tensional stress acting at regional scale, is approximately NNW-SSE oriented (Hippolyte et al. 1994), favouring normal faulting along the roughly E-W seismogenic structure of the north of Mt. Epomeo (Lowrie 2007;Zoback 2010).
As regards the causes of the subsidence, two possible main processes have been suggested in the recent literature. The first process could be associated with the depressurization of a shallow magmatic system (Trasatti et al. 2020), the second one foresees the persistence of a stress field associated with the loading of Mt. Epomeo, and the coupling action of volcano loading and crust rheology De Novellis et al. 2018). For both scenarios, the subsidence of the Mt. Epomeo represents the mechanism that accumulated strain energy along the seismic fault and for which only normal fault mechanism can occur.

Conclusions
In this study, we attempted at reassessing the location and magnitude of two important historical earthquakes of Ischia Island (1881 and 1883) and their possible relation with the recent 21 August 2017 event. We analysed the available macroseismic data and derived quantitative information about the source parameters of the three earthquakes, verifying the similarities of their seismic sources. Although intrinsic approximations and uncertainties are associated with the adopted procedure, the results obtained initially by modelling the data of the 2017 earthquake validates our approach, demonstrating its effectiveness in retrieving consistent information on the source of earthquakes in the northern sector of Ischia Island from macroseismic data. Our results suggest that: • a single seismogenic structure is likely to be responsible for the known destructive earthquakes on the island, located on the northern slope of Mt. Epomeo. The plane has a roughly E-W strike, it is dipping vertically or possibly southwards at high angle, as already evidenced by Acocella et al. (1999)  • the data of both historical and recent earthquakes, joined with geological and geophysical information, indicate a fault dimension of about 2 × 2 km 2 . The maximum magnitude for earthquakes in the northern sector Ischia Island cannot exceed ~ 5.0; • the magnitude of the 1881 and of the 1883 events are evaluated around 4.1 and 4.8, respectively, with the latter value possibly affected by cumulative damage; • the simulated PGA values for the 1881 and 1883 events significantly exceed the PGA values, also referring to the 84th percentile, with probability of exceeding equal to 10% in 50 years reported in the official national seismic hazard map.
We cannot exclude that further seismogenic sources can be activated in the island, although this appears an unlikely scenario, in consideration of the thickness of the fragile crust (which is larger in the northern sector of the island) (Carlino et al. 2006) and of the structural dynamic of Mt. Epomeo. Finally, this work highlights the complexity of the processes leading to the seismic energy accumulation and release in such resurgent calderas and the necessity of refining the national seismic hazard map at a smaller local scale-possibly using a single fault-based approach-and integrating recorded and simulated PGA and PGVs values.

Table 4
Historical seismicity in the island of Ischia as reported by Cubellis and Luongo, (1998) Epicentral intensity and magnitude values from CPTI15