The interseismic velocity field of the central Apennines from a dense GPS network

Since 1999, we have repeatedly surveyed the central Apennines through a dense survey-style geodetic network, the Central Apennines Geodetic Network (CAGeoNet). CAGeoNet consists of 123 benchmarks distributed over an area of ca. 180 km × 130 km, from the Tyrrhenian coast to the Adriatic coast, with an average inter-site distance of 3 km to 5 km. The network is positioned across the main seismogenic structures of the region that are capable of generating destructive earthquakes. Here, we show the horizontal GPS velocity field of both CAGeoNet and continuous GPS stations in this region, as estimated from the position–time series in the time span from 1999 to 2007. We analyzed the data using both the Bernese and GAMIT software, rigorously combining the two solutions to obtain a validated result. Then, we analyzed the strain-rate field, which shows a region of extension along the axis of the Apennine chain, with values from 2 × 10–9 yr–1 to 66·× 10–9 yr–1, and a relative minimum of ca. 20 × 10–9 yr–1 located in the L'Aquila basin area. Our velocity field represents an improved estimation of the ongoing elastic interseismic deformation of the central Apennines, and in particular relating to the area of the L'Aquila earthquake of April 6, 2009.


Introduction
According to current views, the Apennine chain is an arc-shaped, NE-verging belt that is characterized by a complex pattern of thrusts and folds and normal faults.These are related to two superimposed tectonic phases: an upper Miocene-Lower Pleistocene compressional phase that forms NW-SE-trending thrusts and folds; and the subsequent Quaternary extensional phase that forms NW-SE-trending normal faults.These latter are responsible for the formation of large intramontane basins that have been filled by Plio-Quaternary continental sediments (i.e., the L'Aquila, Rieti, Terni, Fucino and Sulmona basins) [Galadini and Messina 1994, Galadini and Galli 2000] (Figure 1).Some studies have explained the change in the tectonic regime as being caused by the flexural retreat through a roll back mechanism of the lithospheric Adriatic plate that dips below the Apennines [Reutter et al. 1980, Boccaletti et al. 1982, Malinverno and Ryan 1986, Royden et al. 1987, Patacca et al. 1990, Doglioni 1991, Doglioni et al. 1994, Frepoli and Amato 1997, Basili and Barba 2007].Other studies have ascribed the change in the tectonic regime as being caused by NE motion of the Adriatic microplate relative to Eurasia, around a rotation pole located in northwest Italy [Anderson and Jackson 1987, Calais et al. 2002, D'Agostino et al. 2005, 2008].At present, the geodetic data show that extensional deformation in the central Apennines is occurring along a narrow belt that is 30 km to 40 km wide [Hunstad et al 2003, Serpelloni et al. 2005, Devoti et al. 2008, Devoti et al. 2011], and which is near the areas where the strongest historical (intensity, ≥XI) and instrumental earthquakes have occurred [Boschi et al. 1998, Selvaggi 1998] (Figure 1).Starting from 1999, a dense survey-mode GPS network, the Central Apennines Geodetic Network (CAGeoNet), was designed and installed in the central Apennines [Anzidei et al. 2005[Anzidei et al. , 2008]].CAGeoNet consists of 123 benchmarks with an average inter-site distance of 3 km to 5 km, which are now surrounded by continuously operating GPS stations (Figure 2) and are positioned across the main active faults, as shown by geological and seismological data [Galadini andGalli 2000, Valensise andPantosti 2001].The high GPS-station density and the quality of the data collected have provided new insights into the presentday deformation of this seismically active area and information that is useful for seismic hazard assessment.
It has been shown that the combination of independent geodetic solutions obtained with different GPS-processing software [Avallone et al. 2010, Devoti et al. 2012] allows eventual systematic errors to be minimized and final velocity solutions to be validated.In the present study, we have estimated the interseismic strain rates from a combination of independent solutions obtained with the Bernese and GAMIT software.This has thus allowed us to investigate the geodetic deformation of the interseismic cycle of the Umbria-Marche Apennines (UMA) and the Lazio-Abruzzo Apennines (LAA).

