Source modelling of the M5–6 Emilia-Romagna, Italy, earthquakes (2012 May 20–29)

Abstract

On 2012 May 20 and 29, two damaging earthquakes with magnitudes M_w 6.1 and 5.9, respectively, struck the Emilia-Romagna region in the sedimentary Po Plain, Northern Italy, causing 26 fatalities, significant damage to historical buildings and substantial impact to the economy of the region. The earthquake sequence included four more aftershocks with M_w ≥ 5.0, all at shallow depths (about 7–9 km), with similar WNW–ESE striking reverse mechanism. The timeline of the sequence suggests significant static stress interaction between the largest events. We perform here a detailed source inversion, first adopting a point source approximation and considering pure double couple and full moment tensor source models. We compare different extended source inversion approaches for the two largest events, and find that the rupture occurred in both cases along a subhorizontal plane, dipping towards SSW. Directivity is well detected for the May 20 main shock, indicating that the rupture propagated unilaterally towards SE. Based on the focal mechanism solution, we further estimate the co-seismic static stress change induced by the May 20 event. By using the rate-and-state model and a Poissonian earthquake occurrence, we infer that the second largest event of May 29 was induced with a probability in the range 0.2–0.4. This suggests that the segment of fault was already prone to rupture. Finally, we estimate peak ground accelerations for the two main events as occurred separately or simultaneously. For the scenario involving hypothetical rupture areas of both main events, we estimate M_w = 6.3 and an increase of ground acceleration by 50 per cent. The approach we propose may help to quantify rapidly which regions are invested by a significant increase of the hazard, bearing the potential for large aftershocks or even a second main shock.

1 INTRODUCTION

On 2012 May 20 at 02:03 UTC an M_w 6.1 earthquake struck the eastern Po Plain (Northern Italy) in the area between Bologna, Ferrara and Modena. It has been followed by an extended aftershock sequence including several M > 5 events, culminating in an M_w 5.9 earthquake on 2012 May 29 at 07:00 UTC. The shallow hypocentral depths of those events (6–10 km) are responsible for the severe damages recorded in and around the towns of Finale Emilia, San Felice and Mirandola. Partial structural collapses, especially on historical buildings and industrial facilities, unfortunately caused 26 fatalities and the evacuation of nearly 15 000 people. The maximal intensities of the 2012 May 20 and the 2012 May 29 earthquakes were I₀ = VIII (MCS), where the intensity of the latter event is conditioned by the pre-damage of the former one. The maximum horizontal acceleration measured in the epicentral area is 0.26 g, whereas the maximum vertical acceleration is 0.31 g (RAN 2012).

The up to 8.5 km thick sedimentary Po Plain foredeep basin is relatively densely populated; industrial activity and land usage there is one of the most intensive in whole of Italy. The Po Plain sedimentary basin is thus very sensitive to ground shaking, resulting in a high vulnerability of structures. The natural seismic hazard in the Po Plain is, therefore, an important factor to understand seismic risk. Until now it has been considered much smaller in comparison to the Apennines and the Southern Alps fold and thrust belt (e.g. Albarello et al.2000), which both border the foredeep basin to the south and north, respectively. The magnitudes of the 2012 May earthquakes were expected at this location (Basili et al.2001 indicated a potential for M_w 6.3 shocks).

Historically and instrumentally recorded earthquakes are reported in the Parametric Catalogue of Italian Earthquakes (CPTI, Rovida et al.2011). In the Po Plain, a diffuse pattern of seismicity is indicated with only a few moderate earthquakes not exceeding M_e 5.6 (equivalent magnitude). The most important seismic event results to be the M_e 5.46 earthquake that occurred on 1570 November 17 near Ferrara, 30 km to the east of the epicentre of the 2012 May 20 main shock. In the western portion of the 2012 ongoing sequence, the CPTI reports damaging earthquakes for the years 1806 (M_e 5.19), 1810 (M_e 5.29) and 1832 (M_e 5.53). In recent times seismic activity was observed especially in the area near Reggio Emilia, 30 km SW of the May 20 main shock, with notable events on 1987 May 2 at 20:43 UTC (M_l 4.6) and on 1996 November 15 at 09:56 UTC (M_l 5.1, M_w 5.4 according to CMT; Selvaggi et al.2001). The lack of important historical earthquakes in the area where the 2012 May events struck resulted in a low hazard estimate, even if the concentration of important industrial facilities implies a high vulnerability and risk (Albarello et al.2000). The earthquakes raise the question whether the seismic hazard of the Po Plain was well estimated. For instance, Burrato et al. (2003) analysed river anomalies in the Po Plain and found geomorphical evidence for active blind thrust faulting at previously unknown locations, thereby indicating that the actual seismic hazard might be larger than what inferred purely from seismicity studies.

On the basis of extensive oil exploration studies, a good knowledge of the upper-crustal structure beneath the Po Plain has been achieved. The formation of the Po Plain is related to the collision between the European and the African plates and represents the foreland of the south verging Southern Alps and the north verging Apennines. The convergence has been accommodated by thrust faults along an almost continuous front, the so-called Northern Apennines Thrust front. According to Pieri & Groppi (1981), Bigi et al. (1989), Scrocca et al. (2007), the outer thrust fronts of the Northern Apennines and the Southern Alps are buried beneath the 7–8.5 km thick syn-orogenic clastic Plio-Quaterny deposits of the Po-Plain (Doglioni 1993). In the southern sector of the Po Plain three main buried arcs have been recognized, that represent the boundary of the Northern Apenninic thrust front. The eastern one is called Ferrara-Romagna arc and is subdivided (from W to E) in the Ferrara-, Romagna- and Adriatic folds (Scrocca et al.2007).

