Articles | Volume 11, issue 5
Research article
04 Sep 2020
Research article |  | 04 Sep 2020

New insights into active tectonics and seismogenic potential of the Italian Southern Alps from vertical geodetic velocities

Letizia Anderlini, Enrico Serpelloni, Cristiano Tolomei, Paolo Marco De Martini, Giuseppe Pezzo, Adriano Gualandi, and Giorgio Spada

This study presents and discusses horizontal and vertical geodetic velocities for a low strain rate region of the south Alpine thrust front in northeastern Italy obtained by integrating GPS, interferometric synthetic aperture radar (InSAR) and leveling data. The area is characterized by the presence of subparallel, south-verging thrusts whose seismogenic potential is still poorly known. Horizontal GPS velocities show that this sector of the eastern Southern Alps is undergoing ∼1 mm a−1 of NW–SE shortening associated with the Adria–Eurasia plate convergence, but the horizontal GPS velocity gradient across the mountain front provides limited constraints on the geometry and slip rate of the several subparallel thrusts. In terms of vertical velocities, the three geodetic methods provide consistent results showing a positive velocity gradient, of  1.5 mm a−1, across the mountain front, which can hardly be explained solely by isostatic processes. We developed an interseismic dislocation model whose geometry is constrained by available subsurface geological reconstructions and instrumental seismicity. While a fraction of the measured uplift can be attributed to glacial and erosional isostatic processes, our results suggest that interseismic strain accumulation at the Montello and the Bassano–Valdobbiadene thrusts it significantly contributing to the measured uplift. The seismogenic potential of the Montello thrust turns out to be smaller than that of the Bassano–Valdobbiadene fault, whose estimated parameters (locking depth equals 9.1 km and slip rate equals 2.1 mm a−1) indicate a structure capable of potentially generating a Mw>6.5 earthquake. These results demonstrate the importance of precise vertical ground velocity data for modeling interseismic strain accumulation in slowly deforming regions where seismological and geomorphological evidence of active tectonics is often scarce or not conclusive.

1 Introduction

Diffuse deformation, slow tectonic rates and long repeat times of large earthquakes make the estimate of seismic hazards at continental plate boundaries a challenging task. This is the case of the Italian Alps, a classic example of a broadly deforming continental collisional belt (Schmid et al., 2004). The present-day convergence between the Adriatic and Eurasian plates is largely accommodated in the eastern Southern Alps (ESA) (e.g., D'Agostino et al., 2005; Cheloni et al., 2014) where the Adriatic lithosphere underthrusts the Alpine mountain belt. This area is characterized by a notable seismic risk due to the high seismic hazard of northeastern Italy (, last access: 31 August 2020), high population density and widespread industrial activities. Geodesy indicates that the highest present-day deformation rates are localized along the mountain front of the Italian ESA (D'Agostino et al., 2005; Bechtold et al., 2009; Cheloni et al., 2014; Serpelloni et al. 2016), which is also the locus of several Mw>6 earthquakes (Rovida et al., 2016). For the Venetian sector of the ESA (between the Schio–Vicenza line and the Cansiglio plateau; see Fig. 1), however, historical and instrumental earthquake records indicate a lower seismic moment release rate with respect to its easternmost (Friuli) sector (Serpelloni et al., 2016). Here, discrepancies between geodetic and seismic deformation rates have been interpreted as being due to partial locking of the ESA thrust front (Cheloni et al., 2014; Serpelloni et al., 2016). However, despite geological or geomorphological evidence of Quaternary deformation (Benedetti et al., 2000; Galadini et al., 2005), the identification of active faults responsible for large past earthquakes is not conclusive. This area is also affected by non-tectonic deformation transients associated with the hydrological cycle in karst regions and precipitations, which are well tracked and described by continuous GPS stations (Devoti et al., 2015; Serpelloni et al., 2018). These deformation transients, which mainly affect the horizontal components of ground displacements, make the estimate of the long-term horizontal tectonic deformation rates more challenging.

Figure 1Seismotectonic map of the study area with the major tectonic lineaments shown in black (from Serpelloni et al., 2016). LM: Lessini Mountains; SV: Schio–Vicenza line; TBC: Thiene–Bassano–Cornuda thrust; MT: Montello thrust; BVT: Bassano–Valdobbiadene thrust; BL: Belluno thrust; VS: Valsugana thrust; CA: Cansiglio plateau; PM: Polcenigo–Maniago thrust; AR: Araba–Ragona thrust; PT: Pozzuolo thrust; MD: Medea thrust; FS: Fella–Sava line; PAF (PG): Periadriatic fault and Pusteria–Gailtal line. Major historical events (Mw>5 are the small squares and Mw>6 are the greater squares) from the parametric catalogue of Italian earthquakes (CPTI) (Rovida et al., 2016) are shown as yellow squares. The blue–purple symbols show instrumental seismicity for the 2000–2017 time span extracted from the National Institute of Oceanography and Applied Geophysics (OGS) bulletins (, last access: 31 August 2020). The focal mechanisms (from Anselmi et al., 2011; Danesi et al., 2015; Serpelloni et al., 2016) are plotted with the same color palette of the OGS instrumental seismicity. The small circles around the geodetic pole of rotation show the motion direction of Adria relative to Eurasia (Serpelloni et al., 2016). Insert map: location of the Euler rotation pole (star), the global navigation satellite system (GNSS) stations (yellow dots) used to define the Adria rotation pole in Serpelloni et al. (2016) and small circles of the Adria–Eurasia rotation pole with respect to the region of interest (black box).

Overall, in the European Alps, vertical GPS velocities correlate with topography (Serpelloni et al., 2013; Sternai et al., 2019) with widespread uplift at rates of up to  2.0–2.5 mm a−1 in the northwestern and central Alps and  1 mm a−1 across a continuous region from the eastern to the southwestern Alps. Proposed mechanisms of uplift include the isostatic response to the last deglaciation (19–11 Ka; Clark et al., 2012), long-term erosion, detachment of the western Alpine slab, and lithospheric and surface deflection due to subalpine asthenospheric convection (Serpelloni et al., 2013; Chery et al., 2016; Nocquet et al., 2016). In the study area, the isostatic response to the last deglaciation and long-term erosion is characterized by long-wavelength submillimeter per annum patterns (e.g., Sternai et al., 2012; Serpelloni et al., 2013), whereas the elastic response to present deglaciation (Barletta et al., 2006) shows localized peaks of uplift (up to the millimeter per annum scale) far from the study area. In the eastern Alps, active shortening contributes to the observed uplift, but the pattern and amount of tectonic uplift due to elastic strain accumulation at faults are still unknown.

In the study area, ∼1 mm a−1 of SE–NW shortening is accommodated in ∼40 km (Serpelloni et al., 2016), with the largest part being accommodated across the southernmost Thiene–Bassano–Cornuda, Montello and Cansiglio thrust faults (Cheloni et al., 2014; Serpelloni et al., 2016). However, horizontal GPS velocities in the Montello and Cansiglio areas appear “noisy”, providing a poorer fit to a single-fault model than in the Friuli sector of the ESA front (Cheloni et al., 2014; Serpelloni et al., 2016). Barba et al. (2013) used a sparse GPS velocity field for the Alps and vertical rates from precise leveling measurements performed by the Istituto Geografico Militare Italiano (IGMI) in the 1952–1984 time interval over a distance of ∼100 km from the Venetian plain to the Southern Alps mountain belt. The leveling measurements show a change from negative to positive vertical rates from the Venetian plain toward the mountains, which has been interpreted by Barba et al. (2013) as being associated with the elastic locking of the Bassano–Valdobbiadene thrust (BVT). However, the leveling line does not cross the Bassano–Valdobbiadene thrust, passing partially through the Montello thrust (MT) and then running along the Fadalto Valley (see Fig. 2), which is interpreted by Serpelloni et al. (2016) as an active left-lateral tectonic feature separating the Carnic from the Dolomitic Alps.

Figure 2Geodetic dataset: the black arrows show the horizontal GPS velocities (with 95 % error ellipses) in a fixed-Adria reference frame, whereas the colored circles show the vertical GPS velocities, where blue colors indicate negative (subsidence) values and red colors indicate positive (uplift) values. The colored diamonds show the vertical ground motion rates estimated from the IGMI leveling data in the 1952–1984 time interval, using the same color scales as the GPS velocities. The blue and red squares indicate the ascending and descending ENVISAT frames, respectively, used in this work. The dashed box shows the A–B cross section shown in Fig. 3.

In this work, we integrate interferometric synthetic aperture radar (InSAR) and GPS ground motion rate measurements across the MT and the BVT, the southernmost portion of the ESA south-verging fold-and-thrust belt (Fig. 1), with the goal of determining a new spatially dense three-dimensional interseismic velocity field, which considers also the available IGMI leveling measurements. We use two InSAR velocity fields obtained from ENVISAT satellite acquisitions along both ascending and descending orbits during 2004–2010 and a GPS velocity field obtained from the analysis of data from continuous and semicontinuous stations between 1998 and 2018 with most of the GPS stations in the study area active after 2005. The datasets and the procedures used to analyze and integrate the geodetic observations are described in Sect. 3. In order to provide insight into the origin of the geodetic uplift, we developed a two-dimensional dislocation model, which jointly inverts GPS, leveling and InSAR velocities along a NNW–SSE oriented profile crossing the MT and BVT, as described in Sect. 4. The results are presented in Sect. 5 and discussed in Sect. 6.

2 Tectonic setting

The Southern Alpine belt, which is part of the Alpine–Carpathian–Dinaridic system (AlCaDi; from Schmidt et al., 2008), is a SSE-vergent orogenic system formed during Neogene times by indenting one continent into another (continent–continent collision). Within this model, the rigid Adriatic microplate has indented obliquely into a softer region along the southern margin of the European plate (Ratschbacher et al., 1991a, b; Frisch et al., 1998). This indenter model is consistent with the overall geometry of the eastern Alps, with the region of maximum horizontal shortening being located north of the indenter front and with long, orogen-parallel structures that fan out toward the east and extrude the orogen eastward toward the Pannonian Basin. The topography of the eastern Alps also reflects the indenter tectonics causing crustal shortening, surface uplift and erosional response (Robl et al., 2017). The present-day convergence between the Adriatic and Eurasian plates is largely accommodated in the ESA where the Adriatic lithosphere underthrusts the Alpine mountain belt. From Serpelloni et al. (2016), the Adriatic microplate rotates counterclockwise at ∼0.3 Myr−1 around a pole of rotation located in the western Po Plain of Italy, implying rates of convergence between Adria and Eurasia between 1.5 and 2.0 mm a−1 in the study region (Fig. 1).

