Assessing Current Coastal Subsidence at Continental Scale: Insights From Europe Using the European Ground Motion Service

Abstract

Beside climate‐change‐induced sea‐level rise (SLR), land subsidence can strongly amplify coastal risk in flood‐prone areas. Mapping and quantifying contemporary vertical land motion (VLM) at continental scales has long been a challenge due to the absence of gridded observational products covering these large domains. Here, we fill this gap by using the new European Ground Motion Service (EGMS) to assess the current state of coastal VLM in Europe. First, we compare the InSAR‐based EGMS Ortho (Level 3) with nearby global navigation satellite systems (GNSS) vertical velocity estimates and show that the geodetic reference frame used to calibrate EGMS strongly influences coastal vertical land velocity estimates at the millimeter per year level and this needs to be considered with caution. After adjusting the EGMS vertical velocity estimates to a more updated and accurate International Terrestrial Reference Frame (ITRF2014), we performed an assessment of VLM in European low elevation coastal flood plains (CFPs). We find that nearly half of the European CFP area is, on average, subsiding at a rate faster than 1 mm/yr. More importantly, we find that urban areas and populations located in the CFP experience a near −1 mm/yr VLM on average (excluding the uplifting Fennoscandia region). For harbors, the average VLM is even larger and increases to −1.5 mm/yr on average. This demonstrates the widespread importance of continental‐scale assessments based on InSAR and GNSS to better identify areas at higher risk from relative SLR due to coastal subsidence.

Introduction

While the understanding and modeling of sea-level rise (SLR) due to ocean density and mass changes have greatly improved over the past few decades, relative SLR contributions due to vertical land motions (VLMs) remain a major source of uncertainty (Fox-Kemper et al., 2021; Slangen et al., 2023). It is critical to downscale global and regional SLR to local relative sea-level change as this is what causes coastal impacts and adaptation needs (Nicholls, Hanson, et al., 2021). In particular, land subsidence amplifies relative SLR, which in turn can strongly exacerbate exposure to coastal flooding, saltwater intrusion, erosion and loss of wetlands, and damage to infrastructure (Cooley et al., 2022; World Bank, 2010).

Nicholls, Lincke, et al. (2021) found that despite an average global relative SLR of 2.6 mm yr⁻¹ over the past two decades, coastal inhabitants are preferentially located in subsiding locations, and experience an average relative SLR up to four times faster (7.8–9.9 mm yr⁻¹). Locally, land subsidence of several meters can be observed over the 20th/early 21st century in major coastal cities founded on coastal lowlands and deltas around the world, especially in Asia such as Tokyo, Shanghai, Bangkok and Jakarta (Cao et al., 2021; Esteban et al., 2020; Nicholls, 2018), and other major cities such as New Orleans (Burkett et al., 2003) and Vancouver (Samsonov et al., 2014). Recently, Wu et al. (2022) analyzed subsidence rates in 99 coastal cities worldwide and found that for most of them, part of the land is subsiding faster than the global mean sea level is rising due to ocean density and mass change.

Although large cities are closely scrutinized as they concentrate most economic assets, subsidence is also typical across large deltaic regions that concentrate about half a billion people worldwide (Giosan et al., 2014; Minderhoud et al., 2019; Tessler et al., 2015). Erban et al. (2014) found that in the Mekong Delta, subsidence rates exceed the climate-induced SLR rate by an order of magnitude, so by the mid-21st century an additional inundation hazard of ∼1 m in portions of the Delta is possible. Recently, by combining groundwater extraction scenarios with climate-change driven SLR and a vertically high resolution digital elevation model of the Mekong Delta, Minderhoud et al. (2020) estimated that ∼40–70% of the Mekong delta could fall below mean sea level by the end of the 21st century depending on the degree of subsidence mitigation. Similarly, the Ganges-Brahmaputra-Meghna Delta experiences ∼2–3 mm/yr subsidence rates that are largest in the most populated area of the delta (Krien et al., 2019). Recently, a study covering the entire US Atlantic coast revealed that a subsidence rate exceeding 3 mm/yr affects most coastal areas, including wetlands, forests, agricultural areas, and developed regions (Ohenhen et al., 2023).

The causes of VLM are multiple (Shirzaei et al., 2021). VLM can be driven by natural processes, such as glacial isostatic adjustment (GIA), which manifests itself particularly in the regions that were covered by or adjacent to the thick ice sheets at the last glacial maximum ∼21,000 years ago like eastern North American and the Baltic region in Europe (Peltier, 2004; Spada, 2017). Tectonics and earthquakes can drive coastal subsidence at active plate margins (Leeper et al., 2017; Milker et al., 2016), so as sediment natural compaction, which is a dominant factor in major sediment depocenters such as deltas (Meckel et al., 2007; Teatini et al., 2011; Törnqvist et al., 2008; Zoccarato et al., 2018) and a well-known subsidence driver of coastal reclaimed land (e.g., Cavalié et al., 2023; Park & Hong, 2021). Finally, significant subsidence can be triggered by anthropogenic activities such as peat oxidation (e.g., Koster et al., 2018) or extraction of subsurface fluids (e.g., Galloway & Burbey, 2011; Hasan et al., 2023; Loschetter et al., 2015).

In Europe, two large subsiding coastal flood plain (CFP) regions have been long recognized and documented: that is, the North Italian coastal plain and the Netherlands. In both regions, the subsidence is caused by combined natural and anthropogenic effects that act at different timescales. In the Northern Italian coastal plain, a long-term component is controlled by tectonics and geodynamics and sedimentation, active on time spans of about 10⁶–10⁴ years (Carminati et al., 2003; Cuffaro et al., 2010; Kent et al., 2002). In a shorter-term (1,000 to 10,000 years), VLM in this region is also influenced by glaciation cycles and GIA processes (Antonioli et al., 2009; Stocchi et al., 2005). The latter results were mainly derived from sedimentological indicators. Over the recent historical period, various techniques based on leveling (e.g., Gatto & Carbognin, 1981), tide gauge analysis (Carbognin et al., 2004), height benchmarking (Zerbini et al., 2017), Synthetic Aperture Radar (SAR) (Tosi et al., 2013) or global navigation satellite systems (GNSS) (Zanchettin et al., 2021), allowed characterizing anthropogenic influence (namely fluids withdrawal) on the sinking of the North Italian Coastal plain, associated with an increasing threat of flooding (Lionello et al., 2021 and reference herein). In the Netherlands, subsidence is the result of combined influence of GIA (e.g., Meijles et al., 2018 and references herein), physical compaction of peat and its oxidation (e.g., Erkens et al., 2016; van Asselen et al., 2009), shrinkage and compression of clay (Fokker et al., 2019; Verberne et al., 2023), and ground fluid extraction related to gas production (Schoonbeek, 1976).

In addition to these two large European subsidence hot-spots, more local-scale ones such as for example, Thessaloniki (Raucoules et al., 2008; Stiros, 2001) have been reported and are associated with a major increase of flooding. Until recently, such coastal subsidence hot-spots could be identified through dedicated local to regional InSAR studies, GNSS and tide gauge measurements or leveling surveys (Poitevin et al., 2019; Raucoules et al., 2008, 2013; Wöppelmann et al., 2013). However, there has been to date no analysis of coastal subsidence at the whole European coast scale. This is the gap we propose to fill in the present study based on the recently released European Ground Motion Service (EGMS) from the Copernicus program (Costantini et al., 2021).

