Gold Open Access: This paper is published under the terms of the CC-BY license.


The eastern flank of the Mount Etna stratovolcano is affected by extension and is slowly sliding eastward into the Ionian Sea. The Pernicana fault system forms the border of the northern part of this sliding area. It consists of three E-W–oriented fault sectors that are seismically active and characterized by earthquakes up to 4.7 in magnitude (M) capable of producing ground rupture and damage located mainly along the western and central sectors, and by continuous creep on the eastern sector. A new topographic study of the central sector of the Pernicana fault system shows an overall bell-shaped profile, with maximum scarp height of 35 m in the center of the sector, and two local minima that are probably due to the complex morphological relation between fault scarp and lava flows. We determined the ages of lava flows cut by the Pernicana fault system at 12 sites using cosmogenic 3He and 40Ar/39Ar techniques in order to determine the recent slip history of the fault. From the displacement-age relations, we estimate an average throw rate of ∼2.5 mm/yr over the last 15 k.y. The slip rate appears to have accelerated during the last 3.5 k.y., with displacement rates of up to ∼15 mm/yr, whereas between 3.5 and 15 ka, the throw rate averaged ∼1 mm/yr. This increase in slip rate resulted in significant changes in seismicity rates, for instance, decreasing the mean recurrence time of M ≥ 4.7 earthquakes from ∼200 to ∼20 yr. Based on empirical relationships, we attribute the variation in seismic activity on the Pernicana fault system to factors intrinsic to the system that are likely related to changes in the volcanic system. These internal factors could be fault interdependencies (such as those across the Taupo Rift, New Zealand) or they could represent interactions among magmatic, tectonic, and gravitational processes (e.g., Kīlauea volcano, Hawaii). Given their effect on earthquake recurrence intervals, these interactions need to be fully assessed in seismic hazard evaluations.


Knowledge of the slip rates of active faults is critical for understanding the earthquake recurrence rate and the tectonic driver(s) of seismicity, and evaluating the rate of interseismic strain accumulation around a locked fault. Fault slip rates are a function of the recurrence rates of characteristic earthquakes that are the starting point for time-dependent seismic hazard assessment (e.g., Field et al., 2015; Peruzza et al., 2011). Numerous studies have compared strain rates recorded over different time scales, for example, by geodesy, seismology (historical or instrumental), and geology (e.g., Friedrich et al., 2003; Pancha et al., 2006; Nicol and Wallace, 2007; Sue et al., 2007; Calais et al., 2008; Rontogianni, 2010; Mouslopoulou et al., 2012). The results suggest that there are cases in which geodetic and seismic/geological measurements are similar (e.g., Hindle et al., 2000; McCaffrey, 2005; Nicol and Wallace, 2007). In other cases, however, significant discrepancies were observed (e.g., Donnellan et al., 1993; Liu et al., 2000; Papanikolaou et al., 2005; Mouslopoulou et al., 2012). These discrepancies reflect limitations and biases in the different measurement techniques used (Gardner et al., 1987; Wang et al., 2011), coupled with other processes, such as fault interactions, magmatic intrusions, rotations around the vertical axis, microseismicity, or changes in the plate-velocity vectors, occurring in addition to fault slip (Friedrich et al., 2003; Nicol et al., 2009; Mouslopoulou et al., 2009).

In general, long-term observation of faults in various tectonic contexts has shown that slip rates vary in space and time and that fault activity is inherently episodic (e.g., Bull et al., 2006; Nicol et al., 2006, 2010; McClymont et al., 2009; Schlagenhauf et al., 2011; Mouslopoulou et al., 2012; Gunderson et al., 2013; Benedetti et al., 2013). Numerical simulations have shown that slip rate variability mainly occurs in response to interactions between adjacent active faults (e.g., Robinson et al., 2009; Cowie et al., 2012; Visini and Pace 2014).

Understanding the temporal variability of fault slips may also provide insights into the associated geodynamic processes. In volcanic areas, active tectonics and seismotectonic approaches are complicated by the interplay of different deformational processes (e.g., Norini and Acocella, 2011). For example, detection and monitoring of gravity-driven volcano deformation are vital for understanding volcanic hazards such as landslides, lateral blasts, and debris avalanches. Deformation has been detected at several large active volcanoes, e.g., Mount Etna (e.g., Solaro et al., 2010), Vesuvius (e.g., Tammaro et al., 2013), and Kilauea (e.g., Bagnardi et al., 2014). However, these systems also exhibit persistent eruptive activity that obscures the gravity-driven signals of ground motion (e.g., Shirzaei et al., 2013). Evaluation of the way in which volcano-tectonic structures and magmatic activity control deformation processes yields fundamental information for evaluating the seismic hazards in volcanic regions. Variability in slip rates and earthquake recurrence could reflect a number of processes, including temporal changes in frictional properties of the faults, variability in the stress interactions between the faults, and/or temporal changes in the strain partitioning to nearby faults and/or from nearby eruptive activity (Nostro et al., 1998; Nicol et al., 2006; Rosenau and Oncken, 2009). Fault rupture can be aided or suppressed by instantaneous changes in stress that impose an irregular timing on the rupture of a specific fault (Stein, 1999; Toda et al., 2012).

