Research ArticleGEOPHYSICS

Seismic signature of active intrusions in mountain chains

See allHide authors and affiliations

Science Advances  03 Jan 2018:
Vol. 4, no. 1, e1701825
DOI: 10.1126/sciadv.1701825


Intrusions are a ubiquitous component of mountain chains and testify to the emplacement of magma at depth. Understanding the emplacement and growth mechanisms of intrusions, such as diapiric or dike-like ascent, is critical to constrain the evolution and structure of the crust. Petrological and geological data allow us to reconstruct magma pathways and long-term magma differentiation and assembly processes. However, our ability to detect and reconstruct the short-term dynamics related to active intrusive episodes in mountain chains is embryonic, lacking recognized geophysical signals. We analyze an anomalously deep seismic sequence (maximum magnitude 5) characterized by low-frequency bursts of earthquakes that occurred in 2013 in the Apennine chain in Italy. We provide seismic evidences of fluid involvement in the earthquake nucleation process and identify a thermal anomaly in aquifers where CO2 of magmatic origin dissolves. We show that the intrusion of dike-like bodies in mountain chains may trigger earthquakes with magnitudes that may be relevant to seismic hazard assessment. These findings provide a new perspective on the emplacement mechanisms of intrusive bodies and the interpretation of the seismicity in mountain chains.


The intrusion of magma and associated fluids into the continental crust plays a major role in controlling the growth, evolution, and composition of the lithosphere (13). The dynamics of the emplacement of intrusive bodies has been investigated by field studies and numerical simulations, and two main mechanisms are proposed (4, 5): (i) magma ascent by fluid-filled cracks (dikes) and (ii) diapiric ascent. It is suggested that, independently from the proposed mechanism, intrusions grow by successive feeding episodes (pulses) (6). However, although processes of magma transfer to the Earth’s surface during volcanic eruptions are relatively well known, those related to the emplacement of intrusive bodies in mountain chains remain enigmatic, because records of geophysical and/or geochemical signals are lacking. In addition, the seismicity of mountain chains is commonly interpreted to be only due to tectonic stress and/or fluid pressure changes (7, 8).

The southern Apennine (SA) fold-and-thrust belt in Italy is associated with the southwestern subduction of the Adriatic plate and separates a western area of diffuse deep CO2 release (2 × 1011 mol year−1) from an eastern, nondegassing foreland region (9). The Apennine chain is presently affected by crustal delamination, with melts from the mantle wedge feeding the Vesuvius, Campi Flegrei, and Ischia active volcanoes, located to the west of SA (Fig. 1A) (10). Mafic intrusions have been drilled in some SA wells located over 100 km away from the volcanoes at depths of 2 to 5 km within the carbonates overlying the crystalline basement, but their age and volume are unknown (11). Earthquakes with magnitude up to 7 affect the first 15 km of the SA crust. These earthquakes and the associated seismogenic normal faults form a northwest-southeast (NW-SE) striking deformation belt overlapping the chain axis (12). The coseismic release of CO2 trapped in the crust through the fault damage zones created by the major earthquakes in SA is thought to play a major role in the temporal and spatial evolution of some seismic sequences (9, 13, 14).

Fig. 1 Seismotectonics of SAs.

(A) The 1990–2016 crustal seismicity with magnitude larger than 2 is shown in yellow (; last accessed, November 2016). The remaining colored dots indicate the most significant seismic sequences in the study area that occurred in the last 20 years. The years 1997–1998 as black dots and 2001 as orange dots indicate the seismicity from the study of Castello et al. (47). The 2005 and 2014–2016 seismicity (green and blue, respectively) are from The 2013 Matese sequence in red is from our study. Moment tensor solutions of events with Mw > 4 are from (last accessed, November 2016). Dashed black lines outline the NW-SE and SW-NE (southwest-northeast) profiles shown in (B) to (D). The white line is the orogenic divide Vezzani et al. (48). Faults are from ITHACA database Comerci et al. (49). Black triangles are the seismic stations. CF, Campi Flegrei; Ve, Vesuvius; Rm, Roccamonfina; Is, Ischia. (B) NW-SE profile showing the seismicity as described above. (C) NW-SE profile showing the 2013 Matese sequence according to the time of occurrence. The shadow gray area shows the inferred aseismic zone, whose size is ~8 km high × 2.5 km wide. (D) SW-NE profile [shown in map view in (A)] in which the seismicity is color-coded as in (B). Crustal section redrawn according to Improta et al. (11): Red lines indicate active normal faults, thrusts are shown as dashed gray lines, and the Tyrrhenian Moho is the dashed blue line. SA units include carbonates and flysch. (E) Histogram of the number of earthquakes as a function of depth. Red bars refer to the 2013 events, whereas gray bars refer to the rest of the seismicity.