This study delivers the first consistent picture of the current state of VLM in low-lying coastal zones prone to marine flooding at the European scale. By doing so, we aim to demonstrate the usefulness of continental scale land motion assessments to better characterize areas at risk from subsidence in coastal zones exposed to SLR and identify challenges. To do so, we use the VLM data set provided by the InSAR-based EGMS service over the 2015–2021 period (Costantini et al., 2021). While our analysis allows the identification and analysis of well-known coastal subsidence hot-spots, it also allows the identification of unrecognized subsiding local scale coastal zones for which flood risk may become a greater concern than expected due to climate change alone. Ultimately, this study supports the development of climate services for coastal adaptation to SLR by providing a first assessment of where subsidence needs to be considered and added to sea-level scenarios in Europe.

The paper is organized as follows. First we describe the data sets used in this study to characterize VLM, coastal zones and associated economic and natural assets. Then, we make a comparison of EGMS VLM estimates with independent GNSS estimates and describe how VLM in CFPs are derived. The results are then presented at various spatial scales in units of CFP but also with respect to exposed assets within these CFPs. We conclude by an extended discussion on the implications of this work for future analyses in Europe and other regions of the world.

Data

In this study, we characterize contemporary VLM velocities within the European CFP using EGMS. EGMS InSAR-based estimates are further compared with GNSS products. Exposed assets and CFP data sets are also described below.

Vertical Land Motion

Three different types of data are used to characterize and quantify VLM: InSAR-based EGMS, GNSS-based velocity estimates and a GIA ensemble model. The terminology of the various vertical land velocity estimates employed hereinafter in this paper is summarized in Table 1.

Table 1 Description of Vertical Land Motion Velocities Estimates Used and Calculated in This Study

Notation	Description	Comments
VEGMS	Vertical velocity estimates from EGMS Ortho product	Raw (or native) vertical velocity data retrieved from EGMS data set. Section 2.1.1 details the data set.
VGNSS	Vertical velocity estimates from GNSS stations	Raw (or native) vertical velocity data retrieved from various GNSS solutions. Section 2.1.2 details the various GNSS data sets used.
VGIA	Vertical velocity estimates of land motions from GIA model	Raw (or native) vertical velocity data retrieved from one GIA model. Section 2.1.3 details the GIA model.
VEGMS/ADJ	Vertical velocity estimates from EGMS adjusted to ITRF2014 reference frame	Post-processed EGMS vertical velocity estimates according to the method described in Section 3.1.
Vresidual	VEGMS/ADJ − VGIA or vertical velocity with GIA contribution removed.	Vertical velocity estimates from EGMS adjusted to ITRF2014 reference frame minus velocity estimates from GIA model.
$$\overline{\mathrm{V}}$$ ($$\widetilde{\mathrm{V}}$$)	Mean (median) spatial velocity estimates computed over more than 1 pixel.	$${\overline{\mathrm{V}}}_{EGMS}$$ indicates for example, the spatial average of raw vertical velocity of all EGMS pixels over a defined area (e.g. one coastal flood plain).

EGMS

EGMS uses Synthetic Aperture Radar Interferometry (InSAR) data derived from Sentinel-1 to detect and measure ground movements across Europe with millimeter precision. The EGMS aims to provide consistent, updated, standardized, harmonized across national borders and reliable information regarding natural and anthropogenic ground motion phenomena over Europe (Crosetto & Solari, 2023). EGMS relies on the advances and improvements of the advanced Differential InSAR (A-DInSAR) data processing and analysis techniques, also referred to as Persistent Scatterer Interferometry (see Crosetto et al., 2016, 2021). EGMS products are being made by the consortium ORIGINAL (OpeRatIonal Ground motion INsar Alliance) comprising four different InSAR Processing Entities using methods proven over the two last decades. All algorithms and methods used to produce EGMS products are detailed in the document Algorithm Theoretical Basis Document (Ferretti et al., 2023).

Although the A-DInSAR technique is powerful to measure localized deformation, it can only provide spatially relative measurements on local spatial scales up to tens of km. On the other hand, GNSS measurements can be considered absolute and accurate over large spatial scales, but are very sparse. In this regard, EGMS exploits the synergy between InSAR and GNSS by calibrating the A-DInSAR product with a low-frequency deformation map obtained from a GNSS gridded model (Crosetto et al., 2020).

EGMS delivers three products that are updated annually:

• The Basic product (or Level 2a) which provides line of sight velocity maps in ascending and descending orbits with annotated geolocation and quality measures per measurement point. Basic products are referred to a local reference point.

• The Calibrated product (or Level 2b) also providing line of sight velocity maps in ascending and descending orbits but referenced to a model derived from GNSS time-series data so that the product is absolute in a geocentric sense (i.e., no longer relative to a local reference point but relative to a geocentric reference frame).

• The Ortho product (or Level 3), constructed from the Calibrated product using trigonometric calculation, which provides horizontal and vertical land velocity components resampled on a 100 m resolution grid and anchored to the geocentric reference frame model.

The EGMS calibration from Level 2b is obtained by referencing the Basic product (Level 2a) to the Augmented EUREF Permanent Network (A-EPND) model derived from GNSS data. The A-EPND gridded model is constructed based on GNSS solution data from the EUREF Densification product (EPND) mainly () and completed by input GNSS data from Nevada Geodetic Laboratory (NGL). The A-EPND model is calculated on a regular 50-km grid covering the whole Europe (∼3,800 grid points) and results in a long-wavelength (>50 km) 3D velocity field in the ETRF2000 realization of the ETRS89 reference frame. The Calibrated product is then obtained from the Basic product by calibrating its lower-frequency phase trends using the A-EPND model. The calibration procedure is based on common existing methodologies (see e.g. Farolfi et al., 2019). More information on the A-EPND model and calibration procedure are provided in the EGMS Algorithm Theoretical Basis Document (Ferretti et al., 2023), the EGMS GNSS Calibration Report (Larsen et al., 2023) and in the Text S1 in Supporting Information S1.

Our study relies on the first update of the Ortho product (Level 3) covering the period 2015–2021 (released in 2022). We focus on the vertical component. The Ortho data set was obtained from the EGMS web portal (), which provides 100 km × 100 km geotiff rasters of the vertical velocity best estimate over the period 2015–2021. The EGMS model applied to estimate the temporal coherence of the measured point and for computing the root mean square error of each displacement time series is a third order polynomial, plus a seasonal component (sinusoidal). Each 100 km × 100 km geotiff raster is accompanied with a tabular file (format.csv) providing, for each pixel (100 m × 100 m) of the raster file, the spatial coordinates, the root mean square error of the time series residuals from the fitting procedure, the mean vertical velocity, acceleration and seasonality, and their standard error, respectively, and the full vertical position time series. In this study, we relied on the contemporary mean vertical velocity estimates by the EGMS fitting model at each pixel (hereinafter V_EGMS, see Table 1). Overall, the Ortho product provides a mean velocity estimate with a precision better than 1 mm/yr (Capes & Passera, 2023; Sala Calero et al., 2023).

To cover the entire European coastline, we downloaded 415 EGMS tiles of 100 km × 100 km spatial extent that we merged into one large pan-European raster file and clipped to the European coastal mask defined within the Dynamic Interactive Vulnerability Assessment (DIVA) modeling framework (see further details in Section 2.2). This leads to more than 10⁶ V_EGMS (i.e., 10⁶ filled pixels) at the scale of Europe.

Global Navigation Satellite System

Contemporary V_EGMS are compared with vertical velocity estimates from continuously operating GNSS stations (hereinafter V_GNSS, see Table 1) in Europe. We used V_GNSS distributed by the SONEL (Système d’Observation du Niveau des Eaux Littorales) web portal (). SONEL serves as the GNSS data assembly center for the Global Sea Level Observing System (GLOSS), which is developed under the auspices of the IOC/UNESCO. This data assembly center, hosted at La Rochelle University in France (Wöppelmann et al., 2021), aims to provide high-quality continuous measurements of sea- and land-levels at the coast from tide gauges (French-only scope) and modern geodetic techniques, especially permanent GNSS stations co-located at or nearby tide gauges (global scope).

