New insights into earthquake precursors from InSAR

We measured ground displacements before and after the 2009 L’Aquila earthquake using multi-temporal InSAR techniques to identify seismic precursor signals. We estimated the ground deformation and its temporal evolution by exploiting a large dataset of SAR imagery that spans seventy-two months before and sixteen months after the mainshock. These satellite data show that up to 15 mm of subsidence occurred beginning three years before the mainshock. This deformation occurred within two Quaternary basins that are located close to the epicentral area and are filled with sediments hosting multi-layer aquifers. After the earthquake, the same basins experienced up to 12 mm of uplift over approximately nine months. Before the earthquake, the rocks at depth dilated, and fractures opened. Consequently, fluids migrated into the dilated volume, thereby lowering the groundwater table in the carbonate hydrostructures and in the hydrologically connected multi-layer aquifers within the basins. This process caused the elastic consolidation of the fine-grained sediments within the basins, resulting in the detected subsidence. After the earthquake, the fractures closed, and the deep fluids were squeezed out. The pre-seismic ground displacements were then recovered because the groundwater table rose and natural recharge of the shallow multi-layer aquifers occurred, which caused the observed uplift.

The identification of earthquake precursors is a key area of investigation in modern seismology. Precursor signals have not been conclusively identified so far, even using modern geodetic techniques (e.g., GPS and SAR interferometry) [1][2][3] . However, the latest satellite missions and processing algorithms may represent progress toward this goal.
We present evidence of ground deformation preceding the 2009 L' Aquila earthquake. We observed this deformation using multi-temporal InSAR techniques, and we propose a plausible causative mechanism.
The April 6, 2009 M w 6.3 L'Aquila earthquake occurred in central Italy. The event was generated by a 13-km-long, NW-striking and SW-dipping normal fault that borders the L' Aquila basin to the northeast [4][5][6][7] . Centroid-moment tensor solutions and the distribution of aftershocks define a fault plane that dips 45° between depths of 1 to 10 km 8,9 . We investigated the seismic cycle associated with the earthquake by applying advanced InSAR techniques to SAR datasets from various satellite missions (RADARSAT-2, Envisat and COSMO-SkyMed) that differ in terms of their wavelengths and spatial resolutions. Specifically, 159 RADARSAT-2 images (75 collected along descending orbits and 84 collected along ascending orbits) that span the six years preceding the mainshock were processed with the SqueeSAR TM software package 10 , and 38 Envisat images collected along descending orbits that cover almost the same temporal interval as the RADARSAT-2 data were processed with the IPTA software package 11 . In addition, 35 COSMO-SkyMed images that were collected along ascending orbits and cover the 16 months following the earthquake were processed using the Persistent Scatterer Pair (PSP) technique 12