Here, we provide seismic and geochemical evidence for active intrusive processes within the crust of SA mountain chain by analyzing an uncommon seismic sequence (Mwmax=5) that occurred on 29 December 2013 in the Matese region (Fig. 1). We analyze the spatiotemporal evolution of the seismic sequence, the attenuation in the crust hosting the hypocenters, and the geochemical composition of the springs in the area. We model the fluid pressure at seismogenic depth and derive the geothermal heat flux. Results identify the geophysical and geochemical signature of magma ascending in the continental crust and reveal the mechanism of pluton emplacement. Our findings open new perspectives on the triggering mechanisms of earthquakes in nonvolcanic areas and on the seismic hazard assessment of mountain chains.


The 2013–2014 Matese seismic sequence

On 29 December 2013 at 17:08 UTC, a moment magnitude (Mw) of 5 earthquake occurred beneath the Matese mountains in SA (Fig. 1A). In the following 50 days, 350 aftershocks were recorded by the seismic network of the Istituto Nazionale di Geofisica and Vulcanologia (; last accessed, November 2016). Fault plane solutions of the mainshock and the largest aftershock (20 January 2014; Mw = 4.2) indicate NW-SE striking normal fault mechanisms with a compensated linear vector dipole component of 14 and 29%, respectively (; last accessed, November 2016; Fig. 1A).

The 2013–2014 Matese sequence was concentrated at depths between 10 and 25 km in the crystalline basement of the chain just above the mantle wedge, whereas the previous seismicity in the same area was generally shallower (Fig. 1, C to E), as is the overall SA seismicity (12). The evolution of the sequence shows that the aftershocks migrate upward and spread southeastward in a few minutes (Figs. 1C and 2A). The hypocenters depict two finger-like clusters, with the two larger shocks located at the bottom of each cluster (black star and black diamond in Fig. 1, B and C). The events that occurred after the larger Mw = 4.2 aftershock cluster in the southernmost cloud. The two finger-like clusters of the Matese sequence surround an ~8-km-high × 2.5-km-long × 1.5-km-wide aseismic zone with a dike-like shape (Fig. 1C). The trend of the hypocentral depth with time in the first 28 days after the mainshock (Fig. 2A) indicates a burst-like rupture process similar to that observed in the fluid injection–induced seismicity (1517). Consequently, we interpret the Matese sequence to have been influenced by fluids during the seismogenic process. With the aim of characterizing the dispersive properties of the rocks, we computed the attenuation 1/Q0 as a function of azimuth and distance for the mainshock using the data set of Convertito et al. (18). Results show that 1/Q0 varies with distance and does not change with azimuth (Fig. 2B). At hypocentral distances shorter than 30 km, the attenuation values are ≥0.0016 with little variation. At larger distances, the attenuation values are mostly lower (<0.0016). This result indicates a higher attenuation zone with a subcircular shape around the hypocentral area, possibly reflecting the presence of fluids and/or hotter material at depth (19). We exclude the interpretation that the attenuation variation with distance is due to horizontal layering, because of the heterogeneous shallow crust in this sector of SA (Fig. 1D) (11). Spectrograms of the mainshock and the larger aftershock at the available recording stations show an overall frequency content well picked below 3 Hz, reaching 6 Hz in a few cases (Fig. 3 and fig. S1, A and B). The same is observed for the later event 20160116 (Fig. 1A for location and fig. S1C). These low-frequency events significantly differ from the slow slip events recorded in other subduction settings, because the latter are deeper and associated to tremor and their hypocenters define low-angle to subhorizontal seismogenic layers (20). Instead, the frequency content of the Matese events is similar to that recorded in volcanic areas, where seismicity is generally associated with fluid-filled cracks (21). Therefore, our data are consistent with the involvement of fluids in the 2013 Matese sequence and suggest tensile-shear ruptures at the boundary of an aseismic zone, possibly representing a magmatic body. Fluids released from this body could trigger the recorded low-frequency events surrounding this zone. To test whether these deep fluids ascend to the surface and to establish their nature, we analyze the composition of the water springs and gas emissions in the Matese area.