The SONEL portal currently distributes V_GNSS at coastal stations worldwide from four different solutions (groups): ULR7A from University of La Rochelle (Gravelle et al., 2023), NGL14 from NGL (Blewitt et al., 2016, 2018), JPL14 from Jet Propulsion Laboratory (Heflin et al., 2020) and GT3 from the German Research Center for Geosciences in Postdam (Männel et al., 2022). The four solutions express their vertical velocity field with respect to the International Terrestrial Reference Frame 2014 (ITRF2014) (Altamimi et al., 2017). In spite of adopting the same International GNSS Service standards and conventions, the different analyst choices (beyond the agreed standards and conventions) at any of the major GNSS processing steps can lead to significant differences in vertical velocity estimates between solutions for a same GNSS station (e.g., Ballu et al., 2019). Here, we compared V_EGMS with V_GNSS from ULR7A and NGL14 solutions as these solutions offered the largest amount of coastal GNSS stations to compare with, and represent two independent GNSS data analysis strategies, the so-called “double difference” and “precise point positioning” in GNSS jargon.

Finally, our analysis of V_GNSS also investigates vertical velocity estimates of the EUREF Permanent GNSS Network (EPN) solution (Bruyninx et al., 2019) as this solution was used for calibrating V_EGMS (Larsen et al., 2023). More specifically, we used the EPN D2150 densification solution (Kenyeres et al., 2019), which retrieved vertical velocities for nearly 2,500 stations in Europe.

GIA

Ongoing changes in the solid Earth caused by the viscous adjustment of the mantle to the substantial mass redistribution that followed the Last Glacial Maximum (named Glacial Isostatic Adjustment; GIA) have an important contribution to VLM in Europe (Peltier & Tushingham, 1991). Typically, the GIA causes a land uplift in the Fennoscandia region (contributing to relative sea-level fall) and a subsidence in its near periphery (i.e., a relative SLR); for example, The Netherlands coast (Vermeersen et al., 2018). Despite the well-established effect of the GIA in Northern Europe, there remain uncertainties on its magnitude and spatial distribution that are reflected by the substantial discrepancies between the output of various existing GIA models (e.g., Bradley et al., 2011, 2023; Kuchar et al., 2012; Lambeck et al., 1998; Meijles et al., 2018; Peltier et al., 2015; Simon et al., 2018).

Here, we rely on the GIA model of Caron et al. (2018), which was built from an ensemble of 128,000 GIA simulations. Ensemble simulations are designed by varying solid-Earth parameters (lithosphere thickness and mantle viscosities) and amplitudes of global deglaciation histories over the past 20,000 years. Each GIA ensemble member provides a consistent set of changes in relative sea level, solid-Earth deformation and changes in equivalent water height and comes with a likelihood that reflects how good the fit is to a data set of vertical GNSS velocities and paleo sea-level records. Therefore, the model allows for a robust quantification of the uncertainties associated with GIA. The spread between the ensemble members depicts the uncertainty in the GIA predictions due to uncertainty in the solid-Earth parameters and the deglaciation history. This GIA model was used to reconstruct and successfully close the basin-mean sea-level budget over the 20th century (Frederikse et al., 2020). The gridded GIA data set of corresponding relative sea level and solid earth deformation mean rates and uncertainty are retrieved from Frederikse, Landerer, and Caron (2019).

Coastal Mask and Flood Plains

Since the study focuses on coastal regions exposed to coastal flooding and SLR, we restrict our analysis to coastal areas and CFPs. The coastal mask and related CFPs were extracted from the DIVA model (Lincke et al., 2022). DIVA has been developed for broad-scale coastal assessment of SLR and used to assess issues such as coastal erosion (Hinkel et al., 2013), coastal flooding (Hinkel et al., 2014), coastal wetland change (Schuerch et al., 2018), coastal migration and adaptation (Lincke & Hinkel, 2021) and subsidence/relative SLR (Nicholls, Lincke, et al., 2021).

More specifically, we used the DIVA coastal mask to construct a V_EGMS coastal raster (i.e., restricting EGMS coverage to the European coastal region). The DIVA coastal mask allows identifying coastal area, which means hydrologically connected to the ocean by an elevation below 20 m or within 1 km distance to the coastline. The resulting V_EGMS coastal raster covers therefore the entire European coastline (i.e., not solely restricted to low-lying CFPs) and is used for V_EGMS/V_GNSS comparisons (see Section 3.1).

Now to focus on flood-prone areas, a CFPs vectorized layer developed within the frame of the Coastal Climate Core Services H2020 project (CoCliCo; ) is used in addition. CFPs are contoured using the 30-m resolution Copernicus Digital Surface Model GLO-30 (European Space Agency Sinergise, 2021), considering a bathtub flooding approach (including hydraulic connections) calculated for a 1-in-100-year event combined with a 2-m SLR. This results, at the scale of Europe, in a total of about 41,000 CFPs spanning a wide range of sizes with a minimum CFP area of 0.01 km² (i.e., equivalent to one EGMS pixel size) and a maximum CFP area of 17.6 × 10³ km². The median CFP area is 0.06 km² while the mean CFP area is 2.5 km²: that is, the distribution of CFPs is strongly positively skewed. The European CFP coverage is shown on Figure 1.

[IMAGE OMITTED. SEE PDF]

Exposure

The pan-European Corine Land Cover (CLC) for the reference year 2018 (CLC2018) is used to characterize land use and land cover within CFPs. The CLC nomenclature consists of 44 land use and land cover classes (Bossard et al., 2000; Feranec et al., 2016). The classes are defined and explained in an illustrated technical guideline (Kosztra et al., 2017). The CLC2018 raster layer is provided at a minimum feature width of 100 m at the scale of Europe, which has the same resolution as the EGMS Level-3 Ortho product. Note that among the 44 land classes, we focus mainly on infrastructures assets (i.e., urban fabrics, networks, commercial units, airports/harbor) and typical coastal natural environments (e.g., marshes).

To analyze the exposed population, we used the population gridded data from Eurostat for the reference year 2021. The total EU population measured on the 1 km² population grid in the census reference year 2021 was 443.2 million residents. The data sets are based on data from the 2021 population and housing census at the level of 1 km² grid cells (Lutz et al., 2019).

Methods and Validation

In this section, we compare vertical land velocities estimates obtained from EGMS (V_EGMS) and GNSS stations (V_GNSS) and we detail the method that has been developed to calculate VLM in CFP. Hereinafter, subsiding and uplifting land motion velocities are displayed in blue and red, respectively.

Comparison Between EGMS and GNSS Vertical Velocities

EGMS Level-3 Ortho is a highly relevant and unique product to characterize the spatial distribution of current VLM as its coverage and resolution is unprecedented at the scale of Europe. This contrasts with GNSS station networks that provide precise and often long-term records (some records are longer than 20 years) of vertical land velocity estimates that are however sparsely distributed spatially. SONEL V_GNSS estimates are independent from EGMS and thus highly relevant to assess V_EGMS estimates.

