Present-day surface deformation of the Alpine region inferred from geodetic techniques

We provide a present-day surface-kinematics model for the Alpine region and surroundings based on a high-level data analysis of about 300 geodetic stations continuously operating over more than 12 years. This model includes a deformation model, a continuous surface-kinematic (velocity) field, and a strain field consistently assessed for the entire Alpine mountain belt. Special care is given to the use of the newest Global Navigation Satellite Systems (GNSS) processing standards to determine high-precision 3-D station coordinates. The coordinate solution refers to the reference frame IGb08, epoch 2010.0. The mean precision of the station positions at the reference epoch is ±1.1 mm in N and E and ±2.3 mm in height. The mean precision of the station velocities is ±0.2 mm a⁻¹ in N and E and ±0.4 mm a⁻¹ in height. The deformation model is derived from the point-wise station velocities using a geodetic least-squares collocation (LSC) approach with empirically determined covariance functions. According to our results, no significant horizontal deformation is detected in the Western Alps, while across the Southern and Eastern Alps the deformation vectors describe a progressive eastward rotation towards Pannonia. This kinematic pattern also makes evident an increasing magnitude of the deformation from 0.1 mm a⁻¹ in the western part of Switzerland up to about 1.3 mm a⁻¹ in the Austrian Alps. The largest shortening is observed along the southern front of the Eastern Alps (in the northern area of the Venetian-Friuli Basin) and in the northern part of the Apennine Peninsula, where rates reach 2 and 3 mm a⁻¹, respectively. The average accuracy of the horizontal deformation model is ±0.2 mm a⁻¹. Regarding the vertical kinematics, our results clearly show an ongoing average uplift rate of 1.8 mm a⁻¹ of the entire mountain chain, with the exception of the southern part of the Western Alps, where no significant uplift (less than 0.5 mm a⁻¹) is detected. The fastest uplift rates (more than 2 mm a⁻¹) occur in the central area of the Western Alps, in the Swiss Alps, and in the Southern Alps in the boundary region between Switzerland, Austria, and Italy. The general uplift observed across the Alpine mountain chain decreases towards the outer regions to stable values between 0.0 and 0.5 mm a⁻¹ and, in some cases, to subsidence like in the Liguro-Provençal and Vienna basins, where vertical rates of −0.8 and −0.3 mm a⁻¹ are observed, respectively. In the surrounding region, three regional subsidence regimes are identified: the Rhône-Bresse Graben with −0.8 mm a⁻¹, the Rhine Graben with −1.3 mm a⁻¹, and the Venetian-Friuli Basin with −1.5 mm a⁻¹. The estimated uncertainty of our vertical motion model across the Alpine mountain belt is about ±0.3 mm a⁻¹. The strain field inferred from the deformation model shows two main contrasting strain regimes: (i) shortening across the south-eastern front of the Alps and the northern part of the Dinarides and (ii) extension in the Apennines. The pattern of the principal strain axes indicates that the compression directions are more or less perpendicular to the thrust belt fronts, reaching maximum values of $\mathrm{20}\times {\mathrm{10}}^{-\mathrm{9}}$ a⁻¹ in the Venetian-Friuli and Po basins. Across the Alpine mountain belt, we observe a slight dilatation regime in the Western Alps, which smoothly changes to a contraction regime in western Austria and southern Germany, reaching maximum shortening values of $\mathrm{6}\times {\mathrm{10}}^{-\mathrm{9}}$ a⁻¹ in north-eastern Austria. The numerical results of this study are available at https://doi.pangaea.de/10.1594/PANGAEA.886889.

1 Introduction

The Alpine orogeny is the result of the continental collision of the African and European plates, more precisely of the push of the African Plate through the Adriatic promontory (also called the Adria or Apulian microplate; Channell et al., 1979) against Central Europe. A wide range of multidisciplinary studies based on geologic, geophysical, and geodetic data have demonstrated that the Alpine orogeny has been active from the Jurassic until the present (e.g.; Dewey et al., 1973, 1989; Channell and Horváth, 1976; Dercourt et al., 1986; Le Pichon et al., 1988; Mueller and Kahle, 1993; Ustaszewski et al., 2008; Handy et al., 2010, 2015). The tectonic evolution and present geodynamic framework have been primarily deduced from seismotectonic syntheses, and it is clear that the Alpine area and, in a more extended sense, the Mediterranean represent a variety of lithosphere blocks of different age, thickness, and rheology, leading to complex kinematic processes, including subduction, back-arc spreading, rifting and thrust, reverse and strike-slip faulting. These processes are superimposed on those phenomena associated with any orogeny; i.e. uplift, deformation, erosion, metamorphism, and foredeep basin changes. During the last three decades, the Global Navigation Satellite Systems (GNSS), like GPS (Global Positioning System) and GLONASS (GLObalnaja NAwigazionnaja Sputnikowaja Sistema), became a fundamental tool to observe, model, and understand present-day kinematic processes, as they provide precise geodetic constraints with a high spatial resolution allowing the improvement of geodynamic models, especially at diffuse plate boundaries like the African–Eurasian boundary.

The use of GPS to detect tectonic surface deformations in Europe started already in the 1990s (e.g. Fejes et al., 1993; Kahle et al., 1994; Chéry et al., 1995; Kaniuth et al., 1995, 1999; van Mierlo et al., 1996; Calais et al., 1998, 2000; Ferhat et al., 1998; Calais, 1999; Grenerczy et al., 2000; Sue et al., 2000). Given the complexity of the deformation zone in the southern part of the continent, scientists concentrated on selected areas to deploy observing stations and to develop local deformation models. Recent studies can be classified following the geographic subdivision of the Alps into the Western, Central, Eastern, and Southern Alps (Fig. 1):

• The Western Alps in south-eastern France, with the boundary region France–Italy and south-western Switzerland; e.g. Tesauro et al. (2006), Delacou et al. (2008), Larroque et al. (2009), Nguyen et al. (2016), and Nocquet et al. (2016).

• The Central and Eastern Alps in Austria and the adjacent regions of Switzerland, Liechtenstein, Germany, Italy, and Slovenia; e.g. Haslinger et al. (2006), Weber et al. (2006), Gosar et al. (2007), Sue et al. (2007), Brückl et al. (2010), and Brockmann et al. (2012).

• The junction zone of the Eastern Alps, the Dinarides, and the Pannonian Basin; e.g. Bada et al. (2007), Bus et al. (2009), Caporali et al. (2009), and Möller et al. (2011).

• the Southern Alps in northern Italy (also considered as a boundary zone for the computation of deformation models in the Adriatic region); e.g. Grenerczy et al. (2005), Grenerczy and Kenyeres (2006), D'Agostino et al. (2005, 2011), Devoti et al. (2008, 2011), Serpelloni et al. (2005, 2013), Cuffaro et al. (2010), and Métois et al. (2015).

