An Integrated Platform for Ground-Motion Mapping, Local to Regional Scale; Examples from SE Europe

: Ground and infrastructure stability are important for our technologically based civilization. Infrastructure projects take into consideration the risk posed by ground displacement (e.g., seismicity, geological conditions and geomorphology). To address this risk, earth scientists and civil engineers employ a range of measurement technologies, such as optical/laser leveling, GNSS and, lately, SAR interferometry. Currently there is a rich source of measurement information provided in various formats that covers most of the industrialized world. Integration of this information becomes an issue that will only increase in importance in the future. This work describes a practical approach to address and validate integrated stability measurements through the development of a platform that could be easily used by a variety of groups, from geoscientists to civil engineers and also private citizens with no training in this ﬁeld. The platform enables quick cross-validation between different data sources, easy detection of critical areas at all scales (from large-scale individual buildings to small-scale tectonics) and can be linked to end-users from various monitoring ﬁelds and countries for automated notiﬁcations. This work is closing the gap between the specialized monitoring work and the general public, delivering the full value of technology for societal beneﬁts in a free and open manner. The platform is calibrated and validated by an application of SAR interferometry data to speciﬁc situations in the general area of the Romanian Carpathians and their foreland. The results demonstrate an interplay between anthropogenically induced changes and high-amplitude active tectono–sedimentary processes creating rapid regional and local topographic variations.


Introduction
Ground and infrastructure stability are important for our technologically based civilization. Infrastructure projects take into consideration the risk posed by ground displacement (e.g., seismicity, geological conditions and geomorphology). To address this risk, earth scientists and civil engineers employ a range of measurement technologies, such as optical/laser levelling, GNSS and, lately, SAR interferometry.
GNSS is a mature technology that was proven beyond doubt to be able to measure 3D topographic changes with millimeter accuracy [1]. GNSS data sets covering the whole world are processed for free and publicly available from recognized groups, such as Nevada Geodetical Laboratory [2], where data from around 19,000 GNSS stations are processed every day. SAR-based persistent scatterers interferometry (PSInSAR) is a newer technique, certified to measure displacements in the radar line of sight (SAR LOS) with accuracies similar to GNSS [3]. The certification process was led by the European Space Agency as part of the Terrafirma [4,5] program designed to validate the PSInSAR techniques and to select a number of service providers, among them, Terrasigna.
Data sources for PSInSAR are satellite SAR imagery from X and C band spaceborne sensors (e.g., TerraSAR-X, Sentinel-1, etc.) Sentinel-1 data, although not of the highest possible spatial resolution, was proven to be of sufficient quality to measure displacements at millimeter scale [3] in the satellite line of sight (LOS) direction. The Sentinel-1 SAR mission, with worldwide periodic coverage and its open data policy, enables wide-scale monitoring, leading to a thriving archive of ground-motion products, many of them freely available. Similar SAR missions, such as L/S band NASA NISAR [6], to be launched in 2023, will follow the ESA Sentinel-1 data acquisition strategy and will complement the existing C-band Sentinel-1 observations in natural areas, potentially doubling the volume of the current datasets.
Existing services, such as EPOS (the European Plate Observing System), were upgraded to include SAR LOS products [7]. New services are providing processing on demand (ESA Geohazards TEP [8]) or Europe-wide ground-motion products (such as Copernicus GMS [9]).
It is a challenge to integrate motion data from different sensors. In [10], GNSS data was integrated with PSInSAR data by projecting the 3D GNSS data into the SAR LOS; although this exercise was important for validating the equivalency of the two measurement techniques, there will never be enough GNSS ground stations able to project the 3D measurements from GNSS into the contrastingly higher density data of regional PSInSAR network.
Furthermore, since GNSS stations are expensive to install and maintain, PSInSAR can assist with the detection and spatial delimitation of homogeneous critical areas (affected by the same displacement process), which in turn will help with the decision of installing the GNSS stations in the apex (or other key areas) of the observed dynamics.
This work describes a practical approach to address and validate integrated stability measurements through the development of a platform that could be easily used by a variety of groups, from geoscientists to civil engineers and also private citizens with no training in this field. The proposed platform is designed to fill the gap between existing data sources and value-added information extracted from displacement measurements. This is demonstrated through the presentation of local-and regional-scale analysis of deformation phenomena in Romania and its neighbors.
The availability of redundant surface motion products (PSInSAR and GNSS) over the same area facilitates an easier cross-validation of the measurements by cross-comparing these products. The comparison takes into account the differences between datasets, such as direction of measurements and use of different reference areas for each dataset. Existing studies can be used to further validate measurements in areas that were previously characterized with sufficient details.
The current work revolves around monitoring wide areas of around 250,000 square kilometers in Europe (Figure 1a,b), spanning multiple countries and geological structures with complex and intricate dynamics (Figure 1c). The challenge consists of detection, extraction and small-to large-scale characterization of critical areas and their eventual interconnection; these areas are often affected by multiple sources of instability that hopefully can be individually distinguished using the multi-scale information (from individual homes to tectonic units). Integrated and cross-validated data displayed in a format that makes possible quick spatial and temporal inspection (see Figure 2) facilitates detection of new critical areas with confidence due to the same techniques having been proven through cross-validation in previously known areas of interest. Additional layers of information are produced to enhance specific motion characteristics in order to facilitate detection and characterization. The temporal displacement patterns are analyzed for information such as: seasonal displacements (periodical dilation/contraction of infrastructure, earth tides), tectonically induced vertical or horizontal motions over large areas, specific seismicaseismic local activation of faults or accelerated displacement (predictor of a disaster event).  Data interpretation requires input from various end-user groups (e.g., geomorphologists, geologists, geophysicists and civil engineers); thus, the platform is designed to be easily accessed and self-sufficient in terms of data inspection and extraction.
Regional tectonic displacements are generally difficult to detect and monitor; stateof-the-art monitoring technology consists of GNSS measurements from a sparse network of stations and heavy spatial interpolation, which may result in significantly contrasting interpretations when movements are in the order of a few mm/year [12]. With this technology it is difficult to determine whether the measured motion is tectonic or local (e.g, due to local ground instability or poor antenna stability) unless measurements show a direct regional consistency at the scale of specific tectonic units. One way to overcome measurement errors is high spatial redundancy; with PSInSAR, millions of measurement points could be available for a single tectonic unit, and any quality issues with a single or a few points would be easily visible. The high density of measurements makes possible the actual delimitation of a tectonic unit by the differential motion at its margins. The effectiveness of this type of analysis is already proven in large, tectonically active induced motions, such as the San Andreas Fault of the western USA or the North Anatolian Fault of Turkey [13]. Furthermore, the high-density data allow the defining of anthropogenically induced topographic changes driven by the exploitation of the subsurface, such as in gas or geothermal reservoirs.
In our study we calibrate and validate the SAR interferometry platform and its correlation with GNSS data with several specific situations in the general region of Romania and neighboring countries, areas in central and eastern Europe centered on the Carpathian Mountains. This region is characterized by an active tectonic regime localized in the SE Carpathians inducing regional vertical and horizontal movements up to 7-8 mm/year [12,14] and associated with large, intermediate-mantle to crustal seismic events and crustal faulting, as detected and characterized by deep geophysical studies [15][16][17]. Other, less well-understood, shallow crustal movements may be driven by salt diapirism in the Transylvania basin or rapid depositional changes driven by sedimentation in the Danube Delta and associated coastal evolution [18][19][20][21]. Local anthropogenic effects induced by gas and water exploitation, rapid re-sedimentation due to coastal protection work or subsidence-induced desiccation of wetlands in the Danube Delta may also have significant effects [22]. Figure 2. PSTool, an online and free-access analysis tool for ground displacement monitoring [23]. Figure 2. PSTool, an online and free-access analysis tool for ground displacement monitoring [23].