Figure 2 shows a comparison between V_GNSS retrieved from 93 GNSS stations distributed across the European coastline and from the ULR7A solution (available on SONEL), and V_EGMS estimates from the EGMS pixels in which the GNSS stations are located. An example of the spatial location of GNSS stations and the EGMS spatial coverage along the Bay of Biscay coastal region is shown on Figure 2b. The scatter plot (Figure 2a) shows a good correlation between V_GNSS and V_EGMS estimates (Pearson's correlation coefficient of 0.94). Similarly, the correlation coefficient between V_GNSS from the NGL14 solution (171 stations, available on SONEL) and V_EGMS is 0.92 (Figure S1 in Supporting Information S1). Note that this analysis was repeated by considering the averaged V_EGMS over pixels located within a 200 m radius around the GNSS stations and the results of the comparison were very similar (Figure S2 in Supporting Information S1).

[IMAGE OMITTED. SEE PDF]

Our comparison hence reveals an overall good agreement between V_EGMS and V_GNSS, however, the linear regression analysis between both estimates suggests that V_EGMS tends to be more negative than V_GNSS as suggested by the intercept of −0.89 mm/yr, and this effect strengthens for large positive land motion velocities (i.e., large ground uplift region). Similar offsets are found when regressing V_EGMS over V_GNSS from the NGL14 solution (Figure S1 in Supporting Information S1). This seeming negative V_EGMS bias against V_GNSS is investigated in more detail below.

As mentioned in 2.1.1 and 2.1.2, V_EGMS are calibrated with respect to EUREF GNSS network (EPN) in the European geodetic reference frame ETRF2000. In this reference frame, the Southern Europe appears to experience an average subsidence of between 0.5 and 1 mm/yr (see EGMS GNSS Calibration report; Larsen et al., 2023). Figure 2b reflects well this effect and shows that most EGMS pixels from the Bay of Biscay coastal region have a V_EGMS comprised between −2 and 0 mm/yr (orange and yellow pixels).

Next, we compare V_GNSS from the EPN D2150 densification solution expressed in the ETRF2000 geodetic reference frame (and used to calibrate V_EGMS) to those of the same GNSS stations but expressed with respect to the ITRF2014 (international geodetic reference frame, common to ULR7A and NGL14 solutions available on SONEL). The results of this comparison are shown in Figure 3. It appears that V_GNSS from EPN D2150 expressed in the ETRF2000 geodetic reference frame are systematically more negative than when expressed in the ITRF2014, and that this effect amplifies toward higher latitudes. This result is consistent with the previous (see Figure 2a) comparative analysis of V_GNSS from ULR7A and NGL14 solutions (expressed in the ITRF2014) and V_EGMS (calibrated in the ETRF2000). Also, the strengthening of the deviation toward large land motion uplift (Figure 2a) is consistent with the amplifying differences in high latitudes of Europe where the land motion uplift caused by GIA is large (Figure 3a).

[IMAGE OMITTED. SEE PDF]

The ITRF is the most precise and stable geocentric reference frame achieved so far by the international community within the International Earth Rotation and Reference Systems framework (see ) and is regularly updated (e.g., ITRF2000, ITRF2005, ITRF2008, ITRF2014 and ITRF2020). Note that the 2000 (2005, etc.) in ITRF2000 underlines the latest geodetic observations considered in the ITRF release. Although not the most recent one, the ITRF2014 still remains a reference frame solution in the data analysis strategies of most groups (e.g., ULR, JPL, NGL) for which highly demanding applications in terms of precision and stability of positioning are required. In this regard, we apply a simple transfer function allowing to adjust V_EGMS calibrated within ETRF2000 into the ITRF2014 (V_EGMS/ADJ hereinafter) to illustrate the importance of the geodetic reference frame choice on EGMS calibration to analyze coastal VLM. The transfer function is shown on Figure 3b, which shows that the difference of vertical velocities estimates between ETRF2000 to ITRF2014 can be adequately modeled by a second order polynomial with respect to the latitude (the determination coefficient of the fit, or R², is 0.99). Note that for the purpose of this study, such a level of precision in the predictive model (root mean square error lower than 0.02 mm/yr) is sufficient.

Assessing VLM in CFPs

Figure 4 illustrates the method that we developed to assess CFP VLM based on contemporary V_EGMS estimates. The Arcachon bay case is used as an example (see also Figure 2b).

[IMAGE OMITTED. SEE PDF]

For each CFP, we first calculate the fraction of spatial area that is covered by EGMS pixels. Figure 4a reveals that the pixel coverage can vary substantially from one CFP to another. For instance, the island located in the middle of the bay (called “l’île aux oiseaux”), which is one entire CFP, is covered by EGMS pixels over less than 20% of its area (lighter gray contour). In contrast, the narrow CFPs located in the northern part of the bay are almost entirely covered by EGMS pixels (more than 80%, black contours). The largest CFP of the bay of Arcachon, along the southern shore, has an EGMS pixel coverage ratio of nearly 50%.

The ratio of EGMS pixels depends on the nature of land coverage and its capacity to reflect the radar signal, but also depends on the area of CFP. Typically, urban areas are very well covered by EGMS pixels in contrast to agricultural land or forest areas. The ratio of EGMS pixels for all 40,746 European CFPs is reported in Table 2. Nearly one quarter of CFPs are empty of EGMS pixels, which however, corresponds to less than 5% of the total CFP coverage in Europe and concerns mostly Norway, Sweden and Finland. Table 2 also shows that CFPs with the lower EGMS pixel coverage ratio are also the most extended ones; they correspond typically to very large low-lying areas that embed various land uses such as infrastructures, networks and urban areas surrounded by agricultural land areas (e.g., in the Netherlands). Therefore, for extended CFPs with a low EGMS coverage ratio, the V_EGMS estimates may not necessarily be representative of the entire CFP behavior, but rather reflect the VLM state of urban assets within those large CFPs.

Table 2 CFPs Coverage by European Ground Motion Service Data (Total; 40,746 CFPs)

EGMS coverage	#CFPs	Total area (103 km2)	Mean CFP area (km2)	Median CFP area (km2)	Max CFP area (km2)
]0; 20]%	2,886	69.0	24	0.5	17.5 × 103
]20; 40]%	3,068	20.6	6.7	0.1	6.0 × 103
]40; 60]%	2,999	5.2	1.7	0.1	2.4 × 103
]60; 80]%	3,900	1.2	0.3	<0.1	154
]80; 100]%	16,613	2.1	0.1	<0.1	50

To obtain an overall view of the VLM velocities in the CFPs, we derive $${\overline{\mathrm{V}}}_{\text{EGMS}}$$ by averaging all V_EGMS within each CFP (Figure 4b). As illustrated on Figure 4b, we find that the largest CFP of the Arcachon bay has a $${\overline{\mathrm{V}}}_{\text{EGMS}}$$ of −2.2 mm/yr (red), hence suggesting a significant subsidence. After adjusting to the ITRF2014 reference frame ($${\overline{\mathrm{V}}}_{\text{EGMS}/\text{ADJ}}$$) (see Figure 3), the VLM velocities of the later CFP is reduced to −1.6 mm/yr (Figure 4c, orange). Note that, for each CFP, we also assess the median (noted $$\widetilde{\mathrm{V}}$$ hereinafter, see also Table 1), the first quartile and the third quartile velocities systematically (see Figure S3 in Supporting Information S1). For most CFPs containing EGMS pixels (i.e., ∼70%), we found that $${\overline{\mathrm{V}}}_{\text{EGMS}}$$ is lower than $${\widetilde{\mathrm{V}}}_{\text{EGMS}}$$; this means that the distribution of V_EGMS within CFPs are generally skewed toward more subsiding values.

Results

Our results are presented in this section. We first analyze the overall state of subsidence in Europe. We then assess in more detail a selection of subsiding coastal hot-spots identified in our data set. Finally, we assess how subsidence may enhance coastal risk in Europe over the 21st Century by combining our EGMS-based data set with assets and exposure inferred from the CLC and population.

State of Vertical Land Motion of European CFPs