Fig. 2 Distance from the 29 December 2013 mainshock and depth versus time distribution and seismic attenuation.

(A) Bottom: Red dots indicate the seismicity that occurred in the first 28 days after the mainshock. The black star is the 29 December 2013 mainshock, whereas the black diamond is the larger aftershock of 20 January 2014. Top: Time-depth distribution of the events. (B) Distribution of the attenuation values (1/Q0) as a function of the station azimuth (bottom) and the hypocentral distance (top) for 34 seismic stations [see the study of Convertito et al. (18) for details]. The mean value of Q0 is 616 ± 128.

Fig. 3 Spectrograms.

Ground velocity spectrograms of the vertical component for the mainshock (20131229) and the larger aftershock (20140120) at seismic stations shown in map view in the left panel at the top. Each spectrogram is computed on a 4-s sliding window with 50% overlap in the frequency range of 0 to 10 Hz and is shown along with the correspondent waveform. Each waveform is normalized to the maximum amplitude.

Water springs and gas emissions

The 2013–2014 epicenters are located within the 812 km2 large Matese aquifer (Fig. 4A), which discharges 22,900 liters s−1 (22). The groundwater circulating in the aquifers dissolves most of the gas coming from depth, because CO2 is relatively soluble in water. Here, we estimate the fraction of deeply derived CO2 entering the aquifers by defining Cext, which is the carbon dissolved in the groundwaters derived from carbon sources external to the aquifers (see Materials and Methods for details). The origin of the carbon dissolved by the groundwaters is investigated by comparing the isotopic composition of Cext with the inverse of the concentration (Fig. 4B). In the diagram, a cluster of points defines a trend characterized by lower Cext content and more negative δ13C values, which represent infiltrating waters (green dots in Fig. 4B). This trend characterizes the groundwater dissolving the carbon from a biogenic source in the soils. When the deep CO2 is added, higher δ13C values and higher concentrations of the dissolved carbon are observed (red dots in Fig. 4B). The isotopic signature of the deep CO2 entering the aquifer, inferred from carbon dissolved in the spring, is very similar to that of the CO2-rich gas emissions located in the area (yellow squares in Fig. 4B). Therefore, we estimate δ13C values close to 0‰ for the deep component dissolved in the nearby springs (red dots in Fig. 4B). Notably, δ13C ~ 0‰ characterizes the gas emissions of the Vesuvius and Roccamonfina magmatic systems, which are located 60 to 70 km away from the Matese aquifers (23, 24). The total mass of deeply derived CO2 transported in solution by the groundwaters, Cdeep, is 380 metric tons day−1 (see Materials and Methods), whereas an amount of ~570 metric tons day−1 is estimated for the entire Matese aquifer, considering that we sampled 66% of the total water discharge. Most of this CO2 (283 metric tons day−1) is discharged by the springs whose recharge area is close to the 2013–2014 earthquakes (samples 1612, 1617, 1618, and 1628 in Fig. 4A). Similar, and also higher, amount of deeply derived CO2 is dissolved by the nearby aquifers (Fig. 4A). These results indicate that deep, plausibly magmatic CO2 reaches the surface in the Matese area.

Fig. 4 Aquifers and isotopic composition.

(A) Location of the springs of higher flow rates discharged by the Matese and nearby aquifers together with the CO2 degassing waters and the emissions of CO2-rich gas. Red dots are the 2013–2014 Matese earthquakes. (B) δCext versus 1/Cext diagram. Gas and water sampling and analysis have been carried out between September 1999 and May 2000.

Geothermal heat flux

