A Combined Velocity Field of the Mediterranean Region

We present a full 3-D velocity field of the Earth’s surface in the Euro-Mediterranean area obtained from a combination of three different velocity solutions computed at the Centro Nazionale Terremoti (CNT) of the Istituto Nazionale di Geofisica e Vulcanologia (INGV). All the publicly available GPS data since 1993, have been fully reprocessed by three different software tools and the final velocity field is estimated combining three independent velocity solutions in a least squares sense. The input velocity solutions are treated as stochastic samples of the true velocity field by loosening the reference frame constraints in the associated variance-covariance matrix. The proposed approach allows for a fast and efficient combination of multi velocity solutions, taking into account the full network covariance, if available. The velocity map for the Euro-Mediterranean region will be updated and released regularly on the web portal of the National GPS Network (http://ring.gm.ingv.it) and made available to the scientific community. Here we show and discuss the data analysis and the combination schemes, and the results of the combined velocity field.


Introduction
For Earth scientists the Mediterranean region represents a unique natural laboratory to test and assess geodynamical theories, being the place where three major continental plates, Eurasia, Africa and Arabia interact in a very complex way, displaying a wide range of crustal deformation patterns.Here several microplates have been described from geological, seismological and geodetic data, among which the major ones are the Anatolian, Aegean, Apulian and Adriatic sub-plates.Thus from a scientific perspective, this region represents a key area for understanding the basic processes of plate tectonics and specifically the interplay of different tectonic styles in a continental collision area.The advent of space geodesy, especially exploiting the increasing number and density of Global Navigation Satellite System (GNSS) networks, allows us to provide an accurate measurement of the 3D Earth's surface motion, revealing the details of the kinematics and strain accumulation rates at different spatial scales.The EUREF permanent network (EPN), represents an important European infrastructure (http://www.epncb.oma.be) that operates and shares GNSS data over a continental scale.The EPN is a voluntary federation of self-funding agencies, research institutions and universities that maintain the Terrestrial Reference System in the European area, realizing and delivering fundamental geodetic products, such as Receiver INdependent EXchange (RINEX) data and position time series, used to build the global International Terrestrial Reference Frame (ITRF).On the other hand, the data provided on a local scale (national and regional) are currently not readily available to the scientific community.Several networks owned by private companies but also by local administrations do not contribute to the data disclosure policy that should sustain and motivate any scientific initiative.For this reason, the European Plate Observing System (EPOS; http://www.epos-ip.org)has initiated a long term project aiming at facilitating the use of integrated data and products in the field of geosciences.As in the past the Mediterranean was the locus of thorough cultural and economical exchanges between diverse peoples, we hope that it still represents a great opportunity for sharing knowledge and awareness.
At present-days, a few thousands of permanent GNSS stations provide unprecedented spatial and temporal coverage of the European deforming plate and its boundaries.This continuous monitoring effort, carried out by various public and private institutions is crucial to understand the large scale kinematics and to shed light on the physics that governs tectonic deformation and seismic, and aseismic faulting.The INGV is the primary Italian research center interested in the collection, management and analysis of GPS measurements.The INGV is archiving all the available GNSS data at national level.Three distinct analysis centers (AC) at INGV process and analyze routinely all the regional GPS data using respectively Bernese, GAMIT/ GLOBK and GIPSY-OASIS scientific software packages.They produce daily position solutions for up to 2000 stations located mostly along the African-Eurasian boundary and on the central and western European continent.
A large amount of work has been recently dedicated to study velocity fields from GNSS data at different geographic scales.In this study we are not willing to present an exhaustive list of such research effort, which is out of the scope of this work, but rather to propose a methodology to obtain rapid and reliable velocity fields on a wide continental scale.Recently, a number of authors published various deformation fields retrieved from GNSS data, on regional or even global scale (e.g.Caporali et al., 2009;Le Pichon and Kreemer, 2010;Pérouse et al., 2012;Nocquet, 2012;Kreemer et al., 2014, and references therein).Their efforts are partly directed to assemble all available geodetic solutions already published and to reduce systematic errors arising from different reference frame adoptions in order to obtain a velocity field as homogeneous as possible.These overviews of crustal deformation models over wide areas are of fundamental importance to test tectonic and geophysical models that govern the plate and plate boundary interplay.Here we present a cross-validated velocity field, based on a complete re-analysis of the whole GPS dataset, focusing on reference frame consistencies, homogeneous modeling over the whole data time-span and evaluating the repeatability of common stations.The purpose of this project is to generate a consistent, combined geodetic velocity model of the Mediterranean area on a regular basis, to offer high-quality geodetic products to a broad community of potential users.This action aims to increase data and products access to the scientific community and promote scientific studies on the deformation processes acting across the Mediterranean basin, but also informing engineers and public policy makers, who may use such results to plan for disaster mitigation and environmental monitoring.

GPS data collection and processing
Many GNSS permanent stations, managed by both scientific and commercial institutions, are available on the Eurasian plate and its boundaries.Although part of them are not specifically devoted for geophysical monitoring (cadastrial, topographic, etc.) and may potentially be of lower quality in terms of monumentation and data flow, their integration has a large potential to improve the resolution of the kinematic patterns of the area.The complete knowledge of metadata is usually not at the level of scientific devoted GNSS stations and needs a distinct analysis and cross-check.
In this study we find that at least 40% of the stations fall in this subsidiary category in which the metadata has to be carefully reconstructed.The data collection rate used in the analysis by all the AC is one sample every 30 seconds.Most of the stations are also streaming data at 1 Hz or even higher sampling rates in real time, but these streams are not processed in this study.
The INGV institute also manages the RING network in Italy (http://ring.gm.ingv.it), a GPS network of about 200 stations that meet specific research criteria based on geographic location and instrumental standardization (monumentation, receiver and antenna type).INGV collects and processes also GPS data from networks not belonging to the International GNSS Service (IGS) or the EUREF Permanent Network (EPN) nor in partnership with the international scientific community, thus expanding and supplementing the GPS database with those regional data whose historical records are often not preserved.Figure 1 shows the distribution of the GPS stations (updated at January 1, 2016) for which INGV collects and archives raw RINEX data.At present, depending on the data availability, more than 1600 European stations are regularly stored every day.These data have been fully re-processed by three AC, in which different approaches and analysis software are adopted (see also Avallone et al., 2010).Figure 2 shows the relative contribution to the combined velocity of the three AC distinguished by different color codes, blue for Bernese, green for GIPSY-OASIS and red for GAMIT/GLOBK processing chain.
The different processing schemes adopted by the three AC, are summarized in Table 1.Essentially two of them (Bernese and GAMIT/GLOBK) are based on the relative positioning concept using phase observables double differencing techniques, whereas the third processing scheme is based on the Precise Point Positioning (PPP) approach using the GISPY-OASIS software.Each AC estimates a variable number of station velocities, ranging from 900 to about 1500 values, with overlapping points.The availability of three independent solutions, in terms of daily stations positions and covariances, secular velocities, seasonal and transient signals, give us the possibility to internally compare and validate the results, with the main goal of assessing the repeatability of the independent velocity estimates.The combined velocity field has to be considered as a "consensus" (of cross-validated) geodetic product, that has been obtained after an iterative trial and error process and a final least squares combination of the overall velocity field.The validation procedure is also a key issue of the INGV effort in the EPOS project.

Details on the Bernese solution
The Bernese Analysis Center (BeAC) uses the Bernese software Ver.5.0 (Beutler et al., 2007), following the recommended guidelines for EUREF Analysis Centers (http://www.epncb.oma.be).Daily coordinate solutions of a network of stations are obtained by means of Ionosphere Free linear combinations of phase observables using the Quasi Ionosphere Free approach to properly solve phase ambiguities (Beutler et al., 2007).
For computational efficiency the full network is divided into sub-networks, each with about 50 or fewer stations.To allow the combination of the sub-networks into a full network daily solution, each sub-network contains a minimal of 11 tie stations.The troposphere modeling consists in the a-priori dry-Niell model fulfilled by the  estimation of zenith delay corrections at 1-hour intervals at each station using the wet-Niell mapping function (Niell, 1996).In addition, one horizontal gradient parameter per day at each site is estimated.Ocean loading is computed using the FES2004 tidal model coefficients (Lyard et al., 2006) provided by the Ocean Tide Loading web service (http://holt.oso.chalmers.se/loading).The GPS orbits and the Earth's orientation parameters are fixed to the combined IGS products (Dow et al., 2009) and an a-priori loose constraint of 10 m is assigned to all site coordinates.The IGS08 absolute antenna phase center correction has been applied to each station receiver antenna.The daily coordinates are thus estimated in a loosely constrained reference frame.In order to express the GPS time series in a unique reference frame, the daily solutions are first projected imposing tight internal constraints (at millimeter level), and then the coordinates are transformed into the IGS realization of the ITRF2008 frame (i.e., IGb08) by a 4-parameter Helmert transformation (translations and scale factor).The regional reference frame transformation uses 45 IGb08 anchor sites located in central Europe.To get rid of common translations of the entire network, the time series are readjusted through a common mode filtering procedure similar to that proposed by Wdowinski et al. (1997).Velocities at GPS stations are estimated by a linear weighted least squares fit of all the coordinate time series simultaneously, using the full daily covariance matrices and modeling secular drifts, episodic offsets and annual sinusoids (Devoti et al., 2014).

Details on the Gamit solution
The Gamit Analysis Center (GaAC) processes double-differenced ionosphere-free GPS carrier phase observations using the GAMIT/GLOBK (Ver.10.4) software, developed by the Massachusetts Institute of Technology (Herring et al., 2015).The estimated daily parameters are site positions and time-variable, piecewise, linear zenith and horizontal gradient tropospheric delay parameters, loosely a-priori constraints are applied to geodetic parameters, and the GPS orbits are fixed to the SOPAC products (http://sopac.ucsd.edu).The whole GPS network is divided into in smaller sub-networks, containing each less than 50 stations, sharing some high quality tie-sites.The ocean loading correction is applied using the FES2004 tidal model (Lyard et al., 2006).The Global Mapping Function (Boehm et al., 2006) is adopted to model both the dry and wet component of the tropospheric delay; pole tide correction is applied accounting for IERS data (pole.usno)(Petit and Luzum, 2010).The IGS08 absolute antenna phase center model for both satellite and ground-based antennas is used.Loosely constrained solutions, in the form of ASCII GA-MIT H-files, are later combined, using the ST_FILTER program of the QOCA software (Dong et al., 1998) with the IGS network solutions available from SOPAC.A global reference frame is realized by minimizing coordinates and velocities of the IGS global core stations, estimating a 7-parameter Helmert transformation (translations, rotations and scale factor) with respect to the IGb08 reference frame.GPS velocities are obtained by fitting a linear trend, annual and semi-annual terms and site specific offsets, assuming a white plus flicker (powerlaw) noise model (see also Serpelloni et al., 2013).The Common Mode Error in the time series is estimated using the Principal Component Analysis (PCA) strategy, following the method proposed by Dong et al. (2006).

Details on the Gipsy solution
The GIPSY-OASIS II software, Ver.6.2 (Zumberge et al., 1997), developed by the Jet Propulsion Laboratory ( JPL) is used to produce Precise Point Positioning (PPP) solutions using ionosphere-free carrier phase and pseudorange observables and using JPL's final fiducial-free GPS orbit products (Bertiger et al., 2010).Tropospheric wet zenith delay and horizontal gradients are estimated as stochastic random-walk parameters every 5 min using the Global Mapping Function (Boehm et al., 2006).Ocean loading is modeled using the FES2004 tidal model coefficients provided by the Ocean Tide Loading web service run by Chalmers University of Technology.The IGS08 absolute antenna phase center variations are used to model the azimuthal and elevation dependence.Station coordinates obtained in the loosely constrained frame of JPL fiducial-free GPS orbits are transformed into the IGS08 reference frame using daily 7-parameter transformations delivered by JPL.In order to reduce the common mode signal, we have specifically developed a terrestrial reference frame (called EU14) suitable for crustal deformation studies in and around that plate following the approach of Blewitt et al. (2013).This frame is defined by 6 Cartesian coordinates and velocities of each of 174 stations selected by specific quality criteria.The EU14 frame is aligned in origin and scale with IGb08.GPS velocities are obtained by fitting a linear trend plus annual and semi-annual terms and site specific offsets at position time-series, assuming a white plus flicker noise stochastic model.

Velocity combination
Prior to the velocity combination phase, all nuisance parameters needed to build the velocity solutions, such as local eccentricities, seasonal variations and position offsets, are solved consistently by each analysis group, in accordance with the station data quality and metadata accuracy.Thus, eventual inconsistencies in the definitions of reference markers and station position eccentricities do not concern the combination process.We consider each velocity field as a sample of the true velocity field while the combined velocity, as the best estimate of the true velocity field.The availability of different samples of the station velocities allows a sort of validation in which the velocity repeatability can be truly assessed.Another important advantage of combining the velocity field follows from the averaging effect, i.e. the combination of velocities obtained with almost independent procedures reduces to a minimum the chance of including biased velocities.
The velocity combination procedure is a generalization of the loosely constraints approach, first proposed by Heflin et al., (1992) and subsequently developed by Davies and Blewitt (2000).Blewitt (1998) evidenced an important property of the standard least squares theory, in which a functional model is fitted to a set of observations affected by a Gaussian noise.He demonstrates that refining the functional model by adding extra unknowns is equivalent, in the limit of unknown apriori information, to an augmentation of the stochastic model i.e. a redefinition of the noise process.In our framework, the functional model is a trivial identity between velocities, but in order to account for reference frame systematic errors, the functional model could be augmented through the estimation of additional reference frame biases.This augmentation is equivalent to enhancing the stochastic model through a loosening transformation of the covariance matrix (sensu Blewitt, 1998), so that it can be conside-red as the implicit estimation of additional parameters but without computing their full covariance.Loosely constrained solutions assign large errors to the implicit parameters and allow treating the differences between the observations as stochastic variables without the need to explicitly estimate them.In our context, loosening the reference frame constraints allows to save computation complexity, by waiving to Helmert transformation parameters and its covariances.
The three dimensional Cartesian velocities are combined in a least squares scheme and by treating the reference frame differences as a stochastic process.For this reason it is important to establish which parameters describe the actual reference frame transformation of our regional network.In global geodetic problems, the transformation of Cartesian coordinates of any point (X) between two terrestrial reference frames (TRF1 and TRF2) is expressed by the well-known Helmert transformation that, in its linearized form, reads where T is the translation vector, s is the scale factor and R is the infinitesimal rotation matrix.All these parameters, including the position coordinates, are time-varying quantities.Then the time derivative of the Helmert transformation is (2) It relates the velocities !X 1 and !X 2 expressed in two different TRFs in the most general way.Depending on the problem we are facing, the two velocity fields are affected by different biases s, R, !T , ! s and !R , some of which possibly negligible or even highly correlated (degenerate case).The selected GPS networks span a geographical region of approximately 30x50 degrees, which is less than 5% of the whole Earth surface, therefore the Mediterranean region could be rightly treated as a small-scale network.In this approximation some parameters may be considered fully correlated and consequently degenerate.For instance R and !R represent respectively the rotation of velocity vectors and the rotation rate of the reference frame (plate-like rotation), thereby their effects may locally be indistinguishable especially if the two rotation axes are orthogonal.The same remarks apply to the translation rate and scale rate factors ( !T and !s ), describing the addition of a constant velocity and a vertically directed velocity, respectively.Locally, the !T vector may mimic a velocity variation in any arbitrary direction.Hence in a confined small region, only three parameters are really independent and are adequate to describe the reference similarity transformation between two velocity fields In this approximation the reference frame transformation is independent from the vector positions, confining the whole combination process in the velocity phase space.This consequence is of fundamental importance, since it makes the velocity combination independent from the knowledge of the station position coordinates.Thus the unknowns that transform the velocity field into different reference frames are three translation rates ( !T ), three rotation angles (R) and one scale factor (s).A further simplification arises if all velocity solutions are obtained using the same datum (IGb08), and if the time series of the stations used to materialize the TRF are sufficiently long (in our case >15 years) so that common differences in the velocities, induced by center of mass motion, are expected to be negligible.In this case the similarity transformation reduces to a four parameters transformation where the rotation matrix (R) and the scale factor (s) describe most of the reference frame biases.
The functional model for the velocity combination is trivial, since each input velocity v s should match the combined velocity v c (4) where A c→s is the reorder matrix that translates the combined velocity vectors to the order of the solution velocity vectors and v is the noise vector.Since the velocity solutions are considered independent among each other, the least squares solution for the combined velocities may be recursively defined as follows where W S is the weight matrix (W s = C s −1 ) of the solution (S), i.e. the inverse of the solution covariance matrix (C s ) and the summation runs over all the solutions.
The combination process consists of two main steps: the stochastic model augmentation, in which rotations and scale uncertainties are increased (i.e.covariance loosening) adopting a diagonal (E) matrix with 10 arc-sec and 10 -4 respectively as assumed loosening parameters.The loosening constraints are in principle arbitrary and should be on the order of the expected systematic differences in order to allow the solutions to rotate and scale by the required amount.The covariance augmentation is provided by the external (E) a-priori covariance matrix that changes the solution covariance matrix (C s ) as follows (6) where the matrix B specifies which linear combinations will be relaxed.As an example, the loosening of a rigid rotation is modeled as usual by the rotation matrix (× Ι = × + Rθ ) , where θ is the vector of the 3 unknown angles and R contains the partials of the vector rotation transformation.In this case the B matrix in eq.6 takes the form of the rotation matrix R (B≡R).The resulting covariance matrix is termed as loosened covariance (Blewitt, 1998) and is associated to the corresponding (unchanged) velocity solution.The second step consists in the least squares estimation of the combined velocity field, where the observations are the velocity solutions with the associated loosened covariances together with a fourth IGS velocity solution, used to establish the ITRF frame.We choose the latest update to the IGS08 solution (ftp://igs-rf.ensg.eu/pub/IGb08/), called IGb08, which contains the best performing IGS stations and through their covariances contribute to realize the ITRF2008 frame.The 59 IGb08 common stations, located on the Eurasian, African and Arabian plates, conveniently define the ITRF2008 reference frame and act as fiducial "anchor" stations providing the datum constraints in the least squares problem.The combination is iterated twice in order to estimate the corresponding solution weighting factors, balancing mutual weights according to each solution chi-squared ( χ 2 ) (Devoti et al., 2010) (7) where Δ are the velocity residuals (v s _ v c ) and W s the weight matrix for each solution (S).
Finally we foresee the possibility of forcing two or more parameters to be estimated together (e.g.tying velocities together).This is achieved using the classical method of Lagrange multipliers (e.g.Arfken et al., 2013), where the least square problem is solved with the equality constraints.
The resulting velocity field includes 1729 stations in the Euro-Mediterranean area from west of the Straits of Gibraltar to east of the Levantine Sea (Fig. 2).The average (median) 1-sigma standard deviation for the combined velocities is respectively 0.3 mm/y for the horizontal components and 0.7 mm/y for the vertical component.The histograms of velocity residuals with respect to the combined solution, in the vertical, east and north components for each input solution, are shown in Figure 3.The weighted root mean square (WRMS) of the residuals ranges from 0.20 to 0.25 mm/y for the three input solutions (see Table 3), whereas the central tendencies (mean) are all within 0.1 mm in modulus except for the north com-ponent of the GAMIT and BERNESE solutions (see Figure 3).The distributions do not differ significantly, although the GAMIT and BERNESE solutions show a slight skewness especially in the northern component showing a mean of -0.12 and 0.17 mm/y respectively.These values, although below the repeatability of the combined velocity, may suggest a slight misalignment of the three reference frames that cannot be accom-  modated by a rigid rotation and scale transformation, the only biases allowed to change in the combination process.At this time we were not able to isolate the input solution (or solutions) that causes such tiny effect, but we will track down the problem in the follow-on combination.The 59 European IGS stations used in this combination show a very little discrepancy with the IGb08 velocities, an overall WRMS of 0.12 mm/y indicates a robust repeatability of these long-lasting stations.Figure 4 shows the histogram of the IGS velocity residuals in the three spatial components.The IGS residuals in the north direction show again a modest bias (0.09 mm/y) with respect to the combined solution.Nevertheless the residuals and biases in the combined velocity field are well below and consistent with the given standard deviations (0.3 and 0.7 mm/y in the horizontal and vertical components respectively).

Technical issues and adopted conventions
A major problem in combining independent velocity solutions is the recognition of the station identity, since no naming convention is guaranteed in an open environment.Worldwide permanent GNSS stations are generally identified by a 4-character ID and, eventually a IERS domes number (9 characters) univocally assigned by the Institut Géographique National (IGN), acting as a central authority to avoid overlays.This procedure has become standard in the IGS community.The GNSS stations must strictly comply to the IGS requirements and should be registered in the IGN data-base which, in turn, is not envisaged as a compelling procedure outside the IGS service.In Italy more than 50 networks contribute to the INGV data archive and are generally not involved in IGS activities.Therefore a thorough regulation of ID uniqueness is difficult to maintain Instead of forcing an a-priori naming convention at the database or data archive level, we decide to adopt an a-posteriori approach based on the assignment of a unique label based on the station positions (i.e.geo-coding).In particular, we choose the GHAM code proposed by Agnew, (2005), to label each GPS station unambiguously.The GHAM code is composed of alternating letters and numbers, providing tags to geographic locations and defining addresses of equal-area cells with arbitrary precision.We choose a 12-character code that corresponds to a cell size of 1.9 m (square root of area), which is sufficiently small to identify a single GNSS antenna installation.The main advantage of this technique is the possibility to automate the site recognition process reducing the amount of knowledge to a minimum (only 12 characters).The geo-code has also an interesting hierarchical sorting property, in that alphabetically sorted codes group stations that would be nearby in space.
For simplicity and traceability of the GNSS stations in the combined solution, all the input 4-character IDs are preserved so that each station can be univocally identified by both the geo-code and the 4-character ID.Thus multiple velocities, referring to the same station (same geo-code) but different 4-character ID, will be estimated as tied velocities and appear independently in the combined solution but having the same velocities and covariances.Forcing two or more velocities to be estimated together (tied velocities) is achieved using the classical method of Lagrange multipliers (e.g.Arfken et al., 2013), where the least squares are solved with the constraints of having equal velocities.
The proposed combination approach provides a computationally efficient algorithm to combine a large number of station velocities in a region of limited size like the Euro-Mediterranean area.It completely neglects the station positions, which may be known only approximately (~1 m) and is able to combine velocities expressed also in different ITRF reference frames sin- ce the systematic differences are treated as stochastic variables Equation ( 6).The time series discontinuities, seasonal variations and local antenna eccentricities are treated in the earlier processing stage and could be sol-ved for independently by each processing center.This may lose some modeling parameters, such as seasonal variations and position offsets, since each time series is processed independently but the main advantage is  a very quick and rather simple combination of GNSS networks addressed to assimilate a large number of a very quick and rather simple combination of GNSS networks addressed to assimilate a large number of crustal velocities in a common reference frame.The method uses the complete input covariances, if available, and definitely provides the complete covariance matrix of the combined velocity field.

Results: Europe-Africa boundary zone deformation
The combined horizontal velocity solution, estimated with respect to the stable Eurasian plate (Figure 5) and the vertical rates in the IGb08 frame (Figure 6), highlight with unprecedented details the 3D kinematics of a large portion of the Euro-Mediterranean region, with dense spatial sampling of crustal deformation across the Mediterranean plate boundary and the most important active fault systems.Although some of the station velocities have been already published elsewhere, and the overall surface kinematics is well known and discussed in several recent papers (e.g.Nocquet, 2012;Serpelloni et al., 2013;Kreemer et al., 2014;Métois et al., 2015, Palano et al., 2015), here the velocity field is represented for the first time at the Eurasian plate scale with homogeneous standards, best available spatial resolution and special care for reference frame stability.Moreover the velocity field is obtained by re-processing the whole data set with different    approaches and combining the velocity fields with the aim of building on a regular basis a high level GPS pro-duct for the scientific community.The combination ensures that solution specific biases, eventually indu- ced by the actual analysis procedure, are averaged out and possibly minimized, in the sense that for a combined solution the chance of having biased velocities is minimized with respect to a single solution.From a user point of view, the availability of standardized regional solutions favors multi-disciplinary access to high-level geodetic products enabling scientific stakeholders unfamiliar with geodetic techniques to benefit of the updated maps of crustal motion in their studies.In particular the wealth of stations allows to accurately map the accumulation of strain and crustal motion over a ~4000 km long plate boundary.The kinematic boundary conditions in the Mediterranean are represented by the ~N-S Africa-Eurasia convergence.The NE-ward to N-ward motion of the Adriatic domain and the W-ward to SW-ward motions of the Anatolian and Aegean plates highlight the kinematics of the most important microplates of the Mediterranean area (McKenzie, 1970).Plates and microplates relative motions are accommodated across the major plate boundary zones, well mapped by the release of seismicity.Deviation from the major plate motions in the Mediterranean implies that additional geodynamic processes are required to explain the observed velocity field.The goal of this work is not that of discussing the tectonic and geodynamic implications of the newly proposed velocity field.In the following we present regional and more detailed velocity maps, including some velocity cross-sections through the major deformation belts of the Mediterranean area.
In the western Mediterranean (Figure 5 -6) the transition from oceanic-oceanic to continent-continent boundary shows a gradual widening of the deformation zone.Around the Alboran Sea, in southern Spain, the western Betics and in northwestern Morocco, the Moroccan Rif show a wide deformation pattern that cannot reflect a simple plate boundary interaction (see e.g.Monna et al., 2014;Chalouan et al., 2014).North Iberia and across the Pyrenees including the Balearic islands show no detectable motion with respect to Eurasia plate, nor significant vertical deformation (Figure 6).North of Sicily (Figure 8), most of the African motion is absorbed offshore in the southern Tyrrhenian thrust system.The Apennines show extension with a strain rate axis at ~90° from the plate convergence vector.The Eastern Alps and the Dinarides accommodate N-NE motion of the Adriatic Sea relative to Eurasia.The Alps and Apennines show present uplift at different rates (1-3 mm/y) whereas subsidence is dominating in the Aeolian islands, NE of Sicily and along the Po plain, in northern Italy (Figure 6).The velocity projections across the Italian peninsula (Figure 9) clearly show how the new combined velocity field well samples the velocity gradients across the major seismically active belts.Profile A-B shows the NS compression in the Southeastern Alps (A-B) and the extensional patterns across the Apennine chain from north (C-D) across the Ferrara Arc facing the Po plain.In particular, this profile shows the crustal extension at a rate of about 3 mm/y across the Apennine belt and the compression at about 2 mm/y towards the Adriatic foreland.Other profiles (E-F, G-H, I-L, M-N) show different rates of extension at 3-4 mm/y along the Apennines from north to south.The profile (O-P) that crosses Sicily from south to north shows a significant extension in correspondence of the Peloritani Mountains, NE of Sicily, and compression close to the northern edge in the Aeolian Islands.
In the eastern Mediterranean the subduction along the Hellenic Arc dominates the tectonic deformation of the whole area (Papazachos, 1988;Wortel et al., 1990).The region between the Dinaric and the Balkan mountains, shows an increasing velocity from the Pannonian basin to the southernmost Macedonia (Figure 10).The profile A-B along the SSE direction shows the progressive increasing of the displacement rate up to 4-6 mm/y in the Halkidiki peninsula.The C-D profile, SSW oriented, depicts the velocity increase across the Isthmus of Corinth and the Peloponnese gaining a rate of about 40 mm/y near the Hellenic trench.
The large number of available GPS stations scattered in different regions of the investigated area, meets spatial densities that allows detailed estimates of the strain rate field, providing important information to improve probabilistic seismic hazard models (e.g.Bird et al., 2015).Under appropriate assumptions that the strain rate can be converted into an estimate of the rate at which strain energy is accumulating, it becomes possible to identify areas where relative changes of strain may correspond to points where energy will be possibly released in future seismic events.

Figure 1 .
Figure 1.History of the number of GNSS stations archived at INGV.The yellow, red and green lines show the evolution of sites number for the Euref, RING and ItalPos networks, respectively.The blue box in the inset show the area for which we provide the combined solution in this work.

Figure 2 .
Figure 2. GPS network arrangement of the combined velocity solution.The following color code of the bullets has been used to highlight the contribution of the three AC: blue (Bernese), green (Gipsy) and red (Gamit).

AFigure 3 .
Figure 3. Histograms of velocity residuals with respect to the combined solution.Top, middle and lower rows are respectively the vertical, east and north components, columns from left to right refer to Gipsy, Gamit and Bernese residuals.

Figure 4 .
Figure 4. Histograms of velocity residuals of only the IGS stations with respect to the combined solution.Top, middle and lower panels show respectively the vertical, east and north components.Mean and standard deviations are indicated in units of mm/y.

Figure 5 .
Figure 5. Map of the horizontal GPS velocities (in the Eurasian-fixed frame) from the combined solution.Error ellipses are not shown here for clarity of the figure.

Figure 6 .
Figure 6.Map of the absolute (i.e.IGb08 frame) vertical GPS velocities of the combined solution.

Figure 7 .
Figure 7. Left: combined Eurasian-fixed horizontal velocities (with 95% error ellipses) for the Iberian region.Right: profile parallel and profile normal velocity components (with error bars showing the 2σ uncertainties) projected along a N10°W cross section from northern Africa to northern Iberia (see the dashed black box in the left panel).The dark grey areas show the average (median) topography in the profile swath, with the light grey and white areas showing the maximum and minimum elevations, respectively.

Figure 9 .
Figure 9. Selected velocity profiles (parallel projections) of the Italian peninsula.A-B profile across the Southeastern Alps, from C-D to P-O profiles along the Apennine chain, from northern Apennines to eastern Sicily.The dark grey areas show the average (median) topography in the profile swath, with the light grey and white areas showing the maximum and minimum elevations, respectively.

Figure 10 .
Figure 10.Left: combined Eurasian-fixed horizontal velocities (with 95% error ellipses) for the Balkans and Aegean region.Right: profile parallel velocity components (with error bars showing the 2σ uncertainties) projected along a N150°E (A-B) and N210E° (C-D) cross sections.The dark grey areas show the average (median) topography in the profile swath, with the light grey and white areas showing the maximum and minimum elevations, respectively.
Figure 7 displays the velocity projections along a N10°W profile across the Strait of Gibraltar and the western Iberian region.The profile parallel components from south to north show a gradual decrease of the velocity projections from 1.5-2 mm/y to ~0 at north of the Betics (about 200 km inland from the Strait) indicating active shortening tectonic process.The profile normal projections show the right-lateral shear between the African plate and the Iberian peninsula with about 4 mm/y accommodated across the Strait of Gibraltar and the Betics, with the velocity components reaching zero at about the same 200km distance.The horizontal velocities also show a clockwise rotation, especially visible in southern Iberia crossing the Strait of Gibraltar from east to west.