The Southern Alps are a typical example of a deformed passive continental margin in a mountain range (Bertotti et al., 1993). Until the Oligocene, this Adriatic domain was the gently deformed retro-wedge hinterland of the Alps, intensively reworked only at its eastern boundary by the Paleogenic Dinaric belt. From the Neogene (∼23 Ma), the southern Alpine fold-and-thrust belt developed and progressively propagated toward the Adriatic foreland, mainly reactivating Mesozoic ( 250–65 Ma) extensional faults (Castellarin et al., 1992). The study area (dashed box in Fig. 1) is located in the Venetian Southern Alps, which are part of the SSE-verging fold-and-thrust belt of the Southern Alps. Here the geometry of the fold-and-thrust belt is that of an imbricate fan with a shortening of 30 km at least (Doglioni, 1992, 1990). The main thrusts are, from the internal parts to the foreland, the Valsugana thrust (VS in Fig. 1), the Belluno thrust (BL in Fig. 1) and the Bassano–Valdobbiadene thrust (BVT in Fig. 1), the latter being associated with a morphological relief of ∼1200 m above the plain (Fig. 3) known as the Pedemountain flexure. The ESA southernmost front is now mainly buried beneath the alluvial deposits of the Venetian plain and sealed by Late Miocene to Quaternary ( 7–2.5 Ma) deposits (Fantoni et al., 2002), and it consists of the Thiene–Bassano–Cornuda thrust (TBC in Fig. 1) and the MT (Benedetti et al., 2000; Fantoni et al., 2002; Galadini et al., 2005; Burrato et al., 2008; Danesi et al., 2015; MT in Fig. 1). The Montello hill (see Fig. 3) is generally interpreted as an actively growing ramp anticline on top of an active, north-dipping thrust that has migrated south of the mountain into the foreland (Benedetti et al., 2000). Benedetti et al. (2000) suggested that the Piave River changed its course because of the growth of the Montello hill at a rate of 0.5 to 1.0 mm a−1 estimated from the ages of the river terraces. Modeling the deformation of the river terraces as a result of motion on the thrust ramp, Benedetti et al. (2000) estimated a constant slip rate of 1.8–2.0 mm a−1.

Figure 3On the left, cross section of GPS velocities, topography, instrumental seismicity and focal mechanisms included in the dashed box of Fig. 2 (and the right panel). Profile-parallel velocity components, in a fixed-Adria reference frame (see Fig. 2), are shown as red circles with 1σ error bars in gray. The dark gray areas show the average (median) topography in the profile swath, with the light gray and white areas showing the maximum and minimum elevations, respectively. The blue–purple circles show instrumental seismicity in the 2012–2017 time interval, taken from Romano et al. (2019), as a function of depth and scaled with magnitude. The gray and black focal mechanisms are from Anselmi et al. (2011) and Danesi et al. (2015), respectively. The continuous and dashed gray lines represent major and secondary faults digitized from the TRANSALP profile interpretation (modified from Fig. 11 of Castellarin et al., 2006). MT signifies the Montello thrust, MBT signifies the Montello backthrust, BV signifies the Bassano–Valdobbiadene thrust, BL signifies the Belluno thrust, and VS signifies the Valsugana thrust. On the right , map showing the instrumental seismicity recorded by the Collalto seismic network (Priolo et al., 2015) taken from Romano et al. (2019), the swath profile (dashed box), and the TRANSALP line (black line).

The ESA mountain front and foreland is the locus of several Mw≥6 earthquakes (green squares in Fig. 1) in the last 2000 years (Rovida et al., 2016), but most of the seismogenic faults associated with these large earthquakes are still debated (Database of Italian Individual Seismogenic Sources, DISS;, last access: 31 August 2020). Nevertheless, there is a general consensus for some of these earthquakes that the fault sources belong to the southernmost thrust fault system, emerging at the boundary between the Venetian plane and the mountain front (Benedetti et al., 2000; Galadini et al., 2005; Burrato et al., 2008; Cheloni et al., 2014). Information from historical sources and modern earthquake catalogues show that at least 12 earthquakes with Mw≥6 occurred in the past 1000 years in the Italian Southern Alps. The last one (Mw 6.4) took place in 1976 in Friuli with several strong aftershocks (5.1ML6.1) occurring in the following months (Slejko et al., 1999; Pondrelli et al., 2001). The largest earthquake in the study area is the 1695 Asolo earthquake, which is estimated, from the interpretation of historical information, to have been a Mw 6.4 event. However, the seismogenic source of this important earthquake is still uncertain. Uncertainties in the location, geometry and seismic potential of active faults in the Venetian Southern Alps persist due to (i) the complex structural framework of the region, inherited by the geodynamic evolution and resulting in deformation distributed over a vast area, (ii) slow deformation rates along mainly blind faults, and (iii) sparse instrumental seismicity causing some areas where strong historical earthquakes occurred to appear presently almost aseismic (Fig. 1).

3 Geodetic observation

3.1 GPS data

We use data from all the continuous GPS (cGPS) stations available for the study region provided by public institutions and private companies, integrated with data from a denser semicontinuous GPS (sGPS) network installed by the National Institute of Geophysics and Volcanology (INGV) in 2009 in the Venetian Southern Alps (Danesi et al., 2015). The solution presented here is part of a wider Euro-Mediterranean geodetic solution which includes ∼3000 cGPS stations that are analyzed in several subnetworks and later combined.

The GPS velocities have been estimated analyzing the position time series, realized in the IGb08 reference frame, of all stations having an observation time span longer than 5 years (in the 1998.0–2017.5 time span to be consistent with IGS08 products and the reference frame) in order to minimize possible biases in the linear trend estimation due to seasonal signals (Blewitt and Lavallée, 2002) and nonseasonal hydrological deformation signals (Serpelloni et al., 2018). The daily position time series were obtained by analyzing raw data with the GAMIT/GLOBK (version 10.60; Herring et al., 2015) and QOCA (, last access: 31 August 2020) software following the three-step procedure described in Serpelloni et al. (2006, 2013). Details on the GPS data processing and postprocessing are provided in the Supplement, and the reader is also referred to Serpelloni et al. (2013, 2016, 2018). The horizontal and vertical GPS velocities are provided in Table S1 of the Supplement. Figure 2 shows the horizontal GPS velocities (black arrows with 95 % confidence error ellipses) in a Adria-fixed reference frame using the rotation pole from Serpelloni et al. (2016) which has been calculated in a consistent IGb08 reference frame. The colored circles and diamonds instead show the vertical GPS and leveling rates, respectively, with white symbols representing stable sites (rates within ±0.5 mm a−1) and blue and red symbols representing points with negative (subsidence) and positive (uplift) rates, respectively.

Horizontal velocities show that the Veneto–Friuli plain is stable to within ±0.3 mm a−1 (see Fig. 2 and Table S1) with respect to the Adriatic microplate, whereas the ESA move toward the SSE, orthogonal to the main thrust front, at rates slower than 2 mm a−1. Figure 3 shows a velocity cross section, with topography and seismicity, obtained by projecting horizontal (Adria-fixed) GPS velocities along a NNW–SSE profile, highlighting the ∼1 mm a−1 shortening accommodated between the Venetian plane and the Dolomites. Consistent with results from Serpelloni et al. (2016), the area where most of the shortening is accommodated, and the horizontal-strain rate is higher, corresponds to the southernmost thrust front (i.e., the MT in this area). Regional kinematic models (Cheloni et al., 2014 and Serpelloni et al., 2016) show that part of the shortening across the southernmost ESA thrust front is likely accommodated by aseismic slip, although the spatial resolution of coupling maps is still poor. The horizontal velocity gradient across the Mt. Grappa–Mt. Cesen anticline is noisy, providing limited constraints on the presence of multiple locked faults across the mountain range.

Importantly, GPS and leveling vertical velocities coherently highlight a change from negative rates in the Venetian plain near the Adriatic coastlines to positive rates (up to ∼1.9 mm a−1) north of the Mt. Grappa–Mt. Cesen anticline (Fig. 2). Vertical velocities from leveling measurements refer to elevation changes measured during two campaigns in 1952 and 1984 by the Istituto Geografico Militare Italiano (see the Supplement for additional information), suggesting that the measured uplift is a steady feature over timescales longer than the GPS observational time span. In Sect. 3.4, this vertical velocity gradient is better detailed and discussed.

3.2 Synthetic aperture radar data

We use synthetic aperture radar (SAR) images (C-band sensor) acquired between 2004 and 2010 by the ENVISAT satellite, operated by the European Space Agency, along both ascending and descending orbits (Fig. 2). The displacement time series and the mean ground velocities in the line of sight (LOS) of the satellite geometry have been obtained by adopting the satellite-based augmentation system (SBAS) multitemporal InSAR technique (Berardino et al., 2002). The InSAR data processing has been performed using the SARscape module of ENVI software (, last access: 31 August 2020) to perform SBAS analysis. The SBAS algorithm is based on a combination of many SAR differential interferograms that are generated by applying constraints on the temporal and perpendicular baselines. The inversion of the interferometric phase, performed with a singular value decomposition (SVD) algorithm, produces a ground displacement time series for each coherent pixel by minimizing possible topographic, atmospheric and orbital artifacts; the achieved accuracy of the displacements can be up to ∼1 mm (Ferretti et al., 2007). In the case of tropospheric delays correlated with topography, we accounted for them as a phase term linearly correlated with elevation. The scaling coefficient is estimated directly from SAR interferograms using a network approach, i.e., exploiting the redundancy of the interferometric pairs to derive a joint estimation for the scaling factor of each interferogram (Elliott et al., 2008; Wang and Wright, 2012; Jolivet et al., 2013). Further details on InSAR data processing and postprocessing are provided in the Supplement.