Results
InSAR results. The RADARSAT-2 and Envisat mean velocity maps show very similar results during the pre-seismic phase; however, no significant spatial pattern of linear deformation trends was observed in these maps ( Figure S1). To identify possible non-linear ground displacement over time, we computed maps of ground acceleration, which was calculated as double the quadratic coefficients of second-order polynomials fitted to the displacement time series of each persistent scatterer (PS). The computed acceleration maps (Figs 1A, S2A and B) reveal noteworthy variations in the deformation velocities measured using the data collected along both ascending and descending orbits. In particular, positive values were observed within the imaged area except in two regions, the Quaternary Preturo and Pizzoli basins. These basins, which are located on the northwestern rim of the coseismic deformation field ( Figure S3), are characterized by generally negative values of up to −3 mm/ yr 2 . Similar patterns can be observed in the data collected along both ascending and descending orbits, demonstrating the existence of predominantly vertical movement. We thus averaged the time series of the PSs within each basin in order to reduce oscillations and noise in the data and to obtain corresponding time series for the In contrast, during the post-seismic phase, the COSMO-SkyMed ascending velocity map (Fig. 1B) shows subsidence associated with tectonic (afterslip) and gravitational phenomena 13,14 over the entire imaged area (negative values), with the exception of the two basins mentioned above. The satellite data indicate progressive uplift within these basins that is characterized by velocities of approximately 5 to 18 mm/yr. This uplift began approximately 8 months after the mainshock and lasted at least 9 months (Fig. 1C).
This analysis emphasizes that the Preturo and Pizzoli basins are the only two regions within the imaged area that subsided before the mainshock and then experienced uplift during the post-seismic phase (Fig. 1). Indeed, the neighboring Quaternary basins ( Fig. 2A,B) show different patterns of displacement during the pre-and post-seismic phases. Namely, during the pre-seismic phase, the averaged time series of the PSs within each basin show very similar seasonal oscillations superimposed on multi-annual ground deformation (Fig. 2C). During the post-seismic phase, all of the neighboring basins show negligible ground deformation, except for the Paganica-Fossa basin, where ground subsidence associated with the afterslip 13 occurs. The Paganica-Fossa basin is located in the hanging wall of the earthquake fault. Geological, hydrogeological and geotechnical setting. This section provides an outline of the geological, hydrogeological and geotechnical features of the Preturo and Pizzoli basins that supports the interpretation of the observed pre-and post-seismic deformation.
The Preturo and Pizzoli plains belong to the intramontane Quaternary basin of L' Aquila, the evolution of which is controlled by two active and SW-dipping normal faults, the Pettino Mt. and Marine Mt. faults. The plains are underlain by continental detrital deposits that unconformably overlie Meso-Cenozoic bedrock (Fig. 3). The latter rocks and the surrounding mountains consist primarily of carbonates, with some terrigenous Meso-Cenozoic lithologies. The Quaternary succession can be divided into three groups of synthems that vary in thickness 15,16 (Fig. 4). Here, we describe these synthems from bottom to top. The basal material within the basins is composed of approximately 50-100 m of coarse-grained upper Piacenzian-Gelasian deposits (slope-derived breccia, debris flow and proximal alluvial fan deposits) that are rarely exposed (layer 3 in Fig. 4B). These deposits are overlain by approximately 50 m of Calabrian meandering fluvial-floodplain-swamp sediments that are characterized by sandy silt and sand beds interlayered with clay and clayey silts (layer 2 in Fig. 4B). Lignite intercalations also occur, and these intercalations become more frequent and thicker to the west. Layers of coarse to medium sand and gravel are present in the middle and lower parts of this sedimentary package. The third synthem consists of upper Pleistocene-Holocene braided fluvial floodplain deposits and slope-derived, clast-supported, medium to coarse gravel beds with primarily sandy matrix material (layer 1 in Fig. 4B).
The Preturo and Pizzoli plains are positioned between two major carbonate aquifers (the Gran Sasso aquifer to the north and the Nuria-Velino Mts. aquifer to the south and west) that differ in terms of their hydrogeologic relationships with the detrital multi-layer aquifers corresponding to the Quaternary basin-filling deposits 17 (Fig. 3). The detrital multi-layer aquifers are fed by groundwater flowing laterally out of the Gran Sasso carbonate aquifer toward the center of the plains and by precipitation falling on the plains. Moreover, the presence of bedrock that is composed of low-permeability Miocene terrigenous and marly lithologies and underlies the multi-layer aquifers in the southern and western parts of the basins indicates that the Preturo and Pizzoli plains represent a minor secondary discharge area for the Nuria-Velino Mts. hydrostructure.
The Gran Sasso and Nuria-Velino Mts. carbonate aquifers are fissured and fault-partitioned media. Infiltration is very high (800 to 1200 mm/yr, which represents more than 60% of total precipitation) 18 . In contrast, runoff is negligible (approximately 5% of total precipitation). The movement of groundwater may be controlled by the presence of faults, which represent discontinuities with high or low hydraulic gradients and low-or high-permeability boundaries. Nevertheless, at the regional scale, the groundwater within the carbonate aquifers may be regarded as continuous. The substantial quantity of groundwater within the carbonate aquifers flows from the core of each aquifer to the borders, where the aquifers are drained by springs that have steady and high discharges (0.1-10 m 3 /s) (Fig. 3). The discharge of groundwater is controlled by the morphology and the altitude of the permeability boundary. The aquifer boundary coincides with the contact with the detrital fill underlying the two plains. The boundary is not totally free of discharge, and groundwater flows into the Quaternary multi-layer aquifer. These background conditions promote the recharge of the Quaternary multi-layer aquifer by the carbonate aquifers.
Regarding the hydrogeologic behavior of the Quaternary deposits, the Calabrian fine-grained sediment (layer 2 in Fig. 4B) corresponds to an aquitard that hydraulically confines or locally semi-confines the underlying upper Piacenzian-Gelasian coarse-grained deposit, whereas the overlying upper Pleistocene-Holocene coarse-grained deposit and the sandy layers at the top of the Calabrian succession correspond to an unconfined aquifer. In spite of this hydrogeologic setting, in which the Calabrian aquitard is enclosed between two aquifers (i.e., the upper Piacenzian-Gelasian deposit and the upper Pleistocene-Holocene coarse-grained deposit, which are found above and below the aquitard, respectively), the piezometric equilibrium of the Quaternary multi-layer aquifer is based on a single water table. Thus, an upward component of groundwater flow into the Calabrian aquitard that occasionally displays artesian behavior cannot be excluded.
Earthquake-induced paleohydrogeological phenomena (sinkholes and liquefaction features) have occurred in the two plains during the middle Pleistocene and in historical times. Together with modern piezometric data, these observations confirm the hydrogeological background described above 19 .
The main geotechnical properties of the Quaternary sedimentary soils are derived from a series of in situ surveys performed in the Preturo basin (Fig. 4A) and laboratory testing 20 . The typical stratigraphy of the area determined from borehole logs (Fig. 4B) shows a sequence of coarse-and fine-grained materials with strata of variable thickness. The results of laboratory tests on undisturbed fine-grained samples taken from the available boreholes show that this material consists primarily of silts ( Figure S4A) (ranging between 20% and 80%), with a slightly lower content of clays (between 20% to 50%) and a small percentage of sand (between 0% and 20%). The consistency limits (i.e., the plasticity index, PI, and the liquid limit, LL) plotted on a Casagrande chart ( Figure S4B) show that these sediments have a homogenous mineralogical composition and that the clay component exhibits medium to high plasticity. The compressibility and swelling indexes (denoted with C C and C s , respectively, in Figure S5C and D) obtained from oedometer tests do not show any significant dependence with depth. Overconsolidation caused by seasonal fluctuations in the water table is present in the upper meters of the basin fill; in Figure S5B, the overconsolidation ratio (OCR) ≫1.
Based on the available results, the fine-grained soils are highly susceptible to consolidation in cases in which the stress field is modified, leading to ground subsidence.

