Measurement(s)ground deformation processTechnology Type(s)GPS navigation systemFactor Type(s)GNSS station • temporal intervalSample Characteristic - LocationAdriatic SeaMachine-accessible metadata file describing the reported data: 10.6084/m9.figshare.13034372
We provide a dataset of 3D coordinate time series of 37 continuous GNSS stations installed for stability monitoring purposes on onshore and offshore industrial settlements along a NW-SE-oriented and ~100-km-wide belt encompassing the eastern Italian coast and the Adriatic Sea. The dataset results from the analysis performed by using different geodetic software (Bernese, GAMIT/GLOBK and GIPSY) and consists of six raw position time series solutions, referred to IGb08 and IGS14 reference frames. Time series analyses and comparisons evidence that the different solutions are consistent between them, despite the use of different software, models, strategy processing and frame realizations. We observe that the offshore stations are subject to significant seasonal oscillations probably due to seasonal environmental loads, seasonal temperature-induced platform deformation and hydrostatic pressure variations. Many stations are characterized by non-linear time series, suggesting a complex interplay between regional (long-term tectonic stress) and local sources of deformation (e.g. reservoirs depletion, sediment compaction). Computed raw time series, logs files, phasor diagrams and time series comparison plots are distributed via PANGAEA (https://www.pangaea.de).
Background & SummaryWe present the results from the analysis of a new dataset of GNSS stations installed on onshore and offshore industrial settlements provided by Eni S.p.A. (https://www.eni.com). This is one of the research activities started by INGV in the framework of the program “CLYPEA: Innovation network for future energy”, as part of the “subsoil deformations” project1. Such a dataset has been collected at continuous GNSS stations installed, respectively closely to 13 onshore infrastructures (e.g. storage centers) located along the Italian Adriatic coast and 24 offshore hydrocarbon production platforms anchored to the seabed in central and northern Adriatic Sea (Fig. 1a).Fig. 1(a) Map distribution of the continuous GNSS stations provided by Eni S.p.A. Legend: 1) GNSS stations installed on storage and treatment centers, 2) GNSS stations installed on offshore platforms, 3) continuous GNSS stations developed and managed by different local, national and international Institutions and Agencies (INGV, ASI, ITALPOS, NETGEO, etc), 4) Oil/Gas prospection and production concession. (b) Temporal raw data availability of the GNSS dataset (see Table 1 for additional details).This dataset improves the current density of GNSS stations along the Adriatic coastal belt, allowing also to capture local and/or regional crustal deformation on a large sector of the northern Adriatic Sea, therefore representing a potentially invaluable dataset, although it involves GNSS stations not realized for geophysical purposes. Moreover, due to its offshore extension, the dataset must be considered as unique since until now, with the exclusion of very few GNSS stations (such as the HARV station installed on the Harvest platform, approximately located 10 km off the coast of central California2), not so extensive offshore datasets have been publicly released.The raw GNSS dataset covers a period ranging from 1998 to the end of 2017, with time series covering time spans from 0.71 to 19.08 years (Table 1 and Fig. 1b). Raw data were provided by Eni S.p.A. in the RINEX format with 30-second sampling rate and containing 24 hours (00:00–23:59 UTC) of continuous tracking of NAVSTAR GPS (NAVigation Satellite Timing And Ranging Global Positioning System) satellites. In addition, Eni S.p.A. provided also a few ancillary information related to adopted receiver and antenna models and their changes over the acquisition time.Table 1GNSS dataset provided by Eni S.p.A. For each station, site id, denomination, notes (P, Platform; STC; Storage and treatment center), coordinates, installation date and duration are reported.SiteNameNoteLong. (°)Lat. (°)Heigth (meters)StartingDuration (years)AGOBAGOSTINO-BP12.47244.55467.132003.7814.22AMEBAMELIA-BP12.66244.40773.261998.9219.08ANBLANNABELLAP13.07944.22976.632007.6010.40ANEBANEMONE-BP12.70544.22963.542007.6010.40ANGAANGELA-ANGELINAP12.34444.39173.341998.9219.08ANLSANNALISAP13.11444.17163.082007.7610.24ANTAANTARESP12.45444.39073.112005.5212.03ANTOANTONELLAP12.77744.21468.062007.6010.37ARIAARIANNA-AP12.62844.30668.842004.4213.58ARMIARMIDAP12.45344.48073.032008.929.07AZABAZALEA-BP12.72044.16778.152005.5312.47BRISBRISIGHELLASTC11.77444.225220.682007.0011.00CERACERVIA-AP12.63944.29467.142003.7814.22CLAECLARA-ESTP14.07243.78065.872007.8510.02CROTCROTONESTC17.10639.10545.782002.7715.23ELEOELEONORAP14.15642.84072.532012.264.99EMILEMILIOP14.24342.93566.182009.198.81FALCFALCONARASTC13.35843.64046.382007.9510.05FANOFANOSTC13.04143.80951.392007.9410.06FAUZFAUZIAP13.55444.05667.002014.743.26FIUNFIUMI UNITISTC12.31644.39743.592002.2015.80FRATFRATELLO CLUSTERP14.16842.61263.692007.000.71GARBGARIBALDI-BP12.53244.48766.392003.7714.23GARCGARIBALDI-CP12.51544.53172.602003.7714.23GRTMGROTTAMMARESTC13.84542.98070.572009.508.50ITGSITALGASSTC13.74443.245105.912009.928.08MANAMANARASTC12.22644.75044.192005.5312.47NAIDNAIDEP12.74544.34362.332005.4912.51NAPANAOMI-PANDORAP12.84744.68968.832002.4615.53PCCPPORTO CORSINI M E CP12.56144.39171.992005.5212.48PCTASPINARONISTC12.26744.49545.162002.2815.72PCWCPORTO CORSINI M W CP12.37344.50972.422003.2714.73PINEPINETOSTC14.01442.64455.972009.508.50RUBIRUBICONESTC12.40944.15047.302007.9510.05SMARSMARLACCASTC12.24744.59444.432002.2815.72TEALTEA-LAVANDA-ARNICAP13.01944.50166.162007.8810.12TREVTREVISOSTC12.24645.66869.712008.049.96Here we provide 3D coordinate time series on a daily basis, as processed by different data analysis centers, providing a robust validation of this particular set of positions time series. Coordinates are consistent with IGb083 and IGS144 reference frames. Both reference frames have origin coinciding with the Earth system center of mass, however, the IGS14 and the IGb08 represent the GNSS realization of the ITRF20145 and ITRF20086 reference frames, respectively. The data files are distributed in text ASCII format (i.e. the “pos” format7, realized by the Plate Boundary Observatory for GNSS time series solutions sharing) via a web-based online repository.MethodsThis section illustrates the analyses carried out by different data processing centers, which analyzed the same dataset by means of different geodetic software and, within the same software adopted different procedures for the reference frames realization. With the aim of producing the most accurate final computation and, at the same time, obtaining a comparable analysis strategy, each data analysis center processed the raw GNSS data by using (i) a double difference approach, applied using the Bernese software, (ii) a double-difference distributed approach, applied by means of the GAMIT/GLOBK software package and (iii) a Precise Point Positioning (PPP) approach, applied by means of GIPSY.The GNSS raw dataset has been processed at four different analysis centers:Istituto Nazionale di Geofisica e Vulcanologia - Sezione di Bologna (hereinafter INGV-BO); this analysis center provided 1 solution processed with GAMIT/GLOBK and QOCA (named IBO_GAMIT).Istituto Nazionale di Geofisica e Vulcanologia - Osservatorio Etneo (hereinafter INGV-OE); this analysis center provided 1 solution processed with GAMIT/GLOBK (named IOE_GAMIT).Istituto Nazionale di Geofisica e Vulcanologia - Osservatorio Nazionale Terremoti (hereinafter INGV-ONT); this analysis center provided 1 solution processed with GIPSY (named ONT_GIPSY) and 1 solution processed with Bernese (named ONT_BERNESE).Dipartimento di Ingegneria Civile, Chimica, Ambientale e dei Materiali - Università di Bologna (hereinafter DICAM-UniBO); this analysis center provided 1 solution processed with GIPSY (named UBO_GIPSY) and 1 solution processed with GAMIT/GLOBK (named UBO_GAMIT).In the following, a brief description of each adopted procedure for obtaining daily positions is provided. Table 2 resumes some relevant characteristics of each processing scheme, including specific models applied in the analysis.Table 2Data processing parameters adopted during the processing of the GNSS dataset provided by Eni S.p.A.SolutionUBO_GAMITUBO_GIPSYIBO_GAMITIOE_GAMITONT_GIPSYONT_BERNESEData processing parametersSoftwareGamit (Ver. 10.61)GipsyX (Ver. rc0.4)Gamit (Ver. 10.70) + QOCAGamit/Globk (Ver. 10.70)Gipsy (ver. 6.3)Bernese (ver 5.0)Processing sampling rate30 sec300 sec30 sec30 sec30 sec30 secElevation mask10°10°10°10°0°10°Antenna phase center modelsigs08.atxigs14_wwww.atxigs08.atxigs14_wwww.atxigs14_wwww.atxigs08.atxGPS orbitsIGS, finalJPL, finalIGS, finalIGS, finalJPL, finalrepro2, IGS, finalOcean tides loading modelFES2004FES2004FES2004FES2004FES2004FES2004Tropospheric mapping functionVMF1VMF1VMF1VMF1VMF1NIELLIonosphere correction2nd order ionosphere correction2nd order ionosphere correction2nd order ionosphere correction2nd order ionosphere correction2nd order ionosphere correctionno 2nd order correctionAmbiguity resolutionYesYes, with wlpb filesYesYesYes, with wlpb filesYesReference frame/system and strategy for alignment of solutionsAlignment strategy7 parameter transformation7 parameter transformation7 parameter transformation7 parameter transformation7 parameter transformation4 parameter transformationReference frame/systemIGb08IGb08IGb08IGb08IGS14IGb08Sinex fileIGb08.snxIGb08.snxIGb08.snxIGb08.snxIGb08.snxReference stationsBUCU, GRAS, GRAZ, MATE, MEDI, NOT1, SOFI, WTZR, ZIMMBUCU, GRAS, GRAZ, MATE, MEDI, NOT1, SOFI, WTZR, ZIMMIGb08 core sitesBRUS, BRUX, COMO, GENO, GRAS,GRAZ, JOZE, MATE, NOT1, PRAT, TORI, ZIMM, ZOUF132 stations used for reference frame alignment45 anchor stations in EuropeRaw data processing by using the Bernese softwareThe Bernese8 GNSS software is a scientific, high-precision, multi-GNSS data processing software developed at the Astronomical Institute of the University of Bern. Data processing with version 5.0 of this software package was performed at INGV-ONT. Daily solutions of station positions were estimated in loosely constraint reference systems. The raw observations were processed forming Ionosphere Free linear combinations and solving for the troposphere biases and phase ambiguities using the Quasi Ionosphere-free approach9. The troposphere modeling consisted in an a priori dry-Niell model fulfilled by the estimation of zenith delay corrections at 1-hour intervals at each site using the wet-Niell mapping function (see Table 2). In addition, one horizontal gradient parameter per day at each site was estimated. Ocean loading was computed using the FES2004 tidal model coefficients as provided by the Ocean Tide Loading provider (http://holt.oso.chalmers.se/loading). The GPS orbits and the Earth’s orientation parameters were fixed to the final IGS products and the site coordinates were constrained to an apriori sigma of 10 m, thus the daily coordinates were estimated in a loosely constrained, unknown reference frame. In order to express the ONT_BERNESE solution in a unique reference frame, the daily covariance was first projected imposing tight internal constraints (at millimeter level), and then the coordinates were transformed into the ITRF2008 reference frame by a 4-parameter Helmert transformation (translations and scale factor). The regional reference frame transformation used 45 anchor sites among the IGb08 stations, located on the Eurasian Plate (see Table 2).Raw data processing by using the GAMIT/GLOBK softwareGAMIT/GLOBK10 is a GNSS analysis package, designed to run under any UNIX operating system and developed at the Massachusetts Institute of Technology (MIT), the Harvard-Smithsonian Center for Astrophysics, Scripps Institution of Oceanography and Australian National University. This software was adopted at DICAM-UniBO (version 10.61) and at INGV-BO and INGV-CT (version 10.7). All the three analysis centers adopted the IGS “Repro2 campaign”11 standards during the raw dataset processing. In this step, INGV-CT and DICAM-UniBO included into the cluster processing some high-quality IGS12 stations (Table 2) in order to improve the overall configuration of the network and to tie the regional stations to an external global reference frame. The INGV-BO solution is part of a continental-scale geodetic analysis, including >3000 continuous GNSS stations13,14. Because of this large number of sites, the GAMIT analysis was performed independently for several sub-networks, each made by <50 stations, with each sub-network sharing a set of high-quality IGS stations, which are used as tie-stations in the combination step.During the processing step, the GAMIT software uses an ionosphere-free linear combination of GNSS phase observables by applying a double differencing technique to eliminate phase biases related to drifts in the satellite and receiver clock oscillators. GPS phase data were weighted according to an elevation-angle-dependent error model10 using an iterative analysis procedure whereby the elevation dependence was determined by the observed scatter of phase residuals. In this analysis the parameters of satellites orbit were fixed to the IGS final products. IGS absolute antenna phase center models (igs08.atx and igs14_wwww.atx available at ftp://ftp.igs.org/pub/station/general/; see also Table 2) for both satellite and ground-based antennas were adopted in order to improve the accuracy of vertical site position component estimations15,16. The first-order ionospheric delay was eliminated by using the ionosphere-free linear combination, while a second-order ionospheric corrections17 were applied using the IONEX files from the Center for Orbit Determination in Europe (CODE). The tropospheric delay was modeled as a piecewise linear model and estimated using the Vienna Mapping Function 1 (VMF118) with a 10° cutoff. The Earth Orientation Parameters (EOP) were tightly constrained to priori values obtained from IERS Bulletin B. The ocean tidal loading was corrected using the FES200419 model. The International Earth Rotation Service (IERS) 2003 model for diurnal and semidiurnal solid Earth tides was also adopted.In a successive step, the GAMIT solutions, in the form of loosely-constrained H-files, were aligned to the IGb08 reference frame. Such an alignment was performed by INGV-CT and DICAM-UniBO through the GLOBK/GLORG10 software package by minimizing the deviations between horizontal positions and velocities of achieved solutions and those available in the IGb08.snx20 reference solutions (Table 2). INGV-BO, instead, performed the alignment to the IGb08 reference frame through the ST_FILTER program of the QOCA21,22 software and by combining all the daily loosely constrained solutions with the global solution of the IGS network made available by MIT (ftp://everest.mit.edu/pub/MIT_GLL/), and simultaneously realizing a global reference frame by applying generalized constraints21. Specifically, the reference frame was defined by minimizing the velocities of the IGS core23 stations, while estimating a seven-parameter transformation with respect to the GNSS realization of the ITRF20086 frame, i.e., the IGb083 reference frame.Raw data processing by using the GIPSY softwareGIPSY is a GNSS-Inferred Positioning System and Orbit Analysis Simulation software package (https://gipsy-oasis.jpl.nasa.gov/) developed by the Jet Propulsion Laboratory. Data processing with this software package was performed at DICAM-UniBO (which used GIPSYX, version rc0.4) and INGV-ONT (which used GIPSY-OASIS II, version 6.3).Raw observations were processed by both analysis centers in a precise point positioning mode applied to ionospheric-free carrier phase and pseudorange data24 and using JPL’s final fiducial-free GNSS orbit products. Ocean loading tidal loading and companion tides, computed using the FES200419 tidal model coefficients (http://holt.oso.chalmers.se/loading) and were applied as a station motion model. The wet zenith troposphere and two gradient parameters were estimated every 5 minutes as a random walk process25 by using the VMF118 with a 10° cutoff. A second-order ionospheric correction26 was applied during the data processing. Station clock errors was treated as a white-noise process. The ambiguity resolution was performed by using the wide lane and phase bias (WLPB) method24, which phase-connects individual stations to IGS stations in common view. Resolving ambiguities reduced significantly the scatter mostly in the east component of time series. Satellite orbits and clock parameters were provided by JPL who determine them in a global fiducial free analysis using a subset of the available IGS core stations as tracking sites.In this step, the ONT_GIPSY solution was aligned to IGS144 by applying a daily seven-parameter Helmert transformation (three rotations, three translations, and a scale component obtained from JPL) specifically calculated using the IGS14 Cartesian coordinates and velocities of 132 stations selected by specific quality criteria (e.g. long observation time, continuity through time, position and velocity constrained with sub-millimeter level accuracy, etc.) and located in and around the Eurasian plate and in the Mediterranean area27. This reference frame realization method applies a continental-scale spatial filter to the station coordinates, leading to a reduction of the common-mode errors28 and to an increase of the signal-to-noise ratio.The UBO_GIPSY solution was aligned to the IGb083 reference frame by applying a daily seven-parameter Helmert transformation which was ad hoc calculated using as reference a regional subset of IGS tracking network (see Table 2). The reference coordinates of this regional subset can be found in the IGb08.snx20 file provided by IGS.Each analysis center, in the end, provided position time-series in a common format (e.g. the UNAVCO PBO “pos” ASCII file format7); the ONT_GIPSY solution was aligned to IGS14 while the other solutions were aligned to IGb08. It must be noted that some few stations show time series with different temporal length because each analysis center requested the RINEX dataset to Eni S.p.A. at different times. Moreover, the ONT_GIPSY solution spans the 2000.0–2017.99 time interval since the limited number of continuous stations doesn’t allow to adequately constraint the adopted continental-scale spatial filter prior to 2000.Data RecordsAll the computed raw time series are distributed via the PANGEA repository29 so that the interested researchers can use them for future studies. The raw time series (in ASCII pos format) were stored in different folders and accordingly renamed, based on the solutions described in paragraph “Raw data processing”.We included in the same repository:29A file, named “ENI_Offsets.json”, containing all manually picked offsets; as mentioned above it can be used as input for later use in the TSAnalyzer software24.All logs files as resulting from the analysis performed with the TSAnalyzer software; the logs files were renamed accordingly to the associated solution. These logs contain all the parameters discussed in section 3.The phasor diagrams in png format.The time series plots (observed, modelled and residual components) in png format; the plots were renamed accordingly to the associated solution.The time series comparison plots in png format.The original raw RINEX data are property of Eni S.p.A. company. Interested users can obtain free access to data under a specific confidential agreement by contacting Eni S.p.A. company (see the Supplementary Information for details).Technical ValidationAs mentioned above, the position time series computed by each analysis center as well as the estimated velocity fields were compared in order to highlight possible differences arising from the use of the different software and from the specific models that were selected and applied in the analysis.Position time series comparisonTo perform a simple comparison between the position time series computed by each analysis center we used the TSAnalyzer30,31 open-source software which reads GNSS position time series with different formats and allows to remove outliers as well as to simultaneously estimate linear and seasonal signals.Roughly speaking, position time series can be modelled as the superposition of three main types of signals: (i) linear long-term deformation, (ii) seasonal signals, (iii) instrumental (e.g., GNSS antenna changes) and tectonic (co-seismic displacement) offsets. Therefore, site motion y(t) for each component can be modelled as:1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y({t}_{i})=a+b{t}_{i}+c\,{\sin }(2\pi {t}_{i})+d\,{\cos }(2\pi {t}_{i})+e\,{\sin }(4\pi {t}_{i})+f\,{\cos }(4\pi {t}_{i})+\mathop{\sum }\limits_{j=1}^{{n}_{g}}\,{g}_{i}H({t}_{i}-{T}_{gj})+{v}_{i}$$\end{document}y(ti)=a+bti+csin(2πti)+dcos(2πti)+esin(4πti)+fcos(4πti)+∑j=1nggiH(ti−Tgj)+viwhere ti for i = 1…N are the daily solution epochs (in units of years). The terms a and b are the site position and the linear rate; coefficients c and d describe the annual periodic motion while e and f stand for the semi-annual motion. The summation term corrects for any number (ng) of offsets, with magnitudes g and epoch Tg using the Heaviside32 step function H and the last term, νi, is the residual.A first visual inspection of all raw coordinates time series reveals that they contain both offsets and outliers. Outliers showing large uncertainties (e.g. blunders) were removed by setting in TSAnalyzer an arbitrary sigma criterion (we set 15, 15 and 20 mm, for North, East and Up components, respectively). Offsets were manually picked and recorded in the “ENI_Offsets.json” file (available on the repository) for later use (see also Table 3). Offsets classified as “tectonic” correspond to co-seismic deformation related to the M >6 earthquakes striking central Italy on 24 August 2016 and on 26 October 201633. Offsets classified as “equipment” correspond to change in antennas while offsets classified as “unknown” are characterized by suspicious motion whose nature should be related to service operations on platform.Table 3Offsets manually picked and recorded in the “ENI_Offsets.json” file which is available on the repository.SiteOffset dateNotesAMEB26/10/2011unknown15/03/2012unknown20/03/2012unknown01/04/2012unknown07/04/2012unknown11/04/2012unknown12/05/2012unknown17/05/2012unknown04/12/2012equipmentELEO30/10/2016tectonicEMIL29/07/2015unknown25/01/2016unknownFALC30/10/2016tectonicGARC06/10/2011unknownGRTM22/06/2010unknown24/08/2016tectonic30/10/2016tectonic14/12/2017unknownITGS24/08/2016tectonic30/10/2016tectonicNAPA01/01/2012unknownPCWC07/10/2009unknown14/07/2011unknownPINE24/08/2016tectonic30/10/2016tectonicSMAR31/08/2011unknownFinally, estimated parameters for each solution are reported as ASCII log files in the online repository. The following analysis is based on the results achieved with this time-series analysis.The Weighted Root Mean Squares (WRMS) values of all time-series were also estimated by using the TSAnalyzer software (see log files on the repository), after filtering outliers and estimating parameters according to Eq. (1). In such a step, each single daily solution was weighted with respect to its formal sigma without taking into account the correlation between the North, East and Up components. Frequency histograms of WRMS values for horizontal and vertical components of each different solution analyzed in this study are reported in Fig. 2, while range, mean and median values are reported in Table 4. All solutions are characterized by very similar values, highlighting similar patterns between the solutions as for instance, (i) the mean and median WRMS values for vertical component are about three times larger than the horizontal ones and (ii) the higher WRMS values are observed always at the same stations (CERA and AMEB; see the related logs file in the repository for additional details). Moreover, computed WRMS values are in the same range of the ones estimated for on-shore stations with comparable observing time windows13,34.Fig. 2Frequency histograms of WRMS values for the North, East and Up components of filtered residual time series.Table 4Range, mean and median values (in mm) for time series residuals WRMS.SolutionNorthEastUpRangeMeanMedianRangeMeanMedianRangeMeanMedianIBO_GAMIT0.94–4.791.641.400.93–3.171.491.393.06–8.014.374.20IOE_GAMIT0.88–4.961.551.280.75–3.261.391.282.85–8.674.123.87ONT_BERNESE1.02–5.212.091.731.21–6.432.081.843.18–9.405.004.55ONT_GIPSY0.89–4.431.581.420.70–6.051.481.272.95–6.954.214.05UBO_GAMIT1.08–4.631.661.441.04–5.181.651.492.91–6.084.083.97UBO_GIPSY0.94–4.541.641.490.80–7.491.631.412.96–7.064.334.27All the examined time series contain significant seasonal signals in both horizontal and vertical directions. Values estimated for each solution show a large agreement between them; annual amplitudes range between 0.5–5.3 mm horizontally and 3–8 mm vertically with higher values observed at sites installed on the offshore platforms (see plots reported in the online repository). Semi-annual amplitudes are usually <0.5 mm horizontally and ~1 mm vertically. Additional features on seasonal signals are provided by considering phasor diagrams of amplitudes and phase signals (see the example reported in Fig. 3). Regarding the annual signals, most of the sites installed on the offshore platforms show maximum amplitudes on January, September and August for North, East and Up components, respectively. Moreover, most of these sites (especially those closely located to the coastal area) are subjected to horizontal oscillations with a prevailing NNW-SSE attitude and to vertical oscillation associated with uplift during summer and subsidence during winter. Sites installed along the coastal area are characterized by maximum amplitudes largely scattered and no clear oscillation patterns can be recognized. Concerning the semiannual signals, most of the sites show maximum amplitudes concentrated during July-September for North component and during April-June for both East and Up components.Fig. 3Example of phasor diagrams of annual (top) and semiannual (bottom) signals for station AGOB (see Table 1). Amplitudes of the estimated sine and cosine parameters (see Eq. (1)) are plotted for the North, East and Up components. The upper-right plot represents the key to correlate the maxima phase direction with the day year period. Phases are referred to January 1 and time increases clockwise.Some position time series (mainly recorded at the stations located on the offshore platforms) exhibit a curved shape, while others show weak post-seismic decay (recorded at few stations located along the coastal area, eastward of the 2016 seismic sequence epicentral area). Both features result in deviations from the linear trend; to quantify such a derivation for each time series, we defined a simple non-linearity index (INL) as the ratio between the standard deviation of the smoothed time series and the standard deviation of its intrinsic noise:2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${I}_{NL}=\frac{\sqrt{{\sum }_{i=1}^{N}\left|{F}_{i}-\bar{F}\right|}}{\sqrt{{\sum }_{i=1}^{N}\left|{R}_{i}-\bar{R}\right|}}$$\end{document}INL=∑i=1NFi−F¯∑i=1NRi−R¯where i = 1…N is the number of data, F and R (with mean values \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bar{F}$$\end{document}F¯ and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bar{R}$$\end{document}R¯) are the smoothed and the residual observations. The smoothed time series should approximate the low frequency non-linear behavior, while the residual time series (computed as the difference between the unsmoothed and the smoothed time series) would contain most of the high frequency noise and should be regarded as a good indicator of the station intrinsic noise. Therefore, by taking into account the time series already de-trended with respect to Eq. (1), we computed i) the smoothed time series by applying a moving average filter with a time span of 1 year, ii) the related residual time series and finally, iii) the INL values (Table 5). Estimations of INL index on a set of high-quality IGS stations (Table 2) highlight that values <0.5 reflect time series with pure linear behavior, therefore in the following we arbitrarily set as “gentle” and “moderately” non-linear all the time series with 0.5 ≤ INL ≤ 1 and INL > 1, respectively (Table 5). Based on this simple definition and excluding FRAT station because of its short time series, 11 stations are generally characterized by linear time series (for all the three components) while the other are characterized by gentle (19 stations) and moderately (6 stations) non-linear time series, at one or more components. Stations with gentle and moderately non-linear time series mainly concentrates offshore on northern Adriatic Sea and onshore, along the coastal belt; the stations with linear time series are mainly located along the coastal onshore belt.Table 5Mean INL index values and associated standard deviation (based on the values computed for all the different time series solutions) for North, East and Up components.SiteNorthEastUpAGOB0.61 ± 0.080.41 ± 0.090.74 ± 0.07AMEB0.19 ± 0.080.55 ± 0.220.45 ± 0.18ANBL0.71 ± 0.110.61 ± 0.130.90 ± 0.12ANEB0.97 ± 0.250.34 ± 0.100.36 ± 0.14ANGA1.04 ± 0.391.09 ± 0.401.61 ± 0.69ANLS0.30 ± 0.050.51 ± 0.200.66 ± 0.21ANTA0.54 ± 0.060.76 ± 0.120.79 ± 0.17ANTO0.46 ± 0.060.39 ± 0.070.33 ± 0.08ARIA1.15 ± 0.130.35 ± 0.100.70 ± 0.10ARMI0.40 ± 0.050.53 ± 0.080.40 ± 0.05AZAB0.34 ± 0.070.33 ± 0.100.49 ± 0.06BRIS0.48 ± 0.190.39 ± 0.140.28 ± 0.04CERA2.48 ± 1.090.72 ± 0.470.50 ± 0.26CLAE1.82 ± 0.111.30 ± 0.280.37 ± 0.03CROT0.98 ± 0.230.47 ± 0.150.51 ± 0.05ELEO0.32 ± 0.080.24 ± 0.070.17 ± 0.05EMIL0.30 ± 0.090.34 ± 0.050.28 ± 0.04FALC0.36 ± 0.050.29 ± 0.030.26 ± 0.05FANO0.54 ± 0.030.48 ± 0.100.28 ± 0.10FAUZ0.10 ± 0.030.24 ± 0.050.32 ± 0.03FIUN1.44 ± 0.263.11 ± 1.111.26 ± 0.09FRAT0.02 ± 0.000.01 ± 0.000.01 ± 0.00GARB0.59 ± 0.140.71 ± 0.180.82 ± 0.07GARC0.58 ± 0.330.52 ± 0.420.67 ± 0.33GRTM0.21 ± 0.050.23 ± 0.030.17 ± 0.04ITGS0.77 ± 0.100.82 ± 0.100.33 ± 0.06MANA0.52 ± 0.090.55 ± 0.130.26 ± 0.04NAID0.40 ± 0.090.45 ± 0.060.53 ± 0.05NAPA0.40 ± 0.110.72 ± 0.140.29 ± 0.07PCCP0.42 ± 0.070.46 ± 0.130.40 ± 0.08PCTA1.16 ± 0.191.08 ± 0.300.79 ± 0.12PCWC0.36 ± 0.110.34 ± 0.070.77 ± 0.23PINE0.28 ± 0.040.23 ± 0.050.22 ± 0.04RUBI0.48 ± 0.050.54 ± 0.060.69 ± 0.07SMAR0.28 ± 0.070.45 ± 0.090.60 ± 0.08TEAL1.14 ± 0.091.10 ± 0.190.71 ± 0.13TREV0.26 ± 0.070.29 ± 0.070.23 ± 0.06As a last step, the time-series were compared between them after correcting the offsets (Table 3) and removing the linear trend as defined in Eq. (1). An example of such a comparison is reported in Fig. 4, while the remaining plots are reported in the online repository. This simple comparison highlights how the main features characterizing the time-series (e.g. seasonal signals, noisier time intervals, short-term transients, etc.) show a good agreement between the different solutions.Fig. 4Example of time series comparison from the different solution described in the main text. Each time series is reported after correcting offsets and removing the linear trend as defined in Eq. (1).Velocity field analysis and comparisonIn the following, we performed some simple comparisons between all the computed linear trend values; however, these analyses (and related results) must be considered with caution because of the non-linear motion detected at many stations.As a first step, we calculated the residual values with respect to the mean ones for each solution and for each component (Table 6). Results are reported as frequency histograms in Fig. 5. Regarding the North component, the IBO_GAMIT, the IOE_GAMIT and the ONT_GIPSY solutions show values up to 0.2 mm smaller than the mean values, while the ONT_BERNESE solution shows values up to 0.4 mm larger than the mean values; the UBO_GAMIT and UBO_GIPSY solutions show values centered (±0.1 mm) on the mean values. Regarding the East component, the IBO_GAMIT, the IOE_GAMIT and the ONT_BERNESE solutions show values up to 0.2 mm larger than the mean values, while the ONT_GIPSY solution shows values up to 0.4 mm smaller than the mean values; the UBO_GAMIT and UBO_GIPSY show again values centered (±0.1 mm) on the mean values. Regarding the Up component, the IBO_GAMIT, the IOE_GAMIT, the UBO_GAMIT and the UBO_GIPSY solutions show prevailing values up to 0.6 mm larger than the mean values, while the ONT_GIPSY and the ONT_BERNESE solutions show prevailing values up to 0.4 mm smaller than the mean values; the overall differences range in the −0.6 to 0.6 mm value interval.Table 6Mean linear rate values and associated standard deviation (in mm) for North, East and Up components. No estimates were performed for the FRAT site because of its short time series.SiteNorthEastUpAGOB16.77 ± 0.1922.62 ± 0.12−7.14 ± 0.42AMEB17.38 ± 0.2420.52 ± 0.18−3.07 ± 0.41ANBL18.32 ± 0.1621.96 ± 0.15−3.74 ± 0.38ANEB17.79 ± 0.0620.98 ± 0.18−4.80 ± 0.40ANGA19.41 ± 0.2222.16 ± 0.14−17.56 ± 0.52ANLS18.12 ± 0.1721.18 ± 0.18−4.23 ± 0.34ANTA18.14 ± 0.1920.74 ± 0.16−4.91 ± 0.35ANTO18.45 ± 0.1621.34 ± 0.16−2.26 ± 0.41ARIA17.84 ± 0.2023.18 ± 0.18−10.79 ± 0.56ARMI17.75 ± 0.1321.94 ± 0.16−4.39 ± 0.41AZAB18.65 ± 0.1822.18 ± 0.15−3.64 ± 0.39BRIS18.96 ± 0.1723.09 ± 0.161.55 ± 0.39CERA18.32 ± 0.1922.74 ± 0.13−11.66 ± 0.36CLAE15.73 ± 0.1521.25 ± 0.15−14.31 ± 0.37CROT17.73 ± 0.1827.24 ± 0.10−0.24 ± 0.60ELEO17.98 ± 0.1623.11 ± 0.100.08 ± 0.52EMIL18.23 ± 0.1023.13 ± 0.21−1.89 ± 0.31FALC18.46 ± 0.1323.00 ± 0.17−0.15 ± 0.38FANO18.91 ± 0.1422.40 ± 0.160.15 ± 0.38FAUZ19.77 ± 0.2020.55 ± 0.26−5.65 ± 0.47FIUN15.56 ± 0.1928.72 ± 0.09−14.18 ± 0.43GARB17.77 ± 0.1820.98 ± 0.14−7.06 ± 0.42GARC17.15 ± 0.0819.73 ± 0.08−8.56 ± 0.28GRTM18.50 ± 0.1023.68 ± 0.15−0.30 ± 0.30ITGS19.07 ± 0.0922.70 ± 0.140.13 ± 0.36MANA17.33 ± 0.1721.59 ± 0.14−2.98 ± 0.46NAID17.82 ± 0.1720.51 ± 0.15−3.53 ± 0.42NAPA17.56 ± 0.2021.45 ± 0.08−1.52 ± 0.34PCCP16.84 ± 0.1823.29 ± 0.15−4.14 ± 0.38PCTA19.65 ± 0.1823.86 ± 0.11−7.57 ± 0.40PCWC17.01 ± 0.2122.11 ± 0.04−5.30 ± 0.51PINE18.64 ± 0.1023.18 ± 0.18−0.74 ± 0.35RUBI19.39 ± 0.1322.13 ± 0.16−2.28 ± 0.36SMAR18.39 ± 0.2120.94 ± 0.05−6.54 ± 0.47TEAL18.21 ± 0.1619.76 ± 0.16−4.55 ± 0.39TREV17.23 ± 0.1620.51 ± 0.15−0.10 ± 0.44As a second step, we computed for each solution the velocity field with respect to a Eurasian35 reference frame; achieved results are reported in Fig. 6. Because of its simple definition, such a rotation affects only the horizontal velocities; therefore, the vertical ones remained constrained in their previous reference frame. Considering the horizontal velocity (Fig. 6), all the solutions show a general agreement in terms of azimuthal pattern and vector magnitudes, however the ONT_BERNESE solution show a systematic counterclockwise rotation with respect to the other ones, probably related to the Helmert constraints imposed for the reference frame transformation (Table 2). The ONT_GIPSY solution, although referred to IGS14, is highly coherent with the other solutions. Most of the stations installed along the onshore coastal area show a prevailing NNW-oriented velocity pattern which is highly coherent with the regional deformation field14,27,36. Stations installed along the Adriatic offshore show a complex deformation field characterized by large variations both in the azimuthal pattern and the vector magnitudes (Fig. 6).Fig. 5Frequency histograms of residual values between computed linear trend values with respect to the mean values (Table 5).Fig. 6Horizontal velocity field comparisons. Velocities are referred to a Eurasian35 reference frame. Lower inset: zoom on Calabria region; upper inset: zoom on the Northern Adriatic region.The vertical velocities are reported in Fig. 7 as averaged values (panel a) and as differences with respect to the average values (panel b). Figure 7a shows that most of the stations installed along the onshore coastal area are characterized by subsidence with rates up to 2 mm/yr in agreement with recent studies8, while all the stations installed on the offshore platforms exhibit a general subsidence with rates up to ~17 mm/yr (see also Table 6). The site-by-site comparisons reported in Fig. 7b confirm the features previously recognized: the ONT_GIPSY and ONT_BERNESE solutions are generally smaller than the mean values, while the other solutions are larger than the mean values. These differences seem related to the adopted reference (the ONT_GIPSY solution is referred to IGS14) and/or to the constraints imposed for the reference frame transformation (the ONT_BERNESE solution adopts a Helmert transformation based only on 4-parameters).Fig. 7(a) mean vertical velocity field 5; (b) comparison of the residual vertical velocity values for all the analyzed solutions.Final considerationsAs mentioned above, computed position time series were analyzed by means of the TSAnalyzer30,31 open-source software in order to estimate the linear and seasonal components as well as the offsets values.We detected some offsets corresponding to equipment changes and to co-seismic deformation. Moreover, some detected offsets, classified as “unknown” in Table 3, are characterized by motion probably related to service operations on platform. All detected offsets have been identified in all solutions and are characterized by very similar values (see related log files on the online repository).Examined position time series contain significant seasonal signals: annual amplitudes range between 0.5–5.3 mm horizontally and 3–8 mm vertically while semi-annual amplitudes are usually < 0.5 mm horizontally and ~1 mm vertically. Higher annual amplitude values were observed only at sites installed on the offshore platforms, with maximum amplitudes on January, September and August for North, East and Up components, respectively. The horizontal oscillations occur with a prevailing NNW-SSE attitude while the vertical oscillation are associated with uplift during summer and subsidence during winter. Conversely, sites installed along the coastal area are characterized by maximum amplitudes largely scattered and no clear oscillation patterns can be recognized. All these observations suggest that offshore stations would be affected by variations in hydrostatic pressure and buoyancy on tubular members of the platforms caused by changes in tide levels. Moreover, since the Adriatic Sea is an almost land-locked basin where atmospheric conditions and vary considerably with the seasons37,38, environmental loads caused by wind (e.g. periodic seiche-like effects), current flows and waves would substantially contribute to the observed seasonal oscillations. Due to the main metallic-fabric of platform, seasonal thermal expansion-contraction cycles of the tiny structures could also contribute to the seasonal oscillations.A large number of stations is characterized also by position time series with non-linear behavior at one or more components, as quantified by the computed INL index (Table 5). Stations with INL values larger than 1 are generally characterized also by high subsidence rates (Table 6 and Fig. 7), clearly suggesting that, the tectonic deformation is superseded by the local sources ones (e.g. reservoir depletion, sediment compaction, etc.). Moreover, these local sources affected also the horizontal deformation rates leading to a complex pattern, strictly depending by the relative position of the station (i.e. the platform) with respect to the sources of deformation. Similar considerations can be done for the stations with 0.5 ≤ INL ≤ 1, where however, the effects of local sources have had a minor influence on station motion. Conversely, stations with INL < 0.5 show a velocity pattern which is highly coherent with the regional deformation field of the investigated area.Based on the analysis and the comparisons performed in this study, the different time series solutions are highly consistent between them despite the use of different software, models, strategy processing and frame realizations. The analyzed dataset represents an invaluable dataset since it allows to improve the current GNSS stations density along the Adriatic coastal belt. In addition, this dataset allowed to discover a complex interplay between regional and local sources of deformation, whose relationships would be addressed in future analysis. Indeed, because of their uniqueness, these data could improve and promote further studies in offshore exploiting contexts including: (i) measurement and modelling of induced subsidence patterns and their spatial and temporal correlation with production data, (ii) estimation of natural and anthropic contributions to the overall ground subsidence, (iii) improvement on the long-term regional crustal deformation on offshore areas, (iv) coast-line dynamic and impact on human activities and natural ecosystems (e.g. coast line setback).Supplementary informationSupplementary Information