For descending and ascending orbits, we analyzed 32 and 45 images, respectively (Table 1), selecting pairs of images and applying constraints on the maximum orbital separation and the temporal distance between the two passages in order to minimize spatial and temporal decorrelation effects. In particular, we chose 450 m and 600 d as the maximum values for the perpendicular and temporal baselines, respectively (Fig. S1 of the Supplement). We use the 90 m Shuttle Radar Topography Mission (available at, last access: 31 August 2020) digital elevation model for the topography subtraction step (Farr et al., 2007). The analysis of both ascending and descending datasets has been performed by applying a multi-looking factor equal to 20 for the azimuth and 4 for the range direction, resulting in a final resolution of 90 m on the ground. In order to avoid unwrapping issues and before the inversion steps, we checked each interferogram, discarding all the pairs showing clear unwrapping errors and keeping the ones with low atmospheric noise. Once the unwrapped phase inversion step is completed, the displacement time series (see the time span reported in Table 1) for each coherent pixel is obtained. Then the mean ground LOS velocities have been estimated by fitting with least-squares a linear trend to the displacement time series of each pixel. We obtain for most of the pixels low values of root mean square error and chi-square (see Figs. S2 and S3 of the Supplement), indicating a mainly linear behavior for most of the dataset. The mean error associated with the linear velocity estimate is ∼1 mm a−1.

Table 1ENVISAT dataset. Ground resolution values are obtained after a multi-looking operation.

Download Print Version | Download XLSX

3.3 GPS and InSAR integration method

The GPS and InSAR velocities have been compared projecting the three components of GPS velocities along the one-dimensional SAR satellite look direction (see Fig. S4). The two datasets show a loose agreement in the velocity pattern with some discrepancies between InSAR and GPS LOS rates (mainly at the border of the frames). These differences may be due to the different reference frame (which consists of a constant offset between the two datasets) and to residual orbital errors and long-wavelength atmospheric delays that are not completely removed from the InSAR mean velocity solution (e.g., Gabriel et al., 1989; Zebker et al., 1994). These last sources of errors are commonly represented by a bidimensional planar ramp in the case of short strips acquisitions, i.e., 1–2 frames (Feng et al., 2012) like in our case, or a quadratic ramp for larger InSAR images (e.g., Lohman and Simons, 2005; Biggs et al., 2007).

The difference between InSAR and GPS velocities are commonly used to correct InSAR velocities and put both geodetic datasets in the same reference frame (e.g., Zebker et al., 1994; Hammond et al., 2010). In particular, we minimize the differences between GPS and InSAR LOS rates by estimating and removing a ramp signal and a constant offset for the whole InSAR dataset. We use 11 and 9 GPS sites for the ascending and descending datasets, respectively, chosen among the sites with longer time spans (>7 years) and uniformly distributed as much as possible over each SAR frame (see Fig. 4). The InSAR velocities, for comparison with GPS LOS rates, are estimated by averaging the values from a certain number of pixels around each of the selected GPS stations and testing several radii of the selection. The optimal radius of selection is defined as the one minimizing the misfit (evaluated in terms of root mean square error, RMSE) between GPS and InSAR LOS rates at the selected stations after the ramp removal (see Fig. S5 in the Supplement for details). Before the correction, the RMSE deviations between GPS and InSAR data for both the ascending and descending datasets were 1.15 and 0.83 mm a−1, respectively. After the correction, the RMSE between GPS and InSAR datasets decreases to 0.55 and 0.63 mm a−1 for the ascending and descending datasets, respectively.

Figure 4Top panels: InSAR line-of-sight (LOS) velocities (after the ramp removal; see Sect. 3.3 for details) for the ascending (left) and descending (right) orbits, with negative (blue) and positive (red) values indicating increasing and decreasing distances between the Earth's surface and the satellite. Colored diamonds indicate the three-dimensional GPS velocities projected along the SAR LOS directions. Bottom panels: InSAR and GPS LOS velocities (red circles with 1σ error bars in gray) along the A–B cross section (dashed box) for the ascending and descending orbits, respectively. The blue line indicates the median value of the InSAR data along the cross section, given the minimum number of 50 SAR pixels in each bin of 1 km, and the light bars represent the data dispersion for each bin.

Figure 4 shows the corrected ascending and descending InSAR and GPS LOS velocities. The SE–NW (A–B in Fig. 4) cross sections show the overall good agreement in the LOS velocities between the two techniques. Both InSAR datasets show a positive velocity gradient of 2–3 mm a−1 in a few tens of kilometers across the BVT and a common area of negative values of 1–2 mm a−1 corresponding to the Venetian plain increasing toward the Adriatic coasts. Moreover, in the ascending dataset, we observe a small area with positive but noisy LOS velocities of ∼1 mm a−1 across the Montello hill that are mostly within the mean uncertainty value of InSAR velocities. Figure 4 shows an area of positive LOS velocities in the northeastern corner of both the ascending and descending orbits south of the Cansiglio plateau, which is more evident in the descending case. This area is located at the edge of both scenes and could be affected by errors associated with the orbital ramp estimation and removal during the processing chain (Berardino et al., 2002). In addition, during the filtering step (Goldstein et al., 1998, Ghulam et al., 2010) the presence of artifacts is possible because the adopted moving window filter size was 64 pixels large close to the frame border where less than 64 pixels remain. However, we cannot exclude that this positive signal may be of tectonic significance, being located near the mountain front where active thrust faults are present (Galadini et al., 2005), but the absolute values could be overestimated (as suggested by the disagreement between InSAR and GPS LOS rates near the global navigation satellite system (GNSS) station PORD, Fig. 2) and cannot be considered reliable. Thus, we do not consider this signal, which requires further investigations, in our analysis. Similarly, we excluded from the analysis the scattered pixels located in the northernmost sector of the ascending frame (Fig. 4) since they might be affected by unwrapping errors.

3.4 Integrated vertical velocity field

In order to estimate the vertical InSAR velocities, we combine the ascending and descending LOS velocities by assuming the N–S component as obtained from an interpolation of the GPS northward velocities in the area covered by the two ENVISAT frames. Given the near-polar sensor orbit direction, in fact, the N–S component of the InSAR velocities cannot be reliably derived. We calculate the East and Up components to solve the following system of equations:


where Asc and Dsc are the ascending and descending LOS rates, North is the N–S interpolated component, det=((e2du2d)-(u2ae2d)) and u2d, e2d, n2d, u2a, e2a, and n2a are the descending and ascending direction cosines, respectively, which account for the LOS angle variability along the respective SAR images swaths (Cianflone et al., 2015).

The vertical GPS and InSAR velocity field is shown in Fig. 5, together with vertical rates from the IGMI leveling measurements. Profile A–B runs normal to the strike of the MT and BVT faults, whereas the C–D profile runs along the leveling route (see the Supplement and Fig. S6 for details about how uncertainties in the leveling rates are obtained). The three geodetic observations, which refer to different temporal intervals (1952–1984 for the leveling, 2004–2010 for SAR and  2005–2017 for GPS data), show good agreement in the vertical ground motion rates across the two profiles. The Venetian plain, near the coast of the northern Adriatic Sea, shows negative rates (e.g., station MSTR) which decrease toward the north, being close to 0 south of the Montello hill across which a submillimeter per annum uplift is apparent from Fig. 5 but not significant considering the uncertainties and dispersion of the geodetic measurements. North of the Montello hill, GPS, InSAR and leveling measurements coherently show a positive vertical velocity gradient, reaching positive rates up to 2 mm a−1 in ∼20 km. However, while this vertical velocity gradient is rather well sampled by leveling data along the Fadalto Valley (profile C–D of Fig. 5), this is not the case across the Mt. Grappa–Mt. Cesen anticline, where a lack of InSAR observations is present due to the highly vegetated area of the Mt. Cesen–Mt. Grappa chain (profile A–B of Fig. 5).

Figure 5On the left, vertical velocity field from InSAR data decomposition (dots), from GPS measurements (circles) and from the IGMI leveling data (diamonds), in which blue colors indicate negative (subsidence) values and red colors indicate positive (uplift) values. On the right, cross sections of vertical velocities along the A–B (top) and the C–D (bottom) profile. The blue lines and red circles are the same as in Fig. 4.

Several studies investigated subsidence processes in the Venetian plain (e.g., Carminati and Di Donato, 1999; Carbognin et al., 2004; Teatini et al., 2005; Bock et al., 2012), which are due to three main causes (of both natural and anthropogenic non-tectonic origin): (1) aquifer compaction after the strong groundwater withdrawal in the second half of the last century (e.g., Gatto and Carbognin, 1981; Carbognin et al., 1995); (2) uncontrolled expansion of coastal settlements and industrial activities (e.g., Tosi et al., 2002); and (3) recent sediment compaction (e.g., Brambati et al., 2003; Fontana et al., 2008). As we can see from the profile A–B of Fig. 5, subsidence rates increase from the center of the plain towards the coasts due to the sum of the aforementioned processes. More localized faster subsidence (>3 mm a−1) is likely due to local anthropogenic processes mainly associated with groundwater exploitation (e.g., the area between stations SUSE and VITT along the profile C–D of Fig. 5). Discrepancies between GPS rates and InSAR or leveling vertical rates (e.g., stations MT02 and SUSE in Fig. 5) are likely due to known problems at GPS sites (e.g., station MT02 has been destroyed and rebuilt twice in a different place) or the onset of recent processes over the time interval of data availability, including anthropogenic processes; e.g., station SUSE is located above a gas storage area (, last access: 31 August 2020), and its data quality and time series are generally poorer than those at other sites. Although basically no InSAR observations are available along the Mt. Grappa–Mt. Cesen chain and on the Cansiglio plateau, the GPS stations there (i.e., MGRD, TAMB, CANV) show vertical positive rates that are less than 0.4 mm a−1. Importantly, these stations are known to be affected by nonseasonal hydrological deformation due to the pressurization of fractured karst aquifers (Devoti et al., 2015; Serpelloni et al., 2018), and their motion rates may be biased by these non-tectonic, nonseasonal, transient processes.

4 Dislocation modeling

In order to interpret the horizontal and vertical velocity gradients across the MT and BVT faults, we develop a two-dimensional fault model constructed from geologic cross sections, seismicity and geophysical prospections (see Fig. 3). The interseismic strain accumulation is represented by slip on buried rectangular dislocations embedded in an elastic homogeneous half-space (Okada, 1985) with a Poisson ratio of 0.25. The length of the dislocations along strike is kept long enough to avoid edge effects, and the width along the dip direction is exaggerated for the deepest fault in order to mimic the far-field long-term convergence rate. This approach has been widely used to simulate interseismic deformation associated with intracontinental thrust faults (e.g., Vergne et al., 2001; Hsu et al., 2003; Grandin et al., 2012; Tsai et al., 2012; Daout et al., 2016).