Online Data Integration Platform
The Synthetic Aperture Radar Persistent Scatterers Online Software Tool (PSTool, [23]) was developed by Terrasigna with the purpose of exploiting the large volume of grounddisplacement data made available with the advent of the ESA Sentinel-1 mission and of the expected similar volume of data to be contributed by the future NASA NISAR [6] and ESA ROSE-L [24] missions.
PSTool is intended to be used by end-users with no knowledge of SAR technology, thus it continuously evolves towards the applications while keeping the platform simple and intuitive.
The PSTool platform can be used to inspect temporal characteristics of the ground motion, select areas of interest and extract temporal averages of displacement rates, export temporal profiles to standard formats for integration in the user's own platforms and upload user-specific layers on top of the displacement information.
One example of the kind of analysis that can be performed with PSTool is illustrated in Figure 2; the user can select an area using an irregular polygon that can be drawn around areas of any shape (in the figure, the cursor in the shape of a blue circle points to the area where the polygon was drawn, and a blow-up image of the polygon is illustrated on the top-right by the arrow), and in the selected area the average LOS displacement time-series will be calculated for all the included PS points. If data from multiple InSAR processing results are available, the averaged time-series and fitted linear trend will be displayed per dataset. In the temporal profiles displayed in Figure 2, for the two measurement layers covering the area of interest, the number of temporal profiles that were averaged is displayed beside every layer name (6 profiles for layer EPOS_36_5 and 31 profiles for layer Track_36c) together with each of the annual displacement rates.
The left column in Figure 2 illustrates the stack of PSInSAR layers available for analysis; each of these layers was obtained from different processing algorithms and using various data sources (X, C and L-band).
For a better interpretation of the results, the user can load custom shapefiles that will be displayed below the PSInSAR results for target identification (in Figure 2 top-right, "Falii_Project_WGS84.shp" represents such a layer).
Moreover, PSInSAR displacement rates or displacement rates plus temporal profiles can be exported in Google Earth or text formats for areas defined by the user.
The architecture of the PSTool application is depicted in Figure 3. The very large number of points in the deformation maps makes it impossible to show them in vector format at each possible zoom level through a web application. Hence, the application uses a map server (GeoServer [25]) and a tile-caching application (Map Proxy [26]) that serve jpeg tiles to the client through a webserver (Apache [27]). The PSTool interface uses OpenLayers [28] for displaying geographical data. The persistent scatterers (geographical coordinates) along with complementary data (deformation profile data, deformation velocity, height, etc.) are stored in a PostgresSQL [29] database with the GIS extension of PostGIS [30]. Data is uploaded in the PostgreSQL database using dedicated scripts for different types of data. The web server uses a MySQL [31] database for storing user rights to accessing maps, user actions, temporary passwords and codes, graphical user interface configurations and other data.
The deformation maps from the PostgreSQL database are periodically synchronized with the MySQL database to ensure proper user-rights management. An administrator web interface is available for managing user rights and point-color maps.
Though a tile caching server is used, the map server is always on, as point data for some high zoom levels will be extracted directly from GeoServer.
All servers (Apache, PostgreSQL, etc.) can run on independent machines to maximize performance, but currently these components are implemented in Docker containers running on the same machine.
The platform also provides programs for deformation-profile post-processing, such as seasonal motion detection/removal and sudden-change detection. At the moment, these are available on demand, only for experts, for off-line processing.
Seasonal motion is detected using harmonic analysis-the deformation profile is resampled to a uniform step, followed by frequency analysis and stop-band filtering the seasonal spectral component (if required). For sudden-change detection, the deformation profile is analyzed using a sliding window of eight measurements. At each step, the local linear deformation trend is computed, the next measurement is predicted and the difference between the predicted and the measured value is computed. If it exceeds an empirically chosen threshold, the profile is classified as exhibiting sudden changes [32].

