Structural setting of the Ischia resurgent caldera (southern Tyrrhenian Sea, Italy) by integrated 3D gravity inversion and geological models Open Access

The reconstruction of the margins and geometry of the resurgent blocks in calderas and the underlying magma bodies is a challenge that, in the scientific literature, has often been faced by following separate paths, starting from surface geological data or from the interpretation of geophysical surveys. This approach has frequently led to disjointed, if not conflicting, interpretations of what can be directly observed, the possible kinematics and dynamics of the uplift, and the characteristics of the magmatic bodies buried at depth. Few attempts at joint interpretation of data from multi-disciplinary studies are known so far, although it has been widely demonstrated that the cross-checking of data from different disciplines and their mutual compatibility can lead to a strongly constrained validation of the interpretative models obtained from them (e.g. for the Neapolitan area: Achauer et al. 2000; Tondi and De Franco 2003, 2006; Di Napoli et al. 2011; Di Vito et al. 2016; Trasatti et al. 2019).

The inversion of gravity data, in particular, has proved to be a very powerful tool for defining the distribution of masses at depth and reconstructing their geometries (Sarkowi and Wibowo 2021a, b; Sarkowi et al. 2021). This approach has also been applied in the Neapolitan areas, through 2–2.5D and 3D inversions and, in particular, also to the island of Ischia (e.g. Berrino et al. 1998, 2008 and references therein; Berrino and Camacho 2008; Paoletti et al. 2009, 2013; Capuano et al. 2013, 2015).

The island of Ischia is one of the most impressive examples of caldera resurgence in the world, with an uplift of about 1000 m, which occurred in less than 50 ka (Selva et al. 2019). The reconstruction of its structure, as well as the kinematics and dynamics of the resurgence process active on the island, have been the subject of numerous studies over the last 30 years (e.g. Selva et al. 2019 and references therein). However, the agreement on a shared interpretative model, which takes into account all available data, is still far from being reached.

The aim of this paper is to define the structural setting of the Ischia resurgent system, from the surface to its deep roots, by a joint interpretation of the results of a 3D inversion of previously collected gravimetric data and all the available geological, geophysical and analogue modelling data.