The deep CO2 of the aquifers (Cdeep) in the Matese area could be supplied by an intrusive body releasing gas (Fig. 4A). To test whether a thermal anomaly related to a magmatic source exists, we refer to sample 1612 (Fig. 4A), which is the main spring (5500 liters s−1) in the area and contributes about 50% to the deeply derived carbon budget. Stable isotopes indicate a meteoric origin for this water (table S1) (23). The 1612 water temperature is not very high (12.4°C), but it is ~2°C higher than that expected for water infiltrating at 1300 to 1400 m of altitude (see Materials and Methods), also considering heating due to gravitational potential energy (GPE) dissipation (25, 26). The geothermal heat flow of the 1612 spring is estimated to be ~242 mW m−2 (table S2), a value evidencing the primary role of aquifers in hiding the real heat flux of mountain and permeable regions, because it is five times higher than the conductive heat flux of the region (40 to 50 mW m−2) (27). This computation indicates that the deeply derived CO2 is associated with ascent of hot fluids, because in the case of normal springs (with no deeply derived CO2), the computed geothermal heating is much lower (<80 mW m−2; see Materials and Methods). Although the heat flow estimation of the 1612 spring is affected by uncertainties that are difficult to quantify because of the low-temperature anomaly of the waters (see Materials and Methods for details), the presence of a thermal anomaly linked to the ascent of CO2-rich fluids is confirmed by the nearby CO2-rich thermal springs with lower flow rate (for example, samples 1617, 1618, and 1628 with temperatures of 20.7°C, 20.6°C, and 22.3°C, respectively; Fig. 4A and table S1). To further investigate the origin of the thermal anomaly, we examine the mixing relations among the waters discharged in this sector of Matese aquifer by plotting the Na content of the waters, Cext, and their enthalpy versus the chlorine content (Fig. 5). Figure 5 shows that sample 1612 (the largest spring with a flow rate of 5500 liters s−1; Fig. 4A) is a mixture between a normal groundwater and a thermal component relatively rich in Cl, Na, and CO2 (that is, Cext). As reported above, the thermal anomaly of sample 1612 (T spring = 12.4°C) is about 2°C. This implies that, in the mixing model shown in Fig. 5, the normal groundwater component of the 1612 mixture should be characterized by a temperature of ~10.4°C, that is, by an enthalpy of 43.7 J g−1. Note that a fluid with such enthalpy and Cl content similar to the normal groundwater (cyan square in Fig. 5C) is a very reasonable cold end-member for the mixture forming the 1612 spring. This computation confirms the reliability of our heat flux estimates.

Fig. 5 Mixing normal and thermal waters.

Variation of the Cl content with Na (A), Cext (B), and enthalpy (C) in the Matese springs and computed mixing lines. The largest spring of the area (sample 1612, 5500 liters s−1) is a mixture of ~15% thermal water (type 1628) and 85% normal cold water; the estimation of the temperature of 10.4°C was obtained on the basis of the deuterium content and is suggested to be a suitable estimate of the cold end-member involved in the mixing that formed spring 1612 [cyan square in (C)]. Spring locations are shown in Fig. 4A.


Our seismological and geochemical analyses and the results of the heat flow modeling show that the 2013–2014 Matese sequence occurred in an SA area characterized by a geothermal anomaly. The seismic events surround a dike-like aseismic zone and show tensile-shear ruptures triggered by a fluid overpressure episode. The high attenuation zone centered on the Matese epicentral area, the dike-like aseismic volume depicted by the hypocenters, the recorded heat flow anomaly, and the release of deep CO2 mirror a condition frequently observed in active volcanic areas (2730). Therefore, the above reported geophysical and geochemical data could be the signature of hot magma feeding a preexisting intrusion and generating a fluid overpressure in the SA crust. The dike-like magmatic reservoir and its fluid-rich aureole account for the seismic attenuation around the subvertical hypocentral cloud. These features cannot be explained by an input of gas alone because the ascent of gas generally produces subvertical seismic swarms and not kilometer-scale attenuation zones in the crust. The unusual distribution of low-frequency events (Fig. 1, B to D) and the burst-like time evolution of the Matese sequence (Fig. 2A) are not unique to SA. These features are also found in Japan at distances greater than 50 km from volcanoes, where seismic bursts are activated by fluid pressure redistribution in the crust (31). This process may explain the Matese seismic sequence, suggesting that the fluid overpressure may be associated with a seismogenic magma “pulse” feeding an intrusion. The roughly NW-SE alignment of the 2013–2014 Matese epicenters and the focal solutions of the main earthquakes indicate NW-SE striking ruptures (Fig. 1A, inset), which occur in response to a NE-SW regional extension dissecting SA. To open cracks and allow fluid migration, the fluid pressure Pf has to be larger than the minimum regional horizontal stress σ3. Assuming a tensile-shear rupture and according to Griffith’s criterion (see Materials and Methods) (32), the relation Pf ≥ σ3 + T, where T [20 MPa (33); see Materials and Methods] is the tensile strength of the rocks, must be satisfied. We assume the pressure source at a depth of 15 km, that is, at the top of the aseismic dike-like zone. For intact basement rocks (see Materials and Methods), we obtain Pf = 308 ± 22 MPa and σ3 = 290 ± 20 MPa at a depth of 15 km. These values indicate that, within the uncertainties related to the values of T and σ3, the condition for crack opening is satisfied (32). The 2013–2014 hypocenters concentrate at depths between 10 and 25 km, within the SA crystalline basement above the mantle wedge and below the Apulian carbonates (Fig. 1D). We suggest that the SA crystalline basement/carbonate interface is a zone of density and/or mechanical discontinuity, where the storage and accumulation of magma are allowed because such an interface may represent the level of magma neutral buoyancy (34) and/or a rigidity barrier within the crust (35). This conclusion is supported by independent seismic and geochemical data (36), showing that, when present in the SA, melts and fluid overpressure zones concentrate below crustal depths of 10 to 15 km at the interface between the crystalline basement and the overlying units.