Herein, we combine a detailed study of the morphology of the most active sector of the Pernicana fault system of Mt. Etna with new age determinations of faulted lava flows in order to construct a model of the spatial and temporal variation in the fault slip rate. This morphochronological approach constrains fault activity over a time period that is longer than can be assessed by instrumental measurements or historical data. Our results may be viewed in the light of fault-based seismic hazard estimates that use individual seismogenic source parameters to evaluate the expected seismicity rates.


The volcanic region of Mt. Etna is characterized by intense ground deformation (Bruno et al., 2012) that is due to the interaction of regional extensional tectonics (Palano, 2015) with flank instability processes and magma intrusion dynamics that affect the volcanic edifice (Azzaro et al., 2013, and references therein). The volcano-tectonic framework is controlled by structural elements mostly located in the eastern sector of the volcano that, although characterized by relevant dynamics expressed by shallow seismicity and surface faulting (both coseismic and creeping), show a different morphology (Fig. 1). The lower southeastern flank is displaced by the Timpe fault system, a suite of right-lateral extensional structures that generally strike north-south. These form several steep escarpments, up to 5–8 km long and up to 200 m high, that offset late Pleistocene to historical-age lava flows. On the basis of the ages of offset geological markers and geodetic measurements, fault slip rates range from 1 to 4 mm/yr (Azzaro et al., 2012).

The Pernicana fault system strikes E-W and dissects the entire northeastern sector of the volcano from the NE Rift to the Ionian Sea (Fig. 1). It represents the northern boundary of the unstable eastern flank of the volcano that is sliding seaward (Lo Giudice and Rasà, 1992; Rust and Neri, 1996; Azzaro et al., 1998). The system has a total length of ∼20 km, and it is formed by two main segments (Fig. 1; Azzaro et al., 2012); the western segment has an individual steep fault escarpment, and the eastern segment has a discontinuous fault zone devoid of morphological expression. In this work, we focused mainly on the western segment of the fault system, which is considered to be the seismically active one, distinguishing it into three sectors (Fig. 2A) that differ in morphology and kinematics along strike. Fault movement is dip slip in the western sector, sinistral oblique slip in the central sector, and purely left lateral in the eastern part (Bonforte et al., 2007). Fault dynamics are characterized by continuous creep movements associated with the eastern flank sliding (Rasà et al., 1996; Azzaro et al., 2001; Acocella and Neri, 2005), punctuated by shallow earthquakes (Mmax 4.7, depth 1–4 km) producing coseismic fault displacements along the western and central sectors of the fault system (Azzaro, 1997, 1999; Obrizzo et al., 2001; Bonforte et al., 2007; Guglielmino et al., 2011). The seismic activity occurs in connection with the intrusion processes that lead to flank eruptions, as well as nonvolcanic causes (Alparone et al., 2013; Bonaccorso et al., 2013).

At the western upslope termination around Piano Provenzana (Fig. 2A), the fault is manifest as a NE-SW–trending morphological flexure, since the faulted lavas are buried under recent and historical lava flows erupted from the NE Rift, to which it appears structurally connected (Azzaro et al., 2012; Tibaldi and Groppelli, 2002; Neri et al., 2004). A few kilometers east, it manifests as a N110°-trending fault scarp, the hanging wall of which hosts sag ponds filled by alluvial deposits (Piano Pernicana and Mandra del Re basins; see Fig. 2A). The fault scarp progressively disappears at Rocca Campana (Fig. 2A), where the western sector of the Pernicana fault system terminates by branching out eastward into a series of splay faults. The eastern sector of the Pernicana fault system develops for 9 km from Rocca Campana to the coast with a general NW-SE trend (Azzaro et al., 2012). It has no morphological manifestations along strike, being instead characterized by hidden features revealed by creep phenomena that have damaged man-made features (Rasà et al., 1996; Azzaro et al., 1998; Neri et al., 2004). An extension of the Pernicana fault system is recognized offshore of Fondachello (Chiocci et al., 2011).

Short-term (annual to decennial) slip rates have been generated from field observations of creep effects (Rasà et al., 1996; Neri et al., 2004) and geodetic measurements (Azzaro et al., 2001; Obrizzo et al., 2001; Bonforte et al., 2011). These mainly concern the eastern sector of the Pernicana fault system and refer to deformation accommodated by a continuous stable sliding, where creep rates range from 20 to 28 mm/yr for the horizontal (sinistral) component of displacement, and 10–20 mm/yr for the vertical component. A slip rate value of 24 mm/yr over the last 1 k.y. has been obtained by paleoseismological trenching (Ferreli et al., 2002; Mandra del Re, see Fig. 2A). On the basis of the fault scarp heights, measured using topographic profiles calculated from digital elevation model (DEM) data, and of the lava flow ages from the geological map (Branca et al., 2011), Azzaro et al. (2012) estimated vertical slip rates ranging from 3.3 to 5.2 mm/yr in the eastern sector of the Pernicana fault system during the Holocene (15–3.9 ka). Despite being the most seismically active part of the Pernicana fault system, the long-term slip rate of the central sector is largely unknown.