We perform an inversion of the geodetic velocities on a predetermined fault geometry (see Fig. 6) composed of a ramp–décollement system pertaining to the Montello thrust (MT) faults, which are constrained by recent instrumental seismicity (Romano et al., 2019; Fig. 3) and focal mechanisms (Anselmi et al., 2011; Danesi et al., 2015), and a fault ramp representing the Bassano–Valdobbiadene thrust (BVT), which is not well depicted by instrumental seismicity but better constrained by geological and geophysical interpretations (e.g., Galadini et al., 2005; Fantoni and Franciosi, 2010; see also Fig. 3), both connecting at depth along with a single deeply rooted thrust extending downdip (Castellarin et al., 2006). We keep the strike fixed perpendicular to the horizontal shortening rate direction (i.e., 65 N; see Fig. 2), the northward dip angles (see Fig. 6 for details) and the fault positions, while the width has been kept constant only for the Montello décollement and the deeper ramp. Assuming this fault geometry, we invert the geodetic velocities for determining two locking-depth (LD) values, which are the depth of the upper edge of the dislocation, indicating the shallow locked seismogenic fault portion for the MT and BVT ramps, and the dip-slip rates of all the fault segments (see Fig. 6).

Figure 6Sketch of the geometries for two-dimensional fault model as described in Sect. 4. LDM (locking depth of Montello fault), LDB (locking depth of Bassano fault), SRMR (slip rate of Montello ramp), SRMD (slip rate of Montello décollement), SRBR (slip rate of Bassano ramp) and SRDD (slip rate of the Deep ramp) are the fault parameters estimated from the data inversion (see Table 2). The dip angles of each fault segment (δ1=45, δ2=5, δ3=25, δ4=22) and the depth of fault junctions (z1=12 km, z2=15.5 km) have been kept fixed.


The inversion method exploits a constrained, nonlinear, derivative-based optimization algorithm (i.e., interior-point; see Byrd et al., 1999; Waltz et al., 2006). It allows us to estimate the optimal parameter solution corresponding to a possible global minimum of the cost function, representing the misfit between the model prediction and the geodetic measurements. These algorithms depend on the gradient and higher-order derivatives in order to guide them through misfit space; thus, they can get trapped in a local minimum (Cervelli et al., 2001) providing the best results when the starting point is near the global minimum. However, in order to ensure that we find a global solution in the inversion, we tested several different initial guesses, finding always the same model estimate. The cost function to be minimized is the weighted residual sum of squares (WRSS):

(3) WRSS = ( d obs - d mod ) T × cov - 1 × ( d obs - d mod ) ,

where cov is the covariance matrix of geodetic data errors and dmod=G(m), where G is the Green's function matrix depending on the fault geometry and slip parameters, m (Okada, 1985). The data covariance matrix cov is computed as follows: covRΣT, where Σ is the diagonal matrix of data uncertainty and R is the data correlation matrix, which is dimensionless and equal to 1 along the diagonal and the off-diagonal elements representing the correlation between each couple of data. Assuming the three geodetic datasets (GPS, InSAR and leveling) are independent among them, the whole covariance matrix is composed of three independent blocks, one for each dataset. The correlation values are nonzero only for the three components of each GPS site, considering the measurements obtained by the GPS stations are uncorrelated among them and for the leveling data, following the approach of Árnadóttir et al. (1992). The InSAR data covariance matrix is instead diagonal with an equal variance of 1 mm2 a−2 for all the pixels. The locking depths have been forced by the inversion algorithm to be within the seismogenic layer, as defined by instrumental seismicity (Fig. 3), while the dip-slip rates have been constrained to be kinematically consistent along the four fault planes. The latter constraint imposes the conservation of the horizontal and vertical motion along the fault junctions by means of the geometric relationships proposed by Daout et al. (2016).

In the inversion we use the horizontal and vertical velocities of 15 GPS sites (highlighted in Table S1 of the Supplement) located within a 35 km wide profile crossing the MT and BVT faults along a −25 N direction (dashed box of Fig. 2, perpendicular to the strike direction of the faults), selecting those stations reliable for site stability. We invert InSAR LOS rates from both ascending and descending datasets (corrected by the planar ramp), exploiting the completeness of information provided by the LOS velocities containing both the horizontal and the vertical contributions.

InSAR velocities within the same swath profile used to select GPS stations have been subsampled in order to reduce the density of pixels while maintaining the first-order ground deformation information (Fig. S7) and allowing for reasonable computational costs. Most of the studies performing a subsampling of InSAR data for geophysical inversions use the quadtree sampling method (e.g., Jónsson et al., 2002; Pedersen et al., 2003; Lohman and Simons, 2005; Jolivet et al., 2012; Maurer and Johnson, 2014) which allows us to reduce the amount of data in order to represent the statistically significant portion of the displacement signal (Jónsson et al., 2002). In our case, with low deformation gradients it is highly risky to apply a subsampling method that depends on the deformation signal itself. For this reason, we apply an alternative method that uniformly reduces the density of pixels, and the specific technical details are provided in the Supplement. The final amount of InSAR data has been reduced to less than 3000 points for each dataset (see Table S3 of the Supplement). The resampled pixels with mean ground velocity less than −0.5 mm a−1 are excluded from the modeling, assuming that subsidence is due to non-tectonic processes that we cannot take into account in this analysis (see Sect. 3.4). The total amount of resampled InSAR data (joining ascending and descending together) used during the inversion amounts to 3115 pixels. We also jointly invert the whole leveling dataset (see Table S2), calculating the associated covariance matrix as described in Árnadóttir et al. (1992) and excluding from the inversion the leveling points with vertical velocities less than −0.5 mm a−1 for the same reason mentioned above.

Given the larger amount of InSAR velocity data compared to the GPS and leveling data, we apply a relative weight, Wsar (varying between 0 and 1), which rescales the InSAR covariance matrix. This approach allows us to calibrate the relative importance of the InSAR dataset with respect to the others with the aim of finding the model that can best reproduce all the geodetic observations. The optimal value of Wsar (equal to 0.46) has been chosen as the “knee point” of the trade-off curves between the WRSS (Eq. 1) of GPS and InSAR data and the WRSS of leveling and InSAR data as the factor of Wsar varies (see Fig. S8 in the Supplement for details). Although the weighting factor affects the InSAR covariance matrix during the inversion, the trade-off curve analysis has been performed considering the original covariance matrix of InSAR data.

5 Results

The estimated model parameters (locking depths and dip-slip rates), with their uncertainties, and residuals in terms of RMSE are reported in Table 2. Figure 7 shows the measured and modeled horizontal and vertical geodetic velocities, as well as the InSAR ascending and descending LOS rates. The error bounds for the estimated optimal fault parameters have been evaluated using the bootstrap method (e.g., Efron and Tibshirani, 1986; Árnadóttir and Segall, 1994) with 1000 random resamplings of data. The bootstrap procedure allows us to estimate confidence intervals of the derived parameters (Segall and Davis, 1997) without making assumptions about the underlying statistics of errors (Amoruso and Crescentini, 2008), reflecting the limitations of the used dataset (Cervelli et al., 2001). Figure S9 in the Supplement shows the frequency distribution of the optimal fault models after the bootstrap procedure, which reports the corresponding 95 percentile confidence bounds (reported in Table 2) by means of the percentile method and the trade-off distributions between fault parameters pairs.

Figure 7Cross sections across the A–B profile in Fig. 2 showing the modeled (black lines) horizontal and vertical velocities, as well as the SAR ascending and descending LOS rates, along with the measured ones. Gray-colored GPS sites have not been used in the inversion unlike the red ones. Green points indicate leveling data and small blue dots represent the subsampled InSAR LOS rates used during the inversion, which are displayed along with the un-subsampled InSAR dataset (light gray points). The bottom panel reports the optimal fault geometry with dip-slip rate and locking depth estimates and 95 percentile confidence bounds (see Sect. 5 for details). BV: Belluno valley.


The inversion results show that the Bassano–Valdobbiadene thrust is characterized by a greater locking depth (LD = 9.1 km) and a faster dip-slip rate (2.1 mm a−1) than the Montello ramp fault (LD = 5.6 km, SR = 0.5 mm a−1). Moreover, the locking depth of the Bassano–Valdobbiadene thrust is better constrained by the data than the one estimated for the Montello ramp, which is weakly constrained because of the noisy submillimeter geodetic signal (close to the geodetic techniques limits) measured above the structure, to which the fault geometry is more sensitive. The trade-off distributions between parameter pairs (Fig. S9) show that the locking depth estimates have no correlation with the other parameters, while the dip-slip rates have strict correlations among them due to the kinematic conservation constraint. The only slip-rate parameter we can consider independent is the Deep ramp slip rate (SRDD in Fig. 6) which is representative of the far-field convergence rate and for which we obtain the widest confidence bounds (Table 2). The obtained results point out that the surface velocity gradient is mostly explained by the BVT, thus suggesting a greater seismogenic potential than that one expected for the Montello thrust fault.

The fault model overall reproduces the observed velocity gradients from all datasets (see Fig. 7), although we find a partial discrepancy between the leveling data and the modeled vertical velocity field (Table 2; Fig. 7). This misfit may be due to the fact that the leveling line from the Montello hill runs along the Fadalto valley (see Fig. 2) and those points might be sensitive to a deformation signal associated with the Cansiglio thrust fault located ∼10 km south of the BVT (Galadini et al., 2005). Regarding the Belluno Valley (see Fig. 7), both GPS and InSAR vertical rates show larger misfits with respect to the fault model. In particular, InSAR data show a steeper velocity gradient than that reproduced by the model, and this may be ascribed to some artifacts not totally estimated and removed during the SAR data processing and is likely due to the correlation between atmosphere and topography gradient (e.g., Elliott et al., 2008; Doin et al., 2009; Jolivet et al., 2012; Shirzaei and Bürgmann, 2012).

Table 2Output model parameters (with 95 percentile confidence intervals) and data residuals in terms of root mean squared errors.

Download Print Version | Download XLSX

6 Discussion

Overall, our results are consistent with the conclusions of Barba et al. (2013) who developed a two-dimensional finite-element model constrained by the leveling vertical rates which shows that the amount of interseismic locking associated with the BVT must exceed that associated with the MT. The shallower locking depth estimated for the Montello ramp is also consistent with other inferences on interseismic fault coupling for the ESA front, constrained by GPS horizontal velocities and indicating the MT as partially elastically locked (Cheloni et al., 2014; Serpelloni et al., 2016).