The CAGeoNet and GPS campaigns
CAGeoNet has provided repeated measurements over the time interval of 1999 to 2007.The surveys were planned taking into account the network grid, the number of stations to be measured simultaneously (up to 11), and the time required to move receivers through the network.Consistent with the logistics, measurements were carried out in approximately the same period of the year, to minimize possible bias due to seasonal variations.Each station was occupied for an average observation window of 48 h, for at least three survey sessions per station, with a sampling rate of 30 s.
Here we discuss the interseismic deformation field that has resulted from the analysis of the velocities obtained from a sub-set of 55 CaGeoNet stations over the time interval of 1999 to 2007.

Data processing and combination procedures
The dataset analyzed (Figure 2) consists of GPS data that were collected on survey-style benchmarks (the CAGeoNet benchmarks) and continuous data that was provided by the continuous GPS networks located in the central Apennines region.The continuous GPS stations belong to different GPS networks: the International Global Navigation Satellite Systems (GNSS) Service (IGS) [http://igscb.jpl.nasa.gov]; the Integrated National GPS Network (Rete Integrata Nazionale GPS; RING) [Avallone et al. 2010]; the Italian Space Agency (Agenzia Spaziale Italiana; ASI) [Vespe et al. 2000]; Leica Geosystems (Italian Positioning Service [ItalPos] network); the Region of Abruzzo; and the Universities of Perugia and L'Aquila.
The GPS data cover the period from 1999 to 2007, and they were arranged into several clusters, each of which shared common fiducial continuous GPS stations that were used as anchor stations in the subsequent combinations.Each cluster was independently processed, and then combined through least-squares combination into a single daily solution.The GPS observations were processed using both the Bernese 5.0 [Beutler et al. 2007] andGAMIT 10.34 [Herring et al. 2006] software.
The Bernese processing was based on the Bernese Processing Engine (BPE) procedure that follows standard analysis for regional networks.The daily station coordinates together with the hourly troposphere parameters were solved using the a-priori Dry Niell troposphere model, with the corrections estimated by the Wet Neill mapping function.The ionosphere was neither estimated nor modeled as we used the L3 (ionosphere free) linear combination of L1 and L2.The a-priori GPS orbits and the Earth orientation parameters were fixed to the precise IGS products.We applied the ocean-loading Finite Element Solution model FES2004, and used the IGS absolute antenna phase-center corrections.The daily solutions were obtained in a loosely constrained reference frame; i.e., all of the a-priori station coordinates were left free to 10 m a-priori sigma.
The GAMIT processing followed the standard procedures for the analysis of regional networks [e.g., McClucsky et al. 2000, Serpelloni et al. 2006], with the application of loose constraints to the geodetic parameters.The GAMIT software used double-differenced, ionosphere-free linear combinations of the L1 and L2 phase observations to generate weighted least-square solutions for each daily session [Schaffrin andBock 1988, Dong andBock 1989].An automatic cleaning algorithm [Herring et al. 2006] was applied to post-fit the residuals, to repair cycle slips, and to remove outliers.The observation weights varied with elevation angle and were derived individually for each station from the scatter of the post-fit residuals that was obtained in a preliminary GAMIT solution.The effects of solid-earth tides, polar motion, and oceanic loading were taken into account according to the International Earth Rotation and Reference Systems Service (IERS)/IGS standard 2003 model [McCarthy and Petit 2004].We applied the FES2004 ocean-loading model and used the IGS absolute antenna phase-center correction table to model the effective receiver and satellite-antenna phase centers.We used orbits provided by the Scrips Orbit Permanent Array Center (SOPAC).The estimated parameters for each daily solution included the three-dimensional Cartesian coordinates for each station, the six orbital elements for each satellite, the Earth orientation parameters (i.e., pole position and rate, UT1 rate) and the integer phase ambiguities, with the application of loose constraints (ca. 10 m) to the a-priori parameters.We also estimated the hourly piecewise linear atmospheric zenith delays at each station, to correct the poorly modeled troposphere, and three east-west and north-south atmospheric gradients per day, to allow for azimuth asymmetry; the associated error covariance matrix was also computed and saved in SINEX format.
Both the analysis procedures (Bernese and GAMIT) produced daily loosely constrained solutions; i.e., free from any a-priori reference frame datum.Coordinates and the com-CENTRAL APENNINE DEFORMATION  plete associated covariance matrices were saved in Solution Independent Exchange (SINEX) format.
The time series of the two solutions were then obtained by applying minimal inner constraints and a four-parameter Helmert transformation, to obtain the coordinates and errors expressed in the IGS05 reference framework (the IGS realization of the ITRF2005 reference framework).Then, we obtained a velocity field for each solution, with the estimation of a linear drift (velocity), annual sinusoid, and occasional offsets due to changes in the station equipment between each time series.

Combined GPS velocity field
The two independent velocity solutions were combined into a unique velocity solution, using a linear least-squares combination approach.The normal matrix was formed from the two independent velocity solutions, and then it was inverted to estimate the unified velocity field of the entire network.As the covariance matrix is usually known separately from a constant multiplier, we also estimated a solution-scale factor together with the combined velocity solution.This ensured that the individual | 2 of each velocity solution was equally balanced (individual solutions do not prevail in the combination process) and the total | 2 is close to unity (realistic errors).The combined solution represents a weighted velocity average that takes into account the correlation ma-trices of the two solutions.The differences between the combined and the individual solutions have low mean values and comparable standard deviations (Figure 3).The values reported in Figure 3 show how the combined solution is placed between the Bernese and GAMIT solutions, and does not give priority to any individual solution.The Bernese solution is slightly more noisy in the east velocity component, while the GAMIT solution is slightly more noisy in the north velocity component.The comparison between the Bernese and GAMIT solutions shows residuals with averages of -0.2 mm yr -1 (Ve component) and 0.3 mm yr -1 (Vn component), and dispersions (at the 1v level) of 0.7 mm yr -1 and 0.9 mm yr -1 for Ve and Vn, respectively.The differences are comparable with the averaged sigma values computed for the permanent stations used in the individual solutions (v E = 0.14 mm; v N = 0.23 mm), which highlights that the solutions are compatible.We carried out statistical screening between the combined solution with respect to each single solution (Bernese and GAMIT) and the single solutions with respect to each other, both for the permanent and nonpermanent stations, with comparisons of the differences obtained between the average values of Ve and Vn.The stations that showed differences greater than the respective 2v values were discarded.
The combined velocity field with respect to a Eurasianfixed reference frame is shown in Figure 4.The velocity components and their uncertainties are reported in Table 1.The Eurasia plate was fixed, minimizing the horizontal velocities of 24 stations located on the stable part of the plate.
The selection of Eurasian stations was statistically inferred using a | 2 test-statistic to select the subset of stations that defined the stable plate [Noquet et al. 2001], starting from the triad of the WSTR, WTZR and ZIMM stations in ITRF2005.
The estimated Euler pole and rotation rate for the Eurasia plate were at 55.85˚N, 95.72˚W and 0.266˚ ± 0.003˚ Myr -1 respectively.
The geodetic strain rate was evaluated by a distanceweighted approach, and computed using all of the stations on a regularly spaced grid, applying the weighting algorithm developed by Shen et al. [1996].The contribution of each station velocity to the strain-rate computed on a given node was down-weighted with the function W = exp(-d 2 /a 2 ), where d is the distance between each node and the stations, and a is the smoothing distance parameter.The algorithm selects the optimal a value from a given a-priori interval, which depends on the spatial distribution of the GPS sites; consequently, strain-rate maps were obtained with spatially variable a.
The second invariant rate was obtained by interpolation of the velocity horizontal components on a 0.1˚ × 0.1˚ regular grid.The smoothing factor down-weights the velocities in the range from 20 km to 100 km, according to the network density (Figure 6).

Results
The combined horizontal velocity field expressed with respect to a fixed Eurasian plate showed: (i) good coherence between the velocities estimated from the CAGeoNet ( and (ii) two different and characteristic main velocity patterns: a NNW-oriented trend in the Tyrrhenian Apennine sector, and a NNE-oriented trend in the Adriatic Apennine sector.A gradual clockwise velocity rotation was clearly evident from W to E, where the velocities were initially NNWoriented and rotated towards the NNE with increasing values (0.9-5.2 mm yr -1 ).This pattern showed an anomaly in the L'Aquila basin, where the vectors were turned in an ca.NS directed, as normal to the main tectonic structures of the Gran Sasso range (Figure 4).
To show the velocity gradients better across the central Apennine chain, we represented the velocity field with respect to a fixed Tyrrhenian coast (Figure 5), and then projected the velocities along two profiles crossing the studied area: ENE-WSW-oriented.The two profiles are parallel to the average direction of the velocity vectors and approximately normal to the main fault systems.Profile (1) was located across the UMA, while profile (2) was located across the LAA sector.The projections contained all of the velocities within the distance of 40 km (profile 1), and 30 km (profile 2).To highlight the velocity gradient, we used a moving-average filter with a 40-km window (Figure 5, gray line).The net extension rates across the two cross-sections were the same, at ca. 2.5 mm yr -1 , but they spread over different distances.Profile (1) showed a velocity variation that was concentrated in a narrow strip of ca.60 km, with a maximum step of 1.5 mm yr -1 at about 30 km, on the western flank of the chain with a strain rate of ca.50•× 10 -9 yr -1 .Profile (2) showed a more irregular velocity variation, with a negative (i.e., shortening component) gradient that roughly corresponded with the L'Aquila basin (between 120 km and 150 km; Figure 5, profile 2), and developed along larger distances (ca. 100 km.) with respect to profile (1), with a lower strain rate of ca.20•× 10 -9 yr -1 .The predominantly extensional deformation was mainly oriented NE-SW, and ranged from 2 ±11•× 10 -9 yr -1 to 66 ±19•×10 -9 yr -1 (Figure 6).
For the UMA, the extensional deformation was distributed over a relatively wide area, and coincided with the culmination of the topographic relief.The distribution of the extension area became narrower and shifted towards the SW, crossing the Ancona-Anzio Line, and had slightly lower extension rates.
These estimations are in general agreement with those obtained by previous geodetic and geologic data for this area.From the geologic data, which included both interseismic and coseismic deformation, Galadini and Galli [2000] obtained an extension rate from 0.7 mm yr -1 to 1.6 mm yr -1 across the main active fault sets recognized in the area.Using fault-slip vectors, Faure Walker et al. [2010] calculated strain rates averaged over 15 kyr of 12•× 10 -9 yr -1 and an extension rate of 1 mm yr -1 over a 160 km × 80 km area, which was consistent with a strain rate ≤38•× 10 -9 yr -1 estimated in 5 km × 80 km boxes that crossed the strike of the central Apennines.From geodetic data, Serpelloni et al. [2005] indicated 31•× 10 -9 yr -1 , while D' Agostino et al. [2008] showed a second invariant band parallel to the chain with values >50•× 10 -9 yr -1 .Devoti et al. [2008Devoti et al. [ , 2011] ] estimated an extension rate of 50•× 10 -9 yr -1 .Our results do not agree with those obtained by Pesci et al. [2010] on a CAGeoNet network sub-set.They found NE shortening corresponding to the UMA, NW shortening in the Liri Valley, and a wide area that was characterized by NE extension in the eastern portion of the LAA.

Figure 2 .
Figure 2. The Central Apennines Geodetic Network (CAGeoNet, triangles) and continuous GPS stations (blue dots) used in the analysis.Inset: Location of the study region within Italy.UMA, Umbria-Marche Apennines; LAA, Lazio-Abruzzo Apennines.

Figure 4 .
Figure 4. GPS combined velocity fields estimated for the time span from 1999 to 2007 with respect to the Eurasian plate.Red arrows, Central Apennines Geodetic Network GPS velocities; blue arrows, continuous GPS velocities.Inset: Location of the study region within Italy.UMA, Umbria-Marche Apennines; LAA, Lazio-Abruzzo Apennines.

Figure 5 .
Figure 5. Velocity projections along the transect directions 1a-1b and 2a-2b.The projections involved vertices at distances of 20 km and 15 km in both perpendicular directions along the transect directions for the profiles 1a-1b and 2a-2b, respectively.Location of the study region within Italy.UMA, Umbria-Marche Apennines; LAA, Lazio-Abruzzo Apennines.

Table 1
(continues on following page).Table1 (continues from previous page).Velocity field dataset of the CAGeoNet and of the surrounding continuous GPS stations used in the analysis.East, North, velocities; sigE, sigma E; sigN, sigma N; smaj-ax, smin-ax, error ellipses; Azim, azimuth.