The height of the scarp is critical for determining the vertical component of the slip rate (throw rate) of an active fault. We measured the height of the fault scarp on the central sector of the Pernicana fault system from 12 high-resolution across-strike topographic profiles (see numbers in Fig. 2B for location) using a Leica total station and global positioning system (GPS) hardware. Heights were measured every 1 m, with an uncertainty of ∼1 cm (Figs. 3–4). Topographic profiles were determined approximately every 500 m, covering a distance of ∼6 km, along the topographically significant part of the fault scarp.

The along-fault profile of the central sector is shown in Figure 3. Five of the topographic profiles are presented in Figure 4 as examples. The height of the fault scarp shows significant variation along strike. The height reaches a maximum (35.5 m) in the central-western sector, in the area around Piano Pernicana (profile 3). Scarp height decreases continuously to the east and west of Piano Pernicana, disappearing at the SP59/IV road, where the fault has mainly strike-slip kinematics (Fig. 5), and in the western sector, where the fault scarp intersects the NE Rift. Local minima in scarp height are present, one to the east and one to the west of Piano Pernicana (profiles 5 and 9 in Fig. 3). These probably resulted from the interaction between the fault scarp and lava flows from the NE crater region. It is likely that the lavas that flowed from south to west have been dammed by, and ponded behind, the E-W–trending fault, modifying the scarp morphology (Azzaro et al., 2012).


By integrating the mapped lava flows (Branca et al., 2011) with detailed field mapping and trace-element geochemistry, we identified five lava flows that are displaced to different degrees in the central sector of the Pernicana fault system: the Piano Pernicana (PP), Pineta di Linguaglossa (PL), Monte Ponte di Ferro (FR), Due Monti (DM), and Millicucco (ML) lavas (Fig. 2B). The major, trace, and rare earth element compositions of these lavas are summarized in Appendix A.1

The oldest lava flows are related to the Mongibello volcano activity (<15 ka; Branca et al., 2011) and crop out along the footwall at the Mandra del Re locality (flows rx, ry, rw on Fig. 2B). These flows are covered by the PP flow, which originated from a fissure along the footwall at the Piano Pernicana locality. This is overlain by the plagioclase-phyric PL flow. These lava flows postdate the rheomorphic lava flows of the final stage of the Ellittico volcano activity (ca. 15 ka; Ragabo member in Fig. 2B) and predate the picritic scoria-fall deposit dispersed on the eastern flank of the volcano (fall stratified pyroclastic deposit [FS layer]) that has been dated by 14C to be 3930 ± 60 yr (Coltelli et al., 2000). Toward the western tip of the Pernicana fault system, the PL flow is covered by the FR flow, which is younger than the fall-pyroclastic lithic-rich deposit (FL) tephra marker bed dated 3150 ± 60 yr (Coltelli et al., 2000; Branca et al., 2011). Two historical lava flows, DM and ML, showing an Early Medieval age (Branca et al., 2015), are displaced between Rocca Pignatello and Rocca Campana (Fig. 2B). Both flows are plagioclase- and pyroxene-phyric ± olivine, and they have an a’a surface morphology, with some portions showing pahoehoe morphology and sparse vegetation cover in places.

Herein, we present new cosmogenic 3He exposure ages for the DM and ML lava flows and a 40Ar/39Ar age for the older PL lava flow. The PP and FR lava flows were unsuitable for dating, though the newly determined age of PL allows new constraints to be determined for the age of the PP lava flow.

DM and ML Lava Flows

Multiple samples of uneroded pahoehoe ropes were collected from the DM and ML lavas. Samples were crushed, rinsed in deionized water, and dry-sieved to isolate the 500–1000 μm size fraction. Pyroxene was purified by magnetic separation, and ∼1 g was handpicked under a binocular microscope. The pyroxene separates were then leached ultrasonically in 10% nitric acid solution to remove adhering groundmass. Individual samples were crushed in a vacuum in an all-metal multisample hydraulic crusher for determination of the magmatic 3He/4He (Stuart et al. 2000). To extract the lattice-hosted helium, crushed and uncrushed grains were melted for 10 min using a 25 W 808 nm diode laser (Foeken et al., 2006). All samples were reheated to ensure total degassing. Replicate values were determined for all samples. The released gases were exposed to two SAES GP50 getters and then a liquid-nitrogen–cooled charcoal trap. Helium isotopes were measured employing a Mass Analyser Products (MAP) 215–50 noble gas mass spectrometer that was calibrated relative to the HESJ standard (He standard of Japan; Matsuda et al., 2002). Fifty calibrations were performed over the period of analysis, with a reproducibility of 0.7% for 3He and varying from 0.3% and 0.9% for 4He (1σ). Analytical procedures were similar to those used in previous studies (Foeken et al., 2009). The 3He and 4He blanks were measured by heating similar masses of degassed pyroxene for the same time as the sample gas extraction (1.9 (±0.4) × 109 atoms of 4He and 1.8 (±0.7) × 104 atoms of 3He). Analytical uncertainties were propagated from uncertainties in the 3He and 4He concentrations, blank corrections, and standard calibration measurements and are presented in Table 1 (1σ).