The Bassano–Valdobbiadene fault is classified as seismogenic and capable of Mw∼6.5 earthquakes in the DISS, but its late Quaternary activity is debated. No major earthquakes are unequivocally associated with this fault (see Fig. 1). The larger seismic event here, the Mw 6.4 1695 Asolo earthquake (see Fig. 1), is located at the foothills of the mountain range delimited by the BVT front but is in general associated with the westernmost portion of the Montello thrust, the Bassano–Corduna segment (Galadini et al., 2005; Burrato et al., 2008). Only recently have local seismic networks (Priolo et al., 2015) and temporary seismic experiments (Anselmi et al., 2011; Danesi et al., 2015) improved the imaging of faults at depth, indicating earthquakes with thrust focal mechanisms for both the Montello and Bassano–Valdobbiadene faults (Fig. 3).

During the 1970s and the 1980s, a large amount of data on the Pliocene–Quaternary fault activity of northeastern Italy have been collected. Castaldini and Panizza (1991) classified the BVT as active (i.e., showing tectonic movements in the middle Pleistocene–Holocene time). They proposed an uplift rate up to 1 mm a−1 based on evidence of deformation in Würmian deposits laying on structures associated with the BVT (Pellegrini and Zanferrari, 1980; Zanferrari et al., 1980), morphotectonic evidence of scarps aligned with the thrust, offsets recorded in karstified surfaces of Mindel–Riss age, and changes in the post-Würmian hydrography (Zanferrari et al., 1980). Galadini et al. (2001) draw an updated summary of the active faults in the eastern Southern Alps; i.e., they are responsible for the displacements of deposits and/or landforms related or subsequent to the last glacial maximum (LGM) set at 26–19 Ka at global scale following Clark et al. (2009) and more specifically at 25–21 Ka for the Alpine area as suggested by Monegato et al. (2017) and references therein. Galadini et al. (2001) produced a map including the main faults whose length is consistent with the occurrence of earthquakes with M≥6.2 (based on the rupture length–magnitude relationship by Wells and Coppersmith 1994) in which the BVT is represented as an unsegmented active fault from the Schio–Vicenza line to the Fadalto line. However, in the most recent literature (e.g., Galadini et al., 2005; Burrato et al., 2009; Moulin and Benedetti, 2018), evidence of recent activity of the BVT appears inconclusive. It is worth considering that the lack of evidence of recent displacements in the mountainous region may have been conditioned by the erosional–depositional evolution of the Alps since the LGM. Intense erosion by the glacial tongues and subsequent rapid deposition of thick successions of alluvial deposits may have hidden the evidence of fault activity, making geomorphological evidence of recent tectonic activity debatable.

On the contrary, the recent deformation of the Montello thrust is testified by several Middle and Upper Pleistocene warped terraces of the Piave River paleo-course flanking the western termination of the fold (Benedetti et al., 2000, and references therein) and by the eastward deflection of the Piave River around the growing anticline. Benedetti et al. (2000) associated with the MT at least three earthquakes with M<6 (778 CE, 1268 and 1859) which have been, however, revised in the Italian macroseismic catalogue (Rovida et al., 2016).

Geological and geomorphological data thus appear to not be fully consistent with the geodetic deformation signals and the interpretation presented in this work. The slip rates estimated for the flat and ramp portions of the Montello thrust are 0.4±0.2 and 0.5±0.3 mm a−1, respectively. These values are smaller than those estimated from modeling for the deformation of the river terraces by Benedetti et al. (2000), who estimated a constant slip rate of 1.8–2.0 mm a−1, values that are, however, incompatible with the actual submillimeter geodetic signal measured above the Montello hill. A possible explanation for the high rate obtained by Benedetti et al. (2000) is related to the new chronologies proposed for the Alpine LGM (as old as 25 Ka; Monegato et al., 2017) and for the Montebelluna megafan (125–30 Ka; Mozzi, 2005; Fontana et al., 2014), much older than the data available in 2000. The estimated slip rate for the BVT (2.1±0.3 mm a−1) is much greater than the one reported in the DISS database (0.3–0.7 mm a−1), which is inferred from regional tectonic and geological interpretations. To our knowledge, there are no other independent estimates of slip rates for this fault, for which conversely our estimates are rather well constrained by a significant geodetic signal. Considering a length of ∼30 km for each fault, a shear modulus of 24 GPa (as a mean value for the whole study region estimated from Vp and Vs values reported in Anselmi et al., 2011), and the estimated slip rates and locking depths (Table 2), the maximum magnitude expected in 1000 years as a recurrence time is Mw=6.9 for the BVT and Mw=6.2 for the MT. These values should be assumed as upper bounds because, given the simple model used, we do not take into account some complexities like the presence of lateral and vertical variations of fault geometry or rheological parameters that may prevent the coseismic rupture propagation and partially coupled fault portions that can reduce the seismogenic potential.

Despite the fault model providing a satisfactory explanation of the geodetic velocity gradients, it is worth considering that the vertical component of motion may be influenced by further ongoing processes on an alpine spatial scale. Indeed, GPS data show that the Alps are undergoing widespread uplift at rates up to  2–2.5 mm a−1 in the northwestern and central Alps and ∼1 mm a−1 across a continuous region from the eastern to the southwestern Alps (Serpelloni et al., 2013; Chen et al., 2018), highlighting a correlation with topography (Serpelloni et al., 2013). Proposed mechanisms of uplift include the isostatic response to the last deglaciation, to the present ice melting and to erosion, as well as the detachment of the western Alpine slab and lithospheric and surface deflection due to mantle convection (see Sternai et al., 2019, for a review). The relative contribution of each proposed mechanism to the geodetically measured uplift is, however, difficult to estimate, but their total contribution is of the order of a few tenths of millimeters per annum (Sternai et al., 2019). Moreover, it is also worth considering that models proposed to quantify these contributions (e.g., Barletta et al., 2006; Spada et al., 2009; Norton and Hampel, 2010; Sternai et al., 2012; Mey et al., 2016) are less reliable at the borders of the Alpine belt (Sternai et al., 2019), where our study area is located, making it even more speculative to apply a plausible “isostatic” correction to the actual uplift rates. The spatial patterns associated with the aforementioned processes have a characteristic wavelength of hundreds of kilometers, and correcting for these contributions would slightly reduce the intensity of uplift rates, leaving the observed velocity gradient mostly unchanged. In order to evaluate the effects of “correcting” the actual vertical velocities for long wavelength isostatic contributions on slip-rate and locking-depth estimates, we remove a linear gradient of the uplift rate along 100 km of distance from the measured geodetic rates. Despite being speculative and based on strong assumptions (see the Supplement for details), we can assume this approach as an upper-bound case to evaluate the impact of correcting or not for long-wavelength uplift signals on the inversion results. After the application of this correction, the locking-depth and slip-rate estimates are slightly reduced (see Fig. S10 in the Supplement) without significant variations, resulting widely within the error bounds of the optimal fault parameters (Table 2) and leading to exactly the same seismogenic potential estimates for the MT and BVT faults.

As these long-wavelength processes could hardly explain the local vertical velocity gradient observed across the Venetian mountain front, a possible source of more localized uplift rate gradient in the study area may come from the contribution of local ice or from present-day glacier shrinkage. At the LGM, the Belluno valley hosted the Piave glacier up to an elevation of ∼1150 m (Pellegrini et al., 2005; Carton et al., 2009). Pellegrini and Zambrano (1979) have estimated a width of ∼15 km for the ice pond, a length of ∼45 km and a maximum thickness of 600–800 m. Given the spatial correlation between the fastest vertical geodetic velocities and the location of the former Piave ice pond, we have employed a simple isostatic model with the aim of providing a first-order estimate of the spatial pattern and range of magnitudes for the isostatic response to the deglaciation of this local ice sheet. The isostatic response has been evaluated using the TABOO software (Spada et al., 2004, 2011) assuming a cap profile for the pond, an equivalent area of 45 km×15 km, a maximum thickness of 800 m (Pellegrini et al., 2005) and an instantaneous deglaciation occurring 10 kyr ago. This is an upper-bound case since recent studies (Pellegrini et al., 2005; Carton et al., 2009) have evaluated the complete disappearance of the glacier having likely occurred before 15 ka. The Earth model includes an elastic crust and a mantle viscosity profile consistent with the ICE-5G(VM2) model of Peltier et al. (2004). Figure S11 in the Supplement shows how the expected vertical velocities due to the Piave ice melting are sensitive to the assumed crustal thickness. Assuming in this region a crustal thickness between 30 and 40 km (Molinari et al., 2015), the model predicts modest glacial isostatic uplift rates (∼0.1 mm a−1).

Although we acknowledge the importance of glacial isostatic contributions to the present-day uplift of the eastern Southern Alps, the vertical velocity gradient identified by terrestrial and space geodetic measurements across the Venetian mountain front is mostly associated with active tectonic processes. Indeed, a simple elastic dislocation model based on seismological, geological and geophysical evidence can satisfactorily explain both the horizontal GPS velocities and the vertical GPS plus InSAR plus leveling rates consistently. However, important constraints from new seismicity data and new surface and subsurface geological observations will be required, together with denser global navigation satellite system data, to better constrain the tectonic rates and seismogenic potential in this area.

7 Conclusions

We study the interseismic deformation at the Adria–Alps boundary in northeastern Italy, focusing on a sector of the Southern Alps fold-and-thrust belt where WSW–ENE oriented subparallel, south-verging thrusts presently accommodate ∼1 mm a−1 of NW–SE shortening. Seismogenic fault databases indicate two possible sources for the larger historical earthquakes of the area: the Bassano–Valdobbiadene thrust (BVT), bordering the mountain front, and the Montello thrust (MT), buried beneath the alluvial Venetian plain. Both GPS and InSAR measurements, consistent with precise leveling measurements, show a positive vertical velocity gradient from the Venetian plain toward the Alps, with the highest rates, of the order of 1.5 mm a−1, located ∼20 km north of the mountain front. The three geodetic measurements span different time periods from 1950 to the present day, suggesting that this uplift pattern is consistent over a longer time interval than the most recent GPS measurements.