The distribution of mean ($$\overline{\mathrm{V}}$$) and median ($$\widetilde{\mathrm{V}}$$) CFP's contemporary VLMs for all of Europe's coast is shown on Figure 5. The distribution is expressed either with respect to the number of CFPs (Figures 5a and 5b), or with respect to the associated CFP area (Figures 5c and 5d). To emphasize the importance of the geodetic reference frame on vertical land velocity estimates, we analyze both the EGMS Ortho product ($${\overline{\mathrm{V}}}_{\text{EGMS}}$$ and $${\widetilde{\mathrm{V}}}_{\text{EGMS}}$$, see also Figure 4b) and the EGMS Ortho product adjusted to ITRF2014 ($${\overline{\mathrm{V}}}_{\text{EGMS}/\text{ADJ}}$$ and $${\widetilde{\mathrm{V}}}_{\text{EGMS}/\text{ADJ}}$$, see also Figure 4d). To quantify how the GIA contributes to VLM in Europe, we also calculate and analyze $${\overline{\mathrm{V}}}_{\text{residual}}$$ and $${\widetilde{\mathrm{V}}}_{\text{residual}}$$; that is, the mean and median EGMS Ortho product adjusted to ITRF2014 to which the GIA effect on solid-Earth vertical deformation has been subtracted. In addition to mean and median estimates, the spatial variability of vertical velocities estimates within each CFP has been analyzed using the first and third quartiles (Figures S3 and S4 in Supporting Information S1). We find that interquartile width for each CFP remains mostly below 1 mm/yr and tends to increase for large CFPs.

[IMAGE OMITTED. SEE PDF]

CFPs that are uplifting on average (i.e., $$\overline{\mathrm{V}}$$ > 0.5 mm/yr) represent a large majority of CFPs in Europe with more than 20,000 CFPs for both $${\overline{\mathrm{V}}}_{\text{EGMS}}$$ and $${\overline{\mathrm{V}}}_{\text{EGMS}/\text{ADJ}}$$ estimates (Figure 5a). The same conclusion is draw from median velocity estimates (Figure 5b). Figures 5c and 5d reveals however that despite the large number of uplifting CFPs, the corresponding spatial extent remains modest as it comprises less than 10% of the total European CFP area both $${\overline{\mathrm{V}}}_{\text{EGMS}}$$ and $${\overline{\mathrm{V}}}_{\text{EGMS}/\text{ADJ}}$$ estimates. This is typically related to small CFPs affected by uplifting effect of the GIA in the Northern European countries. The latter result is also consistent with the fact that when the vertical velocity from GIA model is subtracted from $${\overline{\mathrm{V}}}_{\text{EGMS}/\text{ADJ}}$$ (i.e., $${\overline{\mathrm{V}}}_{\text{residual}}$$ on Figure 5a) or $${\widetilde{\mathrm{V}}}_{\text{EGMS}/\text{ADJ}}$$ (i.e., $${\widetilde{\mathrm{V}}}_{\text{residual}}$$ on Figure 5b), a large amount of uplifting CFPs is moved toward more stable conditions (i.e., range −0.5 to 0.5 mm/yr).

In contrast to uplifting CFPs, the number of CFPs that are on average subsiding at a rate faster than 1 mm/yr is low: 5,000 CFPs for $${\overline{\mathrm{V}}}_{\text{EGMS}}$$ and 2,500 CFPs for $${\overline{\mathrm{V}}}_{\text{EGMS}/\text{ADJ}}$$ (Figure 5a). First, we note that this difference between $${\overline{\mathrm{V}}}_{\text{EGMS}}$$ and $${\overline{\mathrm{V}}}_{\text{EGMS}/\text{ADJ}}$$ is large, which demonstrates that the geodetic reference frame used to express the vertical velocity has a large effect on CFP VLM assessment and needs to be acknowledged. Second, the 2,500 subsiding CFPs (out of ∼30,000 CFPs containing EGMS pixels) determined from $${\overline{\mathrm{V}}}_{\text{EGMS}/\text{ADJ}}$$ represents nearly 50% of the total CFP area in Europe (Figure 5c). Overall, our analysis of $${\overline{\mathrm{V}}}_{\text{EGMS}/\text{ADJ}}$$ suggests that 50% (10%) of the European land area located in CFPs experience an average subsidence stronger than 1 mm/yr (2 mm/yr). We also find that 80% of CFPs that have an average subsidence stronger than 1 mm/yr are covered by 75% of pixels with V_EGMS/ADJ < −0.5 mm/yr (Figure S4 in Supporting Information S1). This means that despite some spatial variability, subsiding CFPs are mostly covered by pixels with negative velocity estimates. Finally, Figures 5c and 5d also reveal that, when subtracting the vertical velocity from the GIA model (see $${\overline{\mathrm{V}}}_{\text{residual}}$$ or $${\widetilde{\mathrm{V}}}_{\text{residual}}$$), more than 60% of the CFP total area becomes stable (range −0.5 to 0.5 mm/yr). This indicates that GIA is a first order contributor to both low-lying coastal land uplift and subsidence.

Next, we examine the distribution of $${\overline{\mathrm{V}}}_{\text{EGMS}/\text{ADJ}}$$ in CFPs per European coastal country. Results are shown in Figure 6 (see also Figure 1 for visualizing European CFPs). Here, vertical range velocity classes are simplified as we assume that an uplifting (subsiding) CFP corresponds to a $${\overline{\mathrm{V}}}_{\text{EGMS}/\text{ADJ}}$$ larger (lower) than 1 mm/yr (−1 mm/yr) and is considered as stable otherwise (i.e., in the range −1 to 1 mm/yr).

[IMAGE OMITTED. SEE PDF]

As expected, CFPs of northern European countries (i.e., Norway, Sweden, Finland and Estonia) are mostly uplifting on average as a consequence of the GIA effect. Just southward of these four uplifting countries, Latvia, Lithuania, Denmark, UK and Ireland show rather stable VLMs, and further south, a widespread subsidence (with exceptions such as Belgium). More specifically, In Germany and the Netherlands, a vast majority (near 75% of the total area) of CFPs are experiencing an average subsidence, which is consistent with the combined effect of GIA and more localized subsidence of other origin (see Introduction). Among the countries of southern Europe with expanded CFPs, Italy, Greece and Romania have more than 75% area of their CFPs that are subsiding on average. Finally, in the Western Europe, more than half of CFPs in France and Portugal are subsiding on average.

The analysis of the distribution of $${\widetilde{\mathrm{V}}}_{\text{EGMS}/\text{ADJ}}$$ in CFPs per European coastal country (i.e., median instead of mean) reveals a slightly different picture (see Figure S5 in Supporting Information S1). While the ratio of uplifting and subsiding CFPs remains nearly unchanged for Fennoscandia and Mediterranean countries, it drastically changes for the Netherlands, Germany and France, where most CFPs now fall into the stable category. This illustrates that the VLM distribution is skewed, which becomes particularly pronounced when CFPs are large (see Discussion section).

Subsiding Hot-Spots Detected in EGMS

In the previous section, EGMS vertical velocities were processed per CFP to provide an overview of the current state of VLM in low-lying vulnerable coastal regions in Europe. Here, we focus on contemporary V_EGMS at a finer scale and illustrate how pixel-scale vertical velocities can characterize subsiding hot-spot regions. Note that V_EGMS/ADJ and V_residual are not investigated in this section as the goal of analysis is a qualitative description of very large and localized land subsidence signals, for which reference frame adjustment and GIA influence are secondary.

We selected four subsidence hot-spots that are illustrated on Figure 7. Two are large scale, well-known and well documented: The Netherlands (see e.g. Erkens et al., 2016), and the northern Italian coastal plain (see e.g. Fiaschi et al., 2018; Tos et al., 2016). The analysis is completed with two more local hot-spots: Thessaloniki in Greece and its surrounding area, reported by Stiros (2001), and the city of Schiavonela in Southern Italy and its periphery, which is here documented as a subsidence hot-spot for the first time.

[IMAGE OMITTED. SEE PDF]