Discussion
Based on the observed InSAR ground deformations and on the regional and local geological, geotechnical and hydrogeologic features, we hypothesize that the most plausible explanation for the observed pre-seismic subsidence and post-seismic uplift is the consolidation and swelling of the fine-grained content of the alluvial deposits that fill the Preturo and Pizzoli basins.
The consolidation would have been caused by the long-term (i.e., not seasonal) lowering of the water table in the multi-layer aquifers in the basins. The contemporaneous, spatially homogeneous subsidence of the two basins is related to a decrease in the aquifer in the Gran Sasso hydrostructure. Indeed, these aquifers are hydrologically connected to the major hydrostructure of the Gran Sasso carbonate range (Fig. 5A); hence, large variations in the former are controlled by large variations in the latter. A lowering of groundwater levels in the Gran Sasso carbonate range was observed before the L' Aquila earthquake 21,22 . This decrease induced the migration of water from the connected multi-layer aquifers in the basins, which led to the observed subsidence (Fig. 5B).
On the other hand, the post-seismic uplift of the basins was induced by rising groundwater levels and the consequent recovery of the elastic portion of the ground settlement (Fig. 5C). Indeed, post-seismic increases in groundwater levels in the Gran Sasso aquifer have also been measured [22][23][24][25] .
Dilatancy theory accounts for the described process of water table variations 23,[26][27][28][29][30] . It implies the formation of cracks due to shearing. Consequently, water diffuses through the rocks within the earthquake source volume (i.e., the focal volume).
In particular, during the pre-seismic phase, shear stresses were applied to the rocks in the focal volume, causing fractures and voids to open and the volume to increase. Hence, water migrated into the focal volume from the surroundings (i.e., the Gran Sasso aquifer), thus lowering the groundwater in the hydrologically connected multi-layer aquifer within the basins. Indeed, a complex sequence of dilatancy and the diffusion of fluid toward the earthquake nucleation area has been inferred from V p /V s variations that began in October 2008 31,32 . Moreover, according to Stacey 33 , a three-year period represents a plausible preparatory phase for an earthquake with a magnitude comparable to the 2009 event.
During the post-seismic phase, the increase in groundwater levels likely resulted from both seasonal recharge of the aquifers (precipitation) and squeezing of the deep groundwater caused by coseismic dislocation, which led to the closure of fractures and voids and the expulsion of water from the focal volume 34,35 .
Although other phenomena can explain the water table oscillations in the two Quaternary basins, no significant industrial or agricultural activities occur in these basins. Thus, anthropogenically driven changes in the groundwater level can be neglected.
The effects of seasonal or multi-annual rainfall variations can be ruled out as well. In fact, we analyzed seasonal and multi-annual precipitation records collected at eight pluviometric stations located within the study area between 2004 and 2010 ( Fig. 2A,B). A comparison of rainfall and ground deformation shows that the ground deformation within the Quaternary basins in the vicinity of Preturo and Pizzoli is clearly correlated with both seasonal and multi-annual rainfall variations (Fig. 2C). Conversely, the ground deformations that occurred within the Preturo and Pizzoli basins display no correlation with any of the multi-annual rainfall variations, only the seasonal ones (Fig. 1C).
Finally, we evaluated the consolidation potential of the fine-grained soils of the Preturo and Pizzoli basins in accordance with the classical principles of Terzaghi 36 , and we compared the numerical results with the InSAR observations.
The key factors affecting the consolidation of soils are the lowering of the water table, which reduces the pore water pressure and increases the overburden effective stress, and the compressibility of the soils at different depths. The magnitude and distribution of subsidence over the studied area should thus reflect the combined action of these two factors.
The analysis is made with reference to the simplified 1D soil profile shown in Fig. 6A. The undisturbed water table level is assumed to be equal to the mean value measured in the available boreholes (approximately 5 meters below the ground level), and the water table drop from the undisturbed condition is varied between 1 and 3 meters. In this calculation, the compressibility and swelling indexes (C C and C s ), the natural unit weight (γ) and the initial void index (e 0 ) have been assigned constant values that are equal to the means of the experimental values shown in Figure S5  The results of the consolidation analysis, which is conducted for different depths (z) and water table drops (1-3 m), are shown in Fig. 6B. The computed settlement values are projected along the line of sight of the RADARSAT-2 satellites to enable comparison with the observations.
The results show that the cumulative InSAR pre-seismic ground displacement range at the time of the April 6 event (gray band in Fig. 6B) is well reproduced by the model when a lowering of the water table between 1 and 3 meters is assumed. The thickness of the fine-grained soils must lie between 10 and 25 meters, in agreement with the retrieved maximum and minimum thicknesses of the fine-grained deposits filling the two Quaternary basins. The similarity between the observed and simulated ground displacements provides a logical explanation for the observed subsidence, which can be considered to be the result of the delayed groundwater flow induced by dilatancy at depth, which caused a drop in the groundwater table. Moreover, the results of the consolidation analysis show that the computed displacements are elastic because the computed effective stress increase caused by the groundwater table lowering is always lower than the preconsolidation stress (Fig. 6C). Consequently, the computed ground displacements are recoverable if the groundwater table rises to its initial level (see the Methods section for details), thereby explaining the amplitude of the uplift measured during the post-seismic phase using the COSMO-SkyMed data (Fig. 1B,C).
In conclusion, this study provides an interpretation of the surface deformation preceding the L' Aquila earthquake in terms of seismic precursors. The InSAR-derived ground accelerations reveal signals that may not be evident in standard mean velocity maps. In seismogenic areas characterized by the presence of Quaternary deposits and favorable hydrogeological settings (i.e., hydraulically connected shallow and deep aquifers), ground acceleration maps can be powerful tools for predicting earthquakes.
Although this study investigates only one earthquake, we believe that it provides a new methodology that can be applied in studies of earthquake precursors and is novel in terms of the technology used and the signals that can be investigated.
The pre-seismic sequence was investigated using RS2 data captured along ascending and descending orbits and Envisat data captured along descending orbits. In contrast, the post-seismic sequence was investigated using CSK data.
The RS2 ascending dataset extends from March 4, 2003 to March 3, 2009. The SAR images were acquired at an incidence angle of 33° and an azimuth of 10°. The RS2 descending dataset extends from April 15, 2003 to March 14, 2009. This dataset was processed with the SqueeSAR technique 10 and includes 84 and 75 SAR images collected along ascending and descending orbits, respectively.
The Envisat data, which were acquired along descending orbits, includes a stack of 50 images collected between November 10, 2002 andMarch 8, 2009. The line-of-sight direction is defined by a 24.5° incidence angle measured relative to the vertical direction and an aspect angle of −12.2 from the east. The dataset was processed with GAMMA IPTA software 11 using a single reference stack.
The CSK dataset was acquired between April 12, 2009 and August 6, 2010 along ascending orbits with an incidence angle of 40°. This dataset consists of 35 SAR images covering an area of approximately 40 × 40 km 2 . The data were processed with the persistent scatterer pair (PSP) technique 12 , which permits estimation of displacement time series for dense sets of locations that correspond to identified persistent scatterers.
The obtained results present a standard deviation of the velocities between 0.6-0.9 mm/year. On the other hand, the standard deviation of the deformation measurement error is on the order of 2 mm and typically ranges between 1 and 3 mm, depending on the quality of the specific PS and the type of ground cover. Consolidation analysis. In the present study, the long-term settlement w caused by groundwater lowering was calculated assuming simplified 1D conditions (equation 1), in accordance with the classical principles of Terzaghi 36 : where H represents the thickness of the compacting layer, and ε z is the vertical strain produced at each depth z by the increase in the effective overburden stress caused by the groundwater head drop. We assumed that the compressibility of the coarse-grained materials is negligible. This assumption was made because the compressibility indexes of silty-clayey and gravelly-sandy soils typically differ by some orders of magnitude. Indeed, whereas the oedometer tests performed on the fine-grained materials yielded compressibility index values ranging between 0.2 and 0.4 (C C in Figure S5), values that are lower by some orders of magnitude are typically found for gravelly soils 37 . Equation 1 implies that horizontal strains are negligible, i.e., the soil deforms under oedometric conditions (ε ε = = 0 x y in Figure S6). This assumption is acceptable if we consider the almost flat geometry of the Preturo and Pizzoli basins (Fig. 4b) and a uniform lowering of the water table. Under these hypotheses, the resulting stress path corresponds to the red curve in Figure S6.
Discretizing the compacting layer into n sublayers with thickness Δ z , equation 1 simplifies to the following: where Δ zi and e 0 are the thickness and the initial void index of the i th layer, respectively, whereas Δ ei is the change in the void index of the i th layer and is given by the following: In equation 3, C C and C S are the compressibility and swelling indexes, respectively; σ′v i 0, is the initial effective vertical stress; σ ∆ ′v i is the effective vertical stress increase; and ′ p c i , is the preconsolidation stress, which is given by the following equation (4): where OCR is the overconsolidation ratio ( Figure S5B). The vertical displacements computed using equations 2 and 3 may or may not be elastic (i.e., recoverable), depending on the relative positions of σ′v i 0, , ′ p c i , and σ ∆ ′v i on the oedometric curve. Initially, the soil has a stress state that is characterized by particular values of the effective vertical stress σ′v i 0, and the void index e 0 (point a in Figure S6). An increase in the vertical effective stress σ ∆ ′v i causes the point to move to the right along the stress path (red line in Figure S6). However, if the initial vertical stresses exceed the preconsolidation stress (σ′ > ′ v p i c i 0, , at point a'), the point moves along a stress path with a slope equal to C C , i.e., the virgin curve ( Figure S6), thus producing elastic (recoverable) and plastic (unrecoverable) deformations. In contrast, if the initial vertical stress is lower than the preconsolidation stress σ′ < ′ v p i c i 0, , at point a), the point moves along a stress path with a slope equal to C S , i.e., the swelling curve ( Figure S6). In such cases, if the increase in the effective stress does not overcome the preconsolidation stress (σ σ ′ + ∆ ′ < ′ v v p i i c i 0, , at point b), the resulting strains are elastic (recoverable). On the other hand, if the increase in effective stresses overcomes the preconsolidation stress (σ σ at point c), the resulting strains will be partially unrecoverable.