We develop an interseismic fault model constrained by InSAR, GPS and leveling rates, which estimates the slip rates and locking depths of the BVT and MT faults, whose geometry is defined by geological and seismological data. This simple model satisfactorily explains both horizontal and vertical geodetic velocities, suggesting that both thrust faults are accumulating elastic strain. The higher locking depth and slip rate estimated for the BVT, which is required to explain the observed uplift rate gradient, point toward a higher seismogenic potential of this fault (compatible with Mwmax>6.5 in 1 kyr), which according to recent results indicates a limited interseismic locking for the MT frontal ramp.

This work demonstrates the important constraints that vertical ground velocity measurements can provide for interseismic deformation studies associated with thrust faults in slowly deforming regions, as already shown in other faster deforming convergent boundaries. More importantly, since several processes contribute to mountain uplift, with the glacial isostatic contribution being the greatest one in the Alps, a more accurate estimate of the contribution of each process at local scale is necessary in order to correct the measured geodetic uplift for its non-tectonic components, guaranteeing a more reliable estimate of the geodetic slip rates. From this point of view, the eastern Southern Alps represent an ideal natural laboratory for further investigations.

Data availability

GPS and leveling velocity data that support the findings of this study can be found in the Supplement. The InSAR line-of-sight velocities obtained in this work are available through the B2SHARE service within the European EUDAT infrastructure at: (Anderlini et al., 2020).


The supplement related to this article is available online at:

Author contributions

ES processed the GPS data, CT and GP processed the InSAR data, and PMDM provided the leveling data. AG wrote the code for the InSAR data subsampling, and LA developed the data analysis, integrated and corrected all the datasets, and carried out the data modeling. GS and LA developed the GIA analysis. LA, ES, CT and GS interpreted the results and wrote the article. All the authors have contributed in the article preparation.

Competing interests

The authors declare that they have no conflict of interest.


This paper does not necessarily represent an official opinion or policy of the Dipartimento di Protezione Civile.


We are grateful to all public and private institutions and companies that make GPS data freely available for scientific applications. In particular, we acknowledge the EPN-EUREF network, Leica-ITALPOS, Topcon-NETGEO, ASI-GEODAF, INGV-RING, InOGS-FREDNET, Rete “Antonio Marussi” in Friuli Venezia Giulia, STPOS and TPOS in Trentino-Alto Adige and Regione Veneto, and the ARPAV-Belluno. Most of the figures have been created using the Generic Mapping Tools (GMT) software (Wessel and Smith, 1998). Mario Anselmi is thanked for fruitful discussions about the modeling and the analysis of results. We thank the reviewers Romain Jolivet, Edi Kissling and Giovanni Monegato for their constructive comments that allowed us to improve this article. We are grateful to the editors Taras Gerya and Charlotte Krawczyk for their editorial support.

Financial support

Letizia Anderlini is supported by INGV-MISE-DGS 2018 project funded by the Italian Ministry of the Economic Development. Giorgio Spada is funded by a FFABR (Finanziamento delle Attività Base di Ricerca) grant from MIUR (Ministero dell'Istruzione, dell'Università e della Ricerca), and by a DiSPeA research grant. This study has been partially supported by the project TRANSIENTI, by the MIUR “Premiale 2014” project, by the Italian Presidenza del Consiglio dei Ministri, Dipartimento di Protezione Civile (project no. S1 2013–2014 INGV-DPC agreement), and by the MIUR “Pianeta Dinamico” institutional project.

Review statement

This paper was edited by Taras Gerya and reviewed by Edi Kissling and Romain Jolivet.


Amoruso, A. and Crescentini, L.: Inversion of synthetic geodetic data for the 1997 Colfiorito events: clues on the effects of layering, assessment of model parameter PDFs, and model selection criteria, Ann. Geophys., 51, 461–475,, 2008. 

Anderlini, L., Serpelloni, E., Tolomei, C., De Martini, P. M., Pezzo, G., Gualandi, A., and Spada, G.: InSAR long-term velocities from Envisat satellites in the Venetian Southern Alps (Italy), EUDAT,, 2020. 

Anselmi, M., Govoni, A., De Gori, P., and Chiarabba, C.: Seismicity and velocity structures along the south-Alpine thrust front of the Venetian Alps (NE-Italy), Tectonophysics, 513, 37–48,, 2011. 

Árnadóttir, T. and Segall, P.: The 1989 Loma Prieta earthquake imaged from inversion of geodetic data, J. Geophys. Res.-Sol. Ea., 99, 21835–21855,, 1994. 

Árnadóttir, T., Segall, P., and Matthews, M.: Resolving the discrepancy between geodetic and seismic fault models for the 1989 Loma Prieta, California, earthquake, B. Seismol. Soc. Am., 82, 2248–2255, 1992. 

Barba, S., Finocchio, D., Sikdar, E., and Burrato, P.: Modelling the interseismic deformation of a thrust system: seismogenic potential of the Southern Alps, Terra Nova, 25, 221–227,, 2013. 

Barletta, V. R., Ferrari, C., Diolaiuti, G., Carnielli, T., Sabadini, R., and Smiraglia, C.: Glacier shrinkage and modeled uplift of the Alps, Geophys. Res. Lett., 33, L14307,, 2006. 

Bechtold, M., Battaglia, M., Tanner, D. C., and Zuliani, D.: Constraints on the active tectonics of the Friuli/NW Slovenia area from CGPS measurements and three-dimensional kinematic modeling, J. Geophys. Res., 114, B03408,, 2009. 

Benedetti, L., Tapponnier, P., King, G. C. P., Meyer, B., and Manighetti, I.: Growth folding and active thrusting in the Montello region, Veneto, northern Italy, J. Geophys. Res.-Sol. Ea., 105, 739–766,, 2000. 

Berardino, P., Fornaro, G., Lanari, R., and Sansosti, E.: A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms, IEEE T. Geosci. Remote Sens., 40, 2375–2383,, 2002. 

Bertotti, G., Picotti, V., Bernoulli, D., and Castellarin, A.: From rifting to drifting: tectonic evolution of the South-Alpine upper crust from the Triassic to the Early Cretaceous, Sediment. Geol., 86, 53–76,, 1993. 

Biggs, J., Wright, T., Lu, Z., and Parsons, B.: Multi-interferogram method for measuring interseismic deformation: Denali Fault, Alaska, Geophysical Journal International, 170, 1165–1179,, 2007. 

Blewitt, G. and Lavallée, D.: Effect of annual signals on geodetic velocity, J. Geophys. Res., 107, 2145,, 2002. 

Bock, Y., Wdowinski, S., Ferretti, A., Novali, F., and Fumagalli, A.: Recent subsidence of the Venice Lagoon from continuous GPS and interferometric synthetic aperture radar, Geochem. Geophy. Geosy., 13, Q03023,, 2012. 

Brambati, A., Carbognin, L., Quaia, T., Teatini, P., and Tosi, L.: The Lagoon of Venice: Geological Setting, Evolution and Land Subsidence, Episodes, 26, 264–268, 2003. 

Burrato, P., Poli, M. E., Vannoli, P., Zanferrari, A., Basili, R., and Galadini, F.: Sources of Mw 5+ earthquakes in northeastern Italy and western Slovenia: An updated view based on geological and seismological evidence, Tectonophysics, 453, 157–176,, 2008. 

Byrd, R. H., Hribar, M. E., and Nocedal, J.: An Interior Point Algorithm for Large-Scale Nonlinear Programming, SIAM J. Optimiz., 9, 877–900, 1999. 

Carbognin, L., Tosi, L., and Teatini, P.: Analysis of actual land subsidence in Venice and its hinterland (Italy), in: Land Subsidence, edited by: Barends, J. F., Brouwer, F. J., and Schröder, F. H., 129–137, A. A. Balkema, Rotterdam, 1995. 

Carbognin, L., Teatini, P., and Tosi, L.: Eustacy and land subsidence in the Venice Lagoon at the beginning of the new millennium, J. Marine Syst., 51, 345–353,, 2004. 

Carminati, E. and Di Donato, G.: Separating natural and anthropogenic vertical movements in fast subsiding areas: The Po plain (N. Italy) case, Geophys. Res. Lett., 26, 2291–2294,, 1999. 

Carminati, E. and Martinelli, G.: Subsidence rates in the Po Plain, northern Italy: the relative impact of natural and anthropogenic causation, Eng. Geol., 66, 241–255,, 2002. 

Carton, A., Bondesan, A., Fontana, A., Meneghel, M., Miola, A., Mozzi, P., Primon, S., and Surian, N.: Geomorphological evolution and sediment transfer in the Piave River watershed (north-eastern Italy) since the LGM, Géomorphologié, 3, 37–58, 2009. 

Castaldini, D. and Panizza, M.: Inventario delle faglie attive tra i fiumi Po e Piave e il lago di Como (Italia Settentrionale), Il Quaternario Italian Journal of Quaternary Sciences, 4, 333–410, 1991. 

Castellarin, A., Cantelli, L., Fesce, A. M., Mercier, J. L., Picotti, V., Pini, G. A., Prosser, G., and Selli, L.: Alpine compressional tectonics in the Southern Alps. Relationships with the N-Apennines, Annales Tectonicae, 6, 62–94, 1992. 

Castellarin, A., Vai, G. B., and Cantelli, L.: The Alpine evolution of the Southern Alps around the Giudicarie faults: A Late Cretaceous to Early Eocene transfer zone, Tectonophysics, 414, 203–223,, 2006. 

Cervelli, P., Murray, M. H., Segall, P., Amelung, F., Aoki, Y., and Kato, T.: Estimating source parameters from deformation data, with an application to the March 1997 earthquake swarm off the Izu Peninsula, Japan, J. Geophys. Res., 106, 11217–11237, 2001. 

Cheloni, D., D'Agostino, N., and Selvaggi, G.: Interseismic coupling, seismic potential, and earthquake recurrence on the southern front of the Eastern Alps (NE Italy), J. Geophys. Res.-Sol. Ea., 119, 4448–4468,, 2014. 

Chen, W., Braitenberg, C., and Serpelloni, E.: Interference of tectonic signals in subsurface hydrologic monitoring through gravity and GPS due to mountain building, Global Planet. Change, 167, 148–159,, 2018. 

Chery, J., Genti, M., and Vernant, P.: Ice cap melting and low-viscosity crustal root explain the narrow geodetic uplift of the Western Alps, Geophys. Res. Lett., 43, 3193–3200,, 2016. 