The depth of historically and instrumentally recorded seismicity in the Po Plain is often under debate. A few hours after the M_w 5.4 Reggio Emilia 1996 earthquake, a local seismic network was deployed to locate about 800 aftershocks with magnitudes 1.6 ≤ M_l ≤ 4.3, sampling a depth range between 12 and 15 km (Selvaggi et al.2001). The authors conclude that basement faults were involved in the deformation process, indicating active compressional tectonics. On the other hand, the relatively large damage and intensities observed for the 2012 May events and the preliminary location of the aftershocks indicate that these earthquakes ruptured at shallower depth. Earthquakes occurring at shallow depths of a sedimentary basin are much more effective in producing strong ground motions. The Reggio Emilia 1996 earthquakes showed a thrust mechanism (Selvaggi et al.2001) and presented the first good evidence of active compressional tectonics in the Po Plain. This is of interest, since the Po Plain lies in a transition zone between the active extensional tectonic regime in the Central and Northern Apennines and the compressional tectonics of the Po Plain Adriatic front. Active blind thrust faults within the sediments have been inferred from seismic reflection profiles (Pieri & Groppi 1981; Pieri 1983), and were indicated by morphological studies, but were rarely evidenced by seismic source mechanisms. The 2012 May earthquakes show a reverse faulting mechanism, indicating a strike direction parallel to the Ferrara–Romagna thrust and fold belt (Fig. 1); the focal mechanism is different to the one of the Reggio Emilia event: strike 94°/217°, dip 54°/53°, rake 132°/47° according to the Global CMT catalogue.

Another striking observation of the 2012 May earthquake sequence is the occurrence of a relatively strong earthquake of M_w 5.9 nine days after the M_w 6.1 main shock. Source locations indicate that the two large events may have ruptured neighbouring patches on the same fault, or along neighbouring faults. Source locations and mechanisms suggest the activation of the Mirandola fault, or a subparallel fault system. This would indicate that the reactivated structure is possibly part of a larger fault system and that a potential of stronger earthquakes might exist. In this paper, we analyse the depth, moment and source mechanism of the main shock and the following largest aftershocks (above magnitude M_w 5.0) by means of regional waveform inversion. One important seismic event has not been included in our analysis: the recording of the M_w 5.1 earthquake, of 2012 May 29 at 11:00:25 UTC, is unfortunately mixed within the coda of the M_w 5.4 shock that occurred 5 min before.

We study the rupture of the main shock by different kinematic rupture inversions and discuss the effect of the depth and rupture directivity on the ground motion. Next, we estimate the probability that the May 29 M_w 5.9 event was induced by co-seismic static stress transfer, and discuss the relation to the spatial-temporal evolution of the aftershock statistics. In particular, we focus on the westernmost part of the sequence to discern whether the increased hazard of this region later on during the seismic sequence could have been foreseen. Finally, we investigate the possibility that the same fault segment could have had the potential of a larger event and the possible scenario that could have followed.

2 POINT SOURCE INVERSION

In the following, we carry out a point source inversion for the five largest events of the seismic sequence (Fig. 1). The inversions rely on the fit of 3-component unrotated (north, east, vertical) full waveform displacement traces and their amplitude spectra, from at least 32 broad-band seismic stations at regional distances (epicentral distances 200–700 km for M_w > 5.5, 200–500 km for M_w < 5.5). An overview of the station distribution is given in Fig. 1. Although the earthquakes occurred in the Po sedimentary basin, the wave travel path to most of the stations employed here samples ordinary crust. Therefore, we assume the IASP91 model (Kennett & Engdahl 1991). Green's functions are computed for the desired range of epicentral distances and for different source depths within the crust, and then stored in a Green's function database to quickly build synthetic seismograms for different source models. We use the Kiwi tools software (Heimann 2010) and a two-step inversion approach, as described in Cesca et al.2010, to derive the best double couple (DC) model. In the first step we fit amplitude spectra in the frequency range 0.01–0.05 Hz (for the two largest events) and 0.02–0.05 Hz (for the remaining 3 events) to obtain fault plane angles, seismic moment and centroid depth. In the second step, displacement traces are fitted to resolve the focal mechanism polarities and to obtain the centroid location and centroid time. DC inversion results are summarized in Fig. 1. In all cases, we find similar focal mechanisms characterized by thrust faulting (rake angles ranging 84°–103°) of common orientation (strike angles 95°–127°) and dipping (dip angles 25°–31° for the less inclined plane). In particular, the two strongest events have almost the same focal mechanism (strike 103º, dip 25°–27°, rake 87°–92°), whereas the other three aftershocks have slightly different orientations (strike 95°–103°). These results show a general agreement with the regional time domain moment tensor inversion results performed, for example, by Pondrelli et al. (2012, QRCMT), Sarao’ & Peruzza (2012) and GFZ solutions. Results by Scognamiglio et al. (2012, TDMT) identify a slightly steeper dip angle, whereas the focal mechanism of the Global CMT solution includes an oblique component. An overview of available focal mechanisms is given in Table 1. Source depths of 7–8 km respectively are found for the two main events, this is slightly deeper than TDMT hypocentral depth estimations (5–6 km). Other solutions (depths 10–15 km) have limited depth resolutions for shallow events. We further apply a bootstrap scheme, based on the repeated focal mechanism inversion using 200 stations configurations (in each inversion, traces are randomly selected with replacement, so that some traces are not considered while others are differently weighted). The bootstrap approach allows us to infer source parameters uncertainties (Table 1). However, it only accounts for perturbation errors related to possible anomalies within the used data set, but it does not consider epistemic errors, which propagate through the model. Therefore, the bootstrap errors should be considered as a underestimation of the true uncertainties which are unknown (see also Custodio et al.2012).

