Confirmation email

85.31.102.204Dear Author: Thank you for publishing with AGU. Your proofs are attached. Following the annotation guidelines and returning the proofs quickly will allow prompt online posting of the final version of your paper. Please note that revised proofs are sent only in the case of We have added several new features to the journal Websites and additional services that we want to bring to your attention. The online version of your paper is enhanced in several ways, and will now include altmetric s," wh ich continua lly trac ks a nd lin ks to mentio ns in news outlets, twitter, blogs, and other social media. A link to scholarly citations is also provided. Papers are also available in a regular PDF and as a ReadCube version, which includes the altmetrics information and dynamic links to references and figures. Please follow our social media sites! Several of our journals have social media sites; these and other AGU sites are indexed here. If!you!have!not!already,!we!encourage!you!to!create!an!ORCID!(Open!Researcher!and! Contributor!ID)!so!that!your!work!can!be!uniquely!tracked.!Your!ORCID!can!be!created!and! linked!through!your!GEMS!account.! Finally,!your!feedback!is!important!to!us.!!If!you!have!questions!or!comments!regarding! your!AGU!Publications experience,!including!information!on!production!and!proofs,!please contact!us!at!publications@agu.org.!!For questions specific to this paper, please contact the production editor at JGRBprod@wiley.com. We will also be contacting you soon after your article is published with an author survey. Please take a few minutes to respond to this online survey; your input is important in improving the overall editorial and production process. Thank!you!again!for!supporting!AGU. extensive corrections upon author request. " Sincerely, Brooks Hanson Director, Publications AGU bhanson@agu.org USING e-ANNOTATION TOOLS FOR ELECTRONIC PROOF CORRECTION Required software to e-Annotate PDFs: Adobe Acrobat Professional or Adobe Reader (version 8.0 or above). (Note that this document uses screenshots from Adobe Reader DC.) The latest version of Acrobat Reader can be downloaded for free at: http://get.adobe.com/reader/ Once you have Acrobat Reader open on your computer, click on the Comment tab (right-hand panel or under the Tools menu). This will open up a ribbon panel at the top of the document. Using a tool will place a comment in the right-hand panel. The tools you will use for annotating your proof are shown below: 1. Replace (Ins) Tool – for replacing text. Strikes a line through text and opens up a text box where replacement text can be entered. How to use it:  Highlight a word or sentence.  Click on .  Type the replacement text into the blue box that appears. 2. Strikethrough (Del) Tool – for deleting text. Strikes a red line through text that is to be deleted. How to use it:  Highlight a word or sentence.  Click on .. 3. Commenting Tool – for highlighting a section to be changed to bold or italic or for general comments. How to use it:   Click on .  Type any instructions regarding the text to be altered into the box that appears. 4. Insert Tool – for inserting missing text at specific points in the text. Use these 2 tools to highlight the text where a comment is then made. How to use it:  Click on .  Click at the point in the proof where the comment should be inserted.  Type the comment into the box that appears. Marks an insertion point in the text and opens up a text box where comments can be entered. Click and drag over the text you need to highlight for the comment you will add.  The text will be struck out in red.  Click on .  Click close to the text you just highlighted. USING e-ANNOTATION TOOLS FOR ELECTRONIC PROOF CORRECTION For further information on how to annotate proofs, click on the Help menu to reveal a list of further options: 5. Attach File Tool – for inserting large amounts of text or replacement figures. Inserts an icon linking to the attached file in the appropriate place in the text. How to use it:  Click on .  Click on the proof to where you'd like the attached file to be linked.  Select the file to be attached from your computer or network.  Select the colour and type of icon that will appear in the proof. Click OK. The attachment appears in the right-hand panel. 6. Add stamp Tool – for approving a proof if no corrections are required. Inserts a selected stamp onto an appropriate place in the proof. How to use it:  Click on .  Select the stamp you want to use. (The Approved stamp is usually available directly in the menu that appears. Others are shown under Dynamic, Sign Here, Standard Business).  Fill in any details and then click on the proof where you'd like the stamp to appear. (Where a proof is to be approved as it is, this would normally be on the first page). 7. Drawing Markups Tools – for drawing shapes, lines, and freeform annotations on proofs and commenting on these marks. Allows shapes, lines, and freeform annotations to be drawn on proofs and for comments to be made on these marks. How to use it:  Click on one of the shapes in the Drawing Markups section.  Click on the proof at the relevant point and draw the selected shape with the cursor.  To add a comment to the drawn shape, right-click on shape and select Open Pop-up Note.  Type any text in the red box that appears. Drawing tools available on comment ribbon THIS IS NOT AN INVOICE JGRB52152 PLEASE COMPLETE THE PUBLICATION FEE CONSENT FORM BELOW AND RETURN TO THE PRODUCTION EDITOR WITH YOUR PROOF CORRECTIONS Please return this completed form and direct any questions to the Wiley Journal Production Editor at JGRBprod@wiley.com. To order OnlineOpen, you must complete the OnlineOpen order form at: https://authorservices.wiley.com/bauthor/onlineopen_order.asp Authors who select OnlineOpen will be charged the standard OnlineOpen fee for your journal, but excess publication fees will still apply, if applicable. If your paper has generated excess publication fees, please complete and return the form below in addition to completing the OnlineOpen order form online (excess fees are billed separately). If you would like to choose OnlineOpen and you have not already submitted your order online, please do so now. YOUR ARTICLE DETAILS Journal: Journal of Geophysical Research: Solid Earth Article: Albano, M., S. Barba, G. Solaro, A. Pepe, C. Bignami, M. Moro, M. Saroli, and S. Stramondo (2017), Aftershocks, groundwater changes and postseismic ground displacements related to pore pressure gradients: Insights from the 2012 Emilia-Romagna earthquake, J. Geophys. Res. Solid Earth, 122, doi:10.1002/2017JB014009. OnlineOpen: No Words: 7,647 Tables: 1 Figures: 8 Total Publishing Units: 24 Journal Base Fee: $1,000 Excess Publishing Units: 0@$125 _______ $0 Publication Fee Total: USD $1,000 An invoice will be mailed to the address you have provided once your edited article publishes online in its final format. Please include on this publication fee form any information that must be included on the invoice. Publication Fees and Length Guidelines: http://publications.agu.org/author-resource-center/ Frequently Asked Billing Questions: http://onlinelibrary.wiley.com/journal/10.1002/(ISSN)2169-8996/homepage/billing_faqs.pdf Purchase Order Instructions: Wiley must be listed as the contractor on purchase orders to prevent delay in processing invoices and payments. Please provide the information requested below. Bill to: Name: ------------------------------------------------------------------------------------------------------------------------------------------------------------ Institution: ------------------------------------------------------------------------------------------------------------------------------------------------------- Address: ---------------------------------------------------------------------------------------------------------------------------------------------------------- Phone: ------------------------------------------------------------- Email: ------------------------------------------------------------------------------------- Signature: --------------------------------------------------------- Date: ------------------------------------------------------------------------------------- Aftershocks, groundwater changes and postseismic ground displacements related to pore pressure gradients: Insights from the 2012 Emilia-Romagna earthquake Matteo Albano Q2 1 , Salvatore Barba1 , Giuseppe Solaro2 , Antonio Pepe2 , Christian Bignami1 , Marco Moro1 , Michele Saroli3,1, and Salvatore Stramondo1 1 Istituto Nazionale di Geofisica e Vulcanologia, Rome, Italy, 2 Istituto per il Rilevamento Elettromagnetico dell'Ambiente, National Research Council (CNR), Naples, Italy, 3 Dipartimento di Ingegneria Civile e Meccanica, Università degli Studi di Cassino e del Lazio meridionale, Cassino, Italy Abstract During the 2012 Emilia-Romagna (Italy) seismic sequence, several time-dependent phenomena occurred, such as changes in the groundwater regime and chemistry, liquefaction, and postseismic ground displacements. Because time-dependent phenomena require time-dependent physical mechanisms, we interpreted such events as the result of the poroelastic response of the crust after the main shock. In our study, we performed a two-dimensional poroelastic numerical analysis calibrated with Cosmo-SkyMed interferometric data and measured piezometric levels in water wells. The simulation results are consistent with the observed postseismic ground displacement and water level changes. The simulations show that crustal volumetric changes induced by poroelastic relaxation and the afterslip along the main shock fault are both required to reproduce the amplitude (approximately 4 cm) and temporal evolution of the observed postseismic uplift. Poroelastic relaxation also affects the aftershock distribution. In fact, the aftershocks are correlated with the postseismic Coulomb stress evolution. In particular, a considerably higher fraction of aftershocks occurs when the evolving poroelastic Coulomb stress is positive. These findings highlight the need to perform calculations that adequately consider the time-dependent poroelastic effect when modeling postseismic scenarios, especially for forecasting the temporal and spatial evolution of stresses after a large earthquake. Failing to do so results in an overestimation of the afterslip and an inaccurate definition of stress and strain in the postseismic phase. 1. Introduction On 20 May 2012, an Mw 5.9 earthquake struck the Emilia-Romagna region (northern Italy) (Figure 1a), causing F1 severe damage and loss of life. The main shock and subsequent aftershocks activated two segments of the Ferrara-Romagna thrust system, which are the Ferrara thrusts (blue lines in Figure 1a) and the Mirandola thrusts (red lines in Figure 1a). Thousands of aftershocks were registered [Govoni et al., 2014], six of them with Mw > 5, with the largest aftershock (Mw 5.7) occurring on 29 May (Figure 1a). In the days following the main shock, the rupture propagated eastward and downdip on the Ferrara thrust system [Govoni et al., 2014] and westward on the adjacent Mirandola thrust, where the Mw 5.7 occurred 9 days later. All aftershocks with Mw > 5 occurred within 15 days of the 20 May event, and more than 70% of the earthquakes took place within 3 months of the main shock [ISIDe Working Group, 2010]. The moment tensor solutions for stronger events display a reverse focal mechanism [Scognamiglio et al., 2012]. After the main shock, several postseismic phenomena occurred. The monitored wells close to the epicentral area registered a sudden increase in water level that was followed by a decay that did not reach the preseismic level during the observation period [Marcaccio and Martinelli, 2012]. Some wells ejected a mixture of water and sand, and more than 700 liquefaction phenomena were observed [Alessio et al., 2013; Chini et al., 2015]. Furthermore, a periodic sampling of both liquid and gaseous phases in surficial aquifers established that deep-seated fluids was mobilized after the main shock and reached surface layers along faults and fractures, thus altering the geochemical composition of fluids [Italiano et al., 2012; Sciarra et al., 2012]. Postseismic deformation has been partly studied, at least for the Mw 5.7 event that occurred on 29 May. Spanning roughly 1 year after the main shock, differential synthetic aperture radar interferometry (DInSAR) ALBANO ET AL. THE 2012 EMILIA-ROMAGNA EARTHQUAKE 1 PUBLICATIONS Journal of Geophysical Research: Solid Earth RESEARCH ARTICLE 10.1002/2017JB014009 Key Points: • The postseismic phenomena that occurred after the 2012 Mw 5.9 northern Italy earthquake depended on the poroelastic response of the crust • Poroelastic relaxation is equally important as the afterslip in terms of the aftershock distribution • DInSAR ground deformation data and in-well water level changes were used to successfully calibrate the poroelastic numerical model Supporting Information: • Supporting Information S1 Correspondence to: M. Albano, matteo.albano@ingv.it Citation: Albano, M., S. Barba, G. Solaro, A. Pepe, C. Bignami, M. Moro, M. Saroli, and S. Stramondo (2017), Aftershocks, groundwater changes and postseismic ground displacements related to pore pressure gradients: Insights from the 2012 Emilia-Romagna earthquake, J. Geophys. Res. Solid Earth, 122, doi:10.1002/2017JB014009. Received 19 JAN 2017 Accepted 12 JUN 2017 Accepted article online 13 JUN 2017 ©2017. American Geophysical Union. All Rights Reserved. Journal Code Article ID Dispatch: 17.06.17 CE: AMI J G R B 5 2 1 5 2 No. of Pages: 17 ME: Revised proofs are sent only in the case of extensive corrections upon request 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 and GPS time series provided postseismic ground displacement data. Cheloni et al. highlighted a rapid postseismic uplift of approximately 2 cm following the Mw 5.7 aftershock on 29 May. These authors assumed that a transient seismic/aseismic slip located on the Mirandola and Ferrara structures caused the postseismic displacement and the aftershock spatial and temporal distributions. Different time-dependent physical mechanisms may account for one or more postseismic events, e.g., the afterslip [Bürgmann et al., 2002], viscoelastic relaxation [Bürgmann and Dresen, 2008], gravitational rebalancing [Albano et al., 2015], and poroelastic rebound [Jónsson et al., 2003]. More than one physical mechanism may develop after a large earthquake. A poroelastic crustal rebound can explain the observed water level changes, the postseismic ground displacements, and the aftershock distribution and decay. Indeed, the crustal poroelastic rebound has been invoked to explain the location and rate of aftershocks [Bosl and Nur, 2002; Shapiro et al., 2003; Antonioli et al., 2005], groundwater modifications [Whitehead et al., 1985; Sneed et al., 2005; Manga and Wang, 2007], and the amplitude and rate of postseismic displacements [Peltzeret al., 1998; Jónsson et al., 2003] after large earthquakes. Figure 1. (a) Tectonic sketch of the study area and earthquakes of Mw > 2 recorded from May 2012 to May 2014. Circles: earthquakes with magnitude Mw > 2. Stars: earthquakes with Mw > 5.0. The red and green stars indicate the events of 20 and 29 May 2012, respectively. The tectonic framework is depicted by the thrusts (solid lines) and back-thrusts (dashed lines) of Ferrara (blue), Mirandola (red), and Pede-Apennines (black) [Boccaletti et al., 2004; Chiarabba et al., 2014]. The purple polygon indicates the footprint of the COSMO-SkyMed (CSK) SAR frame. The dashed black segments indicate the selection limits for earthquakes that are associated with the Ferrara thrusts (see Figure 7). The dotted blue and red boxes identify the Ferrara and Mirandola fault planes obtained from the coseismic inversion modeling of Cheloni et al. . (b) Geological section along the A-A0 profile in Figure 1a [Carminati et al., 2010]. The geologic units are as follows: 1—Quaternary sediments; 2— marine terrigenous deposits (upper to middle Pliocene); 3—evaporitic and terrigenous deposits (upper Miocene to lower Pliocene); 4—marly calcareous and terrigenous levels (Oligocene-Miocene); 5—shallow to deep-water carbonates (Jurassic-Cretaceous) and Triassic evaporites; and 6—crystalline basement. Journal of Geophysical Research: Solid Earth 10.1002/2017JB014009 ALBANO ET AL. THE 2012 EMILIA-ROMAGNA EARTHQUAKE 2 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 The crust of the Earth behaves as a biphasic medium comprising of a solid skeleton and voids filled with fluids [Fyfe, 2012]. At seismogenic depths (8–15 km), such fluids are water or water based. Thus, we can assume that the crust is a fluid-saturated, biphasic medium. When an earthquake affects a fluid-saturated medium, the crustal deformation causes a sudden change in fluid pore pressure with respect to hydrostatic condition. These pore pressures are not in equilibrium, and the fluid then diffuses from overpressurized regions to under pressurized regions until hydrostatic conditions are achieved. The amplitude of the pore pressure change and the duration of the diffusion process obey the poroelastic diffusion laws and depend on the elastic and hydraulic properties of the medium and interstitial fluid [Biot, 1941; Rice and Cleary, 1976]. The slow change in pore pressure caused by fluid diffusion modifies the effective stresses inside the crust of the Earth, i.e., the stresses that are effective in causing displacements [Terzaghi, 1943], according to equation (1) [Wang, 2000]: σ 0 ij ¼ σij þ αpδij (1) where σ 0 ij and σij are the components of the effective and total stresses (positive when the stresses are tensile), p is the pore pressure, δij is the Kronecker delta, and α is the Biot -Willis coefficient. Therefore, further displacements occur at the ground surface and aftershocks occur at depth [Harris, 1998; Bosl and Nur, 2002; Jónsson et al., 2003; LaBonte et al., 2009; Hughes et al., 2010]. Our goal is to investigate the influence of poroelastic rebound on the widespread postseismic phenomena observed after the 20 May 2012 Emilia-Romagna earthquake. To achieve this objective, we calculate the stress perturbation and excess pore pressure produced by the fault dislocation associated with the 20 May event using linear poroelasticity. The poroelasticity theory allows us to couple the rock matrix deformation and fluid flow and simultaneously solve for the stress and pore pressure changes. The modeling results are evaluated with regard to the DInSAR ground displacement data obtained for the 20 May event, the measured coseismic and postseismic piezometric levels in the wells, and aftershock distribution and decay. 2. DInSAR Analysis and Results We studied the coseismic and postseismic displacements associated with the 20 May 2012 event using an unpublished data set consisting of 18 synthetic aperture radar (SAR) images acquired on the descending orbit (side-looking angle of approximately 32° ) of the COSMO-SkyMed (CSK) constellation satellites. The SAR data set, which spans the period from 19 May 2012 to 10 May 2013 (Table S1), allows us to measure the coseismic and 1 year of the postseismic displacement field. The first postseismic acquisition is dated 23 May 2012 (i.e., only 3 days after the main shock), capturing the earliest postseismic displacements. The footprint of the CSK data frames is shown in Figure 1a. Although the SAR images partially cover the main shock area, they represent the best available data set to study the 20 May earthquake and related postseismic phenomena. Beginning with the available SAR data, 55 interferograms were generated (Figure S1) and processed with a well-established small baseline subset (SBAS) algorithm [Berardino et al., 2002] to retrieve a spatially dense map showing the mean displacement velocity of the investigated area and the displacement time series for each coherent distributed target (DT) of the scene. The DInSAR interferograms were flattened using a 90 m Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) of the area [Rosen et al., 2001]. To mitigate the effects of noise, a complex multilook (ML) operation [Rosen et al., 2000] with 10 looks in both the azimuth and range directions and a noise-filtering operation [Goldstein and Werner, 1998] were independently conducted on each interferogram. The interferograms were unwrapped in space and time by applying the method proposed by Pepe and Lanari and subsequently inverted using the SBAS strategy. Residual orbital phase artifacts and topographic mistakes were also corrected [Manzo et al., 2012]. Finally, for each analyzed SAR pixel, the atmospheric phase screen (APS) was estimated and filtered out from the obtained displacement time series based on the assumption that APSs are possibly correlated in space but poorly correlated in time [Ferretti et al., 2001]. The SBAS processing returned 97,902 DT measurement points over an area of approximately 1270 km2 (density of approximately 80 DT/km2 ) (represented by blue dots in Figure S2). The data from DInSAR are herein presented in terms of mean displacement rate maps and time series. DInSAR can only measure the radar line-of-sight (LOS)-projected component of the surface displacements. We retrieved the ground Journal of Geophysical Research: Solid Earth 10.1002/2017JB014009 ALBANO ET AL. THE 2012 EMILIA-ROMAGNA EARTHQUAKE 3 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 displacement pattern by spatially interpolating the displacement values of each DT using the inverse distance to power method. The coseismic LOS displacement pattern on 23 May 2012 (Figure 2a) shows that the ground moved toward F2 the satellite (positive values) over an area encompassing the seismic activity in the hanging wall of the main shock, close to the town of Finale Emilia. A significant component of this motion involves uplift across the Inner and Middle Ferrara thrusts (IF and MF in Figure 1a). Although the obtained displacement lobe partially images the extent of the deformed area, it perfectly overlies the coseismic displacement pattern retrieved by the processing of two Radarsat-2 (RSAT-2) SAR images [Tizzani et al., 2013] (dashed black lines in Figure 2a). However, the CSK coseismic LOS displacement peak at point 1 (Figure 2a) (approximately 16 cm) is lower than the value taken from the coseismic displacement field of the Radarsat-2 data at the same point (17.5 cm) (Figure S3). Given that the acquisition geometry for CSK and RSAT-2 is nearly the same, such a difference is explained by considering that the processed Radarsat-2 pairs are dated 30 April 2012 to 17 June 2012, and the unwrapped displacement field encompasses approximately 28 days of postseismic uplift (from 20 May to 17 June). Figure 2. Ground displacement field encompassing the entire seismic sequence from the CSK SAR data along the descending track (view angle of approximately 32°). (a) Coseismic LOS displacement field associated with the 20 May event (19–23 May). The black dotted lines refer to the coseismic displacement pattern retrieved by the processing of two Radarsat-2 (RSAT-2) SAR images [Tizzani et al., 2013]. (b) Cumulative postseismic LOS displacement field (23 May 2012 to 10 May 2013). The white circles indicate the location of liquefaction phenomena. (c) Mean trend and associated standard deviation of the displacement time series retrieved from a group of DT close to the red, blue, green, and yellow triangles (indicated with the numbers 1, 2, 3, and 4, respectively) in Figure 2b. Journal of Geophysical Research: Solid Earth 10.1002/2017JB014009 ALBANO ET AL. THE 2012 EMILIA-ROMAGNA EARTHQUAKE 4 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 The cumulative postseismic LOS displacements after approximately 1 year following the main shock (10 May 2013) (Figure 2b) encompass the entire postseismic phase and show a complex deformation pattern. The southwestern edge of the CSK SAR data frame is affected by a diffused pattern of subsidence that increases linearly over time (Figure 2b and point 4 in Figure 2c). Such ground displacements are unrelated to the seismic sequence. Indeed, such an area is affected by anthropogenic subsidence (Figure S2) induced by groundwater withdrawal [Stramondo et al., 2007; Modoni et al., 2013]. Proceeding to the NE, the area between S. Agostino and Mirabello villages has undergone postseismic ground subsidence (negative values ieQTbn Figure 2b, representing movements away from the sensor). The LOS displacement time series (point 3 in Figure 2c) show an initial slight postseismic uplift, which is followed by a rapid subsidence that asymptotically decreases over approximately 3 months. Such diffused pattern of ground displacements is ascribable to the widespread liquefaction phenomena observed in the area a few days after the main shock [Alessio et al., 2013] (white circles in Figure 2b). Such phenomena (i.e., sand boils and water leaks from cracks) generally produce ground settlements not only in the affected area but also where the liquefaction phenomena are not recognized [Chini et al., 2015]. Finally, a main lobe of positive LOS displacements (representative of ground movements toward the sensor) with a displacement peak of approximately 3.5 cm is found close to the area that experienced the maximum coseismic uplift (Figure 2b). In this area, the uplift rate is quite fast given that more than 50% of the recorded postseismic displacement occurred in the first 2 to 3 months after the earthquake (points 1 and 2 in Figures 2b and 2c). The spatial correspondence between the aftershock epicenters and the coseismic ground displacement pattern suggests that the observed postseismic displacements have a tectonic origin [Cheloniet al., 2016]. 3. Methods The observed postseismic uplift has been interpreted as the result of the postseismic poroelastic rebound of the crust. Thus, we tested this hypothesis using a fully coupled poroelastic finite element numerical model. Many authors have exploited linear poroelasticity as a useful framework to describe postseismic displacements [Peltzer et al., 1998; Jónsson et al., 2003], aftershock evolution over space and time [Bosl and Nur, 2002], and the change in well water levels [Nespoli et al., 2016]. 3.1. The Poroelastic Model The three-dimensional coupling between fluid flow and medium deformations belongs to the theory of linear poroelasticity [Biot, 1941; Rice and Cleary, 1976; Wang, 2000]. This theory outlines how the stresses are transferred from the solid skeleton to the pore fluid, and vice versa. We exploited the built-in diffusion-structural module of the finite element software package MSC MARC 2015 [MSC Software Corporation, 2015], which allowed us to perform fully coupled poroelastic simulations assuming realistic fault geometries and spatially variable material properties. The constitutive laws that relate strains, fluid mass content per unit volume, stresses, and pore pressure are as follows [Segall, 2005]: 2Gεij ¼ σij ν 1 þ ν σkkδij þ ð Þ 1 2ν α 1 þ ν pδij (2) Δm ¼ m m0 ¼ ð Þ 1 2ν αρ0 2Gð Þ 1 þ ν ∙ σkk þ 3 B p
(3) qi ¼ ρ0 κ η ∂p ∂xj ρ0gδij   (4) Equation (2) relates the stresses and strains in a poroelastic medium, where G is the shear modulus; ν is the drained Poisson ratio; εij and σij are the strain and the stress tensor components, respectively; p is the pore pressure; δij is the Kronecker delta and α is the Biot-Willis coefficient. Equation (3) relates the change in fluid mass per unit volume (Δm) to both the sum of normal stresses (σkk) and pore pressure changes (p). In equation (3), ρ0 is the fluid density, m is the fluid mass content, m0 is the Journal of Geophysical Research: Solid Earth 10.1002/2017JB014009 ALBANO ET AL. THE 2012 EMILIA-ROMAGNA EARTHQUAKE 5 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 fluid mass content for a reference state given by the product of the fluid density ρ0 and the volume of the pore fluid per unit volume of solid (i.e., the porosity n), and B is the Skempton coefficient. Equation (4) is the constitutive law governing the pore fluid diffusion, i.e., the Darcy's law, where qi is the fluid mass flow rate, κ is the permeability, η is the fluid dynamic viscosity, and g is the gravitational acceleration. Equations (3) and (4) are combined with the mass conservation for the infiltrating pore fluid to obtain the diffusion equation. The poroelastic formulation implemented into the MARC code assumes that the solid grains constituting the medium are incompressible, and the Skempton (B) and Biot-Willis (α) coefficients can be written as follows [Rice and Cleary, 1976; Wang, 2000]: B≅ Kf nK þ Kf (5) α≅1 (6) where K and Kf are the bulk moduli of the frame and the pore fluid, respectively. Given the modeling assumptions, a set of eight parameters are required, i.e., the shear modulus G, the drained frame bulk modulus K (or, equivalently, the drained Young modulus E and Poisson's ratio ν), the fluid bulk modulus Kf, the porosity n, the permeability k, and the fluid density ρ0 and dynamic viscosity η. 3.2. Finite Element Model Setup The WNW-ESE alignment of the buried tectonic structures below the Po plain suggests that the fault slip vector is nearly orthogonal to the fault strike line, whereas the along-strike component of displacement is inhibited (i.e., rake ≈ 90°). This hypothesis is confirmed by seismological analyses and both coseismic analytical and numerical fault solutions (dotted blue and red boxes in Figure 1a) [Bignami et al., 2012; Pezzo et al., 2013; Cheloni et al., 2016]. To take advantage of the symmetry in the displacement, the geological section A-A0 (Figure 1b) is modeled to be parallel to the fault slip vector, i.e., orthogonal to the fault strike (Figure 1a). Moreover, a 2-D plane strain finite element model (Figure 3) that extends 70 km horizontally and to a depth of F3 20 km is constructed. The mesh is composed of eight-node, isoparametric quadrilateral elements (20,156 elements). The elements have a variable size, being coarser with increasing depth; the size ranges from approximately 0.35 km at the upper boundary to 0.875 km at the bottom of the model. The mechanical boundary conditions consist of roller support at the sides and bottom. Hydraulic boundary conditions consist of atmospheric pressure constraints at the upper edge, and the lower edge is assumed to be impermeable. The model sides are assigned the hydrostatic pore pressure condition to simulate a flowing boundary. The model geometry is discretized in five strata (Figure 3): a Quaternary layer (layer 1 in Figure 1b), a Pliocene layer (layers 2 and 3 in Figure 1b), an Oligocene-Miocene layer (layer 4 in Figure 1b), a Jurassic-Cretaceous Figure 3. Simplified finite element model of the geological cross section A-A0 of Figure 1b. The model grid extends 70 km horizontally and to depths of 20 km. The red line represents the 20 May 2012 causative fault according to Cheloni et al. . The upper (up) and lower (down) red arrows represent the applied edge loads to simulate the coseismic slip. The black dashed square indicates the area of contour plots of Figures 4, 6, and 7. Journal of Geophysical Research: Solid Earth 10.1002/2017JB014009 ALBANO ET AL. THE 2012 EMILIA-ROMAGNA EARTHQUAKE 6 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 layer (layer 5 in Figure 1b), and a crystalline basement (layer 6 in Figure 1b) [Boccaletti et al., 2004; Carminati and Vadacca, 2010; Carminati et al., 2010; Tizzani et al., 2013]. The 20 May event causative fault corresponds to the Middle Ferrara thrust (MF in Figures 1a and 1b) [Cheloniet al., 2016]. The fault is modeled using a contact interface [Carminati and Vadacca, 2010; Doglioni et al., 2011] dipping 40° and extending approximately 6 km (red line in Figure 3). Nodes are doubled at the interface; thus, the upper and lower parts of the domain can move relative to each other. The contact interface is assumed to be impermeable; thus, fluid cannot flow across it. Although fault planes are generally considered permeable surfaces, observations and experimental data show that as slip increases, faults develop a thick, impermeable gouge layer [Scholz, 1987; Parsons et al., 1999]. Consequently, the modeled fault mostly behaves as a barrier to fluid flow [Piombo et al., 2005]. The slip is simulated by applying edge loads to the top and bottom edges of the fault (up and down arrows in Figure 3). These loads are parallel but have different magnitude and opposite direction. The selection of appropriate elastic and hydraulic properties of the modeled geomaterials is not straightforward. In fact, local geometric and rheological anisotropies, joint density and orientation, and previous stress path significantly influences the mechanical and hydraulic properties of rocks at the scale of tens of meters [Carminati et al., 2010]. Model parameters depend on the size of the domain. In our kilometer-scale model, an equivalent continuum approach is adopted [Sitharam et al., 2001]. It is assumed that each stratum is continuous and isotropic, and we derive the equivalent parameters from the geomechanical properties of both the intact rock and joints. The elastic constants were derived from field and literature data. In particular, drained Poisson's ratios are selected according to literature data [Peltzer et al., 1998; Carminati and Vadacca, 2010; Tizzani et al., 2013]. The stiffness of layer 1 (Quaternary) is constrained with experimental data from a sonic log [Maesano and D'Ambrogi, 2017]. Stiffness for layers 2 to 5 are constrained according to an up-to-date seismic velocity model obtained in the study area after the 2012 seismic sequence [Carannante et al., 2015] (Figure S4). Then, dynamic Young moduli (Table 1) are converted to static Young moduli according to the empirical relationship T1 proposed by Brotons et al. (equation (7)): Est ¼ 11:531∙ρ0:457∙E1:251 dyn (7) where Est is the static Young modulus,ρ is the mass density and Edyn is the dynamic Young modulus (Figure S5). The derived static Young moduli (Table 1) are similar to those used in the same area by other numerical models [Carminati and Vadacca, 2010; Carafa and Barba, 2011]. Q3 Table 1. Elastic, Hydraulic, and State Parameters Adopted in the Numerical Analysisa Material property Parameter Layer 1, Quaternary Layer 2, Pliocene Layer 3, Oligocene-Miocene Layer 4, Jurassic-Cretaceous Layer 5, Crystalline Basement Mass densityb ρ(Kg/m3 ) 2100 2300 2400 2500 2600 Drained Poisson's ratiob ν 0.2 0.28 0.28 0.27 0.26 Drained dynamic Young's modulusc Edyn(Pa) 7.9 × 109 1.3 × 1010 1.8 × 1010 3.7 × 1010 6.7 × 1010 Drained static Young's modulusd Est(Pa) 4.6 × 109 8.5 × 109 1.2 × 1010 2.9 × 1010 6.1 × 1010 Porositye n 0.2 0.03 0.03 0.03 0.02 Permeabilityf k(m2 ) 5.0 × 1014 8.0 × 1014 8.0 × 1014 8.0 × 1014 8.0 × 1015 Fluid bulk modulus Kf (Pa) 2.2 × 109 Fluid dynamic viscosity η(Pa s) 0.001 Fluid density ρ0(Kg/m3 ) 1000 Skempton coefficient B 0.81 0.92 0.89 0.78 0.72 Fault friction coefficient μ 0.08 a These parameters best reproduce the observed pore pressure changes and ground displacements. b Drained Poisson's ratios obtained from literature data [Peltzer et al., 1998; Carminati and Vadacca, 2010; Tizzani et al., 2013]. c Drained dynamic Young's moduli from Carannante et al. and Maesano and D'Ambrogi . d Static Young's moduli are calculated according to equation (7) [Brotons et al., 2016] as a function of the dynamic Young's moduli (Edyn) and the rock mass density (ρ). e Mean value from Styles et al. . f Values from Nespoli et al. (quaternary) and Styles et al. (other layers). Journal of Geophysical Research: Solid Earth 10.1002/2017JB014009 ALBANO ET AL. THE 2012 EMILIA-ROMAGNA EARTHQUAKE 7 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 The hydraulic parameters were selected according to the literature and field data, when available. For the Quaternary layer, an isotropic permeability is assumed because of the large scale of the numerical model. However, according to Lambe and Whitman , the presence of several subhorizontal layers of alternating aquifers and aquitards considerably reduces the equivalent vertical permeability relative to the equivalent horizontal permeability. Therefore, assuming that the fluid flow is mainly vertical, we adopted a permeability (Table 1) close to the smallest value from Nespoli et al., . For the rocky strata, experimental data from diffusivity tests performed at depths between 2.7 and 3.1 km show that the permeability ranges between 4 × 1012 m2 and 3 × 1015 m2 (Figure S6) and that the porosity ranges between 1% and 7% [Laboratorio di monitoraggio Cavone, 2014; Styles et al., 2014]. Such experimental values are assumed to be typical of layers 2, 3, and 4, as these layers have been highly deformed by tectonics. Therefore, the permeability depends more on the degree of fracturing than on the depth. Accordingly, a single value of permeability for layers 2, 3, and 4 has been selected through a parametric analysis to numerically reproduce the amplitude and evolution over time of the observed pore pressures and postseismic ground displacements, and the other parameters have been held constant. The best value (Table 1) falls inside the experimental interval (Figure S6) and is common to other numerical models [LaBonte et al., 2009] and literature data for the upper crust [Kuang and Jiao, 2014]. The results for a range of permeability values will be discussed later. For layer 5, we assumed a permeability one order of magnitude lower than the above strata to account for the dependency of permeability on depth [Manga et al., 2012; Kuang and Jiao, 2014]. The coupling among the solid and fluid phases is simulated through three steps. In the first step, the initial stress and pore pressure distribution are established by applying a gravity load. The pore pressure distribution at the end of this phase increases linearly with depth according to field measurements [Laboratorio di monitoraggio Cavone, 2014]. In the second step, the fault responsible for the 20 May event is allowed to slip instantaneously by imposing two nonuniform edge loads on the upper and lower faces of the fault (Figure 3). No friction is assumed at the contact interface at this stage. The magnitudes of the loads (up = 3.25 × 106 N/m2 and down = 2.6 × 106 N/ m2 ) are calibrated following a trial and error procedure to fit the coseismic CSK LOS displacements (Figure 2a) along section A-A0 , i.e., the ground displacement data available 3 days after the main shock (23 May 2012). In the third step, no further loads are applied. The stress changes caused by the coseismic Δp relaxation are calculated for 2 years after the 2012 main shock. In this phase, a friction law obeying the Coulomb failure criterion has been assumed at the contact interface representing the main fault. The friction coefficient (μ) has been parametrically increased, ranging from an unlocked fault (μ = 0) until reaching to a locked fault (μ = ∞), to fit the observed postseismic displacements. During the postseismic phase, the evolution of the fluid pore pressures and displacements at the ground has been monitored over time. 3.3. Coulomb Stress Changes The change in Coulomb stress (ΔCFS) is determined as follows [Cocco and Rice, 2002]: ΔCFS ¼ Δτ þ μ∙ð Þ¼ Δσn þ Δp Δτ þ μ∙Δσ 0 n (8) where Δτ is the shear stress change (computed in the slip direction), μ is the static friction coefficient [Kinget al., 1994; Harris, 1998], Δσn is the fault-normal total stress change, Δp is the pore pressure change, and Δσ 0 n is the effective fault-normal stress change. The friction coefficient ranges between 0.6 and 0.8 for most rocks [Harris, 1998]. Because crustal rocks are fractured, the lower bound for the friction coefficient is assumed (μ = 0.6) [Piombo et al., 2005], and the Coulomb stresses on preferential fault planes dipping 30° southwest with respect to the horizontal are calculated according to the mean focal mechanisms of the aftershock sequence [ISIDe Working Group, 2010]. Because the numerical model adopts a positive sign convention for tension and a negative convention for compression, a positive ΔCFS indicates fault weakening. 4. Results Here we show the results that best reproduce the measured coseismic and postseismic displacements and the water table changes at wells. The performed parametric analyses will be discussed later. Results are presented in terms of pore pressures and ground deformations in the coseismic and postseismic phases. Journal of Geophysical Research: Solid Earth 10.1002/2017JB014009 ALBANO ET AL. THE 2012 EMILIA-ROMAGNA EARTHQUAKE 8 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 Regarding pore pressures, both positive and negative changes with respect to the hydrostatic values are hereinafter indicated using Δp, and we use "suprahydrostatic" in the case of a pore pressure increase and "subhydrostatic" in the case of a pore pressure decrease. 4.1. Coseismic Δp and Postseismic Evolution Before the main shock, the pore pressure is hydrostatic. In the coseismic phase, the fault dislocation causes shear and volumetric changes inside the medium. In particular, coseismic volumetric changes develop a Δp pattern reaching nearly ±2 MPa close to the fault zone (Figure 4a). These Δp values are transient and gradually F4 dissipate in the postseismic phase because of fluid diffusion. Indeed, the Δp values are 4 times lower after 10 days (Figure 4b) and 10 times lower after 30 days (Figure 4c). The fluid flow rate (black arrows) is higher in the early days after the earthquake (Figure 4b) and gradually decreases over time (Figure 4c). One year later (Figure 4d), the fluid diffusion ceases and the pore pressures correspond to the hydrostatic values. The temporal variation in Δp at depth (points 1 and 2 in Figure 4a) shows that subhydrostatic Δp values develop at point 1 (Figure 4e) after the coseismic slip because of volume dilation. These Δp values Figure 4. Pore pressure changes and decay with respect to the hydrostatic values. (a) Δp pattern after the simulation of the 20 May coseismic slip, (b) 10 days later, (c) 30 days later, and (d) 1 year later. The dashed lines in Figures 4a–4d separate areas with positive and negative Δp. (e) Coseismic Δp peak and postseismic decay at points 1 and 2 (yellow points) in Figures 4a–4d. (f) Coseismic and postseismic Δp values at point 3 compared with the measured water level change at well FE81 (Figure 1a). Journal of Geophysical Research: Solid Earth 10.1002/2017JB014009 ALBANO ET AL. THE 2012 EMILIA-ROMAGNA EARTHQUAKE 9 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 gradually diminish over approximately 3–6 months. At point 2, the coseismic volume compression induces suprahydrostatic Δp values. These Δp vlaues rapidly decay, reaching subhydrostatic values after approximately 70 days because of the arrival of a subhydrostatic fluid pulse from the neighboring regions of the domain. Subhydrostatic Δp values gradually return to hydrostatic values over approximately 3 months. At the surface (point 3 in Figure 4a), coseismic suprahydrostatic Δp values develop (Figure 4f). However, the pressure peak is lower than that at depth (e.g., at point 2). Hydrostatic pore pressures are rapidly restored because the selected point is very close to the ground surface; thus, it is constrained by atmospheric pressure. Moreover, the drainage path is very short. The modeled coseismic and postseismic Δp values are compared with the measured water level change at several wells close to the earthquake area. Indeed, the Quaternary aquifer response to the Emilia seismic event was captured by the monitoring network of the regional agency for environmental protection (Advanced Research Projects Agency) [Marcaccio and Martinelli, 2012]. In particular, three wells (FE80, FE81, and M080 in Figure 1a) recorded a significant water level change after the 20 May event. The water level experienced a sudden increase between 0.2 and 1.6 m, which was followed by a gradual decay that did not recover to the pre-earthquake water level in the observed time interval. The highest peak was observed for the deepest well (M080, depth = 300 m), which was followed by a slower decay with respect to the shallower wells (FE80, depth = 40 m; FE81, depth = 40 m). A smaller peak was also observed after the May 29 event [Marcaccio and Martinelli, 2012]. The comparison has been performed at the FE81 well, as this well is close to the modeled cross section. Wells M080 and FE80 were not considered because earthquake-induced Δp values have a great horizontal variation [Nespoli et al., 2016], and these wells are too far from the modeled 2-D section. The modeled Δp values at point 3 were converted into the equivalent water level change (1 MPa corresponds to a water column height of 101.97 m at 4°C) and compared with the measured water level change (Figure 4f). The observed coseismic peak and postseismic decay are well reproduced by the model, thus validating the Δp pattern shown in Figure 4a. 4.2. Coseismic and Postseismic Deformations In the coseismic phase, the modeled coseismic slip (Figure 5) produces a deformation pattern (Figure 6a) that F5F6 mainly involves a localized area of the medium above the dislocated fault. The displacement pattern and vectors, with a maximum amplitude of approximately 30 cm, highlight the slip along the causative fault and agree with the solution proposed by Tizzani et al. that was obtained from a similar 2-D finite element model, which assumed complete dryness. Our model is completely saturated; thus, the imposed coseismic slip results in the Δp pattern shown in Figure 4a. In the postseismic phase, the gradual dissipation of the coseismic Δp alters the stresses inside the medium, developing further deformations. Indeed, the modeled postseismic deformations 2 years after the 20 May event (Figure 6b) show a further uplift close to the area already affected by the coseismic inflation, which is consistent with the observed postseismic uplift from DInSAR (Figure 2b) and analogous postseismic scenarios [Fattahi et al., 2015; Cheloni et al., 2016]. The modeled coseismic and postseismic deformations are compared with the DInSAR-derived LOS displacements along section A-A0 . In particular, the first (23 May 2012) and last (10 May 2013) available Figure 5. Modeled coseismic (red line) and postseismic (blue lines) slip. The modeled fault is unlocked in the coseismic phase (i.e., μ = 0), while the friction coefficient is parametrically increased in the postseismic phase to fit the observed ground displacements from DInSAR. The plotted postseismic slip is obtained for μ = 0.08 (Table 1). Journal of Geophysical Research: Solid Earth 10.1002/2017JB014009 ALBANO ET AL. THE 2012 EMILIA-ROMAGNA EARTHQUAKE 10