Cianflone, G., Tolomei, C., Brunori, C., and Dominici, R.: InSAR Time Series Analysis of Natural and Anthropogenic Coastal Plain Subsidence: The Case of Sibari (Southern Italy), Remote Sens., 7, 16004–16023,, 2015. 

Clark, P. U., Dyke, A. S., Shakun, J. D., Carlson, A. E., Clark, J., Wohlfarth, B., Mitrovica, J. X., Hostetler, S. W., and McCabe, A. M.: The Last Glacial Maximum, Science, 325, 710–714, 2009. 

Clark, P. U., Shakun, J. D., Baker, P. A., Bartlein, P. J., Brewer, S., Brook, E., Carlson, A. E., Cheng, H., Kaufman, D. S., Liu, Z., Marchitto, T. M., Mix, A. C., Morrill, C., Otto-Bliesner, B. L., Pahnke, K., Russell, J. M., Whitlock, C., Adkins, J. F., Blois, J. L., Clark, J., Colman, S. M., Curry, W. B., Flower, B. P., He, F., Johnson, T. C., Lynch-Stieglitz, J., Markgraf, V., McManus, J., Mitrovica, J. X., Moreno, P. I., and Williams, J. W.: Global climate evolution during the last deglaciation, P. Natl. Acad. Sci. USA, 109, E1134–E1142, 2012. 

D'Agostino, N., Cheloni, D., Mantenuto, S., Selvaggi, G., Michelini, A., and Zuliani, D.: Strain accumulation in the southern Alps (NE Italy) and deformation at the northeastern boundary of Adria observed by CGPS measurements, Geophys. Res. Lett., 32, L19306,, 2005. 

Danesi, S., Pondrelli, S., Salimbeni, S., Cavaliere, A., Serpelloni, E., Danecek, P., Lovati, S., and Massa, M.: Active deformation and seismicity in the Southern Alps (Italy): The Montello hill as a case study, Tectonophysics, 653, 95–108,, 2015. 

Daout, S., Barbot, S., Peltzer, G., Doin, M. P., Liu, Z., and Jolivet, R.: Constraining the kinematics of metropolitan Los Angeles faults with a slip-partitioning model, Geophys. Res. Lett., 43, 11192–11201,, 2016. 

Devoti, R., Zuliani, D., Braitenberg, C., Fabris, P., and Grillo, B.: Hydrologically induced slope deformations detected by GPS and clinometric surveys in the Cansiglio Plateau, southern Alps, Earth Planet. Sc. Lett., 419, 134–142,, 2015. 

Doglioni, C.: Thrust tectonics example from the venetian alps, Studi Geologici Camerti, special volume, 117–129, 1990. 

Doglioni, C.: The Venetian Alps thrust belt, in: Thrust Tectonics, Springer Netherlands, Dordrecht, 319–324,, 1992. 

Doin, M.-P., Lasserre, C., Peltzer, G., Cavalié, O., and Doubre, C.: Corrections of stratified tropospheric delays in SAR interferometry: Validation with global atmospheric models, J. Appl. Geophys., 69, 35–50, 2009. 

Efron, B. and Tibshirani, R.: Bootstrap Methods for Standard Errors, Confidence Intervals, and Other Measures of Statistical Accuracy, Stat. Sci., 1, 54–75,, 1986. 

Elliott, J. R., Biggs, J., Parsons, B., and Wright, T. J.: InSAR slip rate determination on the Altyn Tagh Fault, northern Tibet, in the presence of topographically correlated atmospheric delays, Geophys. Res. Lett., 35, L12309,, 2008. 

Fantoni, R. and Franciosi, R.: Tectono-sedimentary setting of the Po Plain and Adriatic foreland, Rend. Fis. Acc. Lincei, 21, 197–209,, 2010. 

Fantoni, R., Catellani, D., Merlini, S., Rogledi, S., and Venturini, S.: La registrazione degli eventi deformativi cenozoici nell'avampaese veneto-friulano, Mem. Soc. Geol. It., 57, 301–313, 2002. 

Farr, T. G., Rosen, P. A., Caro, E., Crippen, R., Duren, R., Hensley, S., Kobrick, M., Paller, M., Rodriguez, E., Roth, L., Seal, D., Shaffer, S., Shimada, J., Umland, J., Werner, M., Oskin, M., Burbank, D., and Alsdorf, D.: The Shuttle Radar Topography Mission, Rev. Geophys., 45, RG2004,, 2007. 

Feng, G., Ding, X., Li, Z., Mi, J., Zhang, L., and Omura, M.: Calibration of an InSAR-Derived Coseimic Deformation Map Associated With the 2011 Mw-9.0 Tohoku-Oki Earthquake, IEEE Geosci. Remote Sens. Lett., 9, 302–306,, 2012. 

Ferretti, A., Savio, G., Barzaghi, R., Borghi, A., Musazzi, S., Novali, F., Prati, C., and Rocca, F.: Submillimeter Accuracy of InSAR Time Series: Experimental Validation, IEEE T. Geosci. Remote Sens., 45, 1142–1153,, 2007. 

Fontana, A., Mozzi, P., and Bondesan, A.: Alluvial megafans in the Veneto-Friuli Plain: Evidence of aggrading and erosive phases during Late Pleistocene and Holocene, Quaternary Int., 189, 71–90, 2008. 

Fontana, A., Mozzi, P., and Marchetti, M.: Alluvial fans and megafans along the southern side of the Alps, Sediment. Geol., 301, 150–171, 2014. 

Frisch, W., Kuhlemann, J., Dunkl, I., and Brügel, A.: Palinspastic reconstruction and topographic evolution of the Eastern Alps during late Tertiary tectonic extrusion, Tectonophysics, 297, 1–15,, 1998. 

Gabriel, A., Goldstein, R., and Zebker, H.: Mapping Small Elevation Changes Over Large Areas: Differential Radar Interferometry, J. Geophys. Res., 94, 9183–9191,, 1989. 

Galadini, F., Meletti, C., and Vittori, E.: Major active faults in Italy: available surficial data, Neth. J. Geosci., 80, 273–296,, 2001. 

Galadini, F., Poli, M. E., and Zanferrari, A.: Seismogenic sources potentially responsible for earthquakes with M≥6 in the eastern Southern Alps (Thiene-Udine sector, NE Italy), Geophys. J. Int., 161, 739–762,, 2005. 

Gatto, P. and Carbognin, L.: The Lagoon of Venice: Natural environmental trend and man-induced modification, Hydrological Science Bulletin, 26, 379–391, 1981. 

Ghulam, A., Amer, R., and Ripperdan, R.: A filtering approach to improve deformation accuracy using large baseline, low coherence DInSAR phase images, 2010 IEEE International Geoscience and Remote Sensing Symposium, Honolulu, HI, IEEE, 3494–3497,, 2010. 

Goldstein, R. M. and Werner, C. L.: Radar Interferogram Filtering for Geophysical Applications, Geophys. Res. Lett., 25, 4035–4038,, 1998. 

Grandin, R., Doin, M.-P., Bollinger, L., Pinel-Puysségur, B., Ducret, G., Jolivet, R., and Sapkota, S. N.: Long-term growth of the Himalaya inferred from interseismic InSAR measurement, Geology, 40, 1059–1062,, 2012. 

Hammond, W. C., Li, Z., Plag, H. P., Kreemer, C., and Blewitt, G.: Integrated InSAR and GPS studies of crustal deformation in the Western Great Basin, Western United States, Int. A. of the Ph., Rem. Sens. and Spatial Information Sc., Kyoto Japan, XXXVIII, 39–43, 2010. 

Herring, T. A., King, R. W., Floyd, M. A., and McClusky, S. C.: Introduction to GAMIT/GLOBK, Release 10.6, available at: (last access: 31 August 2020), 2015. 

Hsu, Y.-J., Simons, M., Yu, S.-B., Kuo, L.-C., and Chen, H.-Y.: A two-dimensional dislocation model for interseismic deformation of the Taiwan mountain belt, Earth Planet. Sc. Lett., 211, 287–294,, 2003. 

Jolivet, R., Lasserre, C., Doin, M.-P., Guillaso, S., Peltzer, G., Dailu, R., Sun, J., Shen, Z.-K., and Xu, X.: Shallow creep on the Haiyuan Fault (Gansu, China) revealed by SAR Interferometry, J. Geophys. Res., 117, B06401,, 2012. 

Jolivet, R., Lasserre, C., Doin, M. P., Peltzer, G., Avouac, J. P., Sun, J., and Dailu, R.: Spatio-temporal evolution of aseismic slip along the Haiyuan fault, China: Implications for fault frictional properties, Earth Planet. Sc. Lett., 377–378, 23–33,, 2013. 

Jònsson, S., Zebker, H. A., Segall, P., and Amelung, F.: Fault slip distribution of the 1999 Mw 7.1 Hector Mine, California, earthquake, estimated from satellite radar and GPS measurements, B. Seismol. Soc. Am., 92, 1377–1389,, 2002. 

Lohman, R. and Simons, M.: Some thoughts on the use of InSAR data to constrain models of surface deformation: Noise structure and data downsampling, Geochem. Geophy. Geosy., 6, Q01007,, 2005. 

Maurer, J. and Johnson, K.: Fault coupling and potential for earthquakes on the creeping section of the central San Andreas Fault, J. Geophys. Res.-Sol. Ea., 119, 4414–4428,, 2014. 

Mey, J., Scherler, D., Wickert, A. D., Egholm, D. L., Tesauro, M., Schildgen, T. F., and Strecker, M. R.: Glacial isostatic uplift of the European Alps, Nat. Commun., 7, 13382,, 2016. 

Molinari, I., Verbeke, J., Boschi, L., Kissling, E., and Morelli, A.: Italian and Alpine three-dimensional crustal structure imaged by ambient-noise surface-wave dispersion, Geochem. Geophy. Geosy. 16, 4405–4421,, 2015. 

Monegato, G., Scardia, G., Hajdas, I., Rizzini, F., and Piccin, A.: The Alpine LGM in the boreal ice-sheets game, Sci. Rep.-UK, 7, 2078,, 2017. 

Moulin, A. and Benedetti, L.: Fragmentation of the Adriatic Promontory: New Chronological Constraints From Neogene Shortening Rates Across the Southern Alps (NE Italy), Tectonics, 105, 739–721,, 2018. 

Mozzi, P.: Alluvial plain formation during the Late Quaternary between the southern Alpine margin and the Lagoon of Venice (northern Italy), Geografia Fisica e Dinamica Quaternaria, 7, 219–230, 2005. 