The concentration of cosmogenic He (3Hecos) was calculated using the equation from Kurz et al. (1990): 

where 4Hemelt and 3He/4Hemelt are the He concentration and isotopic ratio of the melt steps, respectively, and 3He/4Hecrush is the isotopic ratio of the magmatic He released by crush extraction. Age calculations were computed using the CRONUScalc exposure age online calculator ( This incorporates spallation production for 3He, wherein the modern production rate is calibrated from a combined pyroxene compilation, as described by Goehring et al. (2010). For detailed information on the 3He calibration, see Borchers et al. (2015).

Elevation and latitude scaling of cosmogenic 3He production to the sample locations was calculated using the time- and nuclide-dependent scaling scheme described in Lifton et al. (2014). This model is based on equations from a nuclear physics model incorporating dipole and nondipole magnetic field fluctuations and solar modulation. The scaling factor is based on atmospheric pressure, solar modulation, and cutoff rigidity calculated using trajectory tracing. Nuclide-dependent scaling is implemented by incorporating cross sections for the different reactions (Marrero et al., 2016).

The attenuation of the cosmogenic neutron flux with depth from the surface was corrected for using a 3 cm sample thickness, rock density of 2.7 g/cm3, and attenuation length of 160 g/cm2 (Dunne et al., 1999). The topographic prominence of the sampled sites is not conducive to snow and ash accumulation due to preferential wind erosion. On the basis of the modern snow cover record for Mt. Etna ( and modeled historical ash cover (e.g., Coltelli et al., 2000), we determined that the impact of shielding due to snow and ash covers can be excluded such that no correction for ages would be necessary for burial by snow or ash. The well-preserved lava features of the sampled surfaces suggest that they have not been affected by significant erosion. Neither topographic shielding nor surface dip correction was required. The data are presented in Table 1, in accordance with standard protocols (Dunai and Stuart, 2009).

The 3He/4He ratios of the crushed grains ranged from 5.4 to 7.2Ra, where Ra is the 3He/4He of air (1.384 × 10–6). These are typical of Holocene basalts for Etna (e.g., Martelli et al., 2008; Coulson et al., 2011). The 3He/4He obtained from the melted samples ranged between 7 and 20 times Ra (Table 1). Sample PIGNA yielded the highest magmatic 3He/4He (3He/4Hecrush = 7.2 ± 0.3Ra), which overlaps with the relative 3He/4Hemelt. It is thus not possible to differentiate the cosmogenic 3He from the magmatic 3He component, resulting in large uncertainties in the exposure ages. Weighted averages of cosmogenic ages calculated from the measured cosmogenic He values are 1196 ± 103 yr for the ML lava flow and 1356 ± 294 yr for DM lava flow (Table 1). These results are in agreement with the recent age determinations of the DM and ML lava flows performed by J.C. Tanguy and M. Condomines, as reported in Branca et al. (2015), indicating archaeomagnetic ages of 1516 ± 50 yr and 1316 ± 100 yr, and 226Ra-230Th ages of 1456 +260/–290 yr and 1286 +260/–280 yr, respectively.

PL Lava Flow

The collected basalt samples showed no alteration in hand sample and had phenocrysts ≤1 mm in diameter. Groundmass separates were prepared by standard techniques, including crushing, sieving, dilute nitric acid leaching, and magnetic separation. A final stage of handpicking under a binocular microscope ensured removal of all phenocrysts and alteration phases. Following preparation of the purified separate, 500 mg of material were loaded into a copper packet for irradiation. Samples and neutron flux monitors were placed in foil packets and stacked in quartz tubes. The relative positions of the packets were then carefully measured for later reconstruction of neutron flux gradients. The sample package was irradiated in the Oregon State University reactor, cadmium-shielded facility. Alder Creek sanidine (1.2056 ± 0.0019 Ma; Renne et al., 2011) was used to monitor 39Ar production and establish neutron flux values (J) for the samples. Estimated errors in the neutron flux measurements were within 0.5%. Gas was extracted from samples using an all-metal resistively heated furnace with multistep heating schedules ranging from ∼700 °C to 1500 °C. Liberated argon was then purified of active gases (e.g., CO2, H2O, H2, N2, CH4) using three Zr-Ti-Al getters, one at 25 °C and two at 400 °C. Data were collected on a GVi Instruments ARGUS-5 multicollector mass spectrometer using a variable-sensitivity Faraday collector array in static collection (non-peak-hopping) mode (Sparks et al., 2008; Barfod and Fitton, 2014). Time-intensity data were regressed to t0 with second-order polynomial fits to the data. Mass discrimination was monitored by comparison to running-average values of an air standard. The average total system blanks for furnace extractions, measured between each sample run, were 2 × 10–16 mol for 40Ar, 1 × 10–17 mol for 39Ar, and 2 × 10–17 mol for 36Ar over the temperature range applied for extracting age information. All data were corrected for blank and mass discrimination. Plateau criteria were: (1) the fraction of 39Ar released was ≥50%, (2) the plateau consisted of at least three contiguous steps, and (3) the scatter between the ages of the steps was low, i.e., mean square of weighted deviates (MSWD) close to 1. Isochrons were calculated using only the plateau steps to confirm the composition of the trapped component. A plateau age was accepted if it showed concordance at the 2σ level with the isochron age, had a trapped component indistinguishable from air (298.56 ± 0.31 [1σ]; Lee et al., 2006) at the 2σ level, and met the other three criteria listed.

Three aliquots of the PL flow were dated separately and yielded concordant plateau ages of 9.7 ± 1.6 ka (MSWD = 1.06, p = 0.38, n = 7), 9.5 ± 1.1 ka (MSWD = 0.98, p = 0.44, n = 8), and 10.5 ± 1.3 ka (MSWD = 0.70, p = 0.71, n = 10). All uncertainties were 1σ. The data, plots for individual samples, and the parameters used for the age calculation are presented in Appendix B (see footnote 1). Using only the plateau steps, each experiment yielded age-concordant inverse isochrons and had atmospheric intercepts for the trapped components at the 2σ level. Combining the plateau steps for each experiment and weighting each step (n = 25) by the inverse of variance yielded a “composite” plateau age of 9.9 ± 0.7 ka (MSWD = 0.83, p = 0.70, n = 25) and an associated “composite” isochron of 9.9 ± 2.9 ka (MSWD = 0.87, p = 0.64, n = 25). The composite plateau age is the best estimate of the age of the PL flow. This new result is consistent with the stratigraphic age reported for this lava flow by Branca et al. (2011). Furthermore, it constrains the age of the underlying PP lava flow to the earlier stage of the Mongibello volcano activity, specifically to values between 9.9 ± 0.7 ka and <15 ka.


Slip Rates and Implications

Using the measured vertical displacement of the dated lava flows, the throw rate, and its uncertainties, in the central sector of the Pernicana fault system was determined (Table 2). These findings were integrated with published data from the eastern sector using previously dated lava flows (e.g., Azzaro et al., 2012), geodetic measurements (e.g., Tibaldi and Groppelli, 2002), or permanent scatter analysis (e.g., Bonforte et al., 2011) to improve the space and time scales of observation. In particular, we integrated the new data with three points of measure, one in the Vena-Presa area (the eastern part of the investigated central sector) that gives 100 mm of throw in 60 yr (displacement of a stone concrete wall) and the other two along the Fiumefreddo fault (the eastern segment of the entire Pernicana fault system) that give 1.5 mm/yr of throw rate in a time window of 5 yr (interferometric synthetic aperture radar [InSAR] data in the time window 1995–2000) and ∼1 mm/yr in the interval 15–35 ka (Azzaro et al., 2012). Because a direct comparison of slip rates obtained by different methods through various time intervals may be misleading (Nicol et al., 2009), we applied a statistically derived approach (Gardner et al., 1987; Finnegan et al., 2014; Gallen et al., 2015) to evaluate possible bias in our measurements. Throw data versus measurement intervals were plotted in a log-log graph (Fig. 6). The scaling function was explored using a classical least-squares regression (LSQ) with two weighted regression methods: least absolute residuals (LAR) and bisquare (BIS), in order to minimize the effect of outliers. The LAR method finds a curve that minimizes the absolute difference of the residuals, rather than the squared differences. Therefore, extreme values have a lesser influence on the fit. The BIS method minimizes a weighted sum of squares, where the weight given to each data point depends on how far the point is from the fitted line. Points near the line get full weight; points farther from the line get reduced weight; and points that are farther from the line than would be expected by random chance get zero weight. The BIS weight method is preferred over LAR because it simultaneously seeks to find a curve that fits the bulk of the data using the usual LSQ approach, and it minimizes the effect of outliers. In Figure 6, we report the three different regressions, with the three slope values (LSQ = 0.98, BIS = 0.97, LAR = 0.81). According to the method of Gardner et al. (1987), a regression slope equal or close to 1 indicates the throw rates to be independent of time interval, whereas smaller values might have suggested a dependency by the time interval of observation. Throw rates (TR) were standardized by the scaling function: 

where Δt is the time interval over which the rate is computed, ΔH is the corresponding fault scarp height, and 0.97 is the BIS scaling value (Table 2). We also explored the effect of the LAR scaling value (0.81), but throw rate values of ∼40 mm/yr in the time window of 1000 yr would be needed, as well as erosion of the footwall of the fault by ∼30 m in 1000 yr, and these values are difficult to assume for basalt rocks, or other complex, not recorded, geomorphological process.

To unravel the variation in throw rates along strike, we plotted the integrated and standardized data set in Figure 7. The throw rate of ∼9–10 mm/yr was evaluated at the boundary between the western and central sectors, whereas rates ranging between 2 and 5 mm/yr characterize the westernmost part of the central sector. From the central sector to the easternmost tip, throw rates decrease from ∼15 mm/yr to ∼2 mm/yr at 10–12 km along the strike of the fault system. Because the fault exhibits structural continuity along the strike, and no barriers have been recognized as separate/partition slip distribution, we assumed the throw rate distribution to be consistent with the expected throw rate distribution on an active normal fault segment (Nicol et al., 2005; Manzocchi et al., 2006; Faure-Walker et al., 2009). Assuming this to be the case, variation of throw rates along the strike could be linked to along-strike changes in fault kinematics. In fact, the high vertical displacement rates in the central sector coincide with dominantly normal slip behavior and the locations where the scarp heights are highest (Figs. 3 and 4). The low vertical displacement rates at the eastern end of the fault coincide with dominantly aseismic strike-slip behavior (Fig. 5). Azzaro et al. (2001) previously suggested that the fault near Vena slips at a higher horizontal displacement rate than those in the central and western sectors. Our results are also consistent with the findings of Alparone et al. (2013), who showed the focal solution to mainly be characterized by normal movements in the central sector of the Pernicana fault system and left-lateral kinematics in the eastern sector. Finally, Bonforte et al. (2007) suggested that the along-strike thickness variation of the volcanic pile plays a role in defining the kinematics of the Pernicana fault system. Along the western sector of the Pernicana fault system, the thickness of the volcanic pile ranges from 200 m to 400 m, allowing brittle deformation to dominate. Conversely, the eastern sector of the Pernicana fault system crosses the prominent culmination of the marly-clay basement at Vena. At this location, the volcanic sequence is thinner than in the west (from 0 to 100 m), and brittle deformation is influenced by the plastic behavior of the underlying basement, accommodating the deformation along the southern block with the development of short local structures with different orientations and kinematics, which border the area of minor thickness.

However, our new data, also standardized, support a significant and consistent temporal variation in the throw rate. In fact, the throw along the Pernicana fault system appears to have been at a rate of 2–5 mm/yr integrated over the last 10–15 k.y., but of the order of 6–15 mm/yr since 3 ka (Fig. 8). Assuming the Pernicana fault system to be a continuous structure, throw rates obtained by the lava flows from the same temporal range were interpolated to obtain curves representing the throw rate for along-strike profiles of the Pernicana fault system at two temporal intervals (dashed curves in Fig. 7). From these curves, we infer that: (1) the central sector of the Pernicana fault system has been characterized by a higher throw rate than the western and eastern sectors, at least for the last 15 k.y., and (2) there was a significant change in the throw rate at some time point in the last 3 k.y.

The throw rate history of the central sector of the Pernicana fault system over the last 10 k.y. is shown in Figure 9. The throw rate integrated over the past 3.5 k.y. is 7.5 mm/yr. Had this remained constant over the last 10 k.y., the throw rate profile should have shown a smooth transition toward the western tip (Fig. 7), and the PL lava flow would now be offset ∼75 m (Fig. 9). Conversely, the throw rate profile in Figure 7 shows abrupt change at ∼4 km and ∼6 km from the western termination, which can be explained as a consequence of throw rate variation through time. In order to investigate whether throw rates have changed over the last 10 k.y., we compared the total throw and throw rate profiles (Fig. 9). We assumed a scenario wherein the throw rate was initiated at 10 ka, and then we calculated throw profiles that would develop if the throw rates for the last 1 k.y. had operated for 10 k.y. (see dashed line in Fig. 9). For a constant throw rate, the calculated throw profile should be identical to the measured throw profile. However, we find that the predicted throws are larger than the measured throws. Therefore, the throw rates for the last 1 k.y. cannot have produced the pattern of throws for the last 10 k.y. The throw profile in Figure 7 is, however, consistent with throw rates that have increased through time (Fig. 9).

Therefore, the Pernicana fault system experienced a period of accelerating fault slip beginning at ca. 3.5 ka. This contrasts with the evidence from the eastern sector, which appears to have slipped at a constant throw rate of ∼3 mm/yr over the last 15 k.y. (Fig. 7). The reconstructed slip rate history is broadly consistent with the average rates determined by Azzaro et al. (2013), calculated using the scarp heights from DEM data with the known lava flow ages in the long term, and with geodetic data in the short term.

Possible Causes of Slip Rate Variability

Increased flank instability, perhaps due to the growth of the volcano and/or changes in magma injection rate, may have initiated the observed complex kinematics at 15–3 ka, resulting in compensation of the throw rate deficit with increased horizontal movement. Without a more detailed throw rate history and greater knowledge of the interactions among magmatic, tectonic, and gravitational processes, it is difficult to evaluate the source(s) of throw rate changes. However, it might be feasible to discriminate between internal and external causes based on the time scale of the variability. For example, Naylor and Sinclair (2007) described a simple relationship that can be used to determine if variable slip behavior is caused by external forcing, such as changing surface loads or convergence rates, or processes inherent to the fault network, such as variable fluid pressure or slip partitioning on adjacent structures. They developed a numerical model of an orogenic wedge that allowed the observation of punctuated thrust activity on individual discrete thrusts, despite constant convergence rates and the absence of surface processes. They described the threshold time scale of variability as: 

where TS is the characteristic time scale of internal variability, LF is the cross-sectional length of the fault, and CR is the regional rate (Gunderson et al., 2013). This threshold time scale represents the maximum periodicity at which a fault is expected to display internal variability in the absence of changing external boundary conditions. Variable slip behavior with a periodicity longer than TS requires an external cause, while slip rate variability shorter than TS can result from either external or internal causes.

Alparone et al. (2013) have shown that seismicity linked mainly with the tectonic activity of the Pernicana fault system releases energy to depths of less than 4 km. Assuming LF to be 4 km and average long-term throw rates of 5 and 10 mm/yr (Fig. 7), the characteristic threshold time scale of internal variability (TS) for the Pernicana fault system would be 400–800 k.y. This means that the frequency variability occurring on 10 k.y. time scales in the model slip history for the Pernicana fault system does not need an external cause, and this leads us to conclude that the variability is most likely related to internal fault behavior or the volcanic system, for example, inflation, deflation, and dike injection. While dike intrusions not leading to eruptions, and associated with tectonic and seismic activities, are apparently quite rare in recent times at Etna, geophysical monitoring has revealed evidence of such an event in 1998 (Bonaccorso and Patanè, 2001; Puglisi and Bonforte, 2004). The interaction between gravity instability and magma intrusion is thus widely acknowledged for the eastern flank of the volcano, showing a continuous E-SE sliding (Bonanno et al., 2011; Bonforte et al., 2011, and reference therein; Privitera et al., 2012; Bozzano et al., 2013). The seismic activity along the Pernicana fault system may be another expression of flank unrest (e.g., Acocella and Puglisi, 2013). Causal relationships between Pernicana fault system seismicity and the eruptive activity of Etna have not as yet been well defined (Feuillet et al., 2006; Bonaccorso et al., 2013). Earthquakes or seismic swarms associated with shallow, volcano-tectonic fault systems have been noted to precede volcanic eruptions by days (Azzaro, 1997; Neri et al., 2004) but may also be unrelated (Guglielmino et al., 2011). Bonforte et al. (2011) have, however, shown that the kinematics of the fault can be perturbed by additional stress related to dike intrusions in the upper part of the volcano. This is indicated by an abrupt increase in the slip rate concomitant with eruptive events, highlighting a complex relationship between flank slip and volcanic activity. For example, the Pernicana fault system remained stable during the shallow intrusions feeding the 2004–2005 eruption period (Burton et al., 2005), and minor acceleration was observed during the 2001 flank eruption, which was fed by a deep intrusion located on the upper southern flank of the volcano (Musumeci et al., 2004). This contrasts with high slip rates observed on the Pernicana fault system associated with magma intrusion in the NE Rift during, for example, the onset of the 2002–2003 flank eruption (Branca et al., 2003; Neri et al., 2004). The abundance of explosive volcanism that occurred between 3 and 4 ka in this sector of Mt. Etna (Coltelli et al. 2000; Branca et al., 2011) supports our contention that magma intrusion is a plausible mechanism for the acceleration observed at ca. 3.5 ka.

Implications for Expected Seismicity Rates

Because the seismic moment rate is a function of the slip rate (e.g., Brune 1968), the recognition of throw rate variation in the last 10 k.y., and the complex interaction between volcanic and tectonic processes imply that the release of seismic moment is not constant through time. This has implications for seismic hazard estimations. In fact, with a fault-based approach for evaluating seismic hazard, it will be essential in the future to estimate the rates at which seismicity is released. While remaining conscious that we are not able to separate out the aseismic contribution in the slip rate computations in complex systems such as Etna, which could be important, we tried to evaluate the “maximum” expected seismicity rates using a conservative approach.

Pace et al. (2016) formalized a geometrical and energetic approach in which, using fault geometry and the slip rate, it is possible to evaluate expected rates of seismicity for various magnitude-frequency distributions as well as the maximum magnitude allowed by fault geometry. We used two magnitude-frequency distributions to illustrate the impact of the observed throw rate variation on the estimation of the expected seismicity. Assuming that the central sector of the Pernicana fault system is ∼7 km along strike (Figs. 2 and 3) and that the length along dip is 4 km, we obtained an expected maximum magnitude of 5.1 ± 0.4 (Mw) using the FiSH code (Pace et al., 2016). The link between slip rates and earthquake occurrences is schematically shown by the displacement-time curve in Figure 10. The synthetic slip-time curve was achieved by randomly sampling pairs of slip per event (in the range 0.1–0.5 m) and recurrence intervals (from a few to hundreds of years), and this curve is consistent with historical and instrumental seismic observations (Azzaro, 1997, 1999; Obrizzo et al., 2001; Alparone et al., 2013). This example satisfies the observed throw rate variability (∼7.5 mm/yr and ∼1 mm/yr) and a uniform throw rate of ∼2.5 mm/yr. The displacement-time curve shows the variable earthquake behavior resulting from the complex volcano-tectonic system and is a critical determinant of the expected seismicity rate evaluation. We highlight in the figure an example in which an earthquake clustering produced slip rate variations at wavelengths smaller than the intermediate-term curve. In Figure 11, we show the expected seismicity for seismic release based on throw rates of (1) 2.5 mm/yr (representing the long-term rate), (2) 7.5 mm/yr (representing the intermediate-term rate in the last 1 k.y.), and (3) 1.0 mm/yr (representing a new intermediate-term cycle with a low deformation rate, similar to the period between 10 ka and 3.5 ka). These values were calculated for two magnitude-frequency distributions: a characteristic earthquake distribution (Schwartz and Coppersmith, 1984) with a range of possible magnitudes around a central value (5.1), and a truncated Gutenberg-Richter magnitude-frequency distribution (Kagan, 2002) with magnitudes ranging from 4 to 5.5. These throw rate differences result in large differences in expected seismicity rates (Fig. 11). In the case of a characteristic earthquake distribution, an earthquake with a magnitude of at least 4.7 (the magnitude of the last major earthquake on Pernicana fault system, in 2002) is expected to occur on average once every 95 yr using the long-term average (scenario 1), 32 yr for the short-term rate (scenario 2), and 237 yr for the low-deformation-rate scenario (3). In addition, this model predicts that every 50 yr, these earthquakes are likely to occur 0.5, 1.6, and 0.2 times for cases 1, 2, and 3, respectively. On the assumption of a truncated Gutenberg-Richter magnitude-frequency distribution, one earthquake with a magnitude of at least 4.7 can be expected, on average, every 79 yr, 26 yr, and 198 yr, respectively (Fig. 11). This means 0.6, 1.9, and 0.2 earthquakes, respectively, every 50 yr.


Slip rates on active faults control seismic hazards (Main, 1996; Roberts et al., 2004; Bull et al., 2006; Visini and Pace, 2014) and provide a time scale with which to assess the mechanisms operating during continental deformation (e.g., Cowie et al., 2005). We have used a morphochronological approach to study the slip rate variability of the Pernicana fault system (Mt. Etna) during the last 10 k.y. New age constraints of the lava flows cut by the Pernicana fault system have allowed the long-term vertical slip rate for the last 10 k.y. to be estimated. The data identify an apparent acceleration starting at ca. 3.5 ka. These results, together with the contemporary data, may have implications for fault-based seismic hazard estimations, which use individual seismogenic source parameters to evaluate the expected seismicity rates. Some questions remain, especially as regards the correct evaluation of the causes of slip rate variability, suggesting the need for numerical studies to evaluate how the interaction between volcanic activity and faults influences this variability. However, the methodology used and the results obtained herein hold promise for future studies on other active faults on volcanoes.

Variations in slip rates due to fault interdependencies (e.g., across the Taupo Rift, New Zealand; Nicol et al., 2006) or to interactions among magmatic, tectonic, and gravitational processes (e.g., Kīlauea volcano, Hawaii; Bagnardi et al., 2014) could generate complex fluctuations in the timing and magnitude of earthquakes that need to be considered in earthquake forecast modeling in volcanic contexts.


Holocene slip rate variability along the Pernicana fault system (Mt. Etna, Italy): Evidence from offset lava flows

D. D’Amato, B. Pace, L. Di Nicola, F.M. Stuart, F. Visini, R. Azzaro, S. Branca, and D.N. Barfod

(doi: 10.1130/B31510.1)

The text and references on page 313 have been modified as follows:

However, it might be feasible to discriminate between internal (e.g. variable fluid pressure or interaction with adjacent structures) and external causes (e.g. changes in convergence/extensional rates) based on the time scale of the variability. To this aim, we followed the approach suggested by Naylor and Sinclair (2007), who developed a numerical model of an orogenic wedge, with constant convergence rates and absence of surface processes, but variable slip rate behavior. In particular, it is possible to calculate the threshold time scale of variability as: 

where TS is the characteristic time scale of internal variability, LF is the length of the fault along dip, and CR is the regional rate (Naylor and Sinclair, 2007). In this view, TS represents the threshold of periodicity and if the measured slip rate variability results longer than TS, external causes need to be hypothesized, while slip rate variability shorter than TS can result from either external or internal causes.

We are deeply grateful to M. Condomines and J.C. Tanguy for very helpful comments on earlier versions of the manuscript. We thank B. Goehring, S. Medynski, V. Mouslopoulou, and Associate Editor F.J. Pazzaglia for their detailed and highly instructive reviews.

1GSA Data Repository item 2016306, Appendixes, is available at or by request to
Science Editor: David I. Schofield
Associate Editor: Frank J. Pazzaglia