In the Netherlands (panel 1 of Figure 7), V_EGMS estimates are lower than −1 mm/yr throughout the CFPs (light blue, blue and black pixels). Subsidence tends to diminish landward as revealed by the yellow pixels. The subsidence rate is however not uniformly distributed spatially at the scale of the Netherlands CFPs. In particular, the south-west (around Rotterdam) is characterized by several extended areas where vertical land velocities are lower than −5 mm/yr (black pixels). Subsidence in these areas is consistent with groundwater level lowering in the Holocene that caused oxidation of peat organic matter (Koster et al., 2018). In the north-east, in the area of Groningen, a large subsidence hot-spot, with velocities than can locally be lower than −5 mm/yr (even −10 mm/yr) are measured. The latter has been well documented and attributed to gas extraction (van Thienen-Visser et al., 2016). Between these two distant subsidence hot-spots, the rest of the Netherlands appears to be affected by subsidence, but to a lesser magnitude. This spatial heterogeneity of subsidence in the large CFPs of the Netherlands is very likely one of the primary reasons why the mean and median EGMS velocity over CFPs deviate (see discussion).

The northern Italian coastal plain is also well-known for being a subsidence hot-spot. Not surprisingly V_EGMS confirm this extended subsidence, as shown in panel 2 of Figure 7. In contrast with the Netherlands CFPs case, the subsidence pattern appears to be more uniform in the large-scale northern Italian CFPs with V_EGMS typically in the range −2 to −6 mm/yr. Since the CFPs in Northern Adriatic represent most of the Italian CFPs area and given the rather homogenous distribution of subsidence in them, the median and mean of CFP-scale subsidence estimates are similar (Figure 6 and Figure S5 in Supporting Information S1).

In Southern Italy, our coastal EGMS analysis identified a previously unreported coastal subsiding region in the vicinity of the city of Schiavonea and its fishing harbor (panel 3 of Figure 7). While the corresponding CFPs are not extensive, they experience strong subsidence that is mostly recorded below −5 mm/yr and co-located with exposed agriculture and urban areas. Note that the V_GNSS retrieved from the NGL solution for the sole GNSS station nearby (CORC station) also reveals large negative land movement velocities of nearly −5 mm/yr (calculated over less than 4 years, not shown). This very local example demonstrates the added value of EGMS to identify new and strongly subsiding coastal regions.

The last subsidence hot-spot example is Thessaloniki (panel 4 of Figure 7), which has already been documented (Raucoules et al., 2008; Stiros, 2001). In this case, a strong subsidence (i.e., V_EGMS lower than −5 mm/yr) is measured in two urban plots of the Aksiou delta. The Thessaloniki airport, located on the southern shore of the bay in a CFP is also experiencing strong subsidence with many EGMS pixels that observed V_EGMS lower than −10 mm/yr. In this area again, the subsidence pattern is highly heterogeneous.

In addition to identifying subsidence hot-spots it is also interesting to identify coastal lowlands that have limited subsidence. Both Belgium and the UK have extensive coastal lowlands that have experienced significant loss in land elevation due to drainage and shallow subsidence (de la Vega-Leinert & Nicholls, 2008; Lebbe et al., 2008). For example, the Fens in eastern England has lost up to 5-m elevation since the mid 19th Century (Dawson et al., 2010). Loss of elevation is often assumed to be ongoing. However, in the period of the EGMS observations considered in this paper any subsidence appears negligible except locally (not shown) and contrasts with the observations already discussed in the Netherlands.

Exposed Assets and Population

To assess how subsidence could inflate coastal risk in low-lying coastal areas, we now relate the CFP contemporary VLM to European land cover and population, provided by CLC and Eurostat data sets, respectively (see Sections 2.2 and 2.3). As the purpose here is a quantitative assessment at pan-European scale, the analysis below relies on $${\overline{\mathrm{V}}}_{\text{EGMS}/\text{ADJ}}$$ (i.e., within the ITRF2014 reference frame) rather than $${\overline{\mathrm{V}}}_{\text{EGMS}}$$.

Table 3 relates the land cover to the state of $${\overline{\mathrm{V}}}_{\text{EGMS}/\text{ADJ}}$$ in CFPs in all Europe. The analysis is restricted to a subset of 9 classes (among 44) of CLC those land cover generally include highly exposed anthropogenic (i.e., cities, airports, commercial centers,…) or coastal natural (i.e., marshes) assets. The results reveal that a large majority (more than 70%) of assets within CFPs in Europe are located in CFPs with $${\overline{\mathrm{V}}}_{\text{EGMS}/\text{ADJ}}$$ lower than −0.5 mm/yr. For instance, 92% of the area occupied by transport networks is located in CFPs with $${\overline{\mathrm{V}}}_{\text{EGMS}/\text{ADJ}}$$ lower than −0.5 mm/yr. The ratio of exposed assets is still high when considering CFPs with $${\overline{\mathrm{V}}}_{\text{EGMS}/\text{ADJ}}$$ lower than −1 mm/yr: this is for example, the case for nearly half of discontinuous urban fabric, industries, transport networks or inland marshes located in CFPs. Consistently, half of the European population living in the CFP (which represents ∼16 million inhabitants in the Eurostat data set) is also located in CFPs with $${\overline{\mathrm{V}}}_{\text{EGMS}/\text{ADJ}}$$ lower than −1 mm/yr. Finally, near 10% of the CFPs area covered by exposed assets is located in CFPs with $${\overline{\mathrm{V}}}_{\text{EGMS}/\text{ADJ}}$$ lower than −2 mm/yr.

Although Table 3 provides an overall appraisal of the exposure to subsidence at the CFP scale, the analysis can be further refined by calculating within each CFP, the $${\overline{\mathrm{V}}}_{\text{EGMS}/\text{ADJ}}$$ (i.e., the average) only over pixels associated with each CLC class. The refined analysis is illustrated in Figure 8 for the CLC class “urban fabric” (both continuous and discontinuous combined). Here, for each pixel of CFPs containing the urban fabric class, we select the co-located V_EGMS/ADJ, calculate the average within each CFP (i.e., $${\overline{\mathrm{V}}}_{\text{EGMS}/\text{ADJ}}$$) and draw the corresponding European distribution per CFP (Figure 8). This concerns 4,641 CFPs among ∼41,000. It reveals a bimodal distribution with a first peak slightly negative (near −1 mm/yr) and a second much smaller peak toward large positive values (∼9 mm/yr). This second peak is related to the influence of GIA in Northern Europe. This refined analysis has been repeated for the CLC classes “Harbor” and “Airport” and for the median, first and third quartiles of V_EGMS/ADJ to reflect the spread of the estimates. The main results are summarized in Table 4.

Table 3 Distribution, With Respect to Area, of Different Land Cover Types (e.g., Airports) and Exposure in Coastal Flood Plain as a Function of $${\overline{\mathrm{V}}}_{\text{EGMS}/\text{ADJ}}$$ or $${\widetilde{\mathrm{V}}}_{\text{EGMS}/\text{ADJ}}$$ Ranges of the CFPs