These studies describe in general deformation models of regional scope and cover isolated segments of the Alpine mountain chain and its forelands. The objective of this work is the consistent determination of a continuous horizontal and vertical surface-kinematic field of the Alpine area that provides an integral picture of the ongoing deformation processes in the entire region. A network of about 300 continuously operating GNSS (CO-GNSS) stations with observations collected over 12.4 years is used for the precise determination of station positions and velocities. Based on these results, a continuous kinematic field is derived using a geodetic least-squares collocation (LSC) approach with empirically determined covariance functions. Main results are a deformation model, a continuous surface-kinematic (velocity) field, and a strain field consistently assessed for the entire Alpine mountain belt. The core contribution of this work is the homogeneous analysis of a large number of freely available data from CO-GNSS stations located in the Alpine region, using unified processing standards and a common reference frame for the complete time span covered by the observations. Special care is given to the application of the newest GNSS processing standards like absolute corrections to the GNSS antenna phase centre variations, troposphere delay estimation based on numerical weather models, and atmospheric and oceanic tide loading effects. To ensure consistency, we also use reprocessed GNSS satellite orbits referring to the IGS08/IGb08 reference frame (Rebischung et al., 2012) and including the most recent models adopted by the International GNSS Service (IGS) for orbit determination, such as the albedo model of Rodríguez-Solano et al. (2012), the satellite antenna thrust effect, and the a priori solar radiation pressure model switch-off (Steigenberger et al., 2014). The IGS08/IGb08 reference frame corresponds to a subset of GNSS stations of the ITRF2008 (International Terrestrial Reference Frame 2008; Altamimi et al., 2011), which are used by the IGS as fiducial points for the computation of the satellite ephemeris and the Earth orientation parameters (EOPs).

Figure 1 Simplified tectonic map of the Alps and surrounding areas showing the major tectonic features and faults, modified from Schmid et al. (2004) and Bousquet et al. (2012). Geographic subdivision of the Alps: WA – Western Alps, CA – Central Alps, EA – Eastern Alps, SA – Southern Alps. Tectonic plates/blocks: EU – Eurasia, ADR – Adria, LIG – Liguria, PAN – Pannonia, CS – Corso-Sardinia. Basins: LP – Liguro-Provençal, TS – Tyrrhenian, PO – Po, VF – Venetian-Friuli, MB – Molasse, PB – Pannonian, SB – Styrian, VB – Vienna. Graben: RBG – Rhône-Bresse, RHG – Rhine. Mountains and mountain belts: AP – Apennines, PY – Pyrenees, DI – Dinarides–Hellenides, JM – Jura, MC – Massif Central. Label pr identifies the Periadriatic Line. Topography and bathymetry after the ETOPO1 model (Amante and Eakins, 2009).

In this context, after a brief description of the geological and tectonic framework of the Alpine region, the third section of this paper provides details about the geographical distribution of the GNSS stations, the availability of the GNSS observations, and the data sources. Afterwards, we present in Sect. 4 the GNSS data processing strategy, the geodetic datum realisation, the analysis of the position time series, and the computation of a cumulative (multi-year) solution for the estimation of precise station positions and velocities. Section 5 concentrates on the computation of the kinematic model: correlations between the station velocities and the tectonic settings are analysed, the modelling methodology is presented, and the vertical and horizontal deformation models are briefly described. Finally, the strain field computation is summarised in Sect. 6. The reliability of the results obtained in this study is validated by contrasting them with the deformation patterns and the displacement rates published in the existing literature. With this work, we are providing a present-day surface-kinematic model for the Alpine region based on a high-level GNSS data processing.

Figure 2 CO-GNSS stations processed for the estimation of the surface kinematics in the Alps. The stations of the contributing CO-GNSS arrays are colour coded. Topography and bathymetry after the ETOPO1 model (Amante and Eakins, 2009).

2 Geological and tectonic framework