We propose that the deep magmatic source refilling the Matese intrusion is the SA mantle wedge, where melts accumulate beneath the overriding Apenninic lithosphere (Fig. 1D) (1012, 36). According to the available geochemical data, these melts also feed the deep plumbing systems of the Roccamonfina, Campi Flegrei, and Vesuvius volcanoes, located within 100 km from the Matese area (10).

We show that signals of magmatic intrusions in mountain chains do not qualitatively differ from those observed on active volcanoes. However, in the latter case, earthquakes are generally arranged in swarms or occur uninterruptedly before an eruption at shallow depth, in the upper most levels of the crust (21). The seismic activity during the Matese sequence was not continuous, as expected in a volcanic environment, but of a “burst” type and deeper than that expected in volcanic environments (Fig. 2A). These differences in the seismic signature between active volcanic and intrusive processes allow us to conclude that 2013–2014 Matese earthquakes are related to a pulse of magma within the SA lower crust and do not record a seismically active steady magma accumulation process. The time scale to generate overpressure and trigger seismic events related to the suggested pulse of magma intrusion is on the order of 50 days, that is, the duration of the Matese seismic sequence. Assuming that the dike-like aseismic volume enclosed by the Matese sequence is occupied by a partially melted intrusive body, we estimate that the volume of this body is about 30 km3.

We remark that Apennine seismicity is generally interpreted as due to tectonic stress alone or to the interaction of tectonic stress and crustal overpressurized CO2 reservoirs located in the upper crust (8, 9, 37). Our data show that a deeper seismicity associated with the fluid release from the emplacement of intrusive bodies needs to be also taken into account. As in the case of the Matese sequence, if present, active intrusions can trigger earthquakes deeper than the shallow background seismicity and with magnitude greater than 5. As a result, earthquakes associated with active intrusions should be considered in the evaluation of seismic hazard in mountain chains. This seismicity should also be analyzed in time with the aim to detect a potential recurrence in the magma reservoir assembly process. To identify burst-like seismic sequences characterized by low-frequency content, we suggest searching the seismic record in other mountain belts like the Alpine-Himalayan, North American Cordillera, Andes, and Zagros chains. These seismic signals, if found, could provide important information on the triggering mechanisms of earthquakes related to active intrusion episodes and on the dynamics of magma ascent in the crust.


Seismological analysis