The island of Ischia, located at the northwestern end of the Gulf of Naples in the southern Tyrrhenian Sea (Italy; Fig. 1a, b), is the emerged top of an active volcanic complex rising more than 1000 m above the seafloor (Orsi et al. 1999; Bruno et al. 2002), along the margin of an east–west-trending scarp that borders the Phlegraean Volcanic District to the south (Fig. 1b). This also includes the extinct volcano of the island of Procida and the still active and restless Campi Flegrei caldera (Orsi et al. 1996). The Phlegraean Volcanic District and Mt Vesuvius constitute the Neapolitan volcanic area, which is the active southern termination of the Plio-Quaternary volcanic chain that borders the eastern Tyrrhenian margin (Peccerillo 2005). Tyrrhenian volcanism is related to the extensional tectonic phases that accompanied the anticlockwise rotation of the Italian peninsula, during the complex interaction between the Africa and Eurasian plates, which generated the Apennine thrust-and-fold belt (D'Argenio et al. 1973; Ippolito et al. 1973; Finetti and Morelli 1974; Bartole 1984; Malinverno and Ryan 1986; Cinque et al. 1993; Piochi et al. 2004). Along the Tyrrhenian margin of the Apennine chain, extension was accommodated by the activation of NW–SE normal faults and NE–SW normal to strike-slip transfer fault systems, which dismembered the chain in horst and graben structures and allowed magmas to reach the surface, feeding volcanism (Mariani and Prato 1988; Faccenna et al. 1994; Fusi 1996; Acocella and Funiciello 2006). The Campanian Plain is one of these graben that hosts the Neapolitan volcanic area. It is a NW–SE elongated structural depression, filled by a thick sequence of marine and continental sedimentary deposits, and volcanic-volcanoclastic successions that compensated its subsidence, leading to a complete emersion at around 39 ka (Brocchini et al. 2001; De Vivo et al. 2001). This graben is bordered towards the NW, NE and SE by the Meso-Cenozoic carbonate and terrigenous successions of the Apennine chain, and is subdivided in minor NE–SW-oriented horst-and-graben structures (Carrara et al. 1973; Finetti and Morelli 1974; Fedi and Rapolla 1987; Brancaccio et al. 1991), which lower the carbonates down to 3000–4000 m (Bartole 1984; Moussat et al. 1986; Mariani and Prato 1988; Argnani et al. 1989; Bruno et al. 1998; Florio et al. 1999; Milia and Torrente 1999, 2015; Milia et al. 2003; Caiazzo et al. 2006). Neapolitan volcanoes lie on the structural highs: in particular, the Phlegraean Volcanic District is related to a NE–SW-trending horst-type structure that connects southeastwards, through a graben, to another structural high on which Mt Vesuvius lies (Orsi et al. 1996).

The buried structure of the Campanian Plain has been investigated by many geophysical studies, which allowed us to hypothesize the geometry of the carbonate basement. This has also been outlined through gravity data on land (Oliveri Del Castillo 1966; Carrara et al. 1973; Luongo et al. 1988; Ferri et al. 1990; Cubellis et al. 1995) and offshore (Berrino et al. 1998, 2008; Berrino and Camacho 2008).

The island of Ischia volcanic field covers an area of about 46 km² and is composed of volcanic rocks, landslide deposits and subordinate terrigenous sediments (Fig. 1a), reflecting a complex history of alternating constructive and destructive phases due to interplay among tectonism, volcanism, volcano-tectonism, erosion and sedimentation (de Vita et al. 2006, 2010; Selva et al. 2019). Volcanism dates back to more than 150 ka (Vezzoli 1988; Civetta et al. 1991) and is dominated by the caldera-forming eruptions that generated the so-called Mt Epomeo Green Tuff (MEGT), between 60 and 50 ka (Poli et al. 1987; Vezzoli 1988; Tibaldi and Vezzoli 1998; Brown et al. 2008, 2014; Sbrana and Toccaceli 2011; Sbrana et al. 2018). After the emplacement of the MEGT and the related volcano-tectonic collapse, the caldera depression was the site of submarine sedimentation, with the deposition of landslide deposits and terrigenous sediments (Vezzoli 1988; Barra et al. 1992). Meanwhile, in a not well-defined timespan between 50 and 33 ka, the caldera floor was affected by a resurgence phenomenon (Orsi et al. 1991; Acocella and Funiciello 1999; Molin et al. 2003; Carlino et al. 2006; Marotta and de Vita 2014), which occurred discontinuously through the asymmetrical uplifting of a discrete number of differentially displaced and tilted blocks, likely due to repeated episodes of magma intrusion into the Ischia feeding system (Selva et al. 2019 and references therein). The resurgent area has a polygonal shape, whose geometry is defined by the intersection of regional faults, reactivated during resurgence, and newly formed volcano-tectonic faults (Orsi et al. 1991; Acocella and Funiciello 1999). The most uplifted part of the resurgent area is shown in Figure 1a. The western sector of the resurgent area is bordered by vertical to high-angle (slightly inward-dipping) north–south-, NE–SW- and NW–SE-trending faults, cut by late outward-dipping normal faults due to gravitational readjustment of the slopes (de Vita et al. 2006, 2010; Della Seta et al. 2012; Selva et al. 2019). To the east of the resurgent area, vertical or outward-dipping, north–south, NE–SW and NW–SE normal faults, formed a lowland which is connected westward to the resurgent area of Mt Epomeo through a series of differentially displaced blocks (de Vita et al. 2006, 2010; Della Seta et al. 2012). Towards the NNE the limit of the resurgent area is not well defined, as along the coast beach and shallow-sea fossiliferous deposits are exposed, displaced at different heights above sea-level by east–west, slightly inward-dipping, sub-vertical faults and subordinate NW–SE-trending vertical faults (Alessio et al. 1996; de Vita et al. 2006, 2010; De Novellis et al. 2018; Nappi et al. 2018, 2021; Trasatti et al. 2019; Berrino et al. 2021; Carlino et al. 2022). Resurgence caused the displacement of variably uplifted and inclined blocks, with a general tilt around a NE–SW-oriented horizontal axis, located in the southeastern part of the resurgent area (Selva et al. 2019). The result is an asymmetrical block structure, with a maximum uplifted block in the northwestern part of the resurgent area (Fig. 1a; Rittmann and Gottini 1980; Vezzoli 1988; Orsi et al. 1991; Acocella and Funiciello 1999; Molin et al. 2003; de Vita et al. 2006, 2010; Della Seta et al. 2012).

Volcanism later continued with effusive and explosive, both magmatic and phreatomagmatic, eruptions, concentrated in periods of activity that were separated by centuries to millennia of quiescence, and lasted until the Arso eruption in 1302 CE (de Vita et al. 2010 and references therein).

In the last period of activity, started at around 10 ka, volcanism was very sustained and fed by vents mainly concentrated in the eastern part of the island or along regional fault systems and at the margin of the resurgent block (de Vita et al. 2010; Selva et al. 2019), likely being influenced by the intermittent renewal of resurgence (Orsi et al. 1991; de Vita et al. 2006, 2010; Marotta and de Vita 2014; Galetto et al. 2017).

The structural setting of the Neapolitan volcanic area has been investigated by many geophysical (gravimetric, seismic, aeromagnetic) studies, carried out at different scales since the early 1960s (e.g. Maino et al. 1963; Segre 1967; Cassano and La Torre 1987a, b; Nunziata and Rapolla 1987; Rapolla et al. 1989; Berrino et al. 1998, 2008; Orsi et al. 1999; Tondi and De Franco 2003, 2006; Judenherc and Zollo 2004; Cella et al. 2008; Zollo et al. 2008; Nunziata 2010; Linde et al. 2017).

The results of almost all these studies converge toward an interpretation in which: (a) a sedimentary/crystalline basement, variably displaced by tectonic/volcano-tectonic structures and regional subsidence is present beneath the Campanian Plain; (b) an almost continuous, partially molten body extends underneath the entire Neapolitan volcanic area, whose top is located at a depth of about 8 km; (c) this zone extends over a surface of not less than 30 km², with an approximate thickness of 1 km; (d) the Phlegraean Volcanic District lies on a structural high, which constitutes a second order horst-like ridge included in the Campanian Plain graben.

Geochemical studies corroborated this interpretation, suggesting that a single deep magma reservoir could feed the entire Neapolitan volcanic area. In this reservoir, more or less primitive magmas derived from the mantle, stagnate at a depth of 8–10 km and differentiate, partly assimilating the local continental crust, and then migrate toward shallower crustal levels. Evolved magmas gave rise to the formation of multiple reservoirs, beneath each single volcano of the Neapolitan area, further differentiating and mixing independently, due to the local structural setting and volcanic history (Paoletti et al. 2013 and references therein).

The island of Ischia, as well, has been the object of detailed geophysical studies that investigated its crustal structure, mainly modelling gravimetric and magnetic data or taking into account seismic and geodetic data (Maino and Tribalto 1971; Nunziata and Rapolla 1987; Berrino et al. 1998, 2008; Bruno et al. 2002; Paoletti et al. 2009, 2013, 2017; Fedi et al. 2011; Capuano et al. 2015; Di Giuseppe et al. 2017; Trasatti et al. 2019; Berrino et al. 2021; Galvani et al. 2021).

The results of these studies highlight the possible geometry of the margin of the MEGT-eruption caldera and of the resurgent block, distinguishing the surface structures from those at medium (local) and high (regional) depth. In particular, Maino and Tribalto (1971) evidence a gravimetric maximum to the SW of Mt Epomeo, interpreted as the evidence of magmatic masses solidified at depth, and high gravimetric gradients, corresponding to faults and fractures. Nunziata and Rapolla (1987), on the basis of gravimetric and magnetic surveys, identify differences between short- and long-period anomalies, which would identify the surface distribution of lavas and pyroclastics and crater structures, and a lava substrate with a depth varying from 600 m to 1.4 km, respectively. This lava substrate is partially demagnetized to the SW of Mt Epomeo, due to geothermal fluid circulation, which induces an anomalous thermal state in that sector. Berrino et al. (1998, 2008) defined Ischia as a structure separated from Campi Flegrei and the island of Procida, and interpreted the positive gravity anomaly beneath Ischia as an evidence of the carbonate basement uplift. Particularly, Berrino et al. (2008) also investigated the shallow setting of Ischia, suggesting the presence, above the basement, of a horst-like structure beneath Mt Epomeo. Their model is in agreement with the model of resurgence proposed by Rittmann (1930, 1948) and Orsi et al. (1991).

Orsi et al. (1999) performed 3D modelling of the regional magnetic field of all the Phlegraean Volcanic District, hypothesizing the presence of a partially solidified magma chamber at a depth ranging between 7 and 10 km b.s.l, and a series of shallower bodies, at a variable degree of magnetization beneath the island of Ischia. The magnetometric, side scan sonar and seismic survey by Bruno et al. (2002) showed a greater extension of the Ischia volcanic field in the western offshore of the island, highlighting magmatic intrusions, associated with linear structures, aligned according to the distribution of the main eruptive vents, active between Ischia and the Campi Flegrei caldera. Paoletti et al. (2009, 2013), based on the inversion of aeromagnetic and self-potential data, and the comparison with gravimetric and petrologic data, hypothesized the existence of an east–west elongated trachytic body, with a top located at a depth of about 1200–1750 m. According to the simultaneous presence of a gravimetric maximum and a magnetic minimum in the central part of the island, this body is interpreted as the product of repeated magmatic intrusions at a different degree of cooling, and partially demagnetized in the western underground of Mt Epomeo, probably due to the circulation of geothermal fluids in this sector, or the presence of partially molten-rock bodies in the substrate. Capuano et al. (2015) investigated the shallow crustal structure of Ischia (<2 km b.s.l.) through seismic tomography and gravity data inversion, giving a 3D image of the caldera rim along the perimeter of the island, and evidencing the presence of a high velocity and density area inside the caldera, interpreted as the resurgent block. Paoletti et al. (2017) performed a semblance analysis of the data of a new aeromagnetic survey, compared with the results of a gravimetric survey carried out for the new geological map of Ischia (Fedi et al. 2011). They distinguish a surface structure, dominated by exposed bodies at different density and magnetic intensity (lavas and pyroclastics), and a deep structure, characterized by the presence of a high-density and partially demagnetized buried body, interpreted as the igneous basement, whose demagnetized part (at high temperature) is suggested to be responsible for the resurgence. The magnetotelluric survey carried out by Di Giuseppe et al. (2017) highlighted a large resistivity anomaly below the resurgent structure of Mt Epomeo, interpreted as the apical part of a crystalline magmatic body, intruded below the central part of the island at a depth of about 1 km b.s.l. This body is, in turn, associated with the so-called laccolith of Ischia (Rittmann 1930; Sbrana et al. 2009; Carlino 2012), which produced the bending, fracturing and faulting of the overlying crust. As a result, the uplifted block would have been broken up into minor blocks, with the underlying laccolith possibly developing as a complex structure (Di Giuseppe et al. 2017). The thermo-rheological modelling carried out by Castaldo et al. (2017) allowed the reconstruction of a similar setting of the upper part of the crust beneath Ischia, highlighting individual, shallow magmatic intrusions on the basis of the recorded thermal anomalies. The modelling of the source of the deformation observed at Ischia during the last earthquake of 21 August 2017 was performed by Trasatti et al. (2019) and demonstrated that this kind of deformation, as well as the generalized subsidence of the Mt Epomeo block, is compatible with long-term degassing of a magmatic source located at a depth of about 2 km b.s.l. A similar magmatic source has been hypothesized by Galvani et al. (2021), although they interpreted the related seismicity as mainly reflecting the dynamics of the shallow hydrothermal system, rather than the deflation of the magmatic source.

With the aim of better defining our knowledge on the structural setting of the Ischia resurgent block, and on the nature of the source of this deformation, a multidisciplinary approach was followed, based on comparison between the results of a new 3D inversion of gravity data and all available geological, geophysical and analogical modelling data.

Starting from the Bouguer gravity map of the Neapolitan volcanic area and previous 2.5D modelling (Berrino et al. 1998, 2008), we selected 437 gravity data, both on land and offshore, on and around the island of Ischia (bounded in a blue square on Fig. 2a, b) and carried out a 3D interpretation using the original algorithm proposed by Camacho et al. (2000, 2002) and Berrino and Camacho (2008), respectively named ‘Growth’ and ‘Layers’.

The base of the map is a larger dataset of 2876 gravity values, covering the whole Neapolitan volcanic area and the Gulf of Naples (Fig. 2a). It results from the union of 862 offshore gravity data collected during five cruises carried out from 1988 to 1994 (Berrino et al. 1991, 1998), with about 2000 on-land extant gravity stations made uniform and globally re-analysed (Berrino et al. 1998, 2008 and references therein). Details about the offshore field surveys, instruments, positioning of stations, errors on gravity and depth data, on-land data collection and integration and global data processing are given in Berrino et al. (1998, 2008). Figure 2a shows the Bouguer anomaly map spanning the whole Neapolitan area and drawn with a 2.5 mGal contour interval, taking into account the upper limit of 0.7 mGal of the error affecting the offshore data. It was obtained with reference to the 1980 Ellipsoid (Moritz 1984), using the density values of 2200 kg m⁻³ to reduce the Bouguer and terrain effects, which is the most suitable value for the Neapolitan volcanic areas (Berrino et al. 1998).

In this paper we adapted and merged the basic ideas of two already tested methods for 3D inversion, which detect isolated bodies and layered discontinuities respectively (Camacho et al. 2002; Berrino and Camacho 2008). In the first method (named ‘Growth’), the model is described as being composed of isolated anomalous bodies, which are constructed in a very free growth process as 3D aggregations of cells. In a versatile and non-subjective form, and with few constraints, the process can produce 3D models of the anomalous structures (position, depth, size, shape), which are more valuable if suitable values for the density contrast are previously defined. Conversely, some application problems can arise when the causative structure cannot be clearly associated with isolated bodies. In this case, the inversion model will provide a simplified, rather indicative, solution to the inversion problem that needs further analysis to reach any realistic conclusions. This is the case, for instance, of anomalies due to small distortions of sub-horizontal layers in the underground, like the case of the Neapolitan area, solved by the second algorithm (named ‘Layers’) which describes the underground model introducing sub-horizontal layers, where the irregular discontinuity surfaces are constructed by displacing, step by step (according to a system of connected cells, in a growth process), the original flat mean surface (Berrino and Camacho 2008).

This method was tested at Vesuvius (Berrino and Camacho 2008). Namely, it starts from several horizontal layers, which can be introduced ad hoc. Also, several corresponding density contrasts are previously selected. Then, the algorithm works according to a non-linear explorative approach, to ‘deform’, or better, to ‘shape’, the layers step by step, to finally obtain some irregular shapes that can fit the observed anomaly satisfactorily. Finally, this method allows us to describe the underground model in terms of shallow isolated bodies and sub-horizontal layers in a unique solution. In addition, the sub-horizontal layers can be replaced by a linear or exponential increasing density gradient v. depth.

In short, in order to obtain the 3D model, the subsurface volume close to the surveyed area is dismantled into a 3D partition of small parallelepiped cells. A 3D grid, defining depth of exploration, depth and size of the prisms, is initially automatically defined by the code. But size of the prisms and the exploration depth can be changed by the user. The method considers n gravity stations, located on a rugged topography and which are not necessarily gridded, and their anomalous gravity values (Bouguer anomaly), and initially nh horizontal surfaces with depths and density discontinuities, considering both positive and negative density contrasts.

The 3D model will be described as an aggregation of prismatic cells filled with provided density contrast close to the discontinuity surfaces, thus giving rise to ‘shaped’ layers. When the filled cell is just below the discontinuity surface, this means that there is an intrusion of low density from the upper medium into the lower one. Conversely, when the filled cell is close and above the discontinuity surface, this means there is an intrusion of high density from the lower medium into the upper one.

The gravity attraction Aij at the i-th station due to the j-th prism, per unit density, is computed, as well as the subsequent calculated anomaly values for the resulting model. A regional trend is also simultaneously obtained.

In order to reduce the very large number of degrees of freedom for the model, an exploratory method (in a more restrictive context) is applied not to the global model, but just to every step of its growth process. Under these conditions the exploratory process becomes very effective. This means that the exploration of the possibilities for the entire model is substituted by the exploration of several possibilities of growth (cell by cell) for each step of the surface deformation.

Thus, the prismatic cells that are connected (up or down) to discontinuity surfaces are systematically tested, step by step. Moreover, to solve the non-uniqueness problem, an additional condition of minimization of the model is adopted; i.e. the solution is obtained through a mixed condition between the gravity fitness and the whole anomalous mass quantity.

For one step, some cells are previously filled to modify the geometry of the initial discontinuity surfaces, although not enough to reproduce the anomaly observed. Every cell (up and down) that is connected to every actual discontinuity surface is tested.

This process is repeated successively and stops when a good fit for a final geometry of the discontinuity surfaces and a final regional trend is reached.

A detailed description of the algorithms is given in Camacho et al. (2000, 2002), and Berrino and Camacho (2008) and its last upgrade in Camacho et al. (2021).

The results of the 3D inversion of gravimetric data were subsequently interpreted on the basis of the comparison with surface geological and structural evidence (Vezzoli 1988; Acocella and Funiciello 1999; Marotta 2001; Molin et al. 2003; Chiodini et al. 2004; Tibaldi and Vezzoli 2004; Carlino 2012; Carlino et al. 2006, 2010, 2021; de Vita et al. 2006, 2010; Sbrana and Toccaceli 2011; Della Seta et al. 2012; Sbrana et al. 2018; Nappi et al. 2018, 2021), geophysical and multidisciplinary surveys (Maino et al. 1963; Segre 1967; Maino and Tribalto 1971; Fedi and Rapolla 1987; Nunziata and Rapolla 1987; Rapolla et al. 1989; Orsi et al. 1999; Bruno et al. 2002; Paoletti et al. 2009, 2013, 2017; Di Napoli et al. 2011; Fedi et al. 2011; Capuano et al. 2015; Di Giuseppe et al. 2017; De Novellis et al. 2018; Selva et al. 2019; Trasatti et al. 2019) and the results of analogue-modelling of resurgent calderas (Acocella et al. 2000, 2001, 2004; Acocella 2007; Marotta and de Vita 2014).

To investigate the structure of Mt Epomeo resurgent block and, possibly, the crustal sector hosting the shallow magmatic feeding system of Ischia, we have selected an area of about 15 × 15 km centred on the island (blue square in Fig. 2a, b) in which the selected 437 gravity stations were positioned. Although we aimed to design objective models without any subjective preliminary information or input model, the information provided by the previous 2.5D interpretation and those from old exploration wells (Fig. 2b) (Penta and Conforto 1951; Berrino et al. 1998, 2008; Sbrana and Toccaceli 2011) have been taken into account as a reference, when possible.

We used no initial constraints and therefore the obtained models are totally non-subjective. The 3D model, extended down to a depth of 4 km, has been obtained in a nearly automatic process and was mainly addressed to define the structure of the Mt Epomeo resurgent block and the island crustal structure in the first 3 km from the present sea-level.

The main features of the Bouguer gravity map of the Neapolitan area (Fig. 2a) can be summarized as follows: a strip of maximum gradient runs almost parallel to the Sorrento peninsula and turns towards the NW around the northern limit of the relative minimum of the Sarno Plain. In the southern part of the Gulf of Pozzuoli, a broad gravity minimum separates the relative maximum of the Phlegraean Volcanic District and Somma-Vesuvius. The Campi Flegrei caldera and Somma-Vesuvius are settled at the southern edge of the large gravity minimum of Acerra, to the north of Naples. A well-defined gravity minimum that spans Campi Flegrei caldera is evident. Strong gradients border the island of Ischia, which appears clearly separated by the remnant Phlegraean Volcanic District (Campi Flegrei caldera and Procida island).

The global Bouguer anomaly gravity map has been already interpreted by 2.5D modelling along a series of profiles that provided information about the main volcanic structures, and particularly about the shape and the depth of the sedimentary/crystalline basement. This kind of inversion also allowed the reconstruction of a pseudo-3D pattern of the Meso-Cenozoic sedimentary basement and the delineation of the main tectonic structures underneath the Neapolitan volcanoes (Berrino et al. 1998, 2008). The map of Figure 2a represents the base for the 3D inversion performed in this paper on the area of Ischia (blue rectangle of Fig. 2a).

In the zoom of Figure 2b the Bouguer anomaly has been drawn with a 2.0 mGal contouring interval to obtain a more detailed vision of the island. The Bouguer gravity field on and around Ischia is characterized by an average positive anomaly of about 60 mGal on which a set of local and small anomalies are overlapped. Strong gravity gradients border the island, mainly affecting the southern part. Here, a residual terrain effect, not completely removed because of very rough morphology (Fig. 1), could affect the gravity anomaly field (Berrino et al. 2008).

In our 3D modelling, we first analysed from the Bouguer anomaly the regional trend that is computed simultaneously in the inversion process; it gives an idea of the shape of the deep basement below the island. Figure 3a and b shows the difference between the regional (Fig. 3a) and local (Fig. 3b) anomaly; here the values are indicated in microGal, which is the unit adopted in the inversion process. The regional trend (Fig. 3a) is characterized by a gravity gradient of 2.430 mGal km⁻¹ increasing towards the SW (N14°E), as also indicated by the arrow. Consequently, the local anomaly (Fig. 3b) is generally positive on most of the island and the western offshore, while the remaining part of the offshore (southern, eastern and northern) shows a negative anomaly. The local anomaly (Fig. 3b) ranges from +7.9 to −15.2 mGal and suggests the presence of a large local positive density body, probably a volcanic one, mostly restricted to the emerging island and offshore westwards it. In addition, in order to set Ischia in the more general context of the Neapolitan area, we also extracted a regional west–east profile crossing Ischia (Fig. 4) from the 3D gravity inversion carried out by Berrino and Camacho (2008) on a wide area. The east–west profile crosses the island of Ischia, the Bay of Naples and the Sarno plain (between Mt Vesuvius and the Sorrento Peninsula). Two deep horizontal sections at different depth are displayed (3000 and 6000 m b.s.l.). In these pictures the contribution to the gravity anomaly of bodies with increasing density with depth (from 2400 to 2650 kg m⁻³), shallower in correspondence with the island of Ischia, off the coast of Vesuvius and below the Sorrento Peninsula, is clearly evident.

As we are interested in defining both the distribution at depth of local density bodies and the shape and depth of the basement underneath Ischia, we analysed the observed Bouguer gravity anomalies without removing the regional trend.

The best 3D model (Fig. 4) shows several differences compared with that resulting from the previous 2.5D inversion (Berrino et al. 1998, 2008). Firstly, any introduction of shallow sub-horizontal layers does not give a good fit, from a mathematical point of view and realistic geological models. The best model is given when we introduce a linear density gradient v. depth (of about 100 kg m⁻³ km⁻¹). Another difference with the 2.5D model is the automatically resolved deep density distribution. The density values are lower than that introduced in the 2.5D model. In the present work we have density values ranging from 1700 to 2500 kg m⁻³ against the 2000–2600 kg m⁻³ density range introduced in the previous inversion. Only in the area in which the boreholes are located, the automatically computed density values in our inversion result in the same observed in the cores and, of course, is the same used as constraints in the 2.5D inversion (Berrino et al. 1998, 2008). This is a further confirmation of the goodness of the 3D model.

The distribution of masses at different densities with increasing depth, resulting from the 3D inversion at Ischia is shown in Figure 5 by means of selected horizontal deep slices and some cross-sections through the Mt Epomeo block.

The shallower horizontal slices show several negative density bodies (visible on the −500 and −1000 m slices), distributed around the island, close to the coast, and on land. Some of the negative density bodies on land, reaching the surface from about 1 km depth, are those located in the southern part of the island and in its centre, beneath Mt Epomeo. High-density bodies are scattered throughout the island, in the shallowest horizontal slice, and tend to concentrate in the eastern and western/northwestern sectors of the island at progressively greater depths. The singled out bodies are better defined in the north–south, west–east, NW–SE and SW–NE sections shown in Figure 5, where the shape and the dimensions of the high- and low-density bodies are clearly outlined. The deepest horizontal slice, at 3000 m b.s.l, clearly shows the presence under the whole island, starting from a depth of 2.5 km, of a continuous high-density body, whose density is higher than that of those other bodies, which seem to rise up to the surface in the eastern and western/northwestern parts of the island. The link of these progressively shallower bodies with a common high-density deep root is evident in the west–east and NW–SE profiles (Fig. 5); as well as the central zone in between, filled by less dense rocks. The deepest part of the island is characterized by a very uniform density distribution, visible in the −3000 m slice. This deep high-density body is well shown in the cross-sections of Figure 5, in which it is topped by a flat surface, clearly separated from the overlying layer by a sharp density increase at about 2.5 km b.s.l.

According to previous geophysical and geological studies on the structure of the Neapolitan volcanic area, the regional trend of the Bouguer anomalies in the map of Figure 2a can be interpreted as the product of the differential displacement of the reliefs composed of the Meso-Cenozoic carbonate and terrigenous sedimentary successions, which border the main graben of the Campanian Plain and that are found at variable depths below the plain, being upraised in correspondence of secondary horst-type structures, as for example below Vesuvius and the Phlegraean Volcanic District. In Figure 2b the general southwestward increase in the Bouguer anomaly, from about 40 mGal to around 100 mGal evidences the contribution of the strong regional positive anomaly of the southern Tyrrhenian sea, due to the Moho uplift in this region (Rapolla et al. 1989.). The 2.5D interpretation restricted to the Phlegraean Volcanic District, carried out by Berrino et al. (2008), suggested that Ischia seems to be an isolated block, delimited and displaced by a series of regional-tectonic structures, according to Maino et al. (1963), which postulated a NW–SE-trending regional structure separating the island from the main Phlegraean sector. At a more local scale (Fig. 3a, b) the gravity anomaly obtained in this paper from the Bouguer anomaly map as ‘observed minus regional’ clearly isolated the Ischia gravity features from the remnant Phlegraean Volcanic District, confirming that Ischia might constitute a separated block. This also means that the positive gravity anomaly bordering the island with a strong gradient (Figs 2a, b, 3a & b) is mostly due to the regional field. It could confirm the hypothesis of Berrino et al. (2008) of the uprising basement below the island, likely represented by the Meso-Cenozoic sedimentary succession, displaced during the formation of the Campanian Plain structure (Berrino et al. 2008). The nature of this deep basement is also suggested by the results of the 3D gravity inversion carried out on 6203 gravity stations covering the whole Neapolitan area (ref. Fig. 1a), and extended until about 20 km in depth (Berrino and Camacho 2008). In Figure 4 the regional west–east profile crossing Ischia is shown along with the −2000 m and −6000 m deep horizontal sections. The west–east vertical profile shows that the regional deep basement (green to orange layers) arises under Ischia exactly as under the Sorrento peninsula, where the sedimentary succession is exposed. Moreover, as the previous 2.5D interpretation did not furnish important information about very shallow and isolated bodies, the 3D interpretation performed here has been aimed at better defining the shallow density distribution and trying to discriminate the nature of the uprising basement. We looked for the determination, in a non-subjective way, of a model of the underground density distribution that can reproduce the observed gravity anomaly.

Taking into account the information coming from other geophysical investigations (cf. the section titled ‘Previous geophysical studies’), and mainly the structural setting already furnished by the 2.5D inversions (Berrino et al. 2008), we describe the 3D anomalous density structures mainly by means of shallow isolated bodies and several sub-horizontal discontinuity layers (density discontinuities) with irregular boundaries.

This is a very interesting result since it suggests that the upper part of the medium below the island is not necessary formed by superimposed layer but most probably it is constituted by a continuous medium and by isolated density bodies overlapping a deeper basement as suggested by the regional field trend.

This structural setting is well highlighted by the analysis of the horizontal slices at different depths, shown in Figure 5, and the related density vertical sections that cross Mt Epomeo.

The distribution of bodies with a negative density contrast that mimics a ring around the island, in the −500 m slice, could, at first glance, be interpreted as the expression of the caldera edge of the MEGT, as also dubiously proposed by Paoletti et al. (2013) based on reduced-to-pole aeromagnetic data. However, this interpretation, already discarded by Paoletti et al. (2017), can be definitely ruled out on the basis of the surface geological constraints (Orsi et al. 1991; Tibaldi and Vezzoli 1998; Di Napoli et al. 2011). This feature can be interpreted more effectively if we break it down by analysing the areas where it is observed separately. Along the northern coast, the control exerted by the regional trend of gravimetric anomalies and the contrast with the higher density bodies exposed or buried at shallow depth on land is evident (see Fig. 3a, b); moreover, the areal distribution of the minima seems to trace the surface distribution of low-density bodies (debris avalanche and lahar deposits; Della Seta et al. 2012), also exposed along the western coast. This effect gradually decreases with depth, as seen in the −1000 m slice. The negative density bodies distributed close to the coast in the southern offshore may be an effect of a residual terrain effect, not completely removed because of a very rough morphology (Fig. 1b; Berrino et al. 2008). According to all the surface geological and structural evidence, the possible MEGT caldera rim can be more likely represented by the dashed red line reported in the −500 m horizontal slice of Figure 5. Here the negative density bodies located in the southern part of the island may be the surface effect of thick sequences of loose pyroclastic deposits, related to the post-caldera activity of vents located close to the caldera rim (Vezzoli 1988; de Vita et al. 2010; Sbrana and Toccaceli 2011). Along the western, northern and eastern sectors, the structural margins of the caldera were the site of magmatic intrusions, which in the northern and eastern sectors of the island fed effusive eruptions that emplaced thick lava flows and domes (Vezzoli 1988; de Vita et al. 2010; Sbrana and Toccaceli 2011). These intrusions are marked by the relative positive density bodies that correspond to the exposed lava bodies in the shallowest horizontal slices and tend to concentrate in the eastern and western/northwestern sectors of the island at progressively greater depths (Fig. 5). These progressively deeper bodies are linked to a common root at a depth included between 1 and 2.5 km, characterized by a density of 2300/2350 kg m⁻³, which, considering the lateral variabilities observed, could correspond to the average density of a shallow trachytic reservoir composed of partially molten and high-temperature demagnetized parts, inter-fingered with solid and low-temperature magnetized parts. This interpretation is also supported by aeromagnetic and petrologic studies (Orsi et al. 1999; Chiodini et al. 2004; Di Napoli et al. 2011; Paoletti et al. 2013; Casalini et al. 2017) in which the high geothermal gradient and the hydrothermal fluid circulation of the island are taken into account as well. In the NW–SE and SW–NE cross-sections of Figure 5, similar high density bodies are here interpreted as dykes that intruded tectonic and volcano-tectonic structures in the Ischia eastern and northwestern sectors, feeding the past 10 ka volcanic activity (Orsi et al. 1991; de Vita et al. 2010). These bodies are separated by a central lower-density area, in which the very central part is occupied by a negative density minimum that extends downward to a depth of about 1 km (Fig. 5). In Figure 6 a detail of the north–south and west–east cross-sections is displayed along with the –1500 m horizontal slice, to focus on the central density minimum, which is characterized by an average density of about 1700 kg m⁻³ and thickens westward. It is here interpreted as the MEGT resurgent block, lying on the pre-caldera volcanics, whose average density varies between 2200 and 2250 kg m⁻³, being composed of alternating trachytic lavas and pyroclastics. This pre-caldera volcanic succession is intruded and deformed by more or less evolved trachytic magma bodies (average density between 2300 and 2400 kg m⁻³), which form a complex magmatic reservoir, presently in a variable state of temperature and magnetization (Orsi et al. 1999; Chiodini et al. 2004; Di Napoli et al. 2011; Paoletti et al. 2013; Casalini et al. 2017) at a depth included between 1500 and 2500 m b.s.l. At this depth the density contrast between the magma bodies and the surrounding rocks is almost negligible and approaches a value compatible with that of the Meso-Cenozoic sedimentary basement (c. 2400–2650 kg m⁻³, sandstone, limestone and shale; Fedi and Rapolla 1993). From a depth of about 3000 m the density distribution is almost homogeneous and the contribution to the gravity anomaly is likely due almost exclusively to the sedimentary basement. Density then increases progressively with depth, reaching the values of the crystalline basement at about 6000 m.

The results of 3D gravity inversion carried out in this paper by processing previously acquired Bouguer anomaly data (Tribalto and Maino 1962; Maino et al. 1963; Imbò et al. 1964; Oliveri del Castillo et al. 1964; Maino and Tribalto 1971; Cassano and La Torre 1987a, b; Berrino et al. 1991, 1998, 2008; Italian Geological Survey (SGI), unpublished data) have been interpreted in the light of all the available surface and subsurface geological and structural evidence, geophysical and multidisciplinary surveys and the results of analogue-modelling of resurgent calderas (see the ‘Materials and methods’ section for the reference list). Our interpretation allowed us to outline a possible schematic reconstruction of the geological and structural setting of the crustal sector beneath Ischia, here proposed in Figure 7.

From a depth of about 22 km the Mohorovicic discontinuity separates the crystalline basement from the upper mantle and deepens eastward, as reconstructed in Rapolla et al. (1989) from a very large wavelength gravimetric field. The possible pattern of the Curie isotherm, inferred from magnetic data of the same wavelength, separates the portion of the crystalline basement at a temperature above the Curie point from the underlying magnetic basement (Paoletti et al. 2013). The overlying low density/velocity body postulated by many authors (cf. the ‘Ischia resurgent caldera: state of knowledge’ section for references) is reported in Figure 7 at a depth of about 8 km, according to D'Antonio et al. (2007 and references therein) and Zollo et al. (2008). It is interpreted as a partially molten deep trachybasaltic magma reservoir of regional extent, in which mantle-derived mafic magmas stagnate and differentiate, likely assimilating continental crust (Beccaluva et al. 1991; Rolandi et al. 2003; Civetta et al. 2004; Tonarini et al. 2004; D'Antonio et al. 2007; Di Renzo et al. 2011; Paoletti et al. 2013). From this stagnation level magmas can: (a) rise and reach the surface through regional crustal structures as poorly differentiated melts, directly feeding volcanic eruptions (D'Antonio et al. 2013); (b) rise, reaching shallower multiple reservoirs, where they differentiate before being erupted when conditions are right, or stagnate and crystallize at different depths, possibly remaining in a stationary thermal state for a long time (Casalini et al. 2017).

According to Berrino et al. (1998, 2008) a horst-type structure beneath the Phlegraean Volcanic District raised up the crystalline basement and its Meso-Cenozoic sedimentary cover, so that the top of this succession could have reached a depth of about 1500–2000 m b.s.l., with the top of the crystalline basement at a depth variable between 6 and 8 km (or more; Rapolla et al. 1989) following second-order horst and graben structures.

Repeated magmatic intrusions, occurred after the caldera-forming eruptions of the Mt Epomeo Green Tuff, generated the present asymmetrically-shaped multiple magma reservoir, composed of bodies at different degrees of cooling and magnetization, likely due to the circulation of geothermal fluids, or the presence of partially molten-rocks (Orsi et al. 1999; Paoletti et al. 2009, 2013). This periodically refilled magma reservoir, is considered to be the engine of the resurgence mechanism (Orsi et al. 1991; Tibaldi and Vezzoli 2004; Carlino et al. 2006; de Vita et al. 2006, 2010; Carlino 2012; Marotta and de Vita 2014; Paoletti et al. 2017), which produced the uplifting, fracturing and faulting of the overlying pre-caldera volcanic basement and caldera-filling deposits. Note that, in our reconstruction, the most uplifted sector of the resurgent block, due to the asymmetrically-shaped intrusion, is in correspondence of the thickest overburden above the magma chamber, according to the analogue modelling experiments carried out by Marotta and de Vita (2014). Therefore, rather than a uniformly distributed laccolithic intrusion underneath the resurgent block, geophysical and petrological evidence, along with structural geology observations and the results of the present 3D gravimetric inversion, more likely suggest a composite geometry of the magmatic system that fed the recent volcanic activity at Ischia and conditioned the dynamics and kinematics of resurgence, subsidence and seismicity (Selva et al. 2019; Trasatti et al. 2019).

The authors are very grateful to Rahmat Catur Wibowo and an anonymous reviewer for suggested revisions, and the editorial work by Luca D'Auria, whose suggestions greatly improved the quality of this paper.

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

EM: conceptualization (equal), project administration (lead), supervision (equal), writing – original draft (supporting), writing – review & editing (equal); GB: conceptualization (equal), data curation (lead), formal analysis (equal), software (equal), writing – original draft (equal), writing – review & editing (equal); SDV: conceptualization (lead), formal analysis (equal), visualization (equal), writing – original draft (equal), writing – review & editing (equal); MDV: formal analysis (equal), visualization (supporting), writing – original draft (supporting), writing – review & editing (supporting); AGC: data curation (equal), formal analysis (equal), software (equal), methodology (lead), writing – original draft (supporting), writing – review & editing (supporting).

This work was supported by the INGV project Pianeta Dinamico – Working Earth (CUP 1466 D53J19000170001 – ‘Fondo finalizzato al rilancio degli investimenti delle 1467 amministrazioni centrali dello Stato e allo sviluppo del Paese’, legge 145/2018) – Theme 2 TIFEHO.

Data sharing is not applicable to this article as no datasets were generated during the current study.