GNSS networks for geodynamics in Italy

The recent development of GNSS networks in Italy ma kes it possible to define an increasingly detailed spatial and temporal resolution of the ongoing crus tal deformation and to visualize the complex interplay between different orogens in the Africa-E urasia collision zone. We have analyzed all the available GPS raw data in the time span 1998-2013 o btaining the coordinate time series in a common reference frame. We have finally produced the veloc ity field providing a detailed picture of the kinematics and deformation pattern of the Italian a rea. The horizontal velocities with respect to Eurasia are within 7 mm/yr; vertical rates reach ab out 5 mm/yr along the Apennines and Alpine belts, and -10 mm/yr in the Tyrrhenian backarc area. The d eformation rates are greater than 100 nanostrain/yr in volcanic areas and within 80 nanos train/yr along the Apennines.


Introduction
The Italian peninsula is a rather interesting natural laboratory for geodynamical investigations since its tectonic evolution is driven by the interaction of the African and Eurasian plates.The entire area is characterized by a complex tectonic setting generating slow crustal deformations (at the few mm/yr level) and velocity spatial gradients (strain rates) changing rapidly from point to point.In this perspective, the rapid development of GNSS networks, with relatively low cost and high accuracy positioning, provides great advances in geodynamical studies.
The first attempt to build an Italian nation-wide continuous GPS network was undertaken by the Italian Space Agency (ASI) in the late 1990's.Since then, ASI delivers continuous GPS data from about 30 sites and maintains the regional reference frame in cooperation with the European reference frame consortium (EUREF).In 2001, the Istituto Nazionale di Oceanografia e Geofisica Sperimentale (INOGS) started installing a local GPS network in the Friuli region (northeastern Italy) and in 2004, the Istituto Nazionale di Geofisica e Vulcanologia (INGV) started the construction of the first national GPS network (RING) devoted to geophysical studies (Avallone et al., 2010).
At present, the RING network consists of about 180 stations.More recently an increasing number of permanent GPS sites have been installed by regional administrations and private companies, dedicated mainly to topographic applications and commercial services.These networks have not been conceived to measure long term ground deformations, and the GPS stations have been constructed on a variety of different monument types operating in fairly different environmental conditions, i.e. in the towns or industrial settlements, in the open field or mountain regions.Nevertheless they demonstrate a performance comparable to GPS stations purposely designed for geophysical monitoring, therefore useful to increase the backbone of more reliable geodynamic networks (Devoti et al., in press).
All these datasets are currently archived and processed at INGV providing over 600 raw-data files per day for a mean geometric inter-distance of about 20 km over the whole country, thus realizing a strong and valuable dataset for geodynamical studies.

