Abstract
The Emilia-Romagna seismic sequence in May 2012 was characterized by two mainshocks which were close in time and space. Several authors already modeled the geodetic data in terms of the mechanical interaction of the events in the seismic sequence. Liquefaction has been extensively observed, suggesting an important role of fluids in the sequence. In this work, we focus on the poroelastic effects induced by the two mainshocks. In particular, the target of this work is to model the influence of fluids and pore-pressure changes on surface displacements and on the Coulomb failure function (CFF). The fluid flow and poroelastic modeling was performed in a 3D half-space whose elastic and hydraulic parameters are depth dependent, in accordance with the geology of the Emilia-Romagna subsoil. The model provides both the poroelastic displacements and the pore-pressure changes induced coseismically by the two mainshocks at subsequent periods and their evolution over time. Modeling results are then compared with postseismic InSAR and GPS displacement time series: the InSAR data consist of two SBAS series presented in previous works, while the GPS signal was detected adopting a variational Bayesian independent component analysis (vbICA) method. Thanks to the vbICA, we are able to separate the contribution of afterslip and poroelasticity on the horizontal surface displacements recorded by the GPS stations. The poroelastic GPS component is then compared to the modeled displacements and shown to be mainly due to drainage of the shallowest layers. Our results offer an estimation of the poroelastic effect magnitude that is small but not negligible and mostly confined in the near field of the two mainshocks. We also show that accounting for a 3D fault representation with a nonuniform slip distribution and the elastic-hydraulic layering of the half-space has an important role in the simulation results.
1. Introduction
During the seismic sequence started in May 2012 in the Emilia-Romagna region (northern Italy), a first mainshock ( 6.1, depth 6.3 km) occurred in May 20 at 4:04 AM local time and was followed 9 days later (May 29 at 1:04 AM, local time) by a second mainshock ( 6.0, depth 8 km) about 15 km southwest from the epicenter of the former one (Figure 1). The seismic sequence evolved along an about 50 km long, east-west-oriented thrust belt, including 8 events with [1]. All the events occurred in the alluvial plain of the Po river, in the foredeep of the northern Apennines chain, where the subsoil is composed of Plio-Quaternary sediment layers placed just above the Emilia-Romagna blind thrust system [2–4]. This area is known to be characterized by active tectonic shortening, due to the northward movement of the Apennines at about 2 mm/yr with respect to the stable Adriatic microplate [5].
(a)
(b)
Both the aftershock distribution [1] and the InSAR deformation measurements [3, 6, 7] clearly suggest that the overall seismic sequence involved two different, partially overlapping, fault planes: the May 20 mainshock occurred on the Ferrara fault (the easternmost one), while the May 29 mainshock occurred on the Mirandola fault (the westernmost one). Both faults are mainly east-west oriented, dipping towards south. Several fault geometries, based on different geodetic data inversions, are proposed in the literature (e.g., [3, 6, 8]). According to the model proposed by Pezzo et al. [3] and by Nespoli et al. [7], most of the slip of the May 20 event occurred along the Ferrara fault plane, even though some minor, aseismic slip on the Mirandola fault before the May 27 event is required to explain a tail-shaped deformation detected by the SAR satellites in the western part of the interferogram (Figure 2b in [7]) within the time window May12–27, 2012. Nespoli et al. [7] also show both the ground displacements and the Coulomb failure function (CFF) variations generated by the two mainshocks using a layered half-space. They concluded that the second mainshock has been mechanically promoted by the first event, especially accounting for the CFF change induced by the aseismic slip on the Mirandola fault. Observing the seismicity distribution, they also concluded that the whole sequence was fostered by the two mainshocks.
During the seismic sequence, several fluid-related effects were observed, such as the increase of the level of some water wells [9, 10], episodes of soil liquefaction [11], and the presence of bubbling wells due to the methane seepage [12, 13]. These evidences suggest that the Emilia-Romagna seismic sequence affected the groundwater equilibrium. In fact, the stresses induced by the seismic events on saturated porous rocks can be transferred from the solid grain to the interstitial fluids [14, 15], altering the distribution of the pore pressure within the aquifer and causing water level changes, variations in the discharge rate of streams and springs, and liquefaction. Nespoli et al. [10] simulate the pore-pressure changes induced by the May 20 mainshock in the shallowest layers (≤0.5 km of depth) of the Emilia-Romagna Po Plain. They computed the coseismic stress changes induced by the May 20 mainshock in an elastic half-space, and then following the poroelastic formulation, they calculated the pore-pressure change in undrained conditions [16, 17]. They found that both the magnitude and the time evolution of the water level of three wells close to the May 20 epicenter can be suitably represented as an instantaneous variation of pore pressure and its subsequent recovery. Albano et al. [18] proposed a 2D finite element model in order to simulate the poroelastic response of the May 20 mainshock along a vertical section intercepting the Ferrara fault (Figure 1). In their work, the Ferrara fault was modeled during the postseismic phase as an impermeable contact interface with a uniform 0.08 friction coefficient that can relax poroelastic loading through afterslip. They found that poroelastic readjustment involving afterslip is important in order to reproduce the postseismic displacements measured by InSAR data and also the related CFF changes can be relevant in order to explain a larger percentage of aftershocks than explained when only coseismic CFF changes are considered.
In this work, we aim to study the evolution of pore pressure by using a 3D poroelastic model that extends to the seismogenic depth and that is capable of simulating both the May 20 and May 29 mainshocks. With this purpose, we use the PEGRN-PECMP poroelastic model [19]. Despite this code providing a simpler fluid representation (single-phase and incompressible), than the one for the TOUGH2 numeric simulator [20] used by Nespoli et al. [10], it is capable of representing both the pore-pressure changes induced by the two mainshocks at subsequent times (May 20 and May 29) and the poroelastic evolution of displacements and stresses in the postseismic time interval. The modeled poroelastic response is compared to the postseismic deformation signal measured by GPS and InSAR data. In particular, we focus on the deformation signals extracted from applying a multivariate statistical approach to the GPS ground displacement time series and on displacements extracted from the time series of published InSAR data (obtained with the small baseline subset (SBAS) technique) located in the maximum coseismic displacement area. Due to the sparse temporal resolution of the SBAS time series, especially between the two mainshocks, we do not include them in the multivariate statistical analysis.
1.1. GPS Data
The area affected by the 2012 seismic sequence is monitored by several GPS stations belonging to different continuous and survey-mode networks (see [8], for more details). Here, we have considered the analysis of ground displacements using only data from continuous GPS stations in the May 2012–May 2014 time interval, and we have considered stations within 50 km from the midpoint of the two mainshock epicenters. The displacement time series have been obtained analyzing the raw data following a three-step approach [21], using the GAMIT/GLOBK and QOCA software. Additionally, as in Serpelloni et al. [22], we removed a long-wavelength common mode error component obtained from a principal component analysis, improving the signal-to-noise ratio in the GPS time series. Figure 1 shows the location of the 14 stations considered in this study with respect to the fault position, while Figure s1 (Supplementary Materials) shows the GPS position time series.
The detrended, filtered displacement time series are the input of an independent component analysis (ICA) performed adopting a variational Bayesian approach (vbICA, [23]) and modified in order to take into account missing data [24, 25]. This method, applied to geodetic time series, performs a spatiotemporal separation of the geodetic data into a limited number of signals, subsequently interpreted as the actual physical sources that generated the observed displacements. This method has been successfully used to extract tectonic and nontectonic (e.g., mainly related with hydrology) deformation signals in the GPS time series in different settings [26–29]. In particular, each source has a specific spatial distribution () and follows a specific temporal evolution (). Both the spatial and temporal responses associated with every IC have unit norm, so that a weight coefficient (in mm) is necessary to rescale their contribution in explaining the original data. The displacement time series at a given station can be reconstructed by linearly summing up the contributions from all the ICs, each of which is obtained by multiplying the specific spatial distribution by the associated weight times the temporal evolution. Because of the larger uncertainties and higher noise level of the vertical GPS displacements in the Po Plain area, in this work we focus the analysis on the horizontal components.
The horizontal position time series are well reproduced using three ICs, as asked by the automatic relevance determination criterion adopted to select the number of components. The vbICA results for horizontal components are shown in Figure 2. Due to the low SNR and the low density of stations available, the ICs are noisy. Nonetheless, we can try to give an interpretation based on the signals that we expect to observe. The temporal function of the first independent component (V1, Figure 2(a)) shows a decay with time after the May 20 mainshock. The corresponding spatial distribution (Figure 2(c)) agrees with what we would expect from afterslip occurring on the two faults involved in the seismic sequence, as described by Cheloni et al. [6]. Fitting the temporal dependence of IC1 with an exponential curve, we estimate a decay time of days (Figure 2(b)). Moreover, from the IC1, we can identify a knee point of the curve at about 160 days after the May 20 mainshock that also corresponds to the flattening of the two InSAR displacement time series (Figure 1) and possibly to the end of afterslip. Looking at the remaining two ICs, we do not observe another similar exponentially decaying signal, suggesting that the eventual two afterslip processes on the two faults involved are not separated, likely due to the statistical nonindependence of the two afterslip signals. Even if the characterization of IC2 is not the goal of this study, the source of this signal can be hardly referred to the seismic sequence owing to its spatial pattern, which shows large contributions in the far field. As the second IC exhibits a large time scale, one hypothesis could be that it represents a multiannual signal. Multiannual signals in the GPS displacement time series, in fact, have been recently highlighted in the Apennines and in the Alps [29, 30] and associated with variations of hydrological processes acting at different spatial scales. The temporal function of the third IC shows a rapid increase in the early postseismic stage (first 10 days after the May 20 mainshock). Its spatial response is maximum at MO05, the closest station to the epicentral area. Here, we might expect localized rapid transients associated with fluid flow after a stress change imparted by a mainshock [31], and we interpret this third IC as the poroelastic response at the surface, as shown in the following.
(a)
(b)
(c)
(d)
(e)
1.2. The Poroelastic Model
The PEGRN-PECMP poroelastic model provides numerical solutions, based on the mirror-image technique, for computing earthquake-induced stress changes, displacements, and pore pressure in a multilayered water-saturated half-space. When we refer to pore pressure , we mean the difference with respect to the hydrostatic value. The term “drained limit” is here used when, following a pressure perturbation , the system fully recovers hydrostatic conditions (). The code allows us to represent a horizontal-layered half-space with different elastic and hydraulic parameters. It consists of two separate routines: PEGRN, which calculates the Green functions needed for the strain and stress computations, following user-selected discretization of time and space, and PECMP, which computes the poroelastic variables (such as deformation, stress tensor, pore pressure, and Darcy flow) at receiver locations and at selected time steps. The code solves for the following poroelasticity equations depending on space and time , where and are the Lamé coefficients, is the displacement vector, is the effective stress coefficient, and is the body force per unit volume [19, 32]:
Equation 2 represents the mass conservation equation, where 1/Q is the Biot’s bulk compressibility (the amount of water that can be forced, under pressure, into a constant volume element) and is the Darcy conductivity and it can be expressed as s in terms of the hydraulic diffusivity, , and the specific storage [33]. The fluid volume injection rate, , is here assumed equal to 0. The coefficient (3) is expressed in terms of the drained and undrained Poisson’s ratio ( and ), inserting the Skempton’s coefficient , which represents the pore-pressure change per unit change in confining pressure under undrained conditions.
Once the Green functions are prepared, the PECMP routine allows us to represent an indefinite number of faults that can be discretized with variable slip patches. During the simulation, each fault can be activated, undergoing slip steps, at different time instants. In the postseismic stage, the fault is assumed as locked. This allows us to focus on the effect of pure poroelastic readjustment, excluding concurrent afterslip processes in the model.
In our simulations for both the May 20 and May 29 events, we use the coseismic slip distribution models, the faults geometries, and the elastic layering (Table 1 and Figure 3) proposed by Nespoli et al. [7]. The conversion from diffusivity to permeability is not trivial, but we can estimate the order of magnitude of the permeabilities as , where is the dynamic viscosity of water, is the water density, and is the gravity acceleration. Following the same procedure used by Ingebritsen and Manning [33] and Saar and Manga [34], we estimate permeabilities considering m−1 and and for pure water, varying with depth [34]. The four layers of the poroelastic half-space have different values representing a decrease of permeability over depth (Table 1), as proposed by Ingebritsen and Manning [33]. The chosen shallow permeability values are constrained by evidences obtained by drilling [13, 18, 35]. The resulting dependence of permeability on depth can be expressed , where is the permeability (m2) and is the depth (km). The selected values of diffusivities and permeabilities are consistent with the values proposed by Convertito et al. [36]. We do not have direct evidences of permeability values below few kilometers of depth, but it is reasonable to assume a lowering of permeability over depth, even if we cannot exclude that the fault permeability has locally increased following the earthquake. Furthermore, this profile of permeability is suitable for representing rapid hypocenter migration (see Table 1 in [33]). According to the InSAR data, the tail-shaped slip distribution on the Mirandola fault (Figure 3(a)), described by Pezzo et al. [3] and Nespoli et al. [7], might have occurred between May 20 and May 27. We assume that the time days corresponds to the May 20 mainshock, and we perform two limit case simulations: S1 and S2. In the simulation S1, the tail-shaped slip on the Mirandola fault occurs suddenly at days after the May 20 mainshock, while in the simulation S2, we assume that both the slip on the Ferrara fault and the tail-shaped slip on the Mirandola fault occurred simultaneously at (Table 2). In both S1 and S2, we assumed a Skempton coefficient , which means that we are considering fully saturated rocks in the entire domain. It is worth to notice that some rocks may have even if they are fully saturated [37], anyway considering is equivalent to assuming that the fluid pressure has the maximum possible response to the stress variations, so this assumption defines an upper limit for our model. We will refer to negative stresses as compressional, and we can estimate the undrained coseismic (step-like) pressure change induced by sudden or coseismic slip on a fault at as
We refer and as the undrained pressure change caused by slip on the Ferrara fault (assumed at days) and the tail-shaped slip on the Mirandola fault, respectively, before May 27 (assumed at days in S1 assumed at days in S2). represents the undrained pore-pressure change caused by slip occurred on the Mirandola fault in the time windows May 27–June 4 (assumed at days). The assumed coseismic pressure change history at a certain time can be expressed for S1 as and for S2 as , where is the Heaviside’s function.
2. Modeling Results
In this section, we focus on pressure, ground displacements, and CFF changes over time up to year. The modeled surface poroelastic displacements are then compared with the two InSAR displacement time series obtained at points 1 and 2 (Figure 1) that are located in the middle of the two fault projection where the surface displacement has relative maxima and with the IC3 component (Figure 2).
2.1. Pore Pressure
Figure 4(a) for S1 shows the distribution at the surface, that is, the pressure distribution at (coseismic stage). Red areas represent zones that undergo compression due to a volumetric contraction. In this figure, we also plot the location of fluid-related phenomena such as liquefaction and water level variations in wells. Watching their distribution, we can distinguish concentrations in zones with the greatest pore-pressure change (north of the May 20 epicenter) and just at the edge of the most negative pore-pressure change (east of the May 20 epicenter). Their occurrence is a clear evidence that the first mainshock leads to a perturbation of the preearthquake hydrostatic equilibrium that is mainly manifested in the near field. At the surface, the regions with the greatest positive pressure changes are due to coseismic compression near the projection of the edges of the slip distribution, as it is particularly evident in the case of uniform slip (Figure s2a), while negative pore-pressure changes are due to the coseismic extension created aside thrust slip which reaches the free surface above it (e.g., [38], Figure 8d), and it is reinforced by flexure associated with coseismic uplift detected in the same region (e.g., point 1 in Figure 1). Figure 4(b) shows along the AB vertical section, where the black arrows represent Darcy’s fluid flux in logarithmic scale at . Section AB, which intersects the May 20 mainshock fault (Ferrara), has two negative lobes that are NS oriented, while the positive lobes are mostly distributed in the deep part of the footwall and near the surface below the May 20 epicentral area. In the entire model domain, the pore-pressure variations range from −10 to 30 MPa. As we can expect, the fluid flow follows the pressure gradient directions; thus, they are generally going out from the pressurized zones and ingoing into areas where the pressure is lower than the hydrostatic value. We can distinguish a fluid flow pattern (orange arrows) that overall rises from the overpressurized lobe at depths, bending towards the shallower depressurized lobes with rotation from south to the north in the hanging wall and with a rotation from north to the south in the footwall. As we can see from Figure 4(c), in the CD vertical section, the pore-pressure changes induced by the May 20 event are small, and the fluid flow mainly occurs at shallow depths, where permeability is high. In Figure 5, the pore-pressure evolution is shown for the S2 case. At (Figure 5(a)), we obtain a larger, westward extension of the pore-pressure change at the surface () with respect to S1. At section AB, we can infer that both the tail-shaped slip at days (Figure 5(b) to be compared with Figure 4(b)) and the May 29 coseismic slip that occurred on the Mirandola fault at days (Figure 5(d) to be compared with Figure 5(b)) have negligible effects on the pore pressure. At 1 year after the mainshock, we can see that the pore-pressure changes in the first 3.3 km of depth are significantly decreased (Figure 5(f)) since the complete drainage of the first high permeable layer (above 0.8 km of depth, Table 1) and the almost-total drainage of the second permeable layer (from 0.8 to 3.3 km of depth) are reached. Figure 5(h) shows the pressure evolution at two locations within the high permeability layers, where after 1 year since the mainshock pore pressure is almost vanished in both S1 and S2. Below 3.3 km of depth, the permeability of the modeled rock is lower and the coseismic pore-pressure variation persists. At in the vertical CD section (Figure 5(c)), we can note much higher pore-pressure changes than in the S1 case (Figure 4(c)), as due to the contribution. Starting from 9 days later, the effect is evident in section CD (Figure 5(e)). Even in section CD at depth, the pore-pressure evolution is much slower than in the shallow layers, and at 1 year of simulation (Figure 5(g)), the pore-pressure change magnitude is still high. At from Figure 6(a), we can infer that along longitude, the pore-pressure gradients and consequently the fluid flow have the same magnitude of the ones obtained along latitude (Figure 5) even at seismogenic depth (Figure 6(b)). This means that our model leads to a nonnegligible fluid flow also in the along-strike direction.
2.2. Surface Displacements
We already showed that near the surface, pore-pressure differences between S1 and S2 become negligible after 10 days (Figure 5(h)). In this section, we discuss the surface displacements obtained in the S2 case. The computed May 20 coseismic surface displacements obtained in S2 are shown in Figure 7(a): the vertical deformation pattern (color) is the same as obtained by Nespoli et al. [7], with uplift up to about 20 cm in its maximum and about 16 cm at the point marked with the white cross (point 1). The westernmost positive deformation pattern (around the black cross, point 2) is due to the tail-shaped slip on the Mirandola fault taking place before the second mainshock (Figure 3(a)). The horizontal modeled coseismic displacement pattern (Figure 7(a), black arrows) generally converges toward the largest displacement area, except for the region close to the surface trace of the two faults and ESE along the strike direction, where the horizontal displacement is outgoing from the epicenter zone. After 1 year of simulation (Figure 7(b)), the vertical displacement at point 2 reaches a maximum of 12 cm mostly due to the coseismic effect of the second mainshock.
(a)
(b)
(c)
In order to isolate the influence of the poroelasticity, in Figure 7(c) we show both the vertical (color) and the horizontal (black arrows) modeled postseismic displacement fields obtained subtracting the coseismic offsets of the May 20 and May 29 mainshocks from the total displacements evaluated at the same time. An uplift increase of 1 cm in the area of the maximum coseismic uplift can be noted, while all the surrounding zones undergo less than 1 cm subsidence. In Figure 7(c), the four-petal flower-shaped pattern of the poroelastic postseismic horizontal displacement (black arrows) is very different from the coseismic one, trending normally to the strike out direction from the maximum deformation zone in near field and towards it in the along-strike direction with a magnitude lower than 1 mm at the GPS station locations, except for MO05 station where 5 mm in the SW direction is estimated. The horizontal pattern of the spatial response of IC3 (yellow arrows) is consistent with the modeled poroelastic horizontal deformation pattern, except for the southernmost GPS stations placed near the cities of Bologna and Modena, where IC3 magnitude is larger than the modeled deformation. This zone is well known in the Emilia-Romagna region because it is largely influenced by anthropogenic-driven subsidence, which is due to the massive withdrawal of groundwater for industrial uses [39, 40]. It is thus possible that such nontectonic signal is leaking in the IC3, corrupting its spatial and temporal responses at a certain extent. The vertical poroelastic postseismic responses at the GPS locations are smaller than 4 mm and are within the uncertainty level. Figure s2d shows the poroelastic displacement field obtained after the complete drainage of the half-space. The surface displacement field shown in Figure 7(c) is still far from the complete drained limit (Figure s3a) as it is mainly due to the drainage of the shallowest layers (Figures 5(f) and 5(g)). The poroelastic postseismic displacement due to uniform slip on a thrust fault is presented for comparison in Figure s3. Results of Figure s3b confirm the four-petal flower-shaped pattern of horizontal poroelastic postseismic displacement at 1 year after the mainshock. Even in this case, the surface displacement at 1 year is due to the drainage of the shallowest layers (Figure s2c) and it is far from the fully drained limit (Figure s3c).
Figure 8(a) shows the co- and postseismic displacement evolution modeled along the line of sight (LOS) at points 1 and 2 (see Figure 1). In particular, observed LOS displacements (orange and green dots) are InSAR observations taken, respectively, from Albano et al., [18] and Cheloni et al., [6]. The major modeled poroelastic response leads to a LOS increase of about 1 cm at point 1 (inset of Figure 8) and 0.5 cm at point 2, and it occurs within the first 10 days after the earthquake at both points following the relatively fast poroelastic relaxation of the shallow layers located at less than 3.3 km depth (green and purple lines in Figure 8(a)). After this time, the poroelastic readjustment leads to a slow and small constant decrease of the LOS distance, following the relaxation of the deepest part of the poroelastic medium. As we can infer from the modeled LOS displacement, at point 1 the deformation induced by the poroelastic rebound explains no more than 50% of the total deformation inferred by the interferograms (≈3-4 cm), while at point 2 the measured deformation may be totally imputed to poroelastic effects. In Figure 8(a), we also compare the LOS evolution at point 1 with the temporal evolution of IC3, which was rescaled to represent the point 1 modeled poroelastic LOS displacement in the first 50 days following the May 20 mainshock. We recall that the ICA was performed considering only the horizontal continuous GPS time series. However, comparing it to the InSAR measurements is useful to catch information about the timing of the poroelastic effects in the Emilia-Romagna seismic sequence, since they are supposed to be present with the same time scale in each displacement component. Moreover, at the GPS station MO05, the temporal dependence of the modeled magnitude of poroelastic postseismic horizontal displacement and the low-pass filtered GPS component time series reconstructed from IC3 is similar at least within the first 300 days (Figure 8(b)). In particular, the IC3 rising time (few days after the mainshock) is comparable to the modeled poroelastic signal.
(a)
(b)
2.3. CFF Changes
In order to estimate the fault interaction during the study seismic sequence, we compute the CFF changes on the Mirandola fault plane where the greatest number of seismic events occurred during the time interval between May 29, 2012, and June 28, 2012 ( days, the end of the seismic catalog of relocalized seismicity from [1]). The CFF variation was computed with the following equation: where and are the tangential and the normal stress components on the Mirandola fault plane and is the friction coefficient ( in this work).
For the May 20 coseismic computation on the Mirandola fault plane, is evaluated according to (4). The coseismic CFF changes (Figure s4) are greater, in magnitude, than the ones computed in the corresponding elastic-layered half-space using an apparent friction coefficient 0.4 (Figure S11A in [7]), especially in the deepest and the westernmost parts of the fault. At the May 29 hypocenter location, the coseismic CFF change is here estimated as 5 bar (0.13 bar) in the S2 (S1) case (Figure 9(a)), about 1.5 bar (0.1 bar) greater than in Nespoli et al. [7]. In both S1 and S2 cases, 9 days after the May 20 mainshock, the postseismic poroelastic CFF change at the May 29 mainshock hypocenter (Figure 9(a)) is positive but smaller than 2 kPa; at point 3 (Figure 9(b)), an arbitrarily selected query point within the preseismic slip on the Mirandola fault (Figure 3(a)) is positive but smaller than 5 kPa. In Figure 9(c), at days, we plot the postseismic poroelastic CFF changes together with seismicity [1]. Positive CFF changes larger than 10 kPa occurred within the rupture of the May 29 mainshock where the CFF is severely decreased at days during the coseismic stage (about 1.5 bar, e.g., in the hypocenter, see Figure 9(a)) so that the total (coseismic + postseismic) CFF change is still negative. In order to depict the areas where the poroelastic CFF changes could have contributed to promote the seismic sequence, Figure 9(d) shows at days, only the positive poroelastic occurred where total CFF changes are positive. Most of the seismicity (76%) occurred between days, and days is included in the red areas of Figure 9(d), where the postseismic poroelastic CFF changes are relatively small ranging from 1 to 10 kPa.
3. Discussion and Conclusions
We analyzed the GPS position time series after the May 20 mainshock with the vbICA method, finding three independent signals (Figure 2): the largest one, IC1, is related to afterslip and lasted for about 100–160-day time scale, close to the 180-day estimate of Cheloni et al. [6]; the second one, IC2, interests a wider region and apparently it is not related to the earthquake sequence, while the IC3 is strongly localized in near field and it is characterized by a 10-day rising time.
We estimated the poroelastic effects on the measured deformation field using the coseismic slip models proposed by Nespoli et al. [7] for the May 20 and May 29 mainshocks of the 2012 Emilia seismic sequence. We used the PEGRN-PECMP code [19] that allowed us to model the occurrence of both mainshocks in a layered poroelastic half-space. The elastic properties of the medium are the same as proposed by Nespoli et al. [7], while the permeability values of each layer of our model are computed following the profile suggested by Ingebritsen and Manning [33]. From vertical cross sections (Figure 5), we infer that most of the pore fluid dynamics in our simulation time (1 year) occur in the shallowest layers, that is, above 3.3 km of depth. At larger depths, where the permeability is lower, the fluid flow is not sufficient to drain the rock within the considered time interval.
At the surface, the modeled poroelastic postseismic horizontal displacement reaches the maximum absolute value of 1 cm about 10 days after the May 20 mainshock. This duration is related to the draining time of the shallowest layers of our model, since at this time scale the deeper layers have a second-order influence on surface deformation. Such a short time scale, lower than, for example, the 1-2-month period estimated for the June 2000 Iceland earthquake of 6.5 [31], can be due to the lower magnitude of the mainshocks, the different fault mechanism, and/or the here-assumed higher permeability of the first layers, confirmed by the fast recovery of the three-well water levels that occurred between 1 and 30 days [10]. In order to verify the robustness of our results, we tested two permeability profiles with permeability variations of ±1 order of magnitude with respect to the used values. We do not find substantial differences in the displacement field (Figure s5): after 1 year of simulation, it is mostly the same as the one obtained with the original permeability value (Table 1). This means that permeability changes of ±1 order of magnitude speed up/slow down the poroelastic displacement, without significantly affecting its value at 1 year.
From the poroelastic postseismic vertical deformation computed in scenario S1 (Figure s6), we can exclude that the tail-shaped deformation associated with the aseismic slip inferred on the Mirandola fault in the time between the two mainshocks was due to poroelastic rebound. Indeed, the poroelastic response in that time span leads to subsidence west from the May 20 mainshock, unlike what is observed by geodesy.
After 1 year, the modeled postseismic horizontal poroelastic displacement field (Figure 7(c)) is very different from the coseismic one (Figure 7(a)), resembling the horizontal spatial response associated with the IC3 derived from GPS position time series (Figure 2(e)), which is mainly represented by the SW movement of the MO05 GPS station. For this reason and considering its temporal evolution, we interpreted IC3 as the poroelastic contribution to surface displacement. Due to the short time scale of poroelastic rebound of the shallow layers, this interpretation is mainly supported by GPS data since InSAR time series can be temporally too sparse to sample the related transient displacement signal [41]. According to our model, the poroelastic effects are confined in the near field: they are significant (~5 mm) at the MO05 location and in the locations of the InSAR time series shown in Figure 1 (~10 mm, points 1 and 2 in Figure 1), while at the other GPS stations, they are of the order of 1 mm or smaller. In particular, at point 1, the modeled poroelastic LOS displacement explains about 50% of the total measured LOS displacement (Figure 8(a)), suggesting that here afterslip significantly contributed to the total postseismic displacement. With respect to the model of Albano et al. [18], we obtain a smaller and faster modeled poroelastic LOS displacement evolution at point 1. Possible reasons of this discrepancy are (1) the different assumed configurations used to model the ground displacements (here a 3D one is assumed, with respect to the 2D from previous work), (2) the different permeability profiles, and (3) the fact that we exclude the possibility of poroelastic induced afterslip. Given the different time scale of IC1 and IC3 and the reasonable agreement between IC3 and poroelastic rebound results for a locked fault (Figure 8(b)), we argue that a nonnegligible fraction of postseismic deformation at the surface occurred regardless of the poroelastic process and that even if the permeability layering here adopted is a coarse representation of the Emilia-Romagna subsoil, it is adequate to model both the order of magnitude and the time scale of the poroelastic effects. The low magnitude of surface displacements related to the poroelastic rebound following the Emilia seismic sequence is in agreement with estimates made by Freed [42] for the M6 Parkfield earthquake that occurred on 2004.
Positive postseismic poroelastic CFF changes, spatially correlated with the regions seismically activated, are estimated on the Mirandola fault plane starting from the May 20 mainshock to 30 days of simulations (Figure 9(c)). Nonetheless, the calculated CFF variations are much smaller (~10 kPa) than what was observed in other environments like, for example, in central Oklahoma, where wastewater injection has been imputed as the cause of the observed seismicity and ~0.07 MPa of CFF changes was inferred [43]. In the present model, such small postseismic positive CFF changes often occur in regions that are severely unloaded during the coseismic stage.
At the May 29 hypocenter, we find a 1.5 bar larger coseismic CFF change compared to the results obtained by Nespoli et al. [7] using an apparent friction coefficient model (, ). This result confirms and emphasizes the mechanical interaction of the two faults proposed by Nespoli et al. [7]. The large difference we obtained means that the assumption of the apparent friction model leads to a severe underestimation of the CFF change. At the May 29 hypocenter location, we also find, in the S2 case (Figure 9(a)), that the CFF is further increased by about 1 kPa in the first 9 days due to poroelastic rebound. High-pressurized fluids have been indicated as responsible for aseismic and seismic slip in laboratory experiments (e.g., [44]), midscale experiments (e.g., [45, 46]), and at fault scale level (e.g., [47–49]). In the present model, the estimated pressure variations are relative to hydrostatic conditions. In the depth range 6–9 km, where most of the aftershocks occurred, they can reach several MPa magnitude, as basically due to the coseismic pressure change history , which persists on a one-year time scale (Figure 5(g)). In the same depth range, Pezzo et al. [50] find / variations possibly related to pore-pressure changes. We cannot exclude that the slow fluid flow occurring at hypocentral depth could have played an active role in carrying on the Emilia-Romagna seismic sequence, but given the small postseismic pore-pressure and CFF changes here obtained, we think that the triggering of the Mirandola fault is more due to coseismic pore-pressure and CFF changes than to fluid migration.
This conclusion is supported by the analysis of the spatiotemporal distribution of the seismicity that occurred in the first 30 days (from [1], relocated) or in the first 50 days (INGV catalog, http://cnt.rm.ingv.it/) after the May 20 mainshock (Figure s7). It is possible that the quality of these last localizations does not have the necessary spatial and temporal resolution to detect fluid-driven seismicity; nevertheless, we do not see any clear evidence of fluid-driven hypocenter diffusion both along strike and dip, as observed in previous works (e.g., [49, 51]), but only the sudden westward migration due to the occurrence of the May 29 mainshock.
Finally, we also found that the pore-pressure changes and, consequently, the fluid flow as well have the same magnitude in both the along-strike and along-dip directions, in agreement with previous works that found an important fluid flow along the strike direction, especially if the slip distribution is heterogeneous [52, 53]. This means that, in order to simulate the evolution of pore-pressure induced by earthquakes, a 3D fluid flow modeling can be very important.
Data Availability
The whole data set is available upon request from [email protected]. The InSAR data are available from the the authors of [6] and [18].
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
The authors thank two anonymous referees for appreciating their work and for their useful comments.
Supplementary Materials
GPS time series, additional plots of simulations results, and temporal distribution of seismicity along strike and dip directions of the Ferrara and Mirandola faults. (Supplementary Materials)