We relocated the Matese 2013–2014 seismic sequence using a double difference technique [HypoDD by Waldhauser and Ellsworth (38)]. We first located the earthquakes running HYPOINVERSE location code (39) with P and S readings derived from the Istituto Nazionale di Geofisica e Vulcanologia (INGV) bulletin ( and a one-dimensional P-wave velocity model (11). Second, we applied the HypoDD code, which is based on computing the travel time differences for event pairs at common stations. HypoDD groups the events in clusters of well-linked events (events belong to only one cluster) to have a stable inversion of the differential travel times, producing a final relocated data set that could be significantly smaller than the starting one. In our case, the final double-difference relocated catalog contains 216 events of the 300 in the initial catalog. We finally applied the conjugate gradient least square algorithm with the following parameters: minimum number of observations per event pair is eight (equal to the number of degrees of freedom for an event pair, three spatial and one temporal coordinates for each event pair), and the maximum hypocentral separation allowed between linked events is fixed at 5 km.

To better understand the recorded signals in terms of their frequency content, we computed the spectrograms of the ground velocity for the mainshock (20131229) and the larger aftershock (20140120). We considered a cutoff distance of 60 km using normalized three-component waveforms. Spectrograms were computed on a 4-s sliding window with 50% overlap in the frequency range of 0 to 10 Hz and are shown in Fig. 3 along with the waveforms for the vertical components. Spectrograms for the three components relative to the 20131229 and 20140120 events are shown in fig. S1 (A and B).

To estimate the dispersion of the crustal rocks, we computed the attenuation as the inverse of the quality factor Q0 that does not depend on the frequency. We computed the S-wave displacement spectra, assuming the source model by Boatwright (40) and applying the technique used by Convertito et al. (18). We demonstrated that 1/Q0 does not depend on the azimuth of the stations but only on the distance (Fig. 2B).

Carbon mass balance of the aquifers

The total dissolved inorganic carbon (TDIC) of groundwater from carbonate aquifers is the sum of carbon contributions from the dissolution of carbonate rocks hosting the aquifer (Ccarb) and carbon sources external to the aquifer (Cext). The carbon from the dissolution of atmospheric and soil biogenic CO2 during the rainwater infiltration process (Cinf) and the dissolution of CO2 from deep sources (Cdeep; that is, metamorphism, mantle, and magma) contribute to Cext. The carbon mass balance is described by the equations (41)TDIC =Cext+Ccarb(1)δ13CTDIC×TDIC =δ13Cext×Cext+δ13Ccarb×Ccarb(2)where δ13CTDIC, δ13Ccarb, and δ13Cext are the isotopic composition of TDIC, Ccarb, and Cext, respectively. Equation 2 is valid assuming that isotopic equilibrium exists between all the dissolved carbon species and that no isotopic fractionation occurs during the addition of any carbon sources to the solution (42). Equations 1 and 2 can be applied to the groundwater that has not experienced CO2 degassing and calcite precipitation before the emergence, a condition verified for the Apennine springs (41, 43).

The TDIC of each sample was computed by the PHREEQC code (44) using the field determinations of T, pH and alkalinity, and the major ion concentrations as input data. δ13CTDIC was analytically determined.

The Ccarb was computed from the equationCcarb=[Ca]+[Mg][SO4](3)where concentrations are expressed in molality. Equation 3 considers the dissolution of calcite and dolomite and the presence of gypsum and/or anhydrite in the aquifers. δ13Ccarb was assumed to be equal to the average value of numerous isotopic compositions of Meso-Cenozoic carbonate rocks of the Apennines (+2.21‰) (41).

The Cext and δ13Cext values were computed from Eqs. 1 and 2, whereas the different components of Cext were derived considering the following additional equation Cext=Cinf+Cdeep(4)and by the inspection of the δ13Cext versus Cext diagram (fig. S2). The lowest Cext contents and more negative δ13Cext (derived from organic sources, green dots in fig. S2) characterize the infiltrating water where Cext ~ Cinf. Other samples (CO2-rich waters, red dots in fig. S2) move from the infiltrating water trend to higher Cext values and heavier carbon isotopic compositions because of the addition to infiltrating waters of inorganic carbon from a deep source (Cdeep). The Cdeep content of each CO2-rich water is given by Eq. 4 assuming a Cinf value equal to the mean of the infiltrating waters (Cinf = 2.05 ± 0.4 mmol kg).

Finally, some samples depart from the infiltrating water–deep CO2 mixing trend toward heavier carbon isotopic composition because they are affected by carbon fractionation during water degassing and calcite precipitation (degassed waters, purple dots in fig. S2). The carbon mass balance approach is not applicable in these springs characterized by low flow rates.

Groundwater temperatures and geothermal warming

The geothermal heat flux (Qb) affecting the hydrogeological basin of the springs was computed from the difference between the emergence temperature Ts and the temperature of the recharge water TrT). ΔT (in kelvin) is linked to Qb (in watts per square meter) by the following relation (25)ΔT=Qbρw×Cw×Aq+Δz×gCw(5)where ρw (in kilograms per cubic meter) and Cw (in joules per kilogram per kelvin) are the water density and heat capacity, respectively; A (in square meters) is the extension of the hydrogeological basin; q (in cubic meters per second) is the spring flow rate; Δz (in meters) is the difference between the recharge depth zr and the spring elevation zs; and g (in meters per second squared) is the gravity. In Eq. 5, the term Qbρw×Cw×Aq is the fraction of ΔT related to geothermal heating (Gwarm in table S2), whereas Δz×gCw is the heating caused by the dissipation of GPE (table S2). Equation 5 was solved assuming constant ρw, Cw, and g (1000 kg m−3, 4186 J kg−1 K−1, and 9.81 m s−2, respectively), whereas Ts and elevations zs are from field measurements. Values of A and q were taken from literature (22). The water recharge elevation zr and temperature Tr were estimated on the basis of the δ18O isotopic composition of the springs. Elevations were inferred in the diagram z versus δ18O from the best fitting of 11 meteorological isotopic stations (fig. S3) located in southern Italy (45). Tr values (table S2) were computed by Tr = T0 + zr × ∇Tair considering an average temperature at sea level T0 of 16°C (, the above described estimation of zr, and an air temperature lapse ∇Tair = − 6.5 × 10− 3°C m–1 ( Computations are done for the three springs (table S1) for which oxygen isotopes are available (23, 46). Results indicate high heat fluxes (242 mW m−2) for the CO2-rich sample 1612, whereas much lower values are found for the normal groundwater discharged by springs 1610 and 1615. We emphasized that, due to the high infiltration of recharge waters (28 liters s−1 km−2), the temperature anomalies of the springs are relatively low and this implies some uncertainties in the estimation of Qb.

Computation of the fluid pressure

We assumed an extensional shear rupture mechanism following the Griffith’s criterion and determined the fluid pressure Pf required to activate shear failure and crack opening at a depth of 15 km in the Matese crust. The condition for an extensional-shear failure mode is given by Sibson (32)4T<σ1σ3<5.66T(6)where σ1 = ρgh, T is the tensile strength of the rocks in megapascal, σ1 is the maximum vertical stress, and σ3 is the least compressive stress in a normal stress field like that acting in the Matese area; ρ is the rock density (here assumed to be 2650 kg m−3 from gravity data), h is the depth (in meters), and g is the gravity (9.81 m s−2).

In crustal rocks, T is highly variable ranging between 1 and 20 MPa; however, for crystalline rocks, a reasonable value is 20 MPa (33). Using the selected parameters, σ1 is 390 MPa at a depth of 15 km. Therefore, σ1 − σ3 ~ 100 ± 20 MPa and σ3 = 290 ± 20 MPa. At a depth of 15 km, σ31 = 0.7 ÷ 0.8, a value that is consistent with the one recognized in extensional tectonic regimes like that acting in the Apennines (13).

According to the Griffith’s criterion, the fluid pressure for an extensional shear failure mode is (32)Pf=σ3+8T(σ1σ3)(σ1σ3)216T(7)

Assuming the above selected values of T, σ1, and σ3, we obtained Pf = 308 ± 22 MPa at a depth of 15 km.


Supplementary material for this article is available at

fig. S1. Spectrograms.

fig. S2. δ13Cext versus Cext diagram.

fig. S3. Diagram of the elevation z versus δ18O of rainwater.

table S1. Geochemical composition of the springs.

table S2. Heat balance of the Matese springs.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank A. Esposito, Z. Fu-Guo, L. Improta, G. Milano, P. Persaud, and S. Salvi for useful discussions and suggestions. We also thank the editor G. Gaetani, C. Huber, C. F. Miller, and an anonymous reviewer for their comments, which helped to clarify the manuscript. Funding: This work was supported by INGV. Author contributions: F.D.L., G.V. and N. A. P. conceived the study. F.D.L. relocated the seismicity and provided the spectral analysis. G.C., S.C., and C.C. sampled and analyzed the waters and modeled the heat flow. V.C. and N.A.P. computed the seismic attenuation. C.T. contributed to the first draft of the manuscript. G.V. provided the geological interpretation and modeled the fluid pressure. F.D.L., G.V., G.C. and N. A. P. wrote the manuscript and interpreted the results with inputs from all the authors. Competing interests: The authors declare that they have no competing interests. Data and materials availability: The seismic data used in the paper are available for download at The geochemical data are available in the Supplementary Materials.

Stay Connected to Science Advances

Navigate This Article