Short geodynamic background
The African-Eurasian convergence is a large scale plate collision processes estimated by the NUVEL-1A plate motion model at about 7 mm/yr along the northwestsoutheast direction (DeMets et al., 1994); this process produces in the central Mediterranean a diffuse area of deformation and fragmentation of lithosphere in subplates and includes a varied complexity of tectonic patterns such as subduction, backarc spreading, rifting, thrusting, normal and strike-slip faulting and volcanism.
The Italian area is characterized by a complex tectonic setting where two very different orogens, the Alps and the Apennines, interfere and cause vast areas to deform in a compound way (figure 1).These two orogens belong to larger systems of mountain belts bordering the basins of the Mediterranean region: the Alps-Betics and Dinarides and the Apennines-Maghrebides and Carpathians, respectively double-vergent and single-vergent belts.The Alps are related to the eastsoutheastward subduction of the Eurasian plate underneath the Adriatic plate, whereas the Apennines are the accretionary prism of the westward subduction of the Adriatic plate.
The Apennines subduction retreated eastward generating the arc from the northern Apennines to Sicily, and continuing westward to Morocco along Maghrebides.The eastward migration of the arc was and still is accompanied by the backarc extension with Tyrrhenian basin opening, and has been interpreted as consequence of Apennines subduction rollback rather than relative convergence between Africa and Eurasia.The mechanism of eastward retreating of Adriatic slab has been ascribed to eastward mantle flow (Doglioni, 1991), hypothesis supported by seismological observations of mantle anisotropy (Lucente et al., 2006).The interaction between the northern Apennine retreat and the central-western Alps slab produces bending of lithosphere with consequent subsidence, mostly evident in the Po plain (Carminati et al., 2004).Both the subduction zones experience seismicity, which is more pronounced along the ridge of the Apennines, and in the foothills of the Alps.The Apennines subduction is characterised by compressional seismicity east of the chain, in the frontal accretionary prism, and extensional tectonics to the west, associated to the opening of the Tyrrhenian backarc basin and slab rollback (Pondrelli et al., 2006).
At present, eastern Alps are characterized by active N-S compression with shortening rates of about 2 mm/yr (Devoti et al., 2011;D'Agostino et al., 2005); the whole Apennine belt displays variable deformation styles, in fact the Adriatic slab retreat generates compression along the frontal Apennines accretionary prism, moving towards NE at a rate of 2-3 mm/yr, and extension along the Apennines chain (Serpelloni et al., 2005;Devoti et al., 2011;Bennett et al., 2012).All these mechanisms generate an inhomogeneous deformation belt, characterized by spatial gradients of velocity field, regardless the extensional or compressive tectonic style (Riguzzi et al., 2012).
The analysis of active stress field in Italy using borehole stress measurements, earthquake focal mechanisms and kinematic indicators on faults confirms these observations (Montone et al., 2004).

GNSS networks and data analysis
At INGV we currently archive and process daily data from 21 different Italian GNSS networks.The regional analysis is further augmented with other 50 sites belonging to EUREF and/or IGS networks that are homogeneously distributed in Europe and used for the ITRF2008 reference frame definition (figure 2).The analyzed networks, established for various purposes, using a variety of different monument types and operating in fairly different environmental conditions, cover all the Italian area.
Figure 3 shows all the analyzed GPS stations (black dots) and the color maps the distance from the nearest GPS site: in most sections this distance is between 10 and 20 km and in any case it is nowhere greater than 40 km.In this picture, northwest Italy is slightly sparser populated, since most of the INGV stations were planned to monitor the seismically active Apennines chain, nevertheless the number of new contributing stations is still increasing over time and we foresee a 30% increase of GPS stations in the next two years.The number of analyzed sites has grown during the years in a non-linear fashion, following the deployment of networks and the availability of GPS observations (figure 4).At present we analyze, on average, data from about 600 sites per day, this number is variable because of data gaps due to different causes (real missing data, "bad" data not passing a preliminary quality check).In figure 5 we report, for each network, the ratio between the average number of daily RINEX actually analyzed and the total number of expected daily RINEX, accounting for the network performance in terms of average data presence.We process data using the Bernese GNSS software v. 5.0 (Beutler et al., 2007), following the EUREF Guidelines for the European Permanent Network Analysis Centres (http://www.epncb.oma.be).
The GPS orbits and the Earth's orientation parameters are fixed to the combined IGS products and an apriori loose constraint of 10 m is assigned to all site coordinates.The elevation-dependent phase centre corrections and absolute phase centre calibrations are applied.The troposphere modeling consists in an a priori dry-Niell model fulfilled by the estimation of zenith delay corrections at 1-hour intervals at each site using the wet-Niell mapping function; in addition one horizontal gradient parameter per day at each site is estimated.The ionosphere is not modeled a priori, it is removed by applying the ionosphere-free linear combination of L1 and L2.The ambiguity resolution is based on the QIF baseline-wise analysis.The final network solution is solved with back-substituted ambiguities, if integer; otherwise ambiguities are considered as real valued measurement biases.
Thus the daily GPS solutions are not estimated in a given apriori reference frame but computed in a loosely constrained reference frame.The coordinates are randomly translated or rotated from day-to-day and their covariance matrices have large errors (on the order of meters) as a consequence of the loose constraints applied to the a priori parameters.
To express the coordinate time series in a unique reference frame and to compute the real covariance matrix, we perform two main transformations.First the loose covariance matrix is projected into a well-defined reference frame imposing tight internal constraints (at millimeter level), and then coordinates are transformed into the ITRF2008 by a 4-parameter Helmert transformation (translations and scale factor); the proper set of constraints is driven by the rank deficiency of the normal matrices, as discussed in Devoti et al. (2010).The Helmert transformation uses 45 sites located in central Europe as anchor stations for the regional reference frame realization.
Site velocities are estimated fitting simultaneously a linear drift, episodic offsets and annual sinusoids to all the coordinate time series.Offsets are estimated whenever a change in the GPS equipment induces a significant transient in the time series, whereas seasonal oscillations are accounted for by annual sinusoids.At this stage outliers are rejected whenever the weighted residual exceeds three times the global chi square (χ 2 ).Finally the common mode error signal is filtered out with a procedure similar to that adopted by Wdowinski et al. (1997).
After the common mode filtering the time series scatter decreases by 25-35%, for the overall analyzed network the final median weighted-root-mean-square of the residuals (WRMS) in the vertical and horizontal components are 4.3 and 1.4 mm respectively.In figure 6 we report the mean values of the WRMS for each network, separated for the vertical and horizontal components, accounting for the network performance in terms of time series quality.

Velocity and strain rate fields
The estimated horizontal velocity field with respect to a fixed Eurasian plate of the Italian area is shown in figure 7. The formal velocity errors obtained after the inversion have been rescaled site by site with the corresponding normalized variance factors.Each GPS site spans different life times, so that we show in the figure all those with a minimum observation time of 4 years.The reference Eurasian plate has been realized by minimizing the rigid motion of 15 selected EUREF stations located in stable central Europe.
Figure 8 shows the vertical velocity field with respect to the reference ellipsoid.The vertical velocities are interpolated with The Generic Mapping Tool (Wessel et al., 2013) using a continuous curvature surface gridding algorithm with an adjustable tension factor of T=0.25.The 2D strain rate tensor has been computed from the planar GPS velocities and the associated uncertainties by a distance-weighted approach using all stations on a regularly spaced grid down-weighting velocities with the function W=exp(−d 2 /α 2 ), where d is the distance between each node and the stations and α is the smoothing distance parameter, that depends on the spatial distribution of GPS sites (Shen et al., 1996;Shen et al., 2007).
Figure 9 shows the strain rate second invariant superimposed to the SRTM topography (Rodriguez et al., 2005); it has been computed on a regular grid of 0.1x0.1 degrees, applying a smoothing distance within the distance interval 25-350 km; the second invariant is the scalar accounting for all the 2D strain rate tensor components, thus representing the total amount of deformation rate.

Conclusions
The increased number of GNSS networks allows more and more detailed spatial and temporal resolutions of the ongoing crustal deformation, providing an unprecedented map of intra-plate kinematics of the Italian region.Since the first attempts to measure the convergence rate between Africa and Eurasia at large scale, based on space geodetic methods, see for instance Drewes (1993), current GNSS networks enable us now to study the deformation process at fault scale level (Riguzzi et al., 2013).
At present, the estimated horizontal and vertical velocity fields show interesting features: Apennines and Alpine belts are characterized by general uplift that follows the topographic ridge, whereas the Po Plain shows a gradually increasing subsidence from west to east.The Apennines belt displays significant extension (50-80 nanostrain/yr), while compressive tectonics characterizes northern Sicily, eastern Alps and the northeast front of the northern Apennines (25-50 nanostrain/yr).Second-order deformation patterns, on large scale wavelength (~100 km) have been detected on the accretionary prism of central and southern Apennines that are highly correlated with other geophysical data and related to deep rooted sections (70-100 km), marked by different subduction regimes.Apparently, at this scalelength the observed deformations are governed by the lithosphere as a whole.We interpret these deformations as a result of different subduction mechanisms, such as variations of the subduction rollback velocity affecting different segments of the subduction zone and/or to mantle flows in proximity of the slab edges (Devoti et al., 2011).

Fig. 1 .
Fig. 1.The interaction between the Africa and Eurasia plates generates a diffuse area of deformation and plate fragmentation in the Italian area where Alps and Apennines interact.

Fig. 3 .
Fig. 3. Distance from the nearest GPS site on a regular grid of 0.1x0.1degrees.Black dots are the stations.

Fig. 5 .
Fig. 5. Average percentage of daily RINEX per network.The percentage is given by the number of daily RINEX actually analyzed divided by the total number of daily RINEX expected.The performance of most networks (light-blue bars) is comparable with the one of EUREF/IGS sites (EUR, dark-blue bar).

Fig. 6 .
Fig. 6.Mean WRMS values for the vertical (blue bars) and horizontal components (red bars) of each network.All the networks, except CAT, have WRMS better than EUREF/IGS (EUR, darker bars).

Fig. 7 .
Fig. 7. Horizontal GPS velocity field with respect to Eurasia.Sites with a minimum time span of 4 years are shown.

Fig. 8 .
Fig. 8. Vertical velocity field with respect to the reference ellipsoid (uplift in red, subsidence in blue).Sites with a minimum time span of 4 years are shown.

Fig. 9 .
Fig. 9. Second invariant of the strain rate tensor on a regular grid of 0.1x0.1 degrees, maximum values (> 100 nanostrain/yr) are detected in volcanic areas and fewer along the Apennines.