Present-day uplift of the western Alps

Collisional mountain belts grow as a consequence of continental plate convergence and eventually disappear under the combined effects of gravitational collapse and erosion. Using a decade of GPS data, we show that the western Alps are currently characterized by zero horizontal velocity boundary conditions, offering the opportunity to investigate orogen evolution at the time of cessation of plate convergence. We find no significant horizontal motion within the belt, but GPS and levelling measurements independently show a regional pattern of uplift reaching ~2.5 mm/yr in the northwestern Alps. Unless a low viscosity crustal root under the northwestern Alps locally enhances the vertical response to surface unloading, the summed effects of isostatic responses to erosion and glaciation explain at most 60% of the observed uplift rates. Rock-uplift rates corrected from transient glacial isostatic adjustment contributions likely exceed erosion rates in the northwestern Alps. In the absence of active convergence, the observed surface uplift must result from deep-seated processes.


Data set & processing strategy
The Continuous GPS data used in this study are from the RENAG network (http://webrenag.unice.fr), from the RGP network (http://rgp.ign.fr), from the ORPHEON network (http://reseau-orpheon.fr), from the Piemonte region (http://gnss.regione.piemonte.it), from the RING network from the Istituto Nazionale di Geofisica e Vulcanologia (INGV, http://ring.gm.ingv.it) and from EUREF (http://www.epncb.oma.be). The period 2002.0-2013.5 was analyzed using the GAMIT/GLOBK 10.5 software package [1]. The analysis is divided into 3 steps: (1) raw phase and pseudo-range observations are reduced to produce daily loosely constrained solutions (2) loosely constrained daily solutions are expressed in a consistent reference frame using a Helmert transformation to obtain daily positions time series (3) time series are analyzed in order to detect outliers and offsets. The corrected and cleaned time series are used to assess the noise properties and finally derive the velocities estimates with their associated uncertainties.
For the first step, we use the final combined orbits from the International GNSS Service for Geodynamics [2] derived from the first reprocessing campaign (http://acc.igs.org/reprocess.html). Because the orbits were reprocessed using a consistent strategy for the whole period, they were taken as fixed in the processing. We use the ocean loading correction model FES2004 computed from the tidal hydrodynamic equations and data assimilation [3]. We use the the Vienna Mapping Function (VMF1, [4]) for both hydrostatic and nonhydrostatic models of the tropospheric delay and the global pressure and temperature model. We estimate a zenithal trophospheric delay parameter every 2 hours and two horizontal atmospheric gradient per day as stochastic parameters. We also correct for atmospheric tides and atmospheric loading. The atmospheric loading corrections used together with the VMF1 mapping functions have proven to provide the best results in terms of repeatabilities and noise characteristics [5]. Finally, IGS absolute calibration of Phase Center Variations (PCV) of the antennas have also been used, as recommended by the IGS.
Our second step differs from most used strategies. Instead of expressing our solution in the current ITRF [6] and then try to remove the regional common mode motion, we start by identifying a subset of sites with good repeatabilities and little loss of data over the studied time window. We find that CGPS sites SOPH, MTPL, SJDV, ZIMM and GENO meet such criteria. We then used this subset of sites to define a local reference frame using the following strategy: (1) in a first run, we use their coordinates defined in the ITRF2008 at epoch 2007.0 [6] with 0 velocity as the a priori reference frame and estimate a 7-parameters transformation to derive a first solution (2) then, we estimate the new position at the reference epoch and velocity for the chosen subset of sites, defining the reference frame for the next iteration (3) this process is iterated until estimated positions and velocities agree with the reference frame at the level of 0.1 m for positions and 0.01 mm/yr for velocities. Compared with the classical methodology that uses the raw ITRF as the reference frame, our approach improves the repeatability by a factor of 2, removes the regional common mode signal and improves the noise characteristics. Furthermore, this approach is more rigorous than the classical ITRF/common mode estimation in the sense that the true shape of the network is preserved throughout the procedure because we use only 7-parameters transformations. As a consequence, the obtained results only reflect the change of the network shape through time.
In a third step, we used the approach proposed by [7] and implemented in the CATS software [8] to assess the noise characteristics and estimate offsets, velocity, and amplitude of seasonal variations, together with the associated uncertainties.

Results & precision
GPS results are provided in Supplementary Information Table 1. Extended Data Figure 1 shows the time series for a selection of sites. Long-term repeatabilities on the vertical component are found in the range of 2.0 to 4.5 mm. As for most GPS time series, we find that a combination of white and flicker noise best explains the noise characteristics on the GPS time series [9]. Extended Data Figure 2 shows the histogram of the spectral index κ characterizing the noise when estimated through the Maximum Likehood Estimator embedded in CATS. It shows that −1.6 < κ < −0.4, with a median value of -0.7 for the verticall component. As a result of this analysis, formal errors on velocity components are of the order of a few tenths of mm/yr (0.1-0.6 mm/yr), with the best behaving sites having a formal error of 0.1 mm/yr. Although these formal uncertainties might still be over-optimistic, further external criteria described below indicate that main trends of the vertical velocity field can be trusted at the 0.3 mm/yr level.

Kinematic conditions at the boundaries of the western Alps
We use the best determined sites in our solution to estimate the relative motion between sites located west of the Alpine foreland in France and sites located in the westernmost part of the Po plain in Italy. We calculate the Euler pole for the sites located east of the western Alps with respect to the sites located west of the Alps (Table 3). We do not detect any internal deformation in the western Po plain, as indicated by a wrms of 0.07 mm/yr. The wrms for the Alpine foreland in France is 0.05 mm/yr. The relative pole between the two domains predicts 0.1-0.2 mm/yr of strike-slip motion between the western part of the Adriatic plate and the stable foreland of the western Alps in France ( Figure 1a of the main text). Accounting for the Euler pole uncertainty, we therefore find that 0.3 mm/yr is a strict upper bound of possible motion the western Po plain. Because at this level of precision we might reach the limitation of idealized rigid micro-blocks, we also used the more direct observation of relative velocities between sites west of the Alpine foreland and the western Po plain (Supplementary Table 2). Individual baselines show that the relative motion average over the sites is < 0.1mm/yr. No velocity > 0.25 mm/yr is found between pairs of sites, confirming 0.3 mm/yr as an upper bound for kinematic boundary conditions of the western Alps.

Internal horizontal deformation within the western Alps
We calculate an Euler Pole for the sites located within the Alps and its surroundings and within the Alps only. The first calculation leads to a wrms=0.23 mm/yr and the second to wrms=0. 26 mm/yr, indicating a very small amount of possible horizontal deformation within the western Alps. The largest residuals are found in the northwestern Alps close to the area of maximum uplift. As indicated by Table 4, residuals are found to be statistically highly significant for sites CHTL, JANU, MODA, PUYA, GRAV, LFAZ, VILR, MRGE. Using the sites located in the area of significant uplift, we find that the deformation is primarly extensionnal at a~20 nstrain/yr are in a N120°direction (Table 5) Stable Alpine foreland vertical reference frame As for the horizontal velocities, we want the vertical reference frame to be relevant for the tectonic processes studied here. We therefore discard the use of the ITRF which realizes a center of the mass for the whole Earth frame. Since our study focuses on the Alps dynamics, we choose to define a reference frame outside of the Alps and test its stability. Around the Alps, the sites having the best determined vertical velocities are (MTPL, GENO, SJDV, GRAS). They show a level of internal agreement of 0.13 mm/yr. With respect to this reference frame, all sites surounding the western Alps with a vertical velocity uncertainty σ up < 0.3 mm/yr agree within 0.21 mm/yr (wrms). The vertical rate averaged over (MTPL, GENO, SJDV, GRAS) therefore provide a stable vertical reference frame, with no internal vertical deformation at the precision of the data and making sense from a tectonic point of view. The chosen reference frame certainly does no cause any bias larger than 0.2 mm/yr in the vertical rates seen in the western Alps.

Data set
Levelling data are from the National Geographic Institute of France (IGN, http://www.ign.fr), partly analyzed in [10,11,12]. For each profile i and levelling benchmark j, data were provided as height differences ∆z i j between the two epochs of measurements with respect to a reference benchmark of each profile i following the methodology described in [13]. The time span between two levelling surveys ranges from 58.25 to 112.83 years. Assuming constant velocity, the levelling input data set prior adjustment is ∆v i j = ∆z i j ∆t i , where ∆t i is the time span between two successive surveys of profile i. We used 24 levelling profiles covering the western Alps and its foreland in France (Extended Data Figure 3 & Extended Data Figure 4), including a total of 1655 benchmarks, with an average distance between benckmarks of~2km. Details of the levelling protocol for both periods of surveys are described in [14] and [11].

Combined levelling-GPS adjustment model
We use a combined least-squares adjustment to derive a vertical velocity field including levelling profiles and GPS vertical velocities. The adjustment model relates the vertical velocity v i j of each levelling benchmark j from profile i to the observed velocity difference ∆v i j and nearby GPS vertical velocity where d i (j, 0) is the curved distance between benchmark j and the reference benchmark of profile i. This formula gives an average model error of 0.037 mm/yr/ √ km. Equation 2 is a condition equation constraining the vertical velocities to be equal (within a given uncertainty) at sites common to two levelling profiles i and k. Because we want to avoid any error on a common site to propagate into the adjustment, and because the vertical velocity field is assummed to be smooth at the scale of a few km, we used equation 2 not only at strictly common sites, but at all sites within a certain radius (typically 5 or 10 km), with the weight decreasing with distance d betwwen the two levelling benchmarks: where σ 0 is a reference error value, d jl is the distance between benchmark j from profile k and benchmark l from profile i and d c a critical distance to impose the constraint. In order to choose σ 0 and d c , we checked the WRMS at each area where profiles cross. Equation 3 is a condition equation constraining the velocity to be equal (again, within a given uncertainty) between a levelling site with velocity v i j and the uplift rate of a nearby GPS site v GP S i j . Here the weight of the constraint is given by: where d is now the distance between the levelling benchmark and the nearby GPS site, σ GP S the uncertainty of the GPS vertical velocity, and σ GP S i j a reference level of constraint for levelling-GPS velocities. The chosen adjustment model offers the versatility to test different constraints for levelling-levelling data and GPS-levelling data, and search for outliers. It is further useful to express the levelling results in the GPS reference frame previously defined. Levelling results are often shown with respect to a fixed reference point and the obtained results are therefore highly dependent on the behaviour of this site. Here, equation 3 allows to express levelling results in the same reference (that is average in the least-squares sense) as the GPS results. This is optimal because, on the contrary to levelling results that can suffer from long length scale biases [15], GPS accuracy is largely independent from the scale of the network. Furthermore, we show in the previous paragraph that we could define a tectonically relevant and stable vertical reference frame that we use to express all results.

Levelling adjustment results
We use the previous combination model, with the very weak constraint σ GP S =1 000 mm/yr for GPS-levelling equality condition. This level of constraint is sufficient to ensure that levelling adjusted velocities are expressed in the same reference frame as the GPS results, without inducing any distortion of the intrinsic information provided by the levelling measurements.
We find that σ 0 = 0.2 mm/yr and d c = 4 km provides adjustments with good statistical indicators (histogram of normalized residuals close to the standard normal distribution). Extended Data Figure 5 shows the distribution of residuals. It shows that the average internal consistency among the profiles is~0.2 mm/yr, indicating that this value can be retained as an average indicator of the precision of the levelling results at the 1σ confidence level or~0.4 mm/yr at the 95% confidence level.

Combined GPS-levelling adjustment results -Accuracy estimates
The combined GPS-levelling adjustement has been obtained by progressively increasing the equality constraint of equation 3 by decreasing σ GP S i j . Ideally, if the GPS sites were colocated with levelling benchmarks and GPS and levelling would perfectly agree, we could decrease σ GP S i j to 0. We fixed d c to the average distance of GPS to the closest levelling benchmark (20 km) and find that σ GP S i j = 0.2 mm/yr provides residuals for equation 3 that are in agreement with the GPS uncertainties. Nonetheless, the combined adjustment also indicates areas of significant differences that are further enlighted in Figure 1 of the main text. In the Alpine foreland between longitude 5°E (the Rhone valley) and 6°E, levelling data indicate 0.5-1.0 mm/yr uplift rates that are not seen in the GPS data. To a lesser extent, a similar difference is found in the southern Alps. Finally, a larger subsiding area is found in the Rhone delta. The origin of these differences have yet to be understood. While GPS are located on rocks, levelling data are located along roads, often on sediments. Sediment compaction through time could explain the larger subsidence rates seen in the levelling results, but does not explain the 0.5 mm/yr larger levelling uplift rates recorded in the southern Alps and western pre-Alps. There is no systematic correlation with height that could be explained by a systematic bias in GPS or levelling. Given these results, we chose to interpret only (1) patterns that are consistent between the two techniques; (2) we chose 0.5 mm/yr as a threshold for deciding whether a motion is significant or not.

Map of uplift & models
The map of uplift were simply derived by interpolating the data points using the "surface" algorithm embedded in the GMT software [16](http://gmt.soest.hawaii.edu/). We used an L-curve [17] criterion to choose the tension factor providing the smoothest model explaining the data. The wrms for the chosen model is 0.17 mm/yr for the GPS velocities. For the levelling data, the average overall wrms 0.07 mm/yr. wrms reach 0.15-0.23 mm/yr for profiles 193,174,171,200, were the largest levelling-GPS rate differences are observed.
For ensuring that GIA and erosion physical models and geodetic results are in the same vertical reference frame, the models were corrected by the average uplift rates predicted at the GPS sites used to define the vertical reference frame used throughout the study.