Nocquet, J. M., Sue, C., Walpersdorf, A., Tran, T., Lenôtre, N., Vernant, P., Cushing, M., Jouanne, F., Masson, F., Baize, S., Chery, J., and van der Beek, P. A.: Present-day uplift of the western Alps, Sci. Rep.-UK, 6, 28404,, 2016. 

Norton, K. P. and Hampel, A.: Postglacial rebound promotes glacial re-advances – a case study from the European Alps, Terra Nova, 22, 297–302, 2010. 

Okada, Y.: Surface deformation due to shear and tensile faults in a half-space, B. Seismol. Soc. Am., 75, 1135–1154, 1985. 

Pedersen, R., Jònsson, S., Àrnadòttir, T., Sigmundsson, F., and Feigl, K. L.: Fault slip distribution of two June 2000 Mw 6.5 earthquakes in South Iceland estimated from joint inversion of InSAR and GPS measurements, Earth. Planet. Sc. Lett., 213, 487–502,, 2003. 

Pellegrini, G. B. and Zambrano, R.: Il corso del Piave a Ponte nelle Alpi nel Quaternario, Studi Trentini di Scienze Naturali, LVI, 69–100, 1979. 

Pellegrini, G. B. and Zanferrari, A.: Inquadramento strutturale ed evoluzione neotettonica dell'area compresa nei fogli 23-Belluno, 22-Feltre (p.p.) e 24-Maniago (p.p.), C.N.R. Prog. Final. “Geodinamica”: contributi alla realizzazione della carta neotettonica d'Italia, Pubbl. 356, 459–496, 1980. 

Pellegrini, G. B., Albanese, D., Bertoldi, R., and Surian, N.: La deglaciazione alpina nel Vallone Bellunese, Alpi meridionali orientali, Geografia Fisica e Dinamica Quaternaria, Supplemento 7, 271–280, 2005. 

Peltier, W. R.: Global glacial isostasy and the surface of the ice-age Earth: the ICE-5G (VM2) model and GRACE, Annu. Rev. Earth Planet. Sc., 32, 111–149,, 2004. 

Pondrelli, S., Ekstrom, G., and Morelli, A.: Seismotectonic re-evaluation ofthe 1976 Friuli, Italy, seismic sequence, J. Seismol., 5, 73–83,, 2001. 

Priolo, E., Romanelli, M., Plasencia Linares, M. P., Garbin, M., Peruzza, L., Romano, M. A., Marotta, P., Bernardi, P., Moratto, L., Zuliani, D., and Fabris, P.: Seismic Monitoring of an Underground Natural Gas Storage Facility: The Collalto Seismic Network, Seismol. Res. Lett., 86, 109–123,, 2015. 

Ratschbacher, L., Merle, O., Davy, P., and Cobbold, P.: Lateral extrusion in the eastern Alps, Part 1: Boundary conditions and experiments scaled for gravity, Tectonics, 10, 245–256,, 1991a. 

Ratschbacher, L., Frisch, W., Linzer, H. G., and Merle, O.: Lateral extrusion in the Eastern Alps, part 2: structural analysis, Tectonics, 10, 257–271,, 1991b. 

Robl, J., Heberer, B., Prasicek, G., Neubauer, F., and Hergarten, S.: The topography of a continental indenter: The interplay between crustal deformation, erosion, and base level changes in the eastern Southern Alps, J. Geophys. Res.-Earth, 122, 310–334,, 2017. 

Romano, M. A., Peruzza, L., Garbin, M., Priolo, E., and Picotti, V.: Microseismic Portrait of the Montello Thrust (Southeastern Alps, Italy) from a Dense High-Quality Seismic Network, Seismol. Res. Lett., 90, 1502–1517,, 2019. 

Rovida, A., Locati, M., Camassi, R., Lolli, B., and Gasperini, P. (Eds.): Italian Parametric Earthquake Catalogue (CPTI15), Istituto Nazionale di Geofisica e Vulcanologia (INGV), Rome (Italy),, 2016. 

Schmid, S. F., Genschuh, B., Kissling, E., and Schuster, R.: Tectonic map and overall architecture of the Alpine orogen, Eclogae Geol. Helv., 97, 93–117,, 2004. 

Segall, P. and Davis, J. L.: GPS applications for geodynamics and earthquake studies, Annu. Rev. Earth Planet. Sc., 23, 201–336, 1997. 

Serpelloni, E., Casula, G., Galvani, A., Anzidei, M., and Baldi, P.: Data analysis of permanent GPS networks in Italy and surrounding regions: application of a distributed processing approach, Ann. Geophys.-Italy, 49, 897–928, 2006. 

Serpelloni, E., Faccenna, C., Spada, G., Dong, D., and Williams, S. D. P.: Vertical GPS ground motion rates in the Euro-Mediterranean region: New evidence of velocity gradients at different spatial scales along the Nubia-Eurasia plate boundary, J. Geophys. Res.-Sol. Ea., 118, 6003–6024,, 2013. 

Serpelloni, E., Vannucci, G., Anderlini, L., and Bennett, R. A.: Kinematics, seismotectonics and seismic potential of the eastern sector of the European Alps from GPS and seismic deformation data, Tectonophysics, 688, 157–181,, 2016. 

Serpelloni, E., Pintori, F., Gualandi, A., Scoccimarro, E., Cavaliere, A., Anderlini, L., Belardinelli, M. E., and Todesco, M.: Hydrologically-induced karst deformation: insights from GPS measurements in the Adria-Eurasia plate boundary zone, J. Geophys. Res.-Sol. Ea., 85, 457,, 2018. 

Shirzaei, R. and Bürgmann, R.: Topography correlated atmospheric delay correction in radar interferometry using wavelet transfoms, Geophys. Res. Lett., 39, L01305,, 2012. 

Slejko, D., Neri, G., Orozova, I., Renner, G., and Wyss, M.: Stress field in Friuli (NE Italy) from fault plane solutions of activity following the 1976 main shock, B. Seismol. Soc. Am., 89, 1037–1052, 1999. 

Spada, G., Antonioli, A., Boschi, L., Brandi, V., Cianetti, S., Galvani, G., Giunchi, C., Perniola, B., Piana Agostinetti, N., Piersanti, A., and Stocchi, P.: Modeling Earth's post-glacial rebound, EOS T. Am. Geophys. Un., 85, 62–64,, 2004. 

Spada, G., Stocchi, P., and Colleoni, F.: Glacio–isostatic Adjustment in the Po Plain and in the Northern Adriatic Region, Pure Appl. Geophys., 166, 1303–1318,, 2009. 

Spada, G., Barletta, V. R., Klemann, V., Riva, R. E. M., Martinec, Z., Gasperini, P., Lund, B., Wolf, D., Vermeersen, L. L. A., and King, M. A.: A benchmark study for glacial isostatic adjustment codes, Geophys. J. Int., 185, 106–132,, 2011. 

Sternai, P., Herman, F., Champagnac, J.-D., Fox, M., Salcher, B., and Willett, S. D.: Preglacial topography of the European Alps, Geology, 40, 1067–1070, 2012.  

Sternai, P., Sue, C., Husson, L., Serpelloni, E., Becker, T. W., Willett, S. D., Faccenna, C., Di Giulio, A., Spada, G., Jolivet, L., Valla, P., Petit, C., Nocquet, J.-M., Walpersdorf, A., and Castelltort, S.: Present-day uplift of the European Alps: Evaluating mechanisms and models of their relative contributions, Earth-Sci. Rev., 190, 589–604,, 2019. 

Teatini, P., Tosi, L., Strozzi, T., Carbognin, L., Wegmüller, U., and Rizzetto, F.: Mapping regional land displacements in the Venice coastland by an integrated monitoring system, Remote Sens. Environ., 98, 403–413,, 2005. 

Tosi, L., Carbognin, L., Teatini, P., Strozzi, T., and Wegmüller, U.: Evidence of the present relative land stability of Venice, Italy, from land, sea, and space observations, Geophys. Res. Lett., 29, 1562,, 2002. 

Tsai, M.-C., Yu, S.-B., Hsu, Y.-J., Chen, H.-Y., and Chen, H.-W.: Interseismic crustal deformation of frontal thrust fault system in the Chiayi–Tainan area, Taiwan, Tectonophysics, 554–557, 169–184,, 2012. 

Vergne, J., Cattin, R., and Avouac, J. P.: On the use of dislocations to model interseismic strain and stress build-up at intracontinental thrust faults, Geophys. J. Int., 147, 155–162, 2001. 

Waltz, R. A., Morales, J. L., Nocedal, J., and Orban, D.: An interior algorithm for nonlinear optimization that combines line search and trust region steps, Math. Programm., 107, 391–408, 2006. 

Wang, H. and Wright, T. J.: Satellite geodetic imaging reveals internal deformation of western Tibet, Geophys. Res. Lett., 39, L07303,, 2012. 

Wells, D. L. and Coppersmith, K. J.: New empirical relationships among magnitude, rupture length, rupture width, rupture area, and surface displacement, B. Seismol. Soc. Am., 84, 974–1002, 1994. 

Wessel, P. and Smith, W. H. F.: New, improved version of Generic Mapping Tools released, EOS T. Am. Geophys. Un., 79, p. 579,, 1998. 

Zanferrari, A., Pianetti, F., Mattana, U., Dall'Arche, L., and Tonielli, V.: Evoluzione neotettonica e schema strutturale dell'area compresa nei Fogli 38 – Conegliano, 37 – Bassano del Grappa (p.p.) e 39 – Pordenone (p.p.), in: CNR (1980) – Contributi alla realizzazione della Carta Neotettonica d'Italia, P.F. Geodinamica, Pubbl. 356, 397–435, 1980. 

Zebker, H. A., Rosen, P. A., Goldstein, R. M., Gabriel, A., and Werner, C. L.: On the derivation of coseismic displacement fields using differential radar interferometry: the Landers earthquake, J. Geophys. Res., 99, 19617–19634,, 1994. 

Short summary
The Venetian Southern Alps (Italy) are located in a slowly deforming plate-boundary region where strong earthquakes occurred in the past even if seismological and geomorphological evidence is not conclusive about the specific thrust faults involved. In this study, we integrate and model different geodetic datasets of ground velocity to constrain the seismogenic potential of the studied faults, giving an example of the importance of using vertical geodetic data for seismic hazard estimates.