For example, we obtain for the events of May 20 and 29 that the scalar moment are confined between 1.79 and 2.01×1018 Nm and 8.23–9.31×10¹⁷ Nm, respectively. These values are consistent with the preliminary results of seismic moment estimations published online, for example, by EMSC (http://www.emsc-csem.org/). This is consistent with our experience that moments computed with the Kiwi tools are in agreement with the seismic moments of regional events from EMSC. Magnitude estimations match well also with QRCMT, Sarao’ & Peruzza (2012) and Global CMT solutions (Table 1), whereas the TDMT estimated slightly smaller magnitudes (5.9–5.7). Centroid locations tend to be located towards NE, likely as a consequence of the asymmetric station distribution; therefore, we decided to fix the centroid location for the 2012 May 20 M_w 6.1 event, assuming the GFZ location. The centroid locations of the remaining four events are obtained through the centroid inversion in the time domain, and expressed as relative to the centroid location derived for the main shock. These locations align well along a WNW–ESE direction, with the exception of the 2012 May 29 10:55 UTC event, which is slightly to the south. Centroid depths indicate that the energy was mostly released at shallow depths, with the centroids at about 7–9 km depth, which can partially explain the significant damages reported in the epicentral region. Fig. 2 illustrates the focal mechanism, spectral and waveform comparison for the May 20 M_w 6.1 earthquake.

We also performed a full moment tensor inversion with the Kiwi tools (see Krieger & Heimann 2012; Cesca et al.2013, for details) following a similar multistep inversion approach. DC and full moment tensor are first retrieved by fitting amplitude spectra, whereas polarities are resolved in a second step through a time domain inversion. The full moment tensor solution is finally decomposed into DC, compensated linear vector dipole (CVLD) and isotropic (ISO) components. Results of the full moment tensor inversion are shown in Table 2. Full moment tensor solutions are characterized by minor non-DC components (6–18 per cent) for all events, which suggest that the seismic energy was released by a pure shear failure.

3 FINITE SOURCE INVERSION

The finite source inversion aims at exploring a range of possible extended source models and providing information about the fault plane orientation, the size of the rupture area and the rupture directivity. Our approach makes use of a simplified rupture model. The rupture area is originally circular, and additionally cut at the intersection with the free surface and the base of the crust. Both possible fault planes are tested, as well as different rupture size and nucleations points to account for rupture directivity. The choice of a simplified rupture model allows us to retrieve few important source parameters using a robust and stable inversion approach. The displacement data and the corresponding amplitude spectra used during the point source inversion are now used to derive the extended source parameters. The usage of a common data set and inversion algorithm for both point and extended source inversion allows us to fix the point source parameters (strike, dip, rake, moment and centroid location) during the kinematic inversion, so that only finite source parameters (plane orientation, rupture size, directivity) need to be inverted for. Although the point source inversion was derived in a frequency range between 0.01 and 0.05 Hz, associated with wavelength in between about 65 and 330 km, finite source parameters can only be derived including smaller wavelengths (higher frequencies). Finite source inversions results are only discussed for the two largest events. We apply and compare two independent approaches: in one case we perform an amplitude spectra inversion for several extended source models and find the model which fits better the observations (Cesca et al.2010), whereas in the second case we assume a spatial point source model and detect directivity effects based on the azimuthal pattern of apparent durations of the rupture process (Cesca et al.2011). The first method assumes the eikonal extended source model (Cesca et al.2010; Heimann 2010), which is a simplified finite source model. The source has a circular geometry bounded at the free surface and at the Moho discontinuity. We fix here the rupture velocity (this scales with the S-wave velocity of the model, and has here an average value of 2.24 km s^–1), and test different source geometries: two possible fault plane orientations, source radii in the range 1–15 km and 13 rupture propagation modes, including bilateral and unilateral rupture in different directions. So, though the geometric centre of each tested rupture model is fixed to the centroid location derived in the point source inversion, a number of different hypocentres are tested over the chosen rupture geometry. The L1 norm between observed and synthetic spectra is minimized and the solutions are compared based on their misfit. We repeated the inversion for different frequency ranges, to check the stability of the inversion for different bandpass and to find out a preferred frequency range: the lower frequency corner was fixed at 0.01 Hz, the highest varied between 0.07 and 0.5 Hz. When limiting to frequencies up 0.07 Hz, both finite source models oriented along the Northward and southward planes, perform similarly, as the spectra are best reproduced by a point source model. Differences become more relevant at frequencies up to 0.1 Hz or higher. However, the misfit improvement among models oriented along both possible fault planes does not increase further after 0.2 Hz, indicating that the observed high frequency radiation cannot be better reproduced by any of the two possible models. All inversions are consistent in terms of the identification of the rupture plane orientation and directivity. Based on this findings, we discuss here the results of the 0.01–0.1 Hz inversion. Results of the application of the finite source inversion to real data (Fig. 3 top) indicate for both events that the rupture occurred along the subhorizontal plane, which is dipping about 26° southwards towards the Apennines chain, and the preferred extended source models along these plane fit the amplitude spectra significantly better than models along the alternative steeper plane (L1 norm based misfits are 0.43 against 0.59 for the two planes for the May 20 event, and 0.41 versus 0.53 for the May 29 event). The fault plane probability, based on bootstrap results (200 random station configurations), indicated for both events that the subhorizontal fault geometries are preferred to the subvertical ones. The choice of the causative fault orientation is mostly based on this method, although confirmed from the second approach, which is discussed in the following. The finding of fault planes dipping to the south is confirmed by the known orientation of the buried thrust fault. The hypocentral distribution of the seismic sequence (Scognamiglio et al.2012) confirms the southward dipping plane. Locations are poorly aligned to discuss precisely the dipping of the fault, possibly suggesting a bended plane, steeper at the surface and subhorizontal at about 10–12 km depth. The finding of a shallow dip angle for the southward dipping plane is also supported by moment tensor results for smaller events in the sequence (Sarao’ & Peruzza 2012), where a majority of solutions show similar focal mechanisms. The fits of displacement waveforms and amplitude spectra after the finite source inversion are shown in Fig. 4, for the May 20 event. Note that the algorithm automatically remove few amplitude spectra, those which showed the poorest fits in the point source inversion, before the finite source inversion.

Directivity is well detected for the May 20 main shock, with the rupture propagating towards SE along a rupture length of about 20 km (radius 10 km, depth of the nucleation point 5.7 km); the reliability of this result is further supported by a bootstrap test, which confirms the rupture directivity for all tested station configurations. For the second event, the preferred rupture model is asymmetric, mostly propagating towards south (and down along the dip angle; depth of the nucleation point 6.5 km). The preferred length is also in this case 20 km. Rise times, indicating the duration of energy release for each point source used to describe the extended source model, are fixed to 3.3 s. Rupture times, describing here the time needed for the rupture front to propagate along the whole rupture area, are equal to 5.5 and 4.3 s for the main shocks of May 20 and 29, respectively. The duration of the source time functions result equal to about 9 and 8 s, and comparable to results published on the Geoazúr webpage (about 12 and 9 s), following the SCARDEC inversion approach (Vallée et al.2011). From these estimations and the parameters of the velocity model, we can provide rough values for the average slip, which result equal to 12 and 5 cm for these two events. If the shear modulus of softer sediments were considered, the average slip values would increase.

We employ an additional approach to detect rupture directivity and derive finite source parameters, as described by Cesca et al. (2011). This second method assumes a DC point source representation (focal mechanism, depth and moment are those derived in the DC source inversion) and the inversion is performed for the rise time only, which in this case describes the duration of the rupture process. The inversion is performed in the frequency domain (amplitude spectra inversion), including frequencies up to 0.1 Hz. Inversion tests including higher frequencies led to negligible changes: as for the kinematic source inversion, we tested the directivity inversion performance for different bandpasses with frequencies up to 0.5 Hz, with consistent directivity results. The rise time inversion is repeated for all stations, separately, providing several estimations for the source duration; note that all available components of the seismograms for a station are used together for this inversion. We will refer to ‘apparent’ source duration for the source duration derived by a single station. Stations located in the direction of rupture are expected to show preferably shorter apparent durations, whereas stations in the opposite direction will show preferably longer apparent durations. Cesca et al.2011 showed that, in favourable conditions, also asymmetric and bilateral ruptures may be discriminated with this method. The application of this technique to the two largest events (Fig. 3, bottom) confirms the main results obtained by the previous eikonal source inversion. For the event of 2012 May 20, 02:03 UTC we find a clear directivity, with rupture propagating towards SE. Because of its direction, this directivity effect can only be explained by a rupture occurring along the shallow dipping plane, as found with the previous kinematic inversion. The average rupture duration is about 8 s (comparable with the overall rupture duration derived in the kinematic inversion), and a strong azimuthally dependent directivity effect is determined. The May 29 main shock has an apparent duration in the range 5–8 s, with average duration of 6.5 s. However, the rupture directivity is less well constrained. In any case, it is relevant to note that this result is consistent with the rupture plane and rupture direction inferred by the kinematic inversion. According to Cesca et al. (2011), under the assumption of a simplified purely unilateral linear rupture model, the average rupture duration resembles the sum of rise and rupture time, whereas its variation depends on the wave propagation across the rupture length. Since our inversion is dominated by low-frequency surface waves, this can be likely faster than the rupture velocity. Assuming a rise time of about 3.0 s, as for the kinematic inversion, rupture times result equal to 5.0 and 3.5 s (which match quite well the values of 5.5 and 4.3 s resolved by the kinematic inversion). From the apparent duration curves we can then estimate a source size of about 14 and 8 km (a surface wave velocity of 4.0 km s^–1 was used), which is smaller than the one retrieved by the kinematic inversion (see Fig. 3 and Table 3 for a summary of finite source parameters). A reduction of the rupture size to 14 and 8 km, would lead here for parameters of the IASP91 crustal model to estimate the fault slip equal to 24 and 31 cm.

To check further the stability of our results we tested our inversion approaches (the eikonal source inversion and the quick directivity inversion) in the frequency range up to 0.1 Hz with a synthetic data set generated using a finite source model for the May 20 shock (southward dipping plane with focal mechanism, depth and scalar moment according to our point source solution, rupture propagating unilaterally towards SW along a circular area of 10 km radius, with the assumed rupture velocity and a rise time of 3 s) and a different velocity model (PREM, Dziewonski & Anderson 1981). Synthetic results of the eikonal source inversion correctly identify the target finite source parameters: rupture plane orientation, directivity pattern and size (7.5–10 km radii) of the rupture. Consistent results were also obtained using the alternative directivity approach.

4 STRESS TRANSFER AND EARTHQUAKE TRIGGERING

All results suggest that the two events ruptured within a few days on two neighbouring faults, or even different patches on the same blind low-angle thrust fault, extending WNW–ESE and dipping about 25° to the south. This raises the question whether there was a physical reason for the two earthquakes (or even more of the largest magnitude ones) for being distinct episodes, maybe related to the faults being physically separated, or to a complex fault geometry or rheology, or whether the May 20 and 29 earthquakes could have struck as a single event. In the first hypothesis, we estimate in the following the influence of the 2012 May 20 main shock on the occurrence of the May 29 event and the likelihood that the first event triggered or induced the rupture of the second event; we also analyse how the likelihood of more patches to rupture has been modified by the whole sequence. In the second hypothesis, we estimate what the magnitude of such a cumulative event would have been.

Static stress transfer within the crust is a physical process through which a large earthquake could induce/trigger other earthquakes on nearby favourably oriented seismogenic structures (Harris 1998; Steacy et al.2005). Several studies have looked at the role of static stress transfer with respect to main shock–aftershock sequences, sequences of large earthquakes along major fault systems (Parson et al.2000) or the occurrence of earthquake doublets (see for example Toda & Stein 2003). In particular the latter seems a possible process in play during the ongoing seismic sequence in Emilia. Coulomb stress studies link retrospectively the occurrence of earthquakes to a positive Coulomb stress change caused by a large earthquake or pressurized source in the crust. When areas exposed to a positive Coulomb stress change correspond with the location of large earthquakes occurring shortly thereafter, the earthquakes are often interpreted as triggered. On the other hand, seismicity entirely attributable in terms of loading to a previous earthquake (its aftershocks) or pressurized source (for example the seismicity produced during hydrofracturing) may be called induced (Dahm et al.2013). In this section, we apply a hybrid physical/statistical approach proposed by Passarelli et al. (2013) to the 5 ≤ M ≤ 6 earthquakes of the sequence, with the aim of establishing the likelihood they were induced or triggered by the May 20 event. The approach, here improved to consider an extended target fault geometry, involves calculating the Coulomb stress produced by the main shock on one or more fault planes, calculating the number of expected seismic events from the rate-state model (Dieterich 1994) and then obtaining the probability of occurrence of events overcoming a certain magnitude, with the assumption of a non-homogeneous Poisson earthquake distribution. Finally, the probability that the event was triggered or induced is calculated using the Bayes’ rule, by comparing the probability of occurrence of the target earthquake in two competing scenarios: a scenario in which the target event is entirely due to the regional tectonic loading and a second scenario where the stress induced by the main shock is held responsible for the target event (see Passarelli et al.2013, for more details).

4.1 Application of the probabilistic model to the Emilia 2012 sequence

The Coulomb stress change (ΔCFF = σ_s + μ σ_n, where σ_n and σ_s are the normal (negative when compressive) and shear stresses respectively, and μ is the friction coefficient, which we set equal to 0.8) is calculated for two different geometries of causative and receiver faults, according to the two extended source inversions computed in Section 3. We model the causative and the activated fault as single square edge dislocations (Okada 1992) of area equal to the area of the circular source inversions (Table 2), resulting in squares of side 17.7 and 17.7 km for the causative fault and 12.4 and 7.1 km for the target fault respectively. Strike, dip and rake are chosen according to Table 1. The Coulomb stress is averaged out over square surfaces oriented as, and as large as, the activated fault, and projected onto a horizontal plane at the depth of the May 29 earthquake nucleation point (Fig. 5 A1 and B1 for the two different solutions respectively), as derived in the kinematic inversion. The averaging operation improves the methodology of Passarelli et al. (2013) by allowing us to account for the size of the activated fault and hence the magnitude of the event. The rigidity modulus is taken as 20 GPa, which is compatible with the shear modulus from the velocity model.

We tested the influence of pore pressure in the Coulomb stress calculation but found the changes with respect to the dry formulation to be insignificant. This is mainly due to the minor contribution of normal stresses on the fault plane of the second shock.

Additional parameters needed in the procedure illustrated by Passarelli et al. (2013) are the background seismicity rate, the secular shear stressing rate and a constitutive parameter Aσ for the studied region (see Cocco et al.2010, for a review).

As presented in the introduction, the instrumental catalogue (Iside 2012) lacks events of M > 5 in the area, so that the background rate for M > 5 events is only available from the Italian historical earthquake catalogue. The seismic activity recorded instrumentally in the past 10 years of the area included less than 100 events with maximum magnitude of 4.7 (Iside 2012). In a wider area around the two main shocks (44.75N–44.25N and 11.00E–12.00E, see Fig. 1), the historical catalogue (Rovida et al.2011) includes 4 earthquakes with M > 5, all in the easternmost part of the geological structure, responsible for the ongoing activity. The resulting background rate for the historical quakes with M > 5 is 2.33×10⁻⁷ yr⁻¹ km⁻².

The secular shear-stressing rate can be estimated from the deformation data available for the region. From GPS measurements, the interested area is subject to NNE–SSW shrinkage at a rate of about 3 mm yr^–1 (Bennett et al.2012). We calculate the deformation from the displacement field on a transect of 30 km across the fault resulting in a stressing rate of 2×10⁻⁴ MPa yr^–1 for a rigidity of 20 GPa.

As for Aσ, according to the rate-state earthquake nucleation theory, an independent estimate can be obtained from the stressing rate and the decay time of the aftershock activity to the seismic background. A preliminary estimation of the aftershock decay rate with the available data suggests a decay time of about 10 years yielding to a lower estimate for Aσ of 2×10⁻³ MPa. However, the sequence is still ongoing so that the estimate is not robust yet, and this value seems small if compared with the commonly used value in literature, about 10⁻² MPa. However, the latter value yields a probably unrealistic predicted Omori decay-time of more than 100 years, given the low stressing rate of the region. Given that Aσ is the largest source of uncertainty in the probabilistic methodology as discussed in Passarelli et al. (2013), we decided to calculate the probability of triggering for values of Aσ of 2×10⁻³ MPa and 5×10⁻³ MPa. The latter value of Aσ is chosen since it yields an aftershock decay time of about 50 years, that we see as an upper estimate of such an Omori decay time.

Following Passarelli et al. (2013), we calculate the probability of occurrence of the May 29 event within two competing scenarios. The first scenario is merely Tectonic (T) and the probability of this scenario given the observation of the event, p(T|E), is calculated from the expected number of earthquakes with M > 5 given by the tectonic load only. As for the second competing scenario, which we call first-shock (FS) scenario, its probability conditional to the occurrence of the event, P(FS|E), is inferred by considering the expected number of earthquakes due to the stress induced by the May 20 event on the May 29 event. The uncertainty in the earthquake relative locations is considered by using a spatial Gaussian probability density function with mean equal to the second earthquake epicentre and truncated at 3σ = 5 km. The probability P(FS|E) is obtained by integrating over the probability density function (Passarelli et al.2013) and can be interpreted as the probability of the second shock being induced.

The operative definition of these two scenarios is possible since we can always separate the portion of events due to the tectonic loading and to the stress induced by the first shock within the rate-state framework. Finally, since we do not have any a priori information about the likelihood of each scenario, we set the a priori probabilities for each scenario, that is, p(T) and p(FS), equal to 0.5 as discussed in Passarelli et al. (2013).

4.2 Modelling results

We calculate the daily value of P(FS|E; see Fig. 5 A2 and B2 for the two extended source solutions). For both sets of fault parameters, we find P(FS|E) equal to 0.41 on May 29, see Fig. 5. On the same day, when we increase Aσ to 5×10⁻³ MPa, we find P(FS|E) = 0.24 and P(FS|E) = 0.19 for the first and second set of fault plane parameters, respectively. The values of probability below 50 per cent corroborate the hypothesis that the relatively large May 29 event cannot be fully explained by considering only the stress perturbation of the May 20 earthquake. We conclude that the tectonic stress accumulated on that system of faults must have played a major role in the May 29 earthquake occurrence, or, in other words, that the May 29 event was probably mature on that fault and was not an aftershock of the May 20 event. On the contrary, for one of the largest aftershocks (May, 20 03:02 in the introduction section), occurred on May 20 and only 5 km east than the May 29 event, we find a much higher probability of being induced, that is, 0.99 and 0.98 for the two causative fault parameters, respectively.

In the light of these findings, it can be instructive to evaluate the average shear stress induced by the first event on the fault plane of the second one, to compare it with an estimation of the stress drop associated with the May 29 event. We calculate the static stress drop using the formula Δσ = c μ Δu/L (Lee et al. 2003), where c is a constant set to be equal to one if the characteristic length of the fault, L, is well constrained, μ is the rigidity and Δu is the average slip. By considering the first set of fault parameters we infer a very small value of stress drop Δσ = 6 × 10⁻² MPa. The average shear stress due to the May 20 event on the focal plane of the May 29 one is τ = 8.7×10⁻³ MPa yielding a ratio R₁= τ/Δσ = 0.15. Contrarily, we find a larger stress drop estimate for the second set of fault parameters: Δσ = 9 ×10⁻¹ MPa and τ = 1.5×10⁻² MPa giving a ratio of R₂= τ/Δσ = 0.02. The discrepancy between R₁ and R₂ can be explained simply by the fact that the first set of fault parameters have a larger fault dimension and a smaller average slip, whereas for the second set the fault dimension is smaller and the average slip is higher, ought to the fact that only the seismic moment is well constrained by the waveform inversion. The small values of R₁ and R₂ confirm that the stress induced by the May 20 main shock explains only a small fraction of the total stress released by the May 29 earthquake. Therefore, as already indicated by the estimate of the induced probability, we infer that most of the stress released by the May 29 earthquake was already loaded on the structure. This implies that aftershock prediction models would have not been able to foresee the occurrence of a M = 5.9 shock 9 days after the main shock. However, an extended Coulomb stress model as the one developed in Passarelli et al. (2013) would have highlighted increased probabilities in that area.

5 DISCUSSION

One of the characteristics, observed often during seismic sequences in Italy, is that after the main shock, further earthquakes occur on adjacent faults, showing magnitudes similar to the main shock. The time span between the first and the following shakes can vary significantly, between days (L’Aquila, 2009; Pino & Di Luccio 2009), hours/weeks (Umbria–Marche, 1997; Amato et al.1983) or a few seconds (Irpinia 1980; Deschamps & King 1984). Although in the first two examples the faults are activated at completely different times, in the latter case, the immediate activation of successive faults leads to a sequence of subevents, whose superposition sums up to one large event, in terms of magnitude and damage. Critical for the damage is the time span that elapses between the subevents.

The earthquakes of the Emilia seismic sequence with M_w > 5 occurred on two distinct days and struck two different parts of the Mirandola fault or two subparallel faults. On 2012 May 20 the seismic events of 02:03 (M_w 6.1), 03:02 (M_w 5.2) and 13:18 (M_w 5.1) affected the 20–25 km eastern sector of the fault, whereas the earthquakes of 2012 May 29 occurred on the opposite side. The seismic event of 2012 May 20 was preceded 2 hours before by one M 4.1 earthquake and a few minor events, so that there was no indication for the imminent occurrence of a larger earthquake. This may have been different for the M_w 5.9 event of 2012 May 29. As reported in ‘quasi-real-time’ on the ISIDE bulletin (http://iside.rm.ingv.it), the epicentres of the aftershocks, recorded immediately after M_w 6.1 event of May 20, affected not only the 20–25 km of the eastern sector between Finale Emilia and Ferrara, which already ruptured by the M_w 6.1 event to a major portion, but also an area of similar extent towards west.

One of the main purposes of our contribution is to design a procedure that can be applied quickly, for example, within a few hours, during a seismic sequence, to evaluate where the seismic hazard may have increased. In particular, we have addressed here the question whether the likelihood of a large second event was significant in the Emilia sequence, after the occurrence of the May 20 event, and whether the May 29 event could have been suggested in advance, at least in terms of likelihood of occurrence of a second shock of magnitude similar to the main shock, exclusively on the basis of data of immediate availability (Web-data-server).

Our procedure is articulated in several steps. For the major events we have determined the focal parameters, such as fault plane angles, seismic moment and centroid depth, location and time. We have performed a finite source inversion, to determine the rupture duration, propagation direction (directivity) and average slip and identified on which of the DC-planes the rupture occurred. In the Emilia case, we find a low angle (≈30°) thrust fault of 103° deg strike ruptured the uppermost 10–15 km of the sediments over a length of about 20 km. Theoretical ground motion directivity have been estimated.

On the basis of the kinematic rupture model, we have estimated the earthquake-induced Coulomb stress changes. Finally, by means of the rate-state theory and estimates of the natural tectonic stress rate in the region, we have calculated the probability that earthquakes following the May 20 event of given size and location were purely main shock-induced or were likely driven by pre-existing tectonic stress. The value of probability we obtain for the May 29 event, in the range of 0.2–0.4, is on one hand not large enough for the earthquake to be considered as an aftershock of the May 20 event; on the other hand, from a forward perspective, to obtain such a large probability value for a relatively large magnitude event right after an earthquake, represents a significant warning of increased hazard. In other words, this may serve as an indication that a second main shock may be likely triggered.

The probabilistic approach indicates that the larger ‘aftershock’ of the May 20 main shock, which struck the western part of the segment, is induced with probability of 0.98–0.99. The effect we describe and quantify in a probabilistic way in this paper may be seen in a heuristic, classical way by analysing the magnitude-frequency distribution of early aftershocks. However, such an approach bears inherently the problem of insufficient number of events and therefore an insignificant statistical base.

The evolution of the sequence raises the question whether the fault segments of the M_w 6.1 and the M_w 5.9 event could have ruptured nearly simultaneously, for example, by dynamic triggering, so that the total size of the event would have been larger. Geomorphologic studies allowed the identification of blind active faults associated with subdued topographic expression (Valensise & Pantosti 2001; Burrato et al.2003). One of the recognized seismogenic sources buried below the Po Plain sediments is the so-called Mirandola anticline, which belongs to the Ferrara arc and is considered active during late Pliocene-early Pleistocene times. However, tectonic relative uplift rate of about 0.16 mm y^–1 can still be recognized during the last 125 ky. Horizontal shortening faster than 1 mm y^–1 should be expected in agreement with available GPS data (Scrocca et al.2007). As outlined by Burrato et al. (2003), Boccaletti et al. (2004) and DISS Working Group (2010), the Mirandola structure is considered an active seismogenic fault.

The Mirandola source was identified in the Database of Individual Seismogenic Sources (DISS-3.1.1., Basili et al.2008) as ITIS107 with the following parameters: Location: Lat. 44.8396, Lon. 11.1351, depth: 3.9–7.6 km; fault dimension: length 8.7 km, width 5.8 km; Fault orientation: strike 113°, dip 40°, rake 90°; Magnitude M_w 5.9; Recurrence 900–1800 y.

Previously, in the framework of the Slow Active Faults in Europe(SAFE)-project (Basili et al.2001), the Mirandola source (ITA-107) was already identified as a silent source that has not released any earthquake in historical time, but shows distinct slip rates. The parameters given in the SAFE report for the fault dimension and potential magnitude, are indicated as 18 × 14 km fault dimension in a depth range of 6–12 km, capable to generate a M_w 6.3 earthquake.

How would ground motion have changed if a complex rupture scenario is considered? The results of the point and extended source inversions can be used to model the relative peak ground accelerations (PGA) and its variations in the epicentral region. We consider the preferred source models from the previous inversions and embed them within a 1-D velocity model, derived specifically for the epicentral region by superposing a local crustal model (see electronic supplement) based on CRUST 2.0 (Bassin et al.2000), above the mantle model used for the previous inversion (IASP91, Kennett & Engdahl 1991). The Moho depth has been then fixed at 31 km, a compromise among the outcomes (27–35 km) of different previous studies (Cassinis et al.2003; Di Stefano et al.2009; Piana Agostinetti & Amato 2009); attenuation coefficients have been adapted after Morasca et al. (2010). The finite source models accounts for: (a) a starting circular source model, bounded at depths between the sediments and the Moho, (b) rupture propagation, according to extended source inversion results, (c) rise time fixed to two samples to reproduce high-frequency radiation components and rupture velocity accordingly to the kinematic source inversion approach. We additionally consider the scenario involving a larger, hypothetic rupture spanning the rupture areas of both main events, assuming that stress had been released in a single event. Based on previous results, we choose a rupture size with a maximal length of 36 km and a maximal width of 24 km. The focal mechanism as well as the nucleation point are assumed as for the May 20, event, so that the rupture propagation resembles now the long term evolution of the seismicity towards WNW along the 10 days period including both events, rather than directivities inferred for the largest events. Assuming a similar average slip as for the main event (12 cm), we derive for this rupture scenario a scalar moment of 5.20×10¹⁸ Nm, and thus a magnitude M_w 6.3. This is at the upper bound, but still consistent with the maximum magnitude estimate by Basili et al. (2001).

Maximal peak-to-PGAs (Fig. 6) are computed for a dense grid at the surface, in the epicentral region. Given the sampling of the synthetic accelerograms, they are limited to frequencies up to 2 Hz. PGAs are shown both for the May 20 M_w 6.1 event, the May 29 M_w 5.9 event and the composed rupture scenario (values are normalized with respect to the maximal amplitude of the composed scenario). We limit the discussion to the most general trends of the derived maximal acceleration maps, whose accuracy is limited both, by the assumption of averaged 1-D velocity profiles, and by the choice of a simplified rupture model. The maximal estimated acceleration for the three models, and are mostly localized north of the earthquake centroids, within elongated regions extending WNW–ESE parallel to the strike angles. In the frequency band up to 2 Hz, for which we could compute synthetic maximal accelerations, we observe that the peak acceleration amplitude for the composed scenario with a larger rupture area increases by a factor 1.5, compared to the estimation for the single rupture of the May 20 event. In addition, the comparison of acceleration maps for the modelled rupture and the composed scenario shows that seismic hazard extends to a larger region, prominently towards WNW and ESE.

6 CONCLUSIONS

With the constraint of regional seismic data, we can draw the following conclusions on the 2012 May/June earthquake sequence in the Po Plain, Emilia:

1. We estimate M 6.1 and M 5.9 for the May 20 and May 29, respectively. The strike of the fault planes is inferred to be 103° for both events. Also their dip was probably very similar, as we infer 25° and 27°, respectively. These results are consistent with focal mechanisms solution provided by different institutions (see EMSC-CSEM webpage)

2. We identify the southwards dipping planes as the rupture plane of the May 20 and May 29 events. The dip angles are in agreement with most published moment tensor solutions for the main shocks, and with the average focal mechanisms for major aftershocks. According to Scognamiglio et al. (2012), hypocentres are dipping towards south. Locations are not well aligned but could roughly indicate a bending of the fault, which is steeper at the surface and subhorizontal at the depth of about 10 km.

3. For the May 20 event we estimate a fault length of about 14–20 km and an average slip of about 12–24 cm. For the May 29 fault length and slip have larger uncertainties, we estimate a length of 8–20 km and a slip of 5–31 cm. Ruptures took place in the upper 12 km, in agreement with the depths of located aftershocks, which is mostly confined at shallow depths (Scognamiglio et al.2012).

4. Both events show a clear directivity and strongly asymmetric to unilateral propagation during rupture. The May 20 event nucleated on a NW patch of the fault and the rupture propagated towards SE (downdip and to the east). Our result is not in agreement with results of Piccinini et al. (2012), based on empirical Green's Function and array analysis, who identified two rupture pulses and interpreted them in terms of a rupture directivity towards WSW; the two pulses could also indicate the activation of neighbouring parallel faults and not reflect the main direction of rupture propagating along these rupture planes. The May 29 event nucleated at shallow depth and propagated downdip (or to the south) asymmetrically or even unilaterally. In both cases, we adopted a simplified 1-D model-to-model wave propagation, neglecting lateral structural heterogeneities, path and site effects which may affect finite source inversion results.

5. The May 20 event definitely increased the probability of occurrence of the May 29 event if we assume the fault geometry from the extended source inversions. The May 29 event cannot however be considered fully as an aftershock of the May 20 event because we find a probability in the range 0.2–0.4 that it was induced, and the shear stress induced on the fault plane was not sufficient to cause an event of such magnitude, so that the event must have been more or less mature on the fault. However, it is likely that the action of the May 20 event was to anticipate the occurrence of an event on the Mirandola fault. Also, a probablity in the range of 0.2–0.4 is quite large in terms of hazard estimate during an emergency.

6. We find a very low value for the Aσ rate-and-state parameter, coupled with low stressing rate and very low background. This implies a relatively long time needed to return to background seismicity values. This is in broad agreement with historical accounts for the Ferrara earthquake sequences in the middle age which lasted for years in intensities above the human perception.

7. If the two largest events of the sequence had struck together as one single earthquake, its magnitude would have been about M 6.3.

8. However, this scenario might have not been possible from a physical point of view.

9. The amplitude of the PGA for the composed scenario with a larger rupture area increases by a factor 1.5, compared to the estimation for the single rupture of the May 20 event. The comparison of acceleration maps for the modelled rupture and the composed scenario shows that strong shaking could have affected a larger region,, prominently towards WNW and ESE, if the composed rupture area would have failed at once.

10. The string of several events with magnitudes similar to the main shock highlights a complex fault structure, with neighbouring fault patches slipping in a complex pattern.

Since our study is based on regional seismic data, we do not have the resolution to discriminate the detailed fault geometry or say whether the faults which originated the largest magnitude events were physically separated or not; more information on this issue will come from local seismic data or from deformation studies.

Questions that remain open on this event, or more general questions inspired by our study are:

1. Was there a reason why the two largest events of the sequence struck 9 days apart, or are the faults separated? Why didn't they slip at once?

2. Is there a reason why the directivity of the May 20 event points toward southeast? Was there more stress accumulated on the western side of the fault after the Ferrara earthquakes in the past centuries?

3. Is there a relation in general between the dynamic evolution on the fault plane with the resulting heterogeneous residual strain or stress, and the spatial occurrence of aftershocks or additional main shocks on neighbouring fault patches?

The facilities of GEOFON were used to access seismic waveforms and station metadata. Stations from the following networks have been used: Basilicata University, GEOSCOPE, GEOFON, German Regional Seismic Network, Regional Seismic Network of Northwestern Italy, Hungarian Seismological Network, MEDNET and Austrian Seismic Network. SC received funding from the BMBF project MINE (Programme GEOTECHNOLOGIEN, Grant of project BMBF03G0737); LP, FM and ER received funding from the ERC StG project CCMP-POMPEI. We are grateful to Francesco Mele, Licia Faenza and Enrico Serpelloni for useful advice on the data and the kinematics of the area.

REFERENCES