GNSS Data
As part of a collaboration with the National Institute for Research and Development on Marine Geology and Geo-ecology (GeoEcoMar) and the National Institute for Research and Development of Earth Physics (NIEP), Romania, GNSS data was incorporated in the PSTool platform.
The raw GNSS data was obtained from the GeoPontica permanent GNSS station network-part of the EUXINUS research infrastructure managed by GeoEcoMar and the permanent GNSS station network of NIEP [33], covering the Romanian and Bulgarian Black Sea coast.
The processed GNSS data used in this study was imported from Nevada Geodetical Laboratory (NGL) [2]. Their analyses use data from over 19,000 GNSS locations around the world to estimate the rates, patterns, budgets and sources of the vertical movement of the land surface that is driven by various processes in the Earth, such as plate tectonics, glacial isostatic adjustment, the earthquake cycle, changes in aquifers and hydrological loads. The GNSS data are expressed in ENU local coordinates (east, north, up) to be able to integrate with the PSInSAR results, which are expressed in local SAR LOS (azimuth and incidence angle [34]). Figure 4 illustrates the locations of the NGL-processed GNSS stations available in Romania; any of them can be used as a relative reference to the other GNSS stations. The GNSS stations measure motion in 3D coordinates, and each coordinate can be selected for visualization in PSTool. The data can also be viewed in on-the-fly differential mode (Figure 5a), with the possibility of any of the stations to be used as a reference.
Romania is located between 20-30 • E and 43-48 • N. The maximum deviation between the vertical (up) axes across an area equivalent to Romania is around 8 degrees ( Figure 5b); as an example, if the GNSS reference is chosen in SW Romania (bottom-left corner in Figure 5b) with vertical U SW , and the measuring station is chosen in the NE corner of Romania with vertical U SE , the direction match between the two local N directions is 0.99 (U NW = U SE × cos (8 • ), projection of one direction to the another); thus the differentialdisplacement-rate error is only 1%. Thus, subtracting reference motion from GNSS station motion along each ENU axis is a valid operation as long as the orientation of the ENU coordinates coincide, which is a true assumption for an area the size and location of Romania [35,36].

SAR Data
The total monitored area in this work spans around 250,000 square km, and it was mosaicked from four Sentinel-1 tracks, merging up to four frames for each track. The coverage includes Romania and part of its neighboring countries (Figure 1b), spanning multiple geotectonic units ( Figure 1c) [11,37]. The Copernicus DEM [38] at 30 m resolution was used in order to create a first approximation of the height model.
Around 4000 SAR scenes acquired between 2014-2020 were processed. The InSAR algorithm, described in [39], preserved all the displacement information in order to maximize the chances of detecting slow crustal movements (mm/year). This was performed by extracting the linear deformation information before any kind of spatial or temporal filtering was applied, with the purpose of improving phase statistics. The drawback of the algorithm consists in the challenge of extracting the displacement information from an atmospheric phase that is two orders of magnitude above the displacement information; this requires a careful phase unwrapping check of every residual interferogram, individual correction and displacement model updates.
Ground-motion data produced with the SBAS algorithm [40], provided via the EPOS program (under open license [41]), was used to both complement and cross-validate Terrasigna's PSInSAR product. This algorithm uses multiple-reference, small-baseline interferograms (with negligible topographic phase) in order to maximize temporal coherence and minimize phase. Assuming ground motion does not occur on the short time intervals used to produce interferograms, the main phase component that remains is represented by the atmospheric phase. A spatial filter can thus be applied to reduce the variability of the atmospheric phase, which acts like noise. The ground motion information is then extracted after atmospheric filtering, but the major risk of this approach is the inadequate separation of atmospheric effects from small-scale displacements, leading to the possible filtering of the ground motion on wide areas (e.g., regional motion of tectonic units).
The available (sparse) GNSS stations can be used to validate the 3D motion of homogeneous areas (characterized by similar surface variation dynamics) that are part of a SAR LOS direction of a PSInSAR (dense) network. In [34], the projection formula from GNSS (ENU coordinates) to SAR LOS was derived. With this formula, wherever a GNSS station is available the amplitude and direction of motion vectors from neighboring PSInSAR LOS measurements can be validated (example in Table 1). Conversely, the PSInSAR data could be used to characterize the stability of the area where a GNSS station is installed; GNSS represents a single measurement point, thus it is complicated to figure out if the dynamics measured in that point are local or regional. PSInSAR, with a dense measurement network, can map the dynamics around the GNSS station at various scales, providing insight on what motion the GNSS station actually records.

Methods
In order to extract value-added information, a more automatic post-processing is required. The implemented post-processing approaches have two main scopes: classification based on the dominant motion trend and separation of different causes of ground motion.

Classification Based on the Dominant Motion Trend
Every persistent target is analyzed in time and classified based on the dominant motion (e.g., linear subsidence, sudden changes/fractures in displacement rates or seasonal displacement). The output of this step is represented by separate maps of motion, from which specific motion patterns could become spatially distinct and correlated with areas or structures on the ground. In this way, seasonal displacement could be separated from linear displacements. Moreover, settlement or increased subsidence/uplift events can be detected.
As an example of detecting an area with a variable subsidence rate, Figure 6a illustrates the displacement map of Izmail (Ukraine) and Figure 6c a sample of a relatively linear temporal displacement profile from that area. Figure 6b illustrates an area on the ground affected by a changing subsidence rate, whereas Figure 6d illustrates a corresponding sample of a temporal profile affected by a variable displacement rate (increased subsidence between September 2016-July 2017).
An example of detecting areas affected by seasonal motion is illustrated in Figure 7. The seasonal motion is detected by harmonic analysis of each temporal profile and identifying repetitive patterns with a period of approximately one year. Figure 7a illustrates the whole set of the PSInSAR targets available for monitoring the Sulina, Romania area. By selecting the subset of the PSInSAR targets affected by seasonal motion, the map in Figure 7b was obtained. An example of a temporal profile extracted from this set is illustrated in Figure 7c.
If seasonal motion is not of interest, it can be removed from the temporal profile; in this way, other motion components, such as subsidence patterns, will be more visible. Figure 7d illustrates a temporal profile from which the seasonal motion was removed.

Separation of Different Types/Causes of Ground Motion
On a small-scale map, local and wide-area (including tectonic) motion will always be mixed, making them difficult to distinguish from each other. Thus, criteria of separation of different motion types need to be defined and applied before any interpretation of the results is attempted. In order to distinguish between motion types of different scales, it is important to understand how the analyzed datasets were produced and what they represent. In this way, motion related to surficial processes (e.g., shallow subsidence or ground humidity variation) could be separated from deep, dynamic processes (tectonic). Two criteria are discussed below.
Criterion A: SAR operating frequency used to acquire the data  The depth of the ground penetration of microwaves is proportional with the wavelength [42]. L-band data (JAXA ALOS and future NASA NISAR) tend to remain interferometrically coherent on natural areas; C-band data (ESA ERS/ASAR/Sentinel-1 and CSA Radarsat) remains coherent on infrastructure and some bare-earth areas; X-band data (DLR TerraSAR-X and ASI Cosmo Skymed) tend to remain coherent mostly on infrastructure.
As a result of these differences, it is important to know which targets represent infrastructure and which targets represent natural areas, since the mechanism of their measured motion could vary. The motion of infrastructure anchored on stable bedrock could be related to local geology or to tectonics. The measured displacement of a natural ground surface could be related to deep motion or to surficial motion, such as solifluction (within a few centimeters depth), or could be related even to ground humidity variation, since the characteristic depth from which electromagnetic waves are backscattered depends also on the humidity level of the soil [42].
Criterion B: Processing algorithms used Algorithms using amplitude dispersion and/or spectral phase diversity [39] to detect the candidate measurement points target infrastructure, whereas algorithms employing coherence and spatial multi-looking [40] target both infrastructure and natural areas as long as they remain coherent in time. Figure 8 illustrates the difference in dynamics between infrastructure and natural areas. Figure 8a illustrates the results of PSInSAR analysis where only targets that are persistent for a long period (multiple years) are detected. Figure 8b illustrates the results of lower resolution processing with resolution cells including both infrastructure and natural areas (such as illustrated within the blue circle). Figure 8c,d illustrates the difference in displacement patterns between infrastructure and natural targets in the same area.
Using these criteria, the motion affecting natural areas should be separated from motion affecting infrastructure, since the source of motions could be very different and mixing them could be confusing when interpretation is attempted. Thus, it is important to know what processing algorithms were employed in order to understand what is measured on the ground.  Figure 7b).

Results
The results can be separated into three categories:

Local Scale Results
Local scale results are represented by clearly localized subsidence/uplift areas that can be delimited directly from the PSInSAR maps. In this case, atmospheric effects are negligible and the accuracy of measurements relative to neighboring areas is high. A number of examples are illustrated below, some of them validated and some others representing new information requiring a more detailed analysis by future studies. Due to the large extent of the monitored areas, automatic detection and mapping of critical areas was employed; results were synthesized in a national atlas of ground stability [43].

Salt Exploitations
Provadia area (NE Bulgaria) Located ca. 50 km west of Varna city and the Black Sea coast, the Mirovo salt diapir in the Provadia area has been exploited since the Late Neolithic, Middle and Late Chalcolithic eras, being the oldest salt production center in Europe (5500-4200 BC) and favoring the development of the first prehistoric urban center in Europe (4700-4200 BC) [44]. Industrial exploitation started in the 1950s and is performed using a borehole system injecting fresh water and extracting salt solution [45]. Since the late 1970s, salt exploitation in the Provadia area has been associated with increased seismic activity and large surface subsidence [45].
PSInSAR measurement results for the area indicated by the polygon as illustrated in Figure 9 reveal an average subsidence of −22.6 mm/year for the 2015-2020 timespan in the SAR LOS direction.

Solotvyno (Ukraine)
Similar local land-surface subsidence with a maximum line-of-sight deformation rate of 5 cm/yr associated with sinkhole development in a salt exploitation area was also detected by InSAR methodology at Solotvyno in Ukraine, close to the Romanian border, as reported by Szűcs et al., 2021 [48].
The subsidence effect is also visible from PSTool, serving as cross-validation with the previous study. Figure 10 illustrates with red the subsidence area; a temporal displacement profile selected from the maximum subsidence area indicates a subsidence rate of around 15mm/year. The zero-displacement point in the profile represents the time of the reference scene. Some PSInSAR methods output temporal profiles that intersect the displacement axis (zero-displacement point) at the time of the reference SAR acquisition. This kind of representation is useful because it provides information about the reference date. A different approach (that is also used for GNSS data representation) is to shift the profile in order for the displacement to start from zero; this approach makes sense especially for PSInSAR algorithms that use multiple temporal references. Since multiple datasets can be used in parallel and each one likely starting at a different time, calibrating each dataset to start from zero at the start time does not improve the analysis.

Critical Infrastructure Stability in the City of Cernavoda (Romania)
The Cernavoda nuclear plant is built in the vicinity of the Danube River. The reactors are built on the more stable sandstone bedrock, whereas some other infrastructure is built on the alluvial formation. Subsidence of −3-−8 mm/year (Figure 11a,b) is to be noted in the alluvial Areas 1 and especially 2 (Figure 11c), as elucidated from Sentinel-1 PSInSAR measurements. The subsidence was cross-validated with independent PSInSAR measurements using data from the ESA ASAR mission, as illustrated in Figure 11d and correlated with the geological layers as illustrated in Figure 11c.
A closer look at Area 1 confined within the nuclear plant industrial complex (Figure 12a) reveals an area affected by subsidence of around −3.0~−4.0 mm/year (Figure 12b), and an isolated structure (red dot indicated by black arrow) subsiding with an average of −7.4 mm/year which suffered a 12 mm drop in September 2018 (Figure 12c).

Oil Exploitation in Videle Area (Romania)
In the central part of the Romanian Plain (Southern Romania), PSInSAR measurement results emphasized −5-−7 mm/year subsidence in some areas (e.g., Blejesti and Ciuperceni) situated northwest of the town of Videle, Teleorman County ( Figure 13). In this region an oil field was discovered in 1958, and oil extraction activities began in 1961. The reserves of the oil field are around 300 million barrels, and production averages around 22,000 barrels (3500 m 3 /d) per day [50] Both GNSS and PSInSAR measurements show a subsidence process that affects this structure and the surrounding areas, though at different rates. Figures 14 and 15 illustrate a cross-validation exercise between the SULN GNSS station and PSInSAR measurements. Figure 14a-c illustrates the ENU displacement of SULN relative to the AEGY GNSS reference. Figure 14d illustrates the vertical displacement (U direction) of the SULN GNSS station relative to the IGEO (Chisinau). It can be seen that the SULN GNSS station subsides with 9 mm/year relative to IGEO, but only 5.2 mm/year relative to AEGY; this means the AEGY station itself subsides with 3.8 mm/year relatively to IGEO and proves how important it is to find a common reference when working with wide-area, multiple PSInSAR datasets; with each one having different reference areas, a reference calibration is required before integrating measurements from different sources.
The PSInSAR measurements, with a reference set in Tulcea near the AEGY GNSS station, show a −3 mm/year subsidence of the topographic surface of the SULN GNSS station area in the LOS direction (Figure 15a), as well as a 9 mm subsidence in the immediate surrounding area (Figure 15b).  Projecting the 3D GNSS vector to the Sentinel-1 LOS direction results in a value of −3.91 mm/year subsidence, very close to the value of −3 mm/year measured with PSInSAR (Table 1).
Only 8 km away from the "New Lighthouse", the same −3 mm/year subsidence rate is visible in SAR LOS (Figure 15c) south of Sulina city, whereas in agricultural land a larger subsidence of around 6 mm/year is measured (Figure 15d).
The higher dynamics observed in the area surrounding the lighthouse (but not the lighthouse itself) are related to the hydrotechnical infrastructure or jetty [19,51] (the jetty consists of large boulders grouped together that generate local eddy-like currents with a major role in very fast sediment accumulation). The jetty structure is thus still settling itself on the ground, contributing to the fast compacting of freshly accumulated sediment as a consequence of its weight. The local compaction due to the northern and southern jetties is also visible from the PSInSAR measurements. The rhythmic variations could be due to the impact of storms and/or sea level variations (quite significant at the mouth of the Danube). The two isolated points that are visible in the central-right side (N of the jetties) of Figure 15a,b are related to the shipwreck that was stranded here in the beginning of 2010 slowly drifting SW due to the strong NE storm waves and winds.
In regards to the vicinity of Sulina town, significant subsidence was measured in an anthropogenic area where the city structures are relatively stable, while the agricultural land in the southern part is strongly subsident. The agricultural polder was protected by an embankment between 1970-1990 and is now affected by an accelerated compaction that followed the drying up of the former wetland after the land was reclaimed for agriculture.
Such processes of widespread subsidence in drying peatlands have been quantified in deltaic areas by InSAR data worldwide [52,53].

Wide-Area Displacement Maps
Wide-area displacement maps are usually dominated by a collection of clearly visible (high displacement rate) local motion areas originating at rather shallow depths and possibly masking deep (tectonic), low-amplitude motion.
The Transylvanian Basin is well known for its numerous salt diapirs, whose evolution has been described to be active during post-Miocene geological times, either as exaggerated domes in the north-eastern part of the basin, or fault-controlled salt decollement in its SW sector [20,21,54]. As an example, Figure 16a illustrates the central part of the Transylvanian Basin with numerous subsidence areas identified. Figure 16b illustrates the depths of the salt diapirs (from [21]) overlaid with PSInSAR data. The correlation between diapirs and the local zones of subsidence (Figure 16b) demonstrate that this subsidence takes place in areas where salt is located at 3-4 km depth, and diapirism is in incipient stages (pillows) along their flanks or in the immediate withdrawal synclines. The wavelength of these structures is one order of magnitude higher than that of the observed subsidence, which demonstrates that salt diapirism cannot drive the observed surface subsidence.
Another visual correlation, this time between the subsidence and gas extraction areas (dark green spots) is illustrated in Figure 16c. The correlation is quite evident within the areas outlined with black ellipses.
In Figure 16c, one strong subsidence in the SE part of the city of Targu Mures (blue circle) does not correlate with the local gas extraction sites, but with an underground gas storage facility. Figure 17a illustrates the location of the storage area and of the 18 wells used to inject/extract the gas from the natural storage [55]. Figure 17b illustrates the injection cycles [56] and Figure 17c the PSInSAR-measured temporal behavior of the ground above; a seasonal motion is visible, perhaps related to the injection cycle, although delayed with 2-3 months. This timing is likely driven by the nonlinearly elastic, viscoelastic or non-elastic response of the reservoir overburden, although such a response is highly variable based on the rock rheology and can take place from days to years [57].

Wide-Area Deep Motion
When monitoring large areas, horizontal ground motion becomes important, since such areas could include multiple tectonic units moving in different directions. In [34], the geometry of the SAR measurements is explained; horizontal ground motion towards the nadir of the radar track could be confused with uplift, and motion in the opposite direction (away from the satellite) could be confused with subsidence. In the case of the descending Sentinel-1 track moving approximately NE to SW, subsidence and NW motion could be confused. For an ascending track, subsidence and SE motion could be confused. Thus, GNSS measurements are important to clarify the actual 3D movement vector.
The typical dynamics of tectonic units is within the range of mm/year, which is at both GNSS and PSInSAR error levels [1][2][3][4] (there are a few exceptions, such as in the Anatolia peninsula [13], where the tectonic dynamics are dominant over local dynamics and are clearly visible from direct measurements).

Dobrogea Region (SE Romania)
In Dobrogea and on the Black Sea coast (Romania and Bulgaria), GeoEcoMar and NIEP are operating networks of GNSS stations aiming to detect and monitor vertical and horizontal crustal movements. Data processing results are made available by NIEP [33] and NGL [2] and incorporated into the PSTool platform (Section 2.2) available online and used in this work. The GNSS data indicate very small amplitude horizontal displacements in Dobrogea with respect to the stable Eurasia domain. The resulting vertical movement vectors indicate crustal uplift in the North Dobrogean Orogen, in its prolongation along the North Dobrogean Promontory and in Central Dobrogea.
In order to detect such small-scale, slow motion, which is typically masked by a large number of local areas with higher dynamics and often moving in opposite directions, additional processing is required.
Step 1: determine the dominant motion. In order to determine the dominant ground motion, areas of "subsidence" (or rather, of motion away from the satellite), stability and "uplift" (or rather, of motion towards the satellite) are separated by their motion direction, and their relative dominance is analyzed. Figure 18 (a, total ground motion; b, stable targets; c, subsidence; d, uplift) illustrates an example of how different dynamics are mixing within the same area, making it impossible to extract any smaller-scale motion.  Figure 18a illustrates the total motion (subsidence and uplift) in an area; with this representation it is difficult to separate slow motion at a small scale since multiple pixels with different dynamics are overlapping each other in the figure, with higher amplitude local dynamics dominating the rest. From Figure 18c it can be seen that the spatial density of the subsiding targets is much lower than the spatial density of the 'uplifting' targets ( Figure 18d). Moreover, subsiding targets seem to be grouped in local areas, some of them natural (where it is not clear if ground motion or rather ground moisture variation are measured). From the map in Figure 18d it can be seen that the 'uplift' motion is much more spatially dominant (albeit at lower magnitudes than the spotty subsidence areas) and uniformly distributed, thus it seems to have a more regional characteristic, rather than local. Small-scale motion (regional, tectonic units) is expected to be spatially uniform, with either constant or slow-changing spatial gradients, in agreement with the observed 'uplifting' areas.
As a conclusion, local motion and wide-area motion were separated based on the direction and spatial extent of the motion.
Step two: detect structured dynamics of the ground. By extending this representation to wider areas, certain structured motion becomes apparent. Figure 19a illustrates the subsidence/NW component (negative displacement in radar line-of-sight (LOS) direction) of the measured motion. The spatial distribution pattern of the magnitude of the motion is relatively uniform across the whole region, spanning multiple tectonic units, such as the East European Platform, Dobrogea and the Moesian Platform [37]. No specific structure at tectonic scale can be detected based on its individual motion along this direction. Figure 19b illustrates the uplift/SE component (positive displacement in LOS direction) of the motion; the spatial distribution pattern of this motion is nonhomogeneous. A structured motion is visible in multiple areas delimited by rectangles, which may be interpreted as being related to the south-eastward crustal displacement of the South-East Carpathians Bend Zone and its foreland as emphasized by results of repeated GNSS measurement campaigns, structural, tectonic and exhumation studies [12,17,[58][59][60][61][62]. Therefore, the uplift/SE motion of the SE Carpathians and their foreland, including Dobrogea, is the likely solution.
In Figure 20, the uplifting area delimited by the white rectangle in Figure 19b was visually correlated with the crustal motion map [63]. Using equation (5) in [34] and the input data shown in Table 2 for the time interval 2017-2021, the 3D GNSS trend projected on Sentinel-1 SAR LOS results in a 1.07 mm/year uplift.  Figure 22 illustrates PSInSAR-derived displacement rates around the CRCL GNSS location. Both trends are upward (1.28 mm/year and 1.12 mm/year, respectively), consistent with the GNSS results projected in the PSInSAR LOS direction. Differences between the displacement rates could be related to the measurements having different reference stations (Chisinau Moldova for GNSS, Zoria Ucraine for TRACK_36c PSInSAR and Tulcea Romania for EPOS_36_5 SBAS datasets).

Discussion
The availability of SAR data covering the entire Earth surface and the Sentinel-1 minimum data units covering an area of 250 × 200 square km are leading to a change in the ground-motion measurement paradigm. It is more manageable to start from processing wide areas (Section 2.3) containing a large number of critical spots (known and unknown) rather than starting with local processing for a single critical area. The advantage of the new approach consists in providing a wider image of the dynamics and eventual connections between multiple local critical areas. Moreover, with this approach a number of critical sites that have been previously studied are serving as validation of the measurements (Section 4.1.1). SAR daa is acquired at different operating frequencies. The penetration depth of the microwaves is a function of the microwave's frequency (Section 3.2). C-band microwaves (such as those transmitted by Sentinel-1) tend to be scattered mostly at the surface, whereas L-band microwaves (ALOS PalSAR, the future NASA NISAR and the future ESA ROSE-L) tend to penetrate deeper and scatter at the interaction with the layer underneath or at the interface with high-moisture soil. Therefore, it is very important to know the source of the SAR data and the properties of the ground (infrastructure, bare ground or agricultural) in order to interpret the measurements (e.g., to differentiate between ground motion and variations in soil moisture).
Multiple PSInSAR processing algorithms are currently in use (Section 2.3). Each of them serves a specific purpose. Some algorithms are optimized to extract the exact (nonlinear) motion pattern of local areas; these algorithms make assumptions about the spatial and temporal scale of the motion and filter the data accordingly in order to lower the measurement noise and improve the solution. With these algorithms, small-scale motion (such as tectonics) could be filtered out as atmospheric disturbance.
Other algorithms are optimized to extract the linear component of the displacement rate, which is a parameter of interest for long-term monitoring (such as for the dynamics of the tectonic plates). These algorithms do not produce the exact temporal local motion pattern, but they also do not filter the data before extracting the displacement rate; in this way, all the linear displacement information at any scale will be preserved. When using these algorithms, the unfiltered atmospheric disturbances could influence the linear displacement estimates; to avoid this, processing is performed in small patches for which the atmospheric effects are constant and thus do not influence the displacement information.
It is therefore important to be aware of the type of PSInSAR processing applied to the SAR data before looking at the displacement information.
Ground stability is affected by multiple factors that lead to ground displacement; separation and classification of different motion types (e.g., settlement/accelerated motion, fracture or seasonal variation) helps with identifying the underlying factors contributing to the ground displacement (Section 3.1). Based on temporal characteristics, new layers of data representing classes of motion types can be produced and used to spatially identify areas on the ground affected by specific motion.
At regional scale, multiple sources of motion interfere with each other, making it impossible to separate small-scale from local motion in one map (Section 4.3). Post-processing is required in order to detect wide-area structured motion.
Finally, wide-area maps provide valuable information by detecting critical areas (some of them previously unknown) that require further attention. For those areas, additional processing may be required; this includes local reprocessing of the Sentinel-1 data to densify the monitored network or the use of high-resolution data (TerraSAR-X, CosmoSkyMed, Radarsat-2) to obtain a more detailed map of the motion.
Future developments: Future missions: Sentinel 1 A/B represents only the start of a new era [64], with open access C-band SAR data that will be acquired systematically until at least 2037 [64]; this long-term vision justifies sustained investment in developing the monitoring technology even further. Additionally, complementary L-and S-band missions from ESA and NASA will provide similar data coverage and volume. NISAR [6], to be launched in 2023 and with a minimum lifespan of 3 years, together with ROSE-L [24], to be launched in 2028 and with a lifespan of at least 7.5 years, will provide complete coverage of Earth at L-band for at least the period 2023-2036.
The L-band data will contribute motion information of infrastructure, leading to a densified PSInSAR network. Additionally, L-band monitoring will extend to natural areas [65], providing motion measurements on bare rocks together with ground moisture variation measurements in soil (for agriculture, wetlands, etc.) Future data: A number of programs will provide free ground-motion measurements using PSInSAR and variant techniques. The European Copernicus Ground Motion Service (EGMS) will provide open-access, Europe-wide, ground-motion maps for 2014-2021 produced from Sentinel-1 data [66], possibly with future updates.
The EPOS program [7] also started to make available large area PSInSAR measurements for plate-tectonic motion monitoring under a free and open license [41].
All the current and future open-access data are envisioned to be accessible from the open-access PSTool platform [23] and to use redundancy for quick cross-validation and for extraction of the maximum information in critical areas. Access to platforms hosting national or even continental data (such as EGMS [65] or NASA ASF [66]) will be realized through application programming interface (API) libraries provided by the host platforms developers. Additionally, in-house processed data will fill in missing pieces and will fulfil the requirements specified by end-users.

Future research:
Integration with the GNSS network and the expansion of the scale of the monitored areas will continue. PSInSAR is capable of detecting wide-area motion that could pinpoint interesting locations for new permanent GNSS stations, together with information about locally stable areas suitable for GNSS station foundations. In turn, GNSS can provide the 3D motion vector for the wide-area structures detected with PSInSAR. In gap areas (no infrastructure), radar corner reflectors could be installed in tandem with permanent GNSS stations. The future L-band motion maps are expected to close some of these gaps, thus emphasizing the importance of integrating all sources of information available. The next step in research leading to an advance in the monitoring technology will be to include Artificial Intelligence (AI) in the monitoring process [67]. Correlating previous critical/disaster events measured with PSInSAR with information about external causes (e.g., extreme climate, earthquakes or anthropogenic activities) should lead to a risk analysis and/or the prediction of future disaster events when the conditions are met.

Conclusions
In the last 20 years, the PSInSAR techniques for ground displacement measurements evolved from scientific research to the operational stage. The Sentinel-1 mission filled in the last piece towards the operational service status by providing open access to data fully covering the Earth every 6 and 12 days.
The challenge shifted from local area PSInSAR processing algorithms to wide-area processing and information extraction from the resulting big data.
Increasing the extent of the monitoring areas also introduced challenges related to the dynamics at the scale of tectonic units, as multiple tectonic units can be contained in a single monitored area.
In the current work, a methodology was designed to address these challenges. First of all, a data integration platform (PSTool) was developed (Section 2.1). The purpose of this platform is to bring together ground motion measurements from different sources obtained with various technologies, integrate them in a practical way and make them available to potential end-users to increase their productivity and consequently increase the acceptance of PSInSAR as a ground and infrastructure stability monitoring technique.
GNSS data is harmonized with PSInSAR (Section 2.2) by producing approximate differential GNSS measurements with references located closer to the PSInSAR references. PSInSAR references can also be moved to more convenient areas (more stable or closer to available GNSS stations). Finally, 3D GNSS measurements can be back-projected to the SAR LOS direction to be cross-validated with PSInSAR measurements. Acknowledgments: GeoEcoMar and NIEP GNSS data was made available to the authors through the SETTING National Project (Integrated thematic services in the field of Earth observation -a national platform for innovation), cofinanced from the Regional Development European Fund (FEDR) through the Operational Competitivity Programme 2014-2020 (contract no. 336/390012). We thank NIEP (especially Alexandra Muntean and Eduard Nastase) for their contributions to GNSS data validation and suggestions regarding free data available from NGL and suitable GNSS processing software.

Conflicts of Interest:
The authors declare no conflict of interest.