VLM in CFP exposure	]+∞;−0.5] mm/yr	]−0.5;−∞[ mm/yr	]−1;−∞[ mm/yr	]−2;−∞[ mm/yr
$${\overline{\mathrm{V}}}_{\text{EGMS}/\text{ADJ}}$$	$${\widetilde{\mathrm{V}}}_{\text{EGMS}/\text{ADJ}}$$	$${\overline{\mathrm{V}}}_{\text{EGMS}/\text{ADJ}}$$	$${\widetilde{\mathrm{V}}}_{\text{EGMS}/\text{ADJ}}$$	$${\overline{\mathrm{V}}}_{\text{EGMS}/\text{ADJ}}$$	$${\widetilde{\mathrm{V}}}_{\text{EGMS}/\text{ADJ}}$$	$${\overline{\mathrm{V}}}_{\text{EGMS}/\text{ADJ}}$$	$${\widetilde{\mathrm{V}}}_{\text{EGMS}/\text{ADJ}}$$
Continuous urban fabric	25%	38%	73%	60%	31%	12%	5%	1%
Discontinuous urban fabric	14%	27%	85%	72%	50%	10%	7%	2%
Industrial or commercial units	15%	28%	84%	70%	49%	13%	8%	3%
Road and rail networks and associated land	8%	15%	92%	85%	52%	17%	10%	3%
Port areas	22%	31%	76%	67%	42%	22%	10%	4%
Airports	28%	41%	68%	55%	45%	20%	13%	5%
Inland marshes	15%	22%	81%	74%	66%	20%	17%	2%
Salt marshes	27%	43%	68%	52%	33%	17%	10%	6%
Salinas	12%	35%	79%	56%	44%	24%	14%	12%
Population	22%	27%	74%	69%	49%	16%	7%	3%

[IMAGE OMITTED. SEE PDF]

Table 4 Mean of $${\overline{\mathrm{V}}}_{\text{EGMS}/\text{ADJ}}$$, $${\widetilde{\mathrm{V}}}_{\text{EGMS}/\text{ADJ}}$$, First and Third Quartiles of V_EGMS/ADJ for Urban Fabric, Harbor, and Airport Land Cover

	All Europe	CFP with no uplift GIA
Urban fabric (4641)	No weight	0.57 mm/yr (0.40/0.63/0.81)	−0.57 mm/yr (−0.74/−0.51/−0.35)
Population weighted	−0.52 mm/yr (−0.82/−0.35/−0.04)	−0.86 mm/yr (−1.17/−0.68/−0.35)
Area weighted	−0.48 mm/yr (−0.80/−0.29/0.04)	−0.79 mm/yr (−1.10/−0.59/−0.25)
Harbor (1116)	No weight	−0.51 mm/yr (−0.81/−0.39/−0.09)	−0.95 mm/yr (−1.23/−0.86/−0.58)
Area weighted	−1.01 mm/yr (−1.51/−0.64/−0.11)	−1.46 mm/yr (−1.91/−1.12/−0.60)
Airport (188)	No weight	0.52 mm/yr (0.26/0.66/0.94)	−0.66 mm/yr (−0.94/−0.52/−0.24)
Area weighted	−0.55 mm/yr (−0.93/−0.31/0.11)	−1.07 mm/yr (−1.47/−0.81/−0.37)

The mean of $${\overline{\mathrm{V}}}_{\text{EGMS}/\text{ADJ}}$$ co-located with urban fabric (i.e., mean of the distribution shown on Figure 8) is 0.57 mm/yr when considering the whole Europe (Table 4). If now the mean is weighted by either the population or the area of the corresponding land cover type within each CFP, we find a similar (and consistent) subsiding velocity of urban fabric areas of nearly 0.5 mm/yr. If the CFPs located in a region where GIA triggers a land uplift are discarded, the mean $${\overline{\mathrm{V}}}_{\text{EGMS}/\text{ADJ}}$$ for urban fabric experience stronger subsidence (0.8 mm/yr for area- and population-weighted estimates). In other words, this analysis shows that the European population settled in coastal flood prone areas lives on average where VLM of nearly −1 mm/yr is measured. For airports and harbors, subsidence rates are even stronger (−1 to −1.5 mm/yr). For ports, which are directly connected to the coast and for which sea level is a critical operating parameter, this means that they currently experience an additional relative SLR rate of 1 mm/yr on average. For some ports, this can be even larger. Finally, quartiles estimates display mostly consistent signs with the mean estimates (Table 4). Hence, despite the fact that vertical land velocities estimates vary spatially within CFPs, the mean estimate is relevant to capture the global land motion behavior of CFPs and associated land cover types.

Discussion

Significance of Our Main Findings

In this paper, we present for the first time a pan-European analysis of contemporary VLM and subsidence in low-lying flood prone coastal zones using the EGMS product. First, we perform a comparison between EGMS Ortho (Level 3) vertical velocity estimates and GNSS stations vertical velocity estimates from two different solutions (ULR7A & NGL). This comparison reveals that the geodetic reference frame used to calibrate EGMS can strongly affect the vertical land velocity estimates and needs to be accounted for carefully, especially for coastal applications where the need of sub-mm/yr precision in VLM is required for comparison with SLR changes (as opposed e.g. to landslides or earthquakes that are also characterized by EGMS). As coastal sea-level applications require precision and stability of positioning of the geodetic reference frame, VLM are generally expressed within ITRF (Wöppelmann & Marcos, 2016). We could hence determine and apply a simple polynomial function depending on latitude to transform EGMS VLM velocities from ETRF2000 (used to calibrate EGMS) to ITRF2014 in order to compensate the identified bias. Overall this leads to a reduction of the subsiding trend which amplifies northward. The adjustment we applied is empirical, based on vertical velocity estimates from the same GNSS network of stations (EPN) expressed with respect to these two different reference frames. As shown in Figure 3b, the fit of our transfer function is acceptable given that EGMS pixels are provided with a precision of 1 mm/yr.

Based on our adjusted estimates, we found that 50% of the European CFPs spatial extent is, on average, experiencing subsidence at a rate stronger than 1 mm/yr and 10% are uplifting at a larger rate than 0.5 mm/yr (Figure 4). If the GIA contribution is removed, it reduces the CFP area uplifting to 5% and the CFP area subsiding at rates lower than 1 mm/yr to slightly less than 20% (on average). This shows the substantial first-order effect of GIA in both land uplift and subsidence trends across Europe. However, the detailed results depend on the GIA model employed. Here we used the Caron et al. (2018) model, which has been successfully employed recently to close the basin-scale 20th century sea level budget and which appeared within the confidence intervals of the recent ICE-6G_C(VM5a) GIA model (Frederikse et al., 2020; Peltier et al., 2018). Nonetheless, Vermeersen et al. (2018) found that discrepancies between GIA models obtained from different methods and source data (e.g., data-driven model (Simon et al., 2018), regional GIA model (Bradley et al., 2011) or Ice-sheet history-based 3D-ice sheet simulations (Kuchar et al., 2012)) can significantly affect relative sea level regionally (see e.g. Figure 7 in Vermeersen et al. (2018) or Bradley et al. (2023)). Given the importance of GIA for VLM in Europe, uncertainties related to the GIA model employed need to be analyzed in more depth in future research work.

Comparing Figure 6 (based on mean estimates per CFP) and the Figure S5 in Supporting Information S1 (based on median estimates per CFP), we find consistent results in Northern and Southern European countries, but quite a substantial difference in the Netherlands, Germany and France. Indeed, for the later countries, the median estimates trends are less strong (negatively) than the mean estimates. The panel 1 of Figure 7, which details the Netherlands distribution, shows that within one very expanded CFP, the VLM spatial distribution can be very heterogeneous, suggesting various underlying physical processes (e.g., Erkens et al., 2016). While this heterogeneity of the VLM distribution within one CFP can explain why median and mean estimates are different, it also highlights that for more detailed analysis, these two summary statistics may not always provide enough information to describe the state of VLM in coastal plains, and a complete statistical distribution may be required.

Beside the use of EGMS for the European assessment (Section 4.1), we also demonstrate how to take advantage of the pixel scale estimate (Sections 4.2 and 4.3). We are able to identify a subsiding coastal area (i.e., Schiavonea) that was not previously recognized or considered in previous pan-European coastal impact assessments. The averaged results combined with the CLC data set allow drawing a significant overall behavior. Indeed, we find that urban areas and ∼16 millions residents are experiencing almost a 1 mm/yr subsidence on average (if we discard the uplifting regions due to GIA). For coastal airports and for harbors, the average land motion drops to a 1–1.5 mm/yr subsidence leading to a substantial relative SLR enhancement locally. This may have direct consequences on emerging risks and 21st century SLR adaptation planning and need to be considered to update previous coastal flooding risk assessments (e.g., Hauer et al., 2021).

We also note that while in this work we used EGMS information from a relative SLR perspective, further studies dedicated to local coastal impacts could take further advantage of the gridded VLM information. Indeed, spatial VLM could be used to adjust the elevation of DEMs and DTMs and ultimately to simulate flooding and erosion using coastal impact models (e.g., Miller & Shirzaei, 2021; Minderhoud et al., 2020; Ohenhen et al., 2024).

Broader Implications and Research Needs

As we move from understanding the magnitudes of SLR and the potential for climate mitigation to coastal adaptation, the importance of high quality relative SLR scenarios becomes increasingly apparent (e.g., Bednar-Friedl et al., 2022). Yet, the quality of coastal relative sea-level projections is limited by the lack of information on VLM (Slangen et al., 2023): this is acknowledged in the IPCC AR6, which assigns “low to medium confidence” to VLMs projections presented in the report (Fox-Kemper et al., 2021). In areas affected by subsidence, overlooking VLMs could significantly underestimate the risks and reduce the time available for adaptation as critical thresholds such as seawall heights or the capacity of flood expansion zones in estuaries will be exceeded earlier than projected.

Services such as the EGMS offer an obvious response to this need for improved information on present-day and future VLM. Importantly, such services do not only integrate information from the SAR Sentinel satellites, but also include information from the permanent GNSS networks to link VLMs tiles obtained from SAR interferometry. This, in turn, builds upon three decades of research in space geodesy (e.g., Bruyninx, 2004; Ferretti et al., 2001). For example, our study provides, to our knowledge, a first practical demonstration of the crucial role played by the terrestrial reference frame in which the InSAR results are expressed (e.g., Wöppelmann & Marcos, 2016). The type of analysis illustrated in Figure 3 is crucial to all users of broad scale land motion observations based on InSAR, whether in Europe or in other regions.

Our results show that when developing climate services to assess coastal future flood risks in Europe, estimates of uplift and subsidence were essential. These needs would be stronger in other regions where subsidence is a more important process, most especially east, south-east and south Asia (Nicholls, Lincke, et al., 2021). In these regions, research studies may not be sufficient and operational services similar to EGMS may be useful in order to provide stakeholders with validated and standardized products. We note that such services have emerged nationally before EGMS was opened as a Europe-wide service: the first Ground Motion Service (GMS) covered Italy (Costantini et al., 2017). This was followed by Norway in 2018, Germany in 2019 and then Denmark and the Netherlands. All other countries can consider similar services owing to the almost global coverage of Sentinel-1 constellation of satellite and its past and planned continuity of service.

It is important to note that the subsidence estimates provided here are based on a short period of observations - the five-year EGMS period of 2015–2021. As EGMS is updated on an annual basis, repeating this analysis should allow updating of the VLM estimates provided here and monitoring their continuity in time. This also raises the potential to actively monitor VLM at broad scales, such as Europe, and identify any new or enhanced subsidence hot-spots that might emerge due to changing land use and groundwater exploitation. SAR interferometry, on which EGMS relies, has limitations in vegetated areas due to the lack of backscatterers on the ground. From the perspective of the application presented in this paper, this is not a major limitation as relative SLR projections are most useful in areas with high exposure, that is, in cities and settlements. In such areas InSAR performs well.

Given the huge amount of EGMS data available at the European scale, we had to restrict our analysis to the mean velocity estimates provided at each pixel. For more detailed studies in space (e.g., at the city scale), considering the full time series could be relevant, for example, to test the stationarity of the trends and better quantify their robustness. The latter point also raises the important question of how to develop VLM projections to complement climate-induced sea-level projections: how can we best use EGMS estimates to infer future VLM and especially subsidence? This is a topic for further research.

Conclusions

In this paper we present a pan-European assessment of contemporary VLMs in the coastal flood prone areas based on EGMS InSAR-based products. Our estimates show that near half of the total CFP area is on average subsiding at a rate faster than 1 mm/yr. We find that urban zones and population in flood-prone areas, in European regions not subjected to GIA-induced land uplift (i.e., Fennoscandia), currently experience a subsidence of nearly 1 mm/yr. For harbors and airports located in flood-prone areas, the average coastal subsidence is more than 1 mm/yr. This locally enhances the relative SLR rate by about one third (in comparison with the ongoing climate-induced regional mean sea level in Europe which has risen by an average of 2–4 mm/year over the past 30 years; ) and demonstrates the significance of subsidence in Europe when considering coastal flood risk.

For coastal adaptation practitioners, the next question that arises is to understand to what extent observed VLMs can be extrapolated into the future. The question is important from a policy perspective as adaptation practitioners across the world are increasingly using and refining sea-level scenarios for risk assessment and adaptation planning (Hirschfeld et al., 2023). Another question is how to support decision making at the level of specific infrastructure, as the EGMS shows that some dikes and other coastal infrastructures are subject to subsidence, possibly in part due to their own weight. This will require further research and ensuring the continuity of observations in EGMS to confirm trends and identify potential non-linearities.

Notwithstanding these challenges, our study shows that there is a potential for a service for coastal adaptation practitioners, positioned downstream of EGMS or similar services in other regions worldwide. EGMS is a service of the “land” component of the European Earth Observation program Copernicus, but delivering actionable information to coastal users concerned with relative SLR requires additional analysis, as shown in this paper. Similar services can be developed in other regions of the world owing to the almost global coverage of the world's coasts by the Sentinel 1 constellation of satellites. Given the challenges raised by SLR, it would be important to make sure that such operational services are further tailored to the needs of coastal practitioners concerned with understanding risks and planning for adaptation.

Acknowledgments

This work is funded by CoCliCo and GSEU projects. CoCliCo project has received funding from the European Union's Horizon 2020 research and innovation programme under Grant 101003598. GSEU has received funding from the European Union's Horizon Europe research & innovation programme under Grant 101075609 (HORIZON-CL5-2021-D3-02). We thank the European Copernicus Land Monitoring Service and colleagues who produced and released the EGMS service and CLC. We thank the SONEL platform. We thank researchers Jonathan Chenal, involved in CoCliCo, for the fruitful discussions on geodesy and Aurélie Maspataud for her support as co-leader of the WP5 of the Horizon GSEU program. Alexandra Toimil acknowledges the financial support from the Ministerio de Ciencia e Innovación through the Ramon y Cajal Programme (RYC2021-030873-I with funding from MCIN/AEI and NextGenerationEU/PRTR).

Conflict of Interest

The authors declare no conflicts of interest relevant to this study.

Data Availability Statement

EGMS Ortho products first release (2015–2021) is obtained from European Environment Agency (2023). The ULR-repro3 GNSS products are available from Gravelle et al. (2022). The GNSS NGL14 solution (Blewitt et al., 2018) used in this study is distributed by SONEL at the landing page (). EUREF EPN D2150 densification GNSS solution was obtained from EUREF (). The results from the GIA model were downloaded from Frederikse, Caron, and Landerer (2019). CORINE Land Cover raster is obtained from CLC raster is obtained from European Environment Agency (2020). European population gridded data set from 2021 is obtained from Eurostat (2023). The Copernicus Global Digital Elevation Model at 30-m resolution used to define European coastal flood plains from European Space Agency, Sinergise (2021). Vertical land motion estimates in European coastal flood plains are provided in Thiéblemont (2024).

References

Altamimi, Z., Métivier, L., Rebischung, P., Rouby, H., & Collilieux, X. (2017). ITRF2014 plate motion model. Geophysical Journal International, 209(3), 1906–1912. [DOI: https://dx.doi.org/10.1093/gji/ggx136]