We focus on the estimation of the present-day kinematics in the Alps and surrounding areas based on GNSS measurements gained during the last decade. Previous processes associated with the tectonic evolution and their effects are usually deduced from geophysical modelling, but they are beyond the scope of this work. The main tectonic blocks underlying the Alps and their forelands are the Eurasian Plate in the north and the west, the Adriatic and Ligurian microplates in the south, and the Pannonian fragment in the east (Fig. 1). The boundaries between these units are mainly inferred from (natural and controlled) seismic data and tomographic inversions; see details in Pfiffner et al. (1990), Waldhauser et al. (1998), Doglioni and Carminati (2002), Lippitsch et al. (2003), Kissling et al. (2006), and Brückl et al. (2007). The push of the African Plate against the Eurasian Plate causes an underthrusting motion beneath the Southern Alps as indicated, for instance, by the seismicity in the Friuli area. In addition, the Adriatic promontory rotates anticlockwise, indenting the Eastern Alps and Dinarides in the north (e.g. Nocquet and Calais, 2004; Handy et al., 2015) and subducting beneath the Hellenides in the south (Grenerczy et al., 2005; Bennett et al., 2008). The polarity of subduction changes from one area to the other: while the European lithosphere is subducted beneath the Alps, the African lithosphere is subducted beneath the Apennines and the Dinarides–Hellenides (Argnani, 2009). In this way, Alps, Apennines, and Dinarides–Hellenides are still active orogens, although at different rates (D'Agostino et al., 2005; Devoti et al., 2008; Cuffaro et al., 2010; Cheloni et al., 2014). Each subduction is associated with the usual vertical motions; i.e. subsidence in the foreland basin and uplift in the mountain belt. The continuous curved orogenic belt built by the Alps together with the Carpathians and the Dinarides–Hellenides bifurcates and encircles the Pannonian Basin system (Doglioni, 1987; Royden, 1993; Doglioni and Carminate, 2002; Horváth et al., 2006). South of the Alps are the Po and the Venetian-Friuli basins. As mentioned, these three basins – the Pannonian, Po, and Venetian-Friuli basins – are still active foreland basins; i.e. they continue subsiding and receiving sediments (e.g. Cuffaro et al., 2010; Serpelloni et al., 2006, 2013). The Po Basin is also part of the northern foreland basin of the Apennines (e.g. Laubscher, 1988; Argnani, 2009). In the west of the Apennines and south of the Western Alps, we find the Liguro-Provençal and the Tyrrhenian basins. They are separated by the Corso-Sardinia block (Serpelloni et al., 2006, 2013). The northern foreland basin of the Alps is the Molasse Basin (Ziegler, 1990; Kuhlemann and Kempf, 2002). It runs between the Jura Mountains and the Alps in the north-western part of Switzerland and continues along the Bavarian and Austrian Alps up to Vienna. In contrast to the Pannonian and Po basins, it is not active as a foreland basin; on the contrary, large parts of it are being uplifted and eroded (Becker, 1999). On the west, the Alps are boarded by the Rhône-Bresse Graben. In general, it can be stated that the deformation in the Alpine region is dominated by the anticlockwise motion of the Adria microplate, causing compression in the Eastern Alps, dextral shear in the Central Alps and a very slow deformation in the Western Alps.

3 Distribution of the CO-GNSS stations

The determination of the surface kinematics in the Alpine region is based on the observations collected at 306 CO-GNSS stations (Fig. 2). Stations located at distances up to some thousand kilometres from the Alps are included to improve the geometry of the network and to ensure a reliable geodetic datum realisation. Several sites are part of the global IGS network (Dow et al., 2009), the EUREF Permanent Network (EPN) (Bruyninx et al., 2012, 2013), and the Integriertes Geodätisches Referenznetz Deutschlands (GREF). EUREF stands for the Reference Frame Sub-Commission for Europe of the International Association of Geodesy (IAG). The CO-GNSS station distribution along the Alpine mountain belt is densified with more than 30 GAIN (Geodetic Alpine Integrated Network) stations. Five of these stations were installed and are maintained by DGFI-TUM (Seitz et al., 2014). The GAIN network was established in the frame of the ALPS-GPSQUAKENET project, which was funded by the European Union under the INTERREG IIIB Alpine Space Programme from 2004 to 2007 (Aoudia et al., 2007). Additional data from CO-GNSS sites in Austria were provided by the Space Research Institute (Institut für Weltraumforschung, Graz) of the Austrian Academy of Sciences (Haslinger and Strangle, 2006). Another valuable source of data for this study is the Friuli Regional Deformation Network (FreDNet) established in 2002 (Zuliani et al., 2002; Rossi et al., 2013). While FreDNet densifies the station distribution in the Eastern Alps, the Western Alps are very well covered by CO-GNSS sites of the RENAG (Réseau National GPS, France, Nocquet et al., 2016) and Orpheon networks. Orpheon is a geodetic network providing real time and augmentation services in France. The CO-GNSS network processed in this study does not include stations in Switzerland. However, to get a homogeneous data coverage in the Central Alps, we use IGb08-based station velocities of the Automated GNSS Network for Switzerland (AGNES) provided by Brockmann et al. (2012). Figure 2 shows the station distribution according to the data sources.

The GNSS data cover the time span from 1 January 2004 to 30 May 2016. However, the amount of the available GNSS observations varies over the years (Fig. 3). The period with the maximum number of active stations within the network is between 2012 and 2015, reaching nearly 180 CO-GNSS stations every day. The total number of the daily observation files processed in this study is about 0.57 million.

Figure 3 Amount of daily observation files processed in this study.

4 Analysis of GNSS data and determination of precise station positions and velocities

Effects of ongoing geodynamic processes may be inferred from crustal movements detected by the analysis of repeated (or continuous) geodetic observations. The basic idea is to determine (and compare) the geometry of the network formed by the geodetic stations at certain epochs. The changes in the network geometry over large time spans (e.g. years) are interpreted as a deformation of the network caused by geodynamic processes underlying the studied area. To ensure a high reliability in the determination of the network geometry, the geodetic measurements must be processed over the entire time span following strict standards and procedures. In this study, we therefore apply the same analysis strategy usually applied for the establishment of GNSS-based regional geodetic reference frames like SIRGAS (Sistema de Referencia Geocéntrico para las Américas, e.g. Sánchez et al., 2013; Sánchez and Drewes, 2016) or the EUREF Permanent Network (EPN) as a realisation of the European Terrestrial Reference System (ETRS); e.g. Völksen et al. (2009a, b). The data analysis is based on the conventions outlined by the International Earth Rotation and Reference Systems Service (IERS) for the determination of the ITRF (Petit and Luzum, 2010) and the GNSS-specific guidelines defined by the IGS for the precise processing of the global GNSS network (IGS, 2014).

In a first processing step, free weekly normal equations (NEQ) are computed. “Free” means in this context that the geometry of the GNSS network is fully consistent with the GNSS orbits, but the network's origin and orientation are not constrained yet. They will be constrained in an a posterior processing step, when the geodetic datum for the GNSS network is introduced. The free weekly NEQ covering the complete time span from January 2004 to May 2016 are combined to a multi-year solution with constant velocity for each station. This includes an iterative time series analysis to detect outliers and discontinuities in the weekly station positions. In the multi-year solution, the geodetic datum is realised with respect to the coordinates (positions and velocities) of selected reference stations. The station velocities represent the mean displacement per year of the CO-GNSS stations and they are the primary input data to model the surface kinematics. The following sections describe the key aspects of the GNSS data analysis performed in this study.

4.1 GNSS network weekly solutions

Figure 4 Position time series of the station BOLG located in Bologna, Italy. Discontinuities in 2005 and 2009 (indicated by blue vertical lines) are caused by the poor monument construction of the site. Freezing water in the outer tube incorporating an inner tube carrying the antenna raised the inner tube like a piston. In 2009 the inner tube was put in its original place (Bruni et al., 2014). The blue line in 2012 represents the discontinuity caused by seismic events that occurred on 20 May (M_w 6.0) and 29 May (M_w 5.8) 2012 in northern Italy (11.230^∘ E, 44.890^∘ N and 11.086^∘ E, 44.851^∘ N, respectively). Piecewise sinusoidal lines in light green represent a functional model approximating the seasonal motions detected at the station.

The processing of the daily GNSS observations is accomplished using the double-difference baseline approach and the least-squares adjustment implemented in the Bernese GNSS Software V5.2 (Dach et al., 2015). The main processing characteristics are as follows:

• GPS and GLONASS observations at a 30 s sampling rate are considered.

• The basic observation is the L3 ionosphere-free linear combination.

• The satellite orbits, satellite clock offsets and EOPs are fixed to the values generated by the IGS processing centre CODE (Center for Orbit Determination in Europe, Dach et al., 2017; Steigenberger et al., 2014; http://www.aiub.unibe.ch/download/REPRO_2013, last access: 20 August 2018). To ensure consistency, reprocessed (2004–2013) and final CODE (2014–1016) products based on the reference frame IGS08/IGb08 are used for the analysis.

• The antenna phase centre offsets and direction-dependent phase centre corrections for both satellite transmitting and ground receiving antennas are taken from the model igs08.atx (Schmid et al., 2016; ftp://www.igs.org/pub/station/general/pcv_archive/, last access: 20 August 2018).

• The tropospheric zenith delay is modelled using the global mapping function (GMF, Böhm et al., 2006). Here, the a priori zenith hydrostatic delay (ZPD) values are derived from gridded coefficients based on the global pressure temperature (GPT) model (Böhm et al., 2006) and are refined by computing zenith wet delays (ZWD) in a 1 h interval (also using GMF). In addition to model azimuthal asymmetries, horizontal gradient parameters are estimated using the model described by Chen and Herring (1997).

• Corrections for the solid Earth tide, permanent tide, and solid Earth pole tide are applied as described in Petit and Luzum (2010). The ocean tide loading is estimated with the FES2004 model (Letellier, 2004) and the atmospheric tide loading caused by the semidiurnal constituents S1 and S2 is estimated following the model of van Dam and Ray (2010). The coefficients for the ocean tide loading are provided by Machiel Simon Bos and Hans-Georg Scherneck at http://holt.oso.chalmers.se/loading/ (last access: 20 August 2018). The coefficients for the atmospheric tide loading are provided by Tonie van Dam at http://geophy.uni.lu/ggfc-atmosphere/tide-loading-calculator.html (last access: 20 August 2018).

• Non-tidal loadings such as atmospheric pressure, ocean bottom pressure, or surface hydrology are not reduced.

• The seven daily free normal equations corresponding to a week are combined to generate a free normal equation per week. In this combination, stations with large residuals in any daily normal equation (more than ±10 mm in the horizontal or more than ±15 mm in the vertical) are removed from the corresponding daily equation and the weekly combination is recomputed.

Figure 5 Position time series of the station KARL located in Karlsruhe, Germany. Annual signals are well observed in the north and east components. Equipment changes (indicated by the blue vertical lines) do not seem to affect the time series.

4.2 Determination of precise station positions and velocities

The determination of the station coordinates is performed within a multi-year solution. In a first step, all free weekly NEQ available from January 2004 to May 2016 are combined and solved using the programme ADDNEQ2 of the Bernese GNSS Software V5.2 (Dach et al., 2015). Based on this preliminary solution, in a second step, residual time series for each station are generated to identify outliers or discontinuities that may mislead the computation of the station velocities. In this study, we use as thresholds ±10 mm for the north and east components and ±20 mm for height (about 4-fold the mean RMS). Sporadically occurring outliers are removed from the respective weekly NEQ. If outliers reflect a discontinuity (i.e. they appear in successive NEQ), a new position is estimated for this station. Most of these discontinuities are caused by the replacement of antennas and receivers, either because of defects or to modernise the equipment. In some cases, the discontinuities are caused by local (poor construction of the monument) or seismic effects and they may be followed by changes in the linear movement of the station; thus, new position and velocity values must be estimated for the affected stations (see an example in Fig. 4). Once outliers are removed and discontinuities are identified, a new multi-year solution is computed, new residual time series are generated, and the time series analysis is repeated to identify remaining outliers or discontinuities. This process is iteratively performed until no more outliers or discontinuities remain. The geodetic datum is realised by imposing two minimum constraint conditions with respect to the IGb08 coordinates of the reference stations. The no-net-translation (NNT) condition ensures that the origin [$X=Y=Z=\mathrm{0}$] of our network coincides with the origin of the network built by the reference stations. The no-net-rotation (NNR) condition aligns the orientation of our network to the orientation defined by the network built by the reference stations. In this study, we selected six IGb08 sites (WTZR, ZIMM, GRAS, GRAZ, MEDI, LROC; Fig. 7) as reference stations (see Rebischung et al., 2012; ftp://igs-rf.ign.fr/pub/IGb08/IGb08.snx, last access: 20 August 2018). The main selection criteria is based on the homogenous distribution of the reference stations across the CO-GNSS network, rare incidents of equipment changes, and a nearly complete coverage of the time series from January 2004 to May 2016.

In general, the station position time series describe clear seasonal motions (Fig. 5), which are assumed to be a consequence of non-modelled geophysical loading or local effects superimposed with GNSS inherent systematic errors (e.g. Blewitt et al., 2001; Collilieux et al., 2010, 2012; Drewes et al., 2013; King et al., 2008; Ray et al., 2008). Therefore, to increase the reliability of the constant velocities required for modelling the surface kinematics, stations with time series shorter than 2 years are excluded from further computations. The station position time series of the entire network present an averaged RMS value of about ±1 mm in the horizontal component and ±2.5 mm in the height component (Fig. 6). The final multi-year solution for the CO-GNSS network covering the Alpine region refers to the IGb08, epoch 2010.0. While the averaged RMS precision for the station positions at the reference epoch is ±1.1 mm in N and E and ±2.3 mm in height, the averaged RMS precision for the station velocities is ±0.2 mm a⁻¹ in N and E and ±0.4 mm a⁻¹ in height. Figures 7 and 8 show the horizontal station velocities (after removing the Eurasian Plate motion) and the vertical station velocities, respectively.

Figure 6 Mean RMS values of the station position time series for the period January 2004 to May 2016.

Figure 7 Horizontal station velocities relative to the Eurasian Plate with 95 % error ellipses. Labels identify the reference stations used for the geodetic datum definition: GRAS (Caussols, France), GRAZ (Graz, Austria), LROC (La Rochelle, France), MEDI (Medicina, Italy), WTZR (Bad Kötzting, Germany), and ZIMM (Zimmerwald, Switzerland). Topography and bathymetry after the ETOPO1 model (Amante and Eakins, 2009).

Figure 8 Vertical station velocities: red arrows indicate uplift, blue arrows represent subsidence, and green bars show 95 % confidence. Labels identify the reference stations used for the geodetic datum definition: GRAS (Caussols, France), GRAZ (Graz, Austria), LROC (La Rochelle, France), MEDI (Medicina, Italy), WTZR (Bad Kötzting, Germany), and ZIMM (Zimmerwald, Switzerland). Topography and bathymetry after the ETOPO1 model (Amante and Eakins, 2009).

4.3 Velocity solution

According to the results presented in Fig. 7, we observe in the Alps–Pannonian–Dinarides junction zone horizontal velocities in the range 1–4 mm a⁻¹ relative to the stable part of the Eurasian Plate. These large relative velocities are controlled by the north-eastward motion of the Adria microplate towards the Southern and Central Alps and the Dinarides. Based on our few CO-GNSS sites in the Apennines, we see apparently an extension of 2–4 mm a⁻¹ across these mountains and a 0.5–2 mm a⁻¹ shortening across the southern front of the Eastern Alps. These findings are in agreement with the conclusions presented by Devoti et al. (2017) in a recent study. They infer a crustal extension rate of about 3 mm a⁻¹ across the Apennine belt and a compression of about 2 mm a⁻¹ towards the Adriatic foreland (see Devoti et al., 2017, Fig. 9, profiles A–B and C–D). A progressive eastward rotation of the velocities towards the Pannonian Basin is also detected. The relative velocities of the stations located in the Western Alps and surrounding areas like the Po Basin, the Rhône-Bresse Graben, and southern France are very small (0.1–0.4 mm a⁻¹), allowing us to infer a very slow horizontal deformation in this region. Nocquet et al. (2016) call this a “virtually zero horizontal velocity”. As expected, stations situated in the stable part of the Eurasian Plate show velocities around zero; some exceptions reflect mostly local effects.

While the motion of the Eurasian Plate dominates the horizontal velocities, the vertical velocities (Fig. 8) mirror in general the orogenic processes in the Alpine region: uplift of the mountain chains and subsidence of the active foreland basins. A wide area of the Western Alps presents an uplift rate of approximately 2.2 mm a⁻¹ decreasing westwards. In fact, stations located west of the Alpine foreland show zero vertical rates. The southern part of the Eastern Alps shows a slower vertical uplift rate of about 1.0–1.5 mm a⁻¹, which decreases towards the north. The CO-GNSS stations located in the northern margin of the Alps in Germany and Austria present uplift rates of less than 0.6 mm a⁻¹ decreasing eastwards. Subsidence rates of about 1 mm a⁻¹ are observed in the Rhine and Rhône-Bresse Grabens as well as in the Venetian-Friuli Basin. They increase southwards and reach up to 2 mm a⁻¹ at the coasts in Marseille and Venice, respectively.

In general, it can be stated that the vertical velocities shown in Fig. 8 are a result of the combination of three main mechanisms: (i) uplift caused by the ongoing convergence of Eurasia and Africa; (ii) orogenic gravitational movement due to erosion, collapse, and associated isostatic compensation; and (iii) isostatic rebound caused by the melting of the Alpine ice shields since the last glacial maximum 22 000 years ago, i.e. glacial isostatic adjustment (GIA). The orogenic gravitational movement is much slower than the GIA; however, it is clear that the eroded rock masses are redistributed to the surroundings with lower elevation (delta valleys and sedimentary basins). This may cause a bedrock uplift and surface height lowering in the mountains as well as a bedrock subsidence and a regional surface height increase in the foreland basins. According to Persaud and Pfiffner (2004), Barletta et al. (2006), and Stocchi et al. (2005), the crustal uplift caused by GIA in the Alpine region varies from 0.06 to 0.20 mm a⁻¹. This assessment suggests that the first two mechanisms (plate convergence and orogenic gravitational movement) are the main contributors to the vertical rates obtained in this study.

5 Surface-kinematic modelling based on CO-GNSS network solutions

Based on the precise station velocities estimated previously, the objective of this section is to describe the surface-kinematic modelling performed within this study. We distinguish three main components: a deformation model, a continuous velocity field, and a strain field. A common trend motion is removed from the CO-GNSS station velocities to obtain the so-called residual velocities. The common trend motion represents a large-scale variation, while the residual velocities represent short-scale variations. The common trend motion in the horizontal component is usually well represented by the motion of the tectonic plate underlying the area of study, in this case the Eurasian Plate (see Sect. 5.1). The common trend motion in the vertical component may be assumed to be the average vertical movement of the area of study. In this work, we infer this mean vertical movement from the station velocities directly (see Sect. 5.2 and Appendix A). Once the trend motion is removed from the horizontal and vertical station velocities, the point-wise residual velocities are correlated with the existing tectonic structures (Fig. 1) and then they are interpolated to a regular grid. The interpolation of the residual velocities provides us with the deformation model (Sect. 5.2 and 5.3), which is the basis for the computation of the strain field (Sect. 6). The common trend motion (i.e. the Eurasian Plate motion and the average vertical motion) removed from the initial station velocities is restored to the deformation model to get a continuous velocity field (Sect. 5.4). Various geodetic and geophysical methods based on mathematical interpolation approaches or physical crust models may be applied to represent the continuous surface deformation from point-wise velocities. A frequently used interpolation method is the least-squares collocation (e.g. Moritz, 1973; Drewes, 1978); to include the physical properties of the deforming body, the finite element method is an appropriate tool (e.g. Meissl, 1981; Heidbach and Drewes, 2003). It has been demonstrated that for the sole representation of the Earth surface deformation, the results of both methods are nearly identical (e.g. Drewes and Heidbach, 2005). Therefore, we use here the least-squares collocation method.

5.1 Least-squares collocation (LSC) approach

The primary principle of LSC is the correlation of physical parameters (in this case, station velocities) between adjacent points. The observations are divided into a systematic linearised part (or trend) and two independent random parts: the signal and the observation error (or noise). The parameters describing the systematic component and the stochastically correlated signals are estimated by minimising the noise. The estimated parameters and signals together with a given spatial correlation allow the prediction (interpolation) of “fictitious measurements” at regions where no observations exist. In spatial applications, the correlation function is usually considered as a function dependent on the distance. The so-called correlation distance d should be in agreement with the spacing of the stations and the wavelength of the expected signal. If d is smaller than the mean station spacing, artefacts may appear in the interpolation. If d is too large, a smoothing similar to a filtering effect may occur. According to the station distribution in our network (Fig. 2), the initial maximum distance was set to 100 km for selecting the points of a collocation domain. If not enough points were found, the distance was enlarged until the necessary number of at least four points was available. The grid size for the interpolation was set to 25 km × 25 km. This grid size is chosen as appropriate since it corresponds to the mean station spacing in the most densely covered region within the network (i.e. the French and Swiss Alps). A larger grid size would average out the results and this filtering should be avoided. Appendix A provides a detailed description of the LSC formulation applied in this study.

According to the tectonic plate boundary model PB2002 (Bird, 2003), we assume that the Adria, Pannonia, and Liguria blocks constitute a deformation zone along the southern Eurasian Plate boundary. We do not consider small lithospheric blocks rotating with independent Euler vectors, but a continuous and viscous lithosphere deforming under certain kinematic boundary conditions as suggested, for instance, by Flesch et al. (2000); Vergnolle et al. (2007); or Copley (2008). However, as horizontal and vertical station motions are dominated by different mechanisms (see Sect. 4.3), the LSC procedure is applied separately: as a vector prediction for the horizontal component and as a scalar prediction for the height component. In the empirical procedure, the geocentric Cartesian station positions $\left(X,Y,Z\right)$ and velocities $\left(v_X,v_Y,v_Z\right)$ obtained within the multi-year solution (Sect. 4.2) are transformed together with their precision values into ellipsoidal coordinates $\left(\mathit{\varphi },\mathit{\lambda },h\right)$ and $\left(v_{\mathit{\varphi }},v_{\mathit{\lambda }},v_h\right)$. The Eurasian Plate motion is removed from horizontal velocities (v_ϕ,v_λ) of the stations. The corresponding rotation vector Ω ($\mathrm{260.74}{}^{\circ }\pm \mathrm{0.53}{}^{\circ }$ E, $\mathrm{55.14}{}^{\circ }\pm \mathrm{0.27}{}^{\circ }$ N, 0.2598±0.0011^∘ Ma⁻¹) is inferred from our CO-GNSS stations located on the stable part of the plate following the approach presented by Drewes (1982, 2009). Our values are very similar to the rotation vector derived from the ITRF2008 Ω ($\mathrm{261.15}{}^{\circ }\pm \mathrm{0.85}{}^{\circ }$ E, $\mathrm{55.23}{}^{\circ }\pm \mathrm{0.35}{}^{\circ }$ N, $\mathrm{0.2570}\pm \mathrm{0.0025}{}^{\circ }$ Ma⁻¹) (see Altamimi et al., 2012). The small differences are a consequence of the different station distribution (stations used for the estimation of Ω in this study and in the ITRF2008 are not the same) and the different time span considered for the station velocity estimation (our study includes observations up to May 2016, while the ITRF2008 includes observations up to May 2009; Altamimi et al., 2011).

After removing the Eurasian Plate motion from the station horizontal velocities, the LSC interpolation is applied following the formulation described in Appendix A. In this procedure, outliers with discrepancies larger than 0.5 mm a⁻¹ with respect to the velocities of the surrounding stations are considered to be affected by local effects and they are not included in the collocation.

5.2 Horizontal deformation model

Figure 9 GNSS-inferred horizontal deformation model of the Alpine region. Red arrows represent the predicted surface deformation. Blue arrows represent the input residual velocities. Green dots show the CO-GNSS station distribution. Uncertainty ellipses are separately shown in Fig. 10. Topography and bathymetry after the ETOPO1 model (Amante and Eakins, 2009).

The deformation model presented in Fig. 9 describes the surface kinematics of the Alpine region with respect to the Eurasian Plate. Two distinct kinematic patterns can clearly be seen: the Western Alps as well as the western and northern surrounding areas (Liguro-Provençal Basin, Massif Central mountains, Rhône-Bresse Graben, Jura Mountains, and Molasse Basin) show a small and statistically insignificant motion, indicating that these regions move north-eastward with the Eurasian Plate. This stationary pattern contrasts with the pronounced N- and NE-oriented movement predicted in the Southern and Eastern Alps as well as in the Apennine Peninsula, the Venetian-Friuli Basin, and the northern part of the Dinarides. It is also evident that the magnitude of deformation increases from 0.1 mm a⁻¹ in the western part of Switzerland to the south and to the east, reaching 3.0 mm a⁻¹ in the northern part of the Apennine Peninsula and about 1.5 mm a⁻¹ in the Austrian Alps, respectively. The two kinematic patterns are apparently separated by the so-called Periadriatic Line, a fault system in EW direction along the boundary between the Southern Alps and the Adriatic region (around latitude 46^∘ N, see Fig. 1). The magnitude of the deformation vectors across the westernmost part of the Po Basin and the south of the Western Alps is less than 0.2 mm a⁻¹. The average residual motion with respect to the Western Alps is not larger than 0.1 mm a⁻¹, indicating that presently there is no shortening in this region. These findings are confirmed by Nocquet et al. (2016), who suggest an upper bound of 0.3 mm a⁻¹ (with 95 % confidence) for possible right-lateral strike slip motion in the western Alpine area.

Unlike the Western Alps, where the deformation vectors indicate a very small (insignificant) internal surface deformation, the Central Alps present an increasing north-oriented deformation from 0.2 mm a⁻¹ at the longitude 6^∘ E to 0.6 mm a⁻¹ close to longitude 11^∘ E. This deformation pattern continues over the Southern Alps up to longitude 13^∘ E, where the deformation vectors describe a progressive eastward rotation towards the Pannonian segment, reaching magnitudes up to 1.3 mm a⁻¹ and an orientation of about N20^∘ E near longitude 16^∘ E. These vectors also make evident a shortening of about 2 mm a⁻¹ across the southern front of the Eastern Alps, in the northern area of the Venetian-Friuli Basin. This is in agreement with the results published by Devoti et al. (2017), Métois et al. (2015), and Cheloni et al. (2014). Although our CO-GNSS network presents a sparse station distribution in the northern Apennines and the Dinarides, our prediction supports the conclusions of D'Agostino et al. (2014) and Nocquet (2012) in those regions. In the former case, we identify a motion roughly oriented to the north-east with the largest magnitudes (∼ 3 mm a⁻¹) across the highest topographic features of the Apennines towards the Adriatic Sea. In the opposite direction, towards the Tyrrhenian Sea, the deformation magnitudes decrease to 1.5 mm a⁻¹. According to Mantovani et al. (2015), this may be evidence of a faster motion of the outer Apennine belt with respect to the inner one. In the case of the Dinarides, a shortening of about 3.0 mm a⁻¹ across the Adriatic coastline is found. This magnitude decreases towards the Pannonian Basin to 0.5 mm a⁻¹. In the southern front of the Alps, the orientation of the deformation vectors presents a slightly progressive westward rotation from the area of Venice towards the Po Basin. These findings are quantitatively equivalent to the results presented by Möller et al. (2011) and Devoti et al. (2011, 2017). However, it should be kept in mind that, due to the low number of GNSS stations we processed in the Po Basin, our results are highly influenced by the linear gradient imposed by the LSC approach to the deformation model. To increase the reliability of our model in this particular zone, a major number of GNSS stations covering the Italian Alpine forelands should be considered. Figure 10 shows the uncertainty (95 % confidence) of the horizontal surface deformation model. The largest values (from ±0.5 to ±1.5 mm a⁻¹) occur at those regions with poor station coverage (mainly the Apennines, Dinarides, and easternmost part of the Austrian Alps). In general, the mean accuracy of the model is ±0.2 mm a⁻¹.

Figure 10 Estimated uncertainty (95 % confidence) of the GNSS-based horizontal deformation model of the Alpine region. Green dots show the CO-GNSS station distribution. Topography and bathymetry after the ETOPO1 model (Amante and Eakins, 2009).

5.3 Vertical deformation model

The present-day vertical motion of the Alps and their surrounding areas is displayed in Fig. 11. Our results clearly show an ongoing average uplift of 1.8 mm a⁻¹ of the entire mountain chain, with the exception of the southern part of the Western Alps (6.5^∘ E, 44^∘ N), where no significant uplift (less than 0.5 mm a⁻¹) is detected. The largest uplift rates (more than 2 mm a⁻¹) occur in the central area of the Western Alps (7^∘ E, 45^∘ N), in the Swiss Alps (8^∘ E, 46.5^∘ N), and in the Southern Alps in the boundary region between Switzerland, Austria, and Italy (11^∘ E, 46.6^∘ N). The uplift kinematics we present here for the Western Alps are consistent with the vertical motion patterns identified by Nguyen et al. (2016) and Nocquet et al. (2016). Our vertical rates in the Swiss Alps also agree qualitatively with similar results inferred by Brockmann et al. (2012) using not only GNSS observations, but also levelling data.

Figure 11 GNSS-inferred vertical deformation model of the Alpine region. Green dots show the CO-GNSS station distribution. Topography after the ETOPO1 model (Amante and Eakins, 2009).

In the north-west area, across the Rhône-Bresse and Rhine Grabens, a well-defined regional subsidence pattern with varying rates is observed. Apparently, the largest subsidence of around −1.3 mm a⁻¹ happens in the Rhine Graben (7^∘ E, 48.5^∘ N), while in the Rhône-Bresse Graben (5^∘ E, 46–48^∘ N) the largest rate reaches −0.8 mm a⁻¹. In the complete studied area, the largest subsidence rates are close to −1.5 mm a⁻¹ and are located west of the Venetian-Friuli Basin. We do not detect a significant subsidence in the western part of the Po Basin; this may be a consequence of the poor distribution of our CO-GNSS stations in that region. Actually, Devoti et al. (2017) inferred a mean subsidence rate of about −0.8 mm a⁻¹ in the Po Basin. The maximum (−3 mm a⁻¹) and minimum (−0.5 mm a⁻¹) magnitudes occur near Venice and in the eastern margin of the Western Alps, respectively (see Devoti et al., 2017, Fig. 6). Results of previous studies like Cuffaro et al. (2010) or Serpelloni et al. (2013) suggest subsidence rates up to −11 mm a⁻¹ near Modena, Italy (11^∘ E, 44.6^∘ N). However, we do not detect these extreme vales. This may be explained by the fact that we eliminate stations with strong local effects before the computation of the vertical-kinematics model, while the two studies mentioned base their conclusions on point-wise station velocities. It is well known that GNSS systematic errors and non-modelled local effects accumulate mainly in the estimation of the vertical coordinates. Therefore, we gave special care to identify and to remove those stations reflecting local effects (like water extraction, tunnel construction, and landslides) that may mislead the estimation of vertical motions mirroring mountain building or subsidence processes in the Alpine region.

The general uplift observed across the Alpine mountain chain decreases towards the outer regions to stable values between 0.0 and 0.5 mm a⁻¹ and in some cases to subsidence like in the Liguro-Provençal (5^∘ E, 43.6^∘ N) and Vienna (15.5^∘ E, 48.5^∘) basins, where vertical rates of −0.8 and −0.3 mm a⁻¹ are detected, respectively. Figure 12 displays the estimated uncertainty of the predicted uplift and subsidence rates. The largest values (from ±0.6 to ±0.8 mm a⁻¹) occur in those areas with poor station coverage, mainly in the Apennine Peninsula and the Po Basin. The mean precision of the vertical motion model computed in this study for the Alpine mountain belt is about ±0.3 mm a⁻¹.

Figure 12 Estimated uncertainty (95 % confidence) of the GNSS-based vertical deformation model of the Alpine region. Dark blue dots show the CO-GNSS station distribution. Topography after the ETOPO1 model (Amante and Eakins, 2009).

Figure 13 Velocity profiles in the Western Alps (cross section along longitude 6.5^∘ E, a) and in the Eastern Alps (cross section along longitude 13^∘ E, b). Blue dots with error bars represent the observed station velocities (with 68 % confidence). Red lines represent the north (V_n) and the vertical (V_u) components of the deformation model with their uncertainty (light red stripe). Green line shows the average topography in the profile swath, with light and dark grey shadows showing the maximum and minimum elevations, respectively.

Figure 13 shows two velocity profiles crossing the Alps along longitudes 6.5 and 13^∘ E, respectively. In both cases, the vertical component (V_u) of the surface deformation presents a high correlation with the local topography. The profile in the Western Alps clearly shows the light subsidence (−0.6 mm a⁻¹) close the Liguro-Provençal Basin (near latitude 43^∘ N); the large uplift rates (more than 1.5 mm a⁻¹) in the border area between France and Switzerland (around latitude 45.5^∘ N); and a decreasing uplift rate towards the north up to the Rhine Basin (latitude 48^∘ N), where subsidence is observed. In the Eastern Alps, the gradient of the vertical velocity is stronger in the southern part than in northern area. South of latitude 46^∘ N, the subsidence (up to −1 mm a⁻¹) detected in the Venetian-Friuli Basin can be well observed. North of latitude 46^∘ N, we find uplift rates up to 1.2 mm a⁻¹ in the border area between Italy and Austria (latitude 45^∘ N). These rates diminish up to 0.3 mm a⁻¹ in the border area between Austria and Germany (latitude 48^∘ N). The north component (V_n) of the deformation model along longitude 6.5^∘ suggests a very small (not significant) surface deformation in the Western Alps, while along longitude 13^∘ E it captures the plate boundary region where the Adria microplate collides with the Eurasian Plate. The quite strong velocity gradient of nearly 2 mm a⁻¹ over about 80 km (between latitudes 46 and 47^∘ N) makes evident the NS compression occurring in the southern front of the Southern Alps.

5.4 Continuous velocity field

GNSS techniques are currently the primary tool for precise positioning. In applications of high reliability (accuracy at the 1 cm level or better), the GNSS data processing requires the reference station coordinates to be in the same reference frame and at the same observation epoch as the GNSS satellite orbits. In active seismic regions, strong earthquakes cause large changes in the station positions and velocities of the geodetic reference stations, disabling the use of previous coordinates (see for instance co-seismic and post-seismic effects of the El Maule earthquake in the Latin American reference frame presented by Sánchez and Drewes, 2016). To ensure the long-term stability of geodetic reference frames, the transformation of station positions between different epochs requires the computation of reliable continuous surface deformation (or velocity) models. Based on this kind of models, geodesists monitor the kinematics of reference frames, infer appropriate transformation parameters between pre-seismic and post-seismic (deformed) coordinates (specially in official matters like legal borders, cadastre, and land management), and interpolate surface motions arising from plate tectonics or crustal deformations in areas where no geodetic stations are established. Consequently, after computing the deformation model (Sect. 5.2), the Eurasian Plate motion removed from the CO-GNSS station velocities is restored to the deformation vectors to get a continuous velocity field referred to the ITRF. As shown in Fig. 14, the velocity field mirrors the present-day surface kinematics as a combination of deformation patterns and the Eurasian Plate tectonic motion, which is the dominant kinematic regime in Europe.

Figure 14 Continuous velocity field for the Alpine region referring to the ITRF2008 (a: horizontal velocity field, b: vertical velocity field). Topography and bathymetry after the ETOPO1 model (Amante and Eakins, 2009).

6 Strain field

The strain field represents the zones with compression, extension (or dilatation), and shear strain on the Earth's surface. In this study, the strain rates are computed from the deformation model (Fig. 9) as it provides a better (smoothed) surface coverage than the residual station velocities (Fig. 7). This strategy reduces random errors occurring at individual stations and enhances the major trends on strain patterns. Our computations are based on the infinitesimal strain theory, assuming that the deformation is much smaller (indeed, infinitesimally smaller) than any relevant dimension of the body, in this case the Alpine region. It means very small deformation changes (at the mm a⁻¹ level) over large distances (at the km level) between grid points are presumed. In the empirical procedure, the deformation vectors computed at the grid points are used for the estimation of deformation gradients, and these gradients are used for the computation of the strain components following the equations described by Lambeck (1988, pp. 194–196). Even though the deformation field is available in three components, we concentrate only on the horizontal components. Thus, the obtained strain field is a horizontal projection of the surface deformation. The core contribution of a strain field is the possibility of distinguishing large-scale features (e.g. surface deformation and mountain building) from areas with large strain rates, which are often seismically active.

Figure 15 Strain distribution in the Alpine region. Blue shades represent dilatation; red shades represent compression. Topography after the ETOPO1 model (Amante and Eakins, 2009).

Figure 16 Shear strain distribution in the Alpine region. Blue and red shades represent the left-lateral and right-lateral shear, respectively. For illustration purposes, the left-lateral shear is considered as negative and the right-lateral shear as positive. Topography after the ETOPO1 model (Amante and Eakins, 2009).

Figure 15 shows the principal strain components (compression and dilation) and the first strain invariant (sum of strain components) obtained in this study. Two main contrasting strain regimes are noticed: (i) shortening across the south-eastern front of the Alps and the northern part of the Dinarides and (ii) extension in the Apennines. The pattern of the strain principal axes indicates that the compression directions are more or less perpendicular to the thrust belt fronts (see Fig. 1). This compression presents maximum strain-rate values of about $\mathrm{20}\times {\mathrm{10}}^{-\mathrm{9}}$ a⁻¹ in the Venetian-Friuli and Po basins. These findings are in agreement with previous works published by Serpelloni et al. (2005, 2006), Caporali et al. (2009), Bennett et al. (2012), and Cheloni et al. (2014); in particular, the continuous compression regime we inferred along the Adriatic coastline fits very well with the strain model presented by Carafa and Bird (2016). Nevertheless, the lack of data in the eastern margin of the Adriatic Sea does not permit us to conclude if this regime continues with similar magnitudes to the southern Dinarides. As already stated by D'Agostino et al. (2014), dilatation is dominant in the north-western part of the Apennine Peninsula. Apparently, this dilatation continues to the west across the Po Basin, although at minor rates. However, poor station coverage in this area precludes any definitive conclusions from these observations.

Across the Alpine mountain belt, we observe a slight dilatation regime in the Western Alps, which smoothly changes to a contraction regime in western Austria and southern Germany, reaching maximum shortening values of $\mathrm{6}\times {\mathrm{10}}^{-\mathrm{9}}$ a⁻¹ in the north-eastern part of Austria (eastern margin of the Molasse Basin). The strain rates obtained for the French and Swiss Alps are in the same order of uncertainty of our results; however, they are in agreement with results presented in previous studies like Calais et al. (2002) and Delacou et al. (2008). The magnitude and orientation of the principal strain axes we observe in the boundary region between Switzerland, Austria, and Germany are very similar to the strain-rate patterns described by Tesauro et al. (2006).

The largest shear strain of about $\mathrm{20}\times {\mathrm{10}}^{-\mathrm{9}}$ a⁻¹ is detected in the south-eastern part of the Alps (Fig. 16). The sharp boundaries between the left-lateral and right-lateral shear zones and a quiet area reveal the triple junction between the Eurasian, Adriatic, and Pannonian plates (Möller et al., 2011). According to Brückl and Hammer (2014), the approximate location of the triple junction is 13.7^∘ E, 46.6^∘ N, which is well correlated with the strain model computed in this study. The Western Alps and surroundings do not reveal any specific pattern of the shear strain while the Dinarides show a compressional strain with a right-lateral shear of about 10 to $\mathrm{20}\times {\mathrm{10}}^{-\mathrm{9}}$ a⁻¹. The resulting principal and shear strain fields show that the main tectonic activity is ongoing on the plate boundaries.

7 Data availability

Station positions and velocities as well as velocity and deformation fields computed in the frame of this study are available through the PANGAEA (Data Publisher for Earth and Environmental Science) platform at https://doi.pangaea.de/10.1594/PANGAEA.886889. The dataset collection is composed of the following.

• ALPS2017.SNX: multi-year solution for geocentric Cartesian station positions and velocities referring to the IGb08, epoch 2010.0, in SINEX, the Solution (Software/technique) INdependent EXchange format.

• ALPS2017_XYZ.CRD: geocentric Cartesian station positions.

• ALPS2017_XYZ.VEL: geocentric Cartesian station velocities.

• ALPS2017_NEH.CRD: ellipsoidal station positions.

• ALPS2017_NEH.VEL: ellipsoidal station velocities.

• ALPS2017_EPR.VEL: station velocities relative to the Eurasian Plate.

• ALPS2017_DEF_HZ.GRD: surface deformation model of the Alpine region.

• ALPS2017_DEF_VT.GRD: vertical deformation model of the Alpine region.

• ALPS2017_STR.GRD: strain field of the Alpine region.

• ALPS2017_VEL.GRD: continuous velocity field of the Alpine region.

8 Conclusions

In this study, we provide a present-day surface-kinematics model homogeneously assessed for the whole Alpine mountain chain. It is based on a high-level GNSS data analysis covering 12.4 years of continuous geodetic observations. Together with a coordinate solution for more than 300 CO-GNSS stations, our results comprise a deformation model, a continuous velocity field, and a strain field covering the area between 4–16^∘ E and 43–50^∘ N at a spatial interval of 25 km × 25 km. The average uncertainty of our model is ±0.2 mm a⁻¹ in the horizontal component and ±0.3 mm a⁻¹ in the height component. Although the station coverage is heterogeneous, our model describes coherent patterns of ongoing uplift processes and horizontal deformations with respect to the Eurasian Plate. The determination of the constant velocities, which form the basis for the deformation model and strain field, was particularly complex, because each little discontinuity had to be revealed, analysed, and taken into account in determining the constant velocity of each site. Because the geodetic network consists of many permanently operating GNSS stations, which have to fulfil different purposes, the analysis was difficult and required a lot of manual intervention. Some sites fulfil geodetic services requiring frequent updates in the set-up, which in turn cause inconsistencies in the time series that had to be recovered. Not all of these changes were well documented. The GNSS data processed in this study are unique as they cover the entire Alpine range and the associated foreland. The comparison of our derived deformation patterns with the results of other studies, whose research areas were mostly limited to smaller regions of the Alps, has proven the high reliability of our model. With this data set, a deformation model of the Alps derived from a uniform analysis is now available.

Appendix A: Least-squares collocation equations

The basic least-squares collocation formula is given by (e.g. Moritz, 1973; Drewes, 1978)

$$\begin{array}{}\text{(A1)}& {\displaystyle }{\mathit{v}}_{\mathrm{pred}}={\mathbf{C}}_{\mathrm{new}}^T{\left({\mathbf{C}}_{\mathrm{obs}}+{\mathbf{C}}_{\mathrm{nn}}\right)}^{-\mathrm{1}}{\mathit{v}}_{\mathrm{obs}}.\end{array}$$

v_obs contains the velocities observed at the geodetic stations (Figs. 7 and 8). v_pred contains the velocities to be predicted at the continuous grid of 25 km × 25 km (Figs. 9 and 11). C_obs is the correlation matrix between observed velocities. C_new is the correlation matrix between predicted and observed velocities. C_nn is the noise covariance matrix; i.e. it contains the uncertainty of the station velocities obtained within the multi-year solution (error ellipses and error bars in Figs. 7 and 8).

The empiric correlation between the observed velocities may be written as

$$\begin{array}{}\text{(A2)}& {\displaystyle }{\mathbf{C}}_{\mathrm{obs}}\left(d_{ik}\right)=E\left\{v_i\cdot v_k\right\},\end{array}$$

where E is the expectation; i,k represent the geodetic stations; and (d_i,k) is the distance between the station i and k. The C_obs values are classified in Δ_j class intervals with j=1…10. The cross-covariance C_obs(Δ_j) and auto-covariance ${\mathbf{C}}_{\mathrm{obs}}\left(d=\mathrm{0}\right)=C_{\mathrm{0}}$ are determined for each interval Δ_j:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathbf{C}}_{\mathrm{obs}}\left({\mathrm{\Delta }}_j\right)={\displaystyle \frac{\mathrm{1}}{n_j}}\sum _{i<k}^j{v}_i\cdot v_k;\\ \text{(A3)}& {\displaystyle }& {\displaystyle }{\mathbf{C}}_{\mathrm{obs}}\left(d=\mathrm{0}\right)=C_{\mathrm{0}}={\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n{v}_i^{\mathrm{2}},\end{array}$$

where n stands for the number of stations available at the domain defined by d, while n_j represents the number of stations available at each class interval Δ_j. After estimating the discrete empirical covariances with Eq. (A3), they are approximated by a continuous function C(d_ik): in this case, the exponential function

$$\begin{array}{}\text{(A4)}& {\displaystyle }C\left(d_{ik}\right)=a{e}^{-b\cdot d_{ik}}.\end{array}$$

The function parameters a and b are estimated by a least-squares adjustment. C_obs is symmetrical and its main diagonal (i=k) contains the values C₀. The elements of C_new are computed using the same Eq. (A4) as a function of the distance between the grid node to be interpolated and the geodetic stations. The precision ${\mathit{\sigma }}_{\mathrm{new}}^{\mathrm{2}}$ of the predicted values is assessed with (Moritz, 1973)

$$\begin{array}{}\text{(A5)}& {\displaystyle }{\mathit{\sigma }}_{\mathrm{new}}^{\mathrm{2}}=C_{\mathrm{0}}-{\mathbf{C}}_{\mathrm{new}}^T({\mathbf{C}}_{\mathrm{obs}}+{\mathbf{C}}_{\mathrm{nn}}{)}^{-\mathrm{1}}{\mathbf{C}}_{\mathrm{new}}.\end{array}$$

To fulfil the statistical condition of stationarity, the correlation function Eq. (A4) has to be identical for all points located in the d region; i.e. the exponential function must be homogeneous over the entire domain of known and predicted points. A consequence of the isotropy and stationarity conditions is that any trend (e.g. plate tectonics) in the given set of point parameters has to be removed before starting the collocation procedure. As horizontal and vertical station motions are dominated by different mechanisms (see Sect. 4.3), it is necessary to apply the LSC procedure separately. For the horizontal model, the north-east movement of the Eurasian Plate defines the trend motion to be removed from the station velocities. For the vertical model, the trend motion is given by the mean value of the vertical station velocities available at each d region.

As mentioned in Sect. 5, the deformation model is computed at a regular grid of 25 km spacing interval. This grid size is chosen as appropriate since it corresponds to the mean station spacing (about 20 to 30 km) in the most densely covered region within the network (the French and Swiss Alps, see Fig. 2). For the horizontal component of the deformation model, the minimum correlation distance d is defined to be 100 km. This radius provides an overlap between adjacent grid points and allows a certain degree of smoothing. For areas with poor station coverage, the correlation distance is extended until at least four stations are available for the LSC. For the vertical component, a larger correlation distance (up to 300 km) is needed as the number of rejected outliers is higher for the vertical station velocities than for the horizontal ones.