GOM20: A Stable Geodetic Reference Frame for Subsidence, Faulting, and Sea-Level Rise Studies along the Coast of the Gulf of Mexico

We have established a stable regional geodetic reference frame using long-history (13.5 years on average) observations from 55 continuously operated Global Navigation Satellite System (GNSS) stations adjacent to the Gulf of Mexico (GOM). The regional reference frame, designated as GOM20, is aligned in origin and scale with the International GNSS Reference Frame 2014 (IGS14). The primary product from this study is the seven-parameters for transforming the Earth-Centered-Earth-Fixed (ECEF) Cartesian coordinates from IGS14 to GOM20. The frame stability of GOM20 is approximately 0.3 mm/year in the horizontal directions and 0.5 mm/year in the vertical direction. The regional reference frame can be confidently used for the time window from the 1990s to 2030 without causing positional errors larger than the accuracy of 24-h static GNSS measurements. Applications of GOM20 in delineating rapid urban subsidence, coastal subsidence and faulting, and sea-level rise are demonstrated in this article. According to this study, subsidence faster than 2 cm/year is ongoing in several major cities in central Mexico, with the most rapid subsidence reaching to 27 cm/year in Mexico City; a large portion of the Texas and Louisiana coasts are subsiding at 3 to 6.5 mm/year; the average sea-level-rise rate (with respect to GOM20) along the Gulf coast is 2.6 mm/year with a 95% confidence interval of ±1 mm/year during the past five decades. GOM20 provides a consistent platform to integrate ground deformational observations from different remote sensing techniques (e.g., GPS, InSAR, LiDAR, UAV-Photogrammetry) and ground surveys (e.g., tide gauge, leveling surveying) into a unified geodetic reference frame and enables multidisciplinary and cross-disciplinary research.


Introduction
The Gulf of Mexico (GOM) is located at the south end of the North American tectonic plate and north of the Caribbean plate. This region is bounded on the northeast, north and northwest by the Gulf Coast of the United States, on the southwest and south by Mexico, and on the southeast by Cuba   [13]. The topographic and bathymetric datasets are from the GEBCO_2019 Grid (https://www.gebco.net).

Continuous GNSS Stations within the GOM Region
Due to the inherent benefits of obtaining high-precision daily-positions using GNSS, and thereby obtaining high-accuracy displacement time series and site velocities, a dense network of Continuously Operating Reference Stations (CORS) has been established in the GOM region by the geodesy research and land surveying communities. As shown in Figure 1, there are over 1000 long-history CORS in the GOM region that each has accumulated over 1000 observational files (one file per day) spanning over at least three years as of the end of 2019. A detailed analysis of the observational history of these GNSS stations is depicted in Figure 2. Over 500 CORS have a history of over eight years and have more than 2400 daily observational files. These GNSS stations were installed and operated by the following research and surveying communities: the National Geodetic Survey (NGS) (http://geodesy.noaa.gov/CORS), The number of permanent GNSS stations is continuously increasing in this region. These longterm permanent GNSS stations have provided fundamental datasets for delineating ground deformation associated with faulting and subsidence within the GOM region. The majority of GNSS antennas utilized in this study are mounted on building-based monuments. Recent investigations have confirmed that long-term (e.g., >3 years) building-based GNSS monuments could provide reliable site velocity estimates for ground deformation studies provided that the buildings are stable [14,15].
This study applied the precise point positioning (PPP) method employed by the GipsyX (Version 1.2) software package (https://gipsy-oasis.jpl.nasa.gov) for calculating single-receive phaseambiguity-fixed daily solutions with respect to IGS14. The methodology and algorithms used for the PPP solutions are described in Zumberge et al. [16] and Bertiger et al. [17]. The Nevada Geodetic Laboratory at the University of Nevada, Reno, provides daily PPP positions from over 17,000 GNSS stations to the public [18]. The Geodetic Facility for the Advancement of Geoscience (GAGE), operated by UNAVCO, also provides daily PPP solutions of numerous GNSS stations to the public [19]. The PPP method has attracted broad interests in ground deformation and structural health monitoring [20][21][22][23].

Realization of GOM20
GNSS positions are initially provided as a set of coordinates with respect to a global geodetic reference frame, such as IGS14, which is based on the International Terrestrial Reference Frame 2014 (ITRF14) [24]. In general, a global geodetic reference frame is a no-net-rotation (NNR) reference frame that is realized with an approach of minimizing the overall horizontal movements of a group of selected reference stations distributed worldwide [12]. As a result, GNSS-derived site movements with respect to a global reference frame are often dominated by factors such as the long-term drift and rotation of the tectonic plate on which the site is located and other secular movements [25]. In general, site velocities with respect to a global reference frame are difficult to interpret visually from a regional or local-scale geophysical perspective. Localized and/or temporal ground deformation, such as subsidence and fault creeping, could be obscured or biased by those common motions. A regional-or local-scale reference frame that fixed on the stable portion of a tectonic plate is needed when local-scale ground deformation is of interest to researchers.
There has been a conscious need for a consistent and stable reference frame in the GOM region since the early 2000s due to the rapid increase of publicly available GNSS data for the research community as well as the broad extent of land subsidence within the GOM region. During the past

Realization of GOM20
GNSS positions are initially provided as a set of coordinates with respect to a global geodetic reference frame, such as IGS14, which is based on the International Terrestrial Reference Frame 2014 (ITRF14) [24]. In general, a global geodetic reference frame is a no-net-rotation (NNR) reference frame that is realized with an approach of minimizing the overall horizontal movements of a group of selected reference stations distributed worldwide [12]. As a result, GNSS-derived site movements with respect to a global reference frame are often dominated by factors such as the long-term drift and rotation of the tectonic plate on which the site is located and other secular movements [25]. In general, site velocities with respect to a global reference frame are difficult to interpret visually from a regional or local-scale geophysical perspective. Localized and/or temporal ground deformation, such as subsidence and fault creeping, could be obscured or biased by those common motions. A regional-or local-scale reference frame that fixed on the stable portion of a tectonic plate is needed when local-scale ground deformation is of interest to researchers.
There has been a conscious need for a consistent and stable reference frame in the GOM region since the early 2000s due to the rapid increase of publicly available GNSS data for the research community Remote Sens. 2020, 12, 350 5 of 29 as well as the broad extent of land subsidence within the GOM region. During the past two decades, Interferometric Synthetic Aperture Radar (InSAR), light detection and range (LiDAR), and Unmanned Aerial Vehicle (UAV) photogrammetry techniques have also been frequently employed for mapping land subsidence, faulting, and coastal erosion within the GOM region [8,[26][27][28][29][30]. Applications of remote sensing techniques further prompted the need for a unified reference frame to align measurements from different remote sensing systems.

Reference Stations
The accuracy of GNSS-derived long-term site velocities does not solely rely on equipment (antennas and receivers), but largely depends on the available regional geodetic infrastructure, which can be defined by two fundamental components: a dense continuous long-history GPS network and a stable regional reference frame [11]. The former is referred to as the hardware component, and the latter is referred to as the firmware component of a regional geodetic infrastructure. In the surveying and geodesy community, a regional reference frame is often developed through a simultaneous Helmert transformation from a well-established and broadly used global reference frame. A group of common points (reference stations), are used to link these two reference frames. The key principle for selecting a group of reference stations is that these reference stations are not expected to move relative to each other.
The selection of reference stations is critically crucial for realizing a stable reference frame. A reference station that is not locally stable will degrade the accuracy (stability) of the reference frame. Unfortunately, there are no fixed criteria for selecting reference stations. In practice, the selection is largely based on the availability of long-term GNSS stations in the study area. A stable reference frame is a secular reference frame that applies a linear motion model predicting station positions within a temporal window of interest. Therefore, the linearity of the displacement time series, which indicates how well their daily positions can be described by a constant velocity, is the key criterion for selecting reference stations. Additionally, geographic distribution, time span, data continuity, and site stability are considered for selecting reference stations. The primary criteria for selecting reference stations for realizing GOM20 include: (1) distance within 500 km to the GOM coastline; (2) having over 2400 daily files and spanning over at least eight years; (3) active through 2019 with a continuous segment of minimum 5 years; (4) outside of known subsiding or fault-creeping areas.
The vertical displacement time series with respect to IGS14 provides a quick clue to exclude those potential subsiding stations. According to our previous investigations [31,32], within the southern portion of the North American plate, a GNSS site that is not affected by localized subsidence retains a vertical site velocity well below 3 mm/year with respect to IGS global reference frames (e.g., IGS08, IGS14). In other words, if a GNSS station depicts a steady downward movement larger than 3 mm/year with respect to IGS14, this site can be confidently described as experiencing localized subsidence. Applying the four aforementioned selection criteria resulted in approximately 200 remaining candidate frame stations. Further manual screening weeded out about 100 stations with visible issues, such as frequent data gaps, significant positional shifts, and unusually large seasonal variations.
A seven-parameter transformation method was employed to realize an interim reference frame with these 100 candidate reference stations. The process for realizing a regional reference frame will be addressed in the next section. The three-component velocities of these reference stations with respect to the interim reference frame were calculated. Ideally, a stable site retains near-zero (e.g., sub-millimeter per year) three-component site velocities with respect to the stable reference frame. Reference stations with a horizontal velocity larger than 2 mm/year and a vertical velocity larger than 3 mm/year with respect to the interim reference frame were excluded as outliers. The first iteration removed 24 stations. The procedure was executed several times, with each cycle removing a few outlier stations with iteratively tighter velocity thresholds. The final threshold was set to 1 mm/year in both horizontal and vertical directions. Finally, by careful consideration of the overall geographic coverage of the remaining candidate stations, fifty-five stations were selected as reference stations for realizing GOM20. The locations of these reference stations and their site velocities with respect to IGS14, North American Datum 1983 (NAD83), and GOM20 are depicted in Figure 3.
Remote Sens. 2020, 12, 350 6 of 29 realizing GOM20. The locations of these reference stations and their site velocities with respect to IGS14, North American Datum 1983 (NAD83), and GOM20 are depicted in Figure 3.  Table A1.
We also conducted several tests by adding to or removing a few stations from the group of 55 reference stations. We found that the changes provided little to no impact on the final displacement time series of coastal GNSS stations. In theory, three common stations will be able to provide the seven parameters for reference frame transformation [31]. However, more reference stations improve the reliability of the reference frame transformations and provide higher confidence in assessing the stability of the reference frame [33]. The definition of frame stability will be addressed in the next session. To take advantage of the dense and long-history CORS infrastructure that is available within the GOM region, we use 55 reference stations to realize the regional reference frame. Note that the redundancy of reference stations does not add any burden to GOM20 users. Users do not need any reference stations in obtaining GOM20 coordinates, and they even do not need to worry about the future status of these reference stations.
The detailed information of these 55 stations is listed in Appendix Table A1. The vast majority of reference stations span over 15 years. The average observational history is 13.5 years, with about 330 daily files (11 months) per year. The station with the shortest observational history is MSEV, with a history of 7.8 years as the end of 2019. MSEV is located in Ellisville, Mississippi. There are not so many qualified reference stations in the state of Mississippi compared to other nearby states. The three-component displacement time series of MSEV retain high linearity with little gaps. Therefore, MSEV is selected as a reference station for realizing GOM20.

Coordinate Transformation from IGS14 to GOM20
In the surveying and geodesy research community, the Helmert transformation method is used to transform Cartesian coordinates from a well-established global reference frame to a regional or local reference frame. There are two approaches for transforming the positional time series between  Table A1. We also conducted several tests by adding to or removing a few stations from the group of 55 reference stations. We found that the changes provided little to no impact on the final displacement time series of coastal GNSS stations. In theory, three common stations will be able to provide the seven parameters for reference frame transformation [31]. However, more reference stations improve the reliability of the reference frame transformations and provide higher confidence in assessing the stability of the reference frame [33]. The definition of frame stability will be addressed in the next session. To take advantage of the dense and long-history CORS infrastructure that is available within the GOM region, we use 55 reference stations to realize the regional reference frame. Note that the redundancy of reference stations does not add any burden to GOM20 users. Users do not need any reference stations in obtaining GOM20 coordinates, and they even do not need to worry about the future status of these reference stations.
The detailed information of these 55 stations is listed in Appendix A Table A1. The vast majority of reference stations span over 15 years. The average observational history is 13.5 years, with about 330 daily files (11 months) per year. The station with the shortest observational history is MSEV, with a history of 7.8 years as the end of 2019. MSEV is located in Ellisville, Mississippi. There are not so many qualified reference stations in the state of Mississippi compared to other nearby states. The three-component displacement time series of MSEV retain high linearity with little gaps. Therefore, MSEV is selected as a reference station for realizing GOM20.

Coordinate Transformation from IGS14 to GOM20
In the surveying and geodesy research community, the Helmert transformation method is used to transform Cartesian coordinates from a well-established global reference frame to a regional or local reference frame. There are two approaches for transforming the positional time series between two frames: a daily 7-parameter method and a total 7-parameter method [31]. The daily 7-parameter method is often used for transforming GNSS-derived positional coordinates between two reference frames day by day. For examples, UNAVCO employs the daily 7-parameter transformation method in producing daily solutions for the Plate Boundary Observatory (PBO) network [19]; the Nevada Geodetic Laboratory at the University of Nevada, Reno employs the daily 7-parameter transformation method in producing daily solutions with respect to the North American Reference Frame (NA12) [34]. However, calculating daily transformation parameters can be a burden for most practical users. During the past years, we have developed a total 7-parameter method for transforming Earth-Centered-Earth-Fixed (ECEF) Cartesian (XYZ) coordinate time series from a global reference frame to a regional reference frame. The details of the method are introduced in our recent publications for establishing regional reference frames in North China [35], Houston metropolitan region [36], and the Caribbean [33].
GOM20 is aligned in origin and scale with IGS14. The PPP solutions with respect to IGS14 are presented in an ECEF-XYZ coordinate system. The ECEF-XYZ coordinates with respect to GOM20 can be obtained by the following formula, Equation (1): IGS14 , and Z(t) IGS14 are the ECEF-XYZ coordinates (at epoch t) of a site with respect to IGS14; X(t) GOM20 , Y(t) GOM20 , and Z(t) GOM20 are the ECEF-XYZ coordinates of the site with respect to GOM20. The units of these XYZ coordinates are meters. t 0 is the epoch that aligns the coordinates with respect to these two reference frames, which is fixed at epoch 2015.0 (year) for GOM20. A site retains identical XYZ coordinates with respect to IGS14 and GOM20 at epoch 2015.0. T x , T y , T z , R x , R y , and R z are constant parameters indicating the rates (one-time derivates) of translational shifts and rotations along the x, y, an z coordinate axes between two reference frames. These seven parameters: t 0 , T x , T y , T z , R x , R y , and R z , are listed in Table 1. To study ground deformation at the Earth's surface, firstly, the ECEF-XYZ coordinates are transformed to GOM20 from IGS14; secondly, the geocentric XYZ coordinates are converted to a geodetic orthogonal curvilinear coordinate system (longitude, latitude, and ellipsoid height) referencing to the GRS80 ellipsoid; thirdly, the longitude and latitude coordinates are projected to a two-dimensional local horizontal plane for tracking superficial ground deformation in the north-south (NS) and east-west (EW) directions at each site. The change of ellipsoid heights over time is used to measure the up-down Remote Sens. 2020, 12, 350 8 of 29 (UD) movement (subsidence or uplift). In practice, the vertical displacements (subsidence or uplift) derived from the ellipsoid heights and orthometric heights are the same [37].

Stability of GOM20
There is no widely accepted definition for the stability or accuracy of reference frames. These two terms are often applied interchangeably. The basic rule is that a stable site should retain near-zero velocities over time in all three directions (NS, EW, and UD) with respect to the stable reference frame. The term accuracy emphasizes the difference between the real site velocities and the supposed zero velocities at stable sites, while the stability emphasizes the accumulation of the positional shift over time. The stability of a static reference frames degrades over time since a linear movement (over time) is assumed in realizing the stable frame transformation. Eventually, the positional error becomes too imprecise with respect to user expectations. The stability of a reference frame determines its ability to extrapolate station coordinates accurately into the past and the future beyond the frame range [34]. In practice, the stability or accuracy of a regional reference frame is often evaluated by the root-mean-square (RMS) of the velocities of all reference stations with respect to the reference frame [33]. In essence, the lifetime of a static reference frame depends on the rigidity of the crustal block that the network of reference stations covers. The reference station network of GOM20 indicates a relative horizontal RMS accuracy of 0.3 mm/year and a vertical RMS accuracy of 0.5 mm/year (see Appendix A Table A1), suggesting that the GOM region is rigid at a half-millimeter per year level. Nevertheless, users should exercise appropriate cautions in interpolating site velocities at a sub-millimeter per year level for long-term ground deformation studies.
According to our previous investigations, the RMS-accuracy (repeatability) of the PPP daily positions is approximately 2 to 4 mm in the horizontal directions and 6 to 8 mm in the vertical direction within the GOM region [32,38]. As mentioned earlier, the stability of GOM20 is at a level of 0.3 mm/year in the horizontal directions and 0.5 mm/year in the vertical direction. Thus, GOM20 could result in accumulated positional-errors of 4.5 mm in the horizontal directions and 7.5 mm in the vertical direction across a 15-year window, which are still comparable with the RMS-accuracy of the daily GNSS measurements (PPP solutions). GOM20 is aligned to IGS14 at 2015.0. Hence, the potential positional errors 15 years before and after epoch 2015.0 (2000-2030) are still within the accuracy of the daily GNSS measurements. Though GPS datasets as early as 1990 are available within the GOM region, only a few stations have data before 1995, and over 98% GNSS stations were installed after 2000 (see Figure 2). In general, the daily positions derived from GPS data in the 1990s retain a poorer positional accuracy than the daily positions since 2000 because of fewer available satellites and reference frame stations available in the 1990s [38,39]. Thus, the accumulated positional errors associated with GOM20 may still lie within the accuracy of daily GPS measurements, even for the data before 2000. Accordingly, GOM20 can be confidently used for a time window from the 1990s to 2030 without causing positional errors larger than the accuracy of daily GNSS measurements.
In order to verify the stability of GOM20, we checked the positional time series of two long-history LKHU is located at Lake Houston, a northern suburb of Houston. The antenna pole of LKHU is also anchored at the bottom of the Evangeline aquifer, approximately 591 m below the local land surface. Detailed information about the two GNSS stations is presented in our previous publications [40,41]. The Evangeline aquifer is the major aquifer providing groundwater for industry and public use in the greater Houston area [42,43]. According to our previous investigations, present sediment compactions at the two sites occurred within the top portion of the Evangeline aquifer and its overlay aquifer, Chicot aquifer [41]. There are no reported ongoing faulting activities adjacent to the two locations. Thus, the two sites should be stable in all three directions with respect to GOM20. The plots in Figure 4 indicate that there are no steady displacement trends at both sites over 27 years. That means both sites are stable over time with respect to GOM20. The 27-year displacement time series verify that GOM20 provides a reliable reference in assessing the long-term site stability within the GOM region.

Applications of GOM20
GOM20 provides a coherent reference frame for aligning ground deformation within the GOM region over time and space. Figure 5 depicts horizontal and vertical velocity vectors (with respect to GOM20) derived from GNSS observations within the GOM region. To minimize the temporal variations of vertical site velocities, we only depict the velocity vectors at stations that span at least three years and have at least 1000 observational files within the period from 2010 to 2019. The velocities are calculated with the Median Interannual Difference Adjusted for Skewness (MIDAS) method [44]. MIDAS breaks displacement time series into 1-year-long (365 days) segments and selects the median slope of the selected segments as the slope (trend) of the whole time series. We compared the site velocities calculated by MIDAS and the conventional least-squares method. It is found that the site velocities (trends) obtained by the two methods agree well in general and the MIDAS method does a better job in minimizing the effects of outliers, step discontinuities, and seasonal movements superimposed into GNSS time series.
MIDAS computes a velocity uncertainty that is based on the distribution of sampled slopes of 1year-long displacement time series. For the 500 stations shown in Figure 5, the velocity uncertainty varies from ±0.2 to ±0.6 mm/year in the horizontal directions and ±0.3 to ±1.0 mm/year in the vertical direction. The MIDAS velocity uncertainty approximates to one standard deviation of the selected slopes, which indicates a 68% confidence interval (CI) of the velocity estimate. The uncertainty is inversely correlated with the time span of GNSS datasets. A longer observational history and fewer data gaps often result in a smaller uncertainty of the velocity estimate. In general, the uncertainty is resulted by noises superimposed into the GNSS-derived displacement time series. The noise sources are dominated by seasonal ground deformation due to changes in groundwater levels and atmospheric loading, monument instability (random walk noise), and errors in satellite orbits and GNSS equipment (receivers, antennas). The noises can be best described by a combination of white noise and flicker noise [45].
The velocity vectors illustrated in Figure 5 indicate that the GOM region is overall stable except

Applications of GOM20
GOM20 provides a coherent reference frame for aligning ground deformation within the GOM region over time and space. Figure 5 depicts horizontal and vertical velocity vectors (with respect to GOM20) derived from GNSS observations within the GOM region. To minimize the temporal variations of vertical site velocities, we only depict the velocity vectors at stations that span at least three years and have at least 1000 observational files within the period from 2010 to 2019. The velocities are calculated with the Median Interannual Difference Adjusted for Skewness (MIDAS) method [44]. MIDAS breaks displacement time series into 1-year-long (365 days) segments and selects the median slope of the selected segments as the slope (trend) of the whole time series. We compared the site velocities calculated by MIDAS and the conventional least-squares method. It is found that the site velocities (trends) obtained by the two methods agree well in general and the MIDAS method does a better job in minimizing the effects of outliers, step discontinuities, and seasonal movements superimposed into GNSS time series.
MIDAS computes a velocity uncertainty that is based on the distribution of sampled slopes of 1-year-long displacement time series. For the 500 stations shown in Figure 5, the velocity uncertainty varies from ±0.2 to ±0.6 mm/year in the horizontal directions and ±0.3 to ±1.0 mm/year in the vertical direction. The MIDAS velocity uncertainty approximates to one standard deviation of the selected slopes, which indicates a 68% confidence interval (CI) of the velocity estimate. The uncertainty is inversely correlated with the time span of GNSS datasets. A longer observational history and fewer data gaps often result in a smaller uncertainty of the velocity estimate. In general, the uncertainty is resulted by noises superimposed into the GNSS-derived displacement time series. The noise sources are dominated by seasonal ground deformation due to changes in groundwater levels and atmospheric loading, monument instability (random walk noise), and errors in satellite orbits and GNSS equipment (receivers, antennas). The noises can be best described by a combination of white noise and flicker noise [45].
The velocity vectors illustrated in Figure 5 indicate that the GOM region is overall stable except the southwest region, which is affected by the subducting of the Rivera and Cocos Plates beneath the North American Plate along the Middle America subduction zone ( Figure 1). Rivera and Cocos Plates are relatively young (10-25 Ma) oceanic plates; both plates are subducting along the Middle American Trench at different convergence rates [46]. The GNSS-derived velocity vectors along the western coast of Central Mexico vividly indicate that the convergence rate increases progressively to the southeast along the Middle American Trench. Driven by the subduction of the Cocos Plate, the western coast of Central Mexico retains steady movements towards the northeast (1 to 3 cm/year) and downward

Urban Subsidence
Among these 1000 continuous GNSS stations plotted in Figure 1, seven stations recorded subsidence rates over 2 cm/year, which are located in major cities of central Mexico ( Figure 6). Over twenty stations in the greater Houston area recorded subsidence rates between 1 cm/year and 2 cm/year. GNSS-derived land subsidence and faulting activities in the greater Houston area have been investigated by Harris-Galveston Subsidence District [47], US Geological Survey [43], and the University of Houston [39,[48][49]. Therefore, this study will not explore the details of subsidence and faulting activities within the greater Houston area.

Urban Subsidence
Among these 1000 continuous GNSS stations plotted in Figure 1, seven stations recorded subsidence rates over 2 cm/year, which are located in major cities of central Mexico ( Figure 6). Over twenty stations in the greater Houston area recorded subsidence rates between 1 cm/year and 2 cm/year. GNSS-derived land subsidence and faulting activities in the greater Houston area have been investigated by Harris-Galveston Subsidence District [47], US Geological Survey [43], and the University of Houston [39,48,49]. Therefore, this study will not explore the details of subsidence and faulting activities within the greater Houston area.
These seven GNSS stations experienced rapid subsidence (>2 cm/year) are MMX1 (subsidence rate: 27.1 cm/year), INEG (subsidence rate: 4.0 cm/year), UAGU (subsidence rate: 7.0 cm/year), CECM (subsidence rate: 6.3 cm/year), UNTO (subsidence rate: 5.5 cm/year), UIRA (subsidence rate: 6.7 cm/year), and TNZA (subsidence rate: 7.5 cm/year). They are located in urban areas that rely mainly on groundwater to sustain agriculture, industrial, and domestic water needs [50]. The GPS-derived subsidence time series of these stations are depicted in  UAGU and INEG are two GNSS stations located at the northern and southern Aguascalientes, the capital of the Mexico state Aguascalientes, 430 km northwest of Mexico City. Aguascalientes is situated over a vast unconfined aquifer system encompassing three Mexico states. The growing population and increased agricultural and industrial activities have resulted in intensive groundwater extraction since the 1970s [51]. The distance between UAGU and INEG is 7.6 km ( Figure 6). The subsidence rate at INEG was 8.  (Figure 7b). The subsidence time series at UAGU and INEG depict extraordinary temporal and spatial variations of urban subsidence, which have also been observed in many other sinking cities [2,32,39]. In general, rapid subsidence in urban areas is primarily caused by groundwater pumping; the rate of subsidence varies over time and space in response to the variations of pumping amount, groundwater levels, and geologic and hydrologic conditions of local aquifers.  UAGU and INEG are two GNSS stations located at the northern and southern Aguascalientes, the capital of the Mexico state Aguascalientes, 430 km northwest of Mexico City. Aguascalientes is situated over a vast unconfined aquifer system encompassing three Mexico states. The growing population and increased agricultural and industrial activities have resulted in intensive groundwater extraction since the 1970s [51]. The distance between UAGU and INEG is 7.6 km (Figure Land subsidence in major cities in central Mexico has been frequently investigated using InSAR techniques [50,52,53]. Since InSAR employs the difference in the carrier signal phase between Synthetic Aperture Radar (SAR) images to recover relative displacements over space, there is a benefit to tying InSAR-derived displacements to absolute positions with respect to GOM20, which is aligned to the global reference frame IGS14 at epoch 2015.0. The GNSS-derived positional time series with respect to GOM20 will provide ground-truth to calibrate and constrain InSAR-derived subsidence measurements over time and space [54].

Coastal Subsidence and Faulting
Subsidence has been widely recognized as a driver of geomorphic change in the northern Gulf of Mexico [2,55,56]. However, there is considerable disagreement over the subsidence rates and the interpreted temporal and spatial variability in these rates. The vertical velocity vectors depicted in Figure 5 indicate that vertical land deformation rates along the Florida portion of the GOM coast are within ±1 mm/year with respect to GOM20.   Figure 9 illustrates the subsidence time series of these seven coastal GNSS stations on the Mississippi River delta. Overall, these seven sites retain steady subsidence over time, which are remarkably different from urban subsidence dominated by groundwater withdrawal, as shown in Figure 7b. Modern subsidence in the Louisiana coastal area has been described as a result of (1) fluid withdrawal (groundwater, oil and gas) (up to a few centimeters per year) [2], (2) deep-seated tectonic movements (fault processes, up to a few millimeters per year) [6], (3) Holocene sediment compaction (up to a few millimeters per year) [58,59], (4) sediment loading (below 0.5 mm/year) [60][61][62], (5) glacial isostatic adjustment (below 1 mm/year) [63], and (6) surface water drainage and management [64]. It is worth pointing out that these processes are not fully independent from one another.
Dakka et al. [6] suggested that the steady subsidence in the Mississippi River delta area is likely due to deep-seated faulting. They proposed that the Mississippi delta area is located on the hanging wall of a large listric normal fault system.   There is a clear trend that the subsidence rate generally decreases westward, eastward, and landward with respect to the Mississippi River delta. From the inland to the front of the delta, subsidence rate increases from 2 to 4 mm/year (MGW3: 2.0 ± 1.2 mm/year; AME4: 2.8 ± 1.1 mm/year; LAGM: 3.1 ± 0.8 mm/year; HOUM: 3.8 ± 0.6 mm/year; BVHS: 3.9 ± 0.6 mm/year) to 6.5 mm/year (LMCN: 6.3 ± 0.6 mm/year; GRIS: 6.5 ± 0.5 mm/year). Figure 9 illustrates the subsidence time series of these seven coastal GNSS stations on the Mississippi River delta. Overall, these seven sites retain steady subsidence over time, which are remarkably different from urban subsidence dominated by groundwater withdrawal, as shown in Figure 7b. Modern subsidence in the Louisiana coastal area has been described as a result of (1) fluid withdrawal (groundwater, oil and gas) (up to a few centimeters per year) [2], (2) deep-seated tectonic movements (fault processes, up to a few millimeters per year) [6], (3) Holocene sediment compaction (up to a few millimeters per year) [58,59], (4) sediment loading (below 0.5 mm/year) [60][61][62], (5) glacial isostatic adjustment (below 1 mm/year) [63], and (6) surface water drainage and management [64]. It is worth pointing out that these processes are not fully independent from one another.  Dakka et al. [6] suggested that the steady subsidence in the Mississippi River delta area is likely due to deep-seated faulting. They proposed that the Mississippi delta area is located on the hanging wall of a large listric normal fault system. The distance between these two stations is approximately 2 km. Based on our experience with urban faults in Houston [9], the highly coherent and steady horizontal movements observed at the two sites are likely associated with local faulting activities.

Delineating the Absolute Sea-level Rise Rate along the GOM Coast
Global sea level has been rising over the past century, and the rising rate has increased in recent decades due to global warming [65]. Sea level plays a critical role in wetland loss, coastal flooding, shoreline erosion, and flooding hazards from storms. In the long term, rising sea level also creates stress on groundwater systems and coastal ecosystems. As the sea-level rises, saltwater intrudes into freshwater aquifers, many of which sustain coastal ecosystems and provide water for municipal and agriculture. The Louisiana and Texas coastal areas are currently experiencing subsidence varying from 4 to 6.5 mm/year, as depicted in Figure 8.
Sea level is primarily measured using tide gauge stations along the coasts and satellite laser altimeters in the open ocean. The Center for Operational Oceanographic Products and Services (COOPS) at the National Ocean Service, National Oceanic and Atmospheric Administration (NOAA) of the US operates 25 long-history (>30 years) tide gauge stations along the GOM coast (Figure 11), which are a part of the National Water Level Observation Network (NWLON). The network operates 210 long-term, continuously operating water level stations throughout the US and its territories. NOAA/COOPS also provides routinely-updated analyses of the long-term trends and variability from tide gauge stations operated by international agencies. There are four long-history (>30 years) tide gauge stations on the Mexico portion of the GOM coast and one tide gauge station on the coast of Cabo San Antonio, Cuba. In total, there are 30 long-history tide gauge stations along the GOM coast ( Figure 11). Long-term tide gauge records used for this study are downloaded from the NOAA COOPS (https://tidesandcurrents.noaa.gov/sltrends/sltrends.html). These measurements are monthly averages, which removes the effect of higher frequency phenomena to compute an accurate linear sea-level trend.
which are a part of the National Water Level Observation Network (NWLON). The network operates 210 long-term, continuously operating water level stations throughout the US and its territories. NOAA/COOPS also provides routinely-updated analyses of the long-term trends and variability from tide gauge stations operated by international agencies. There are four long-history (> 30 years) tide gauge stations on the Mexico portion of the GOM coast and one tide gauge station on the coast of Cabo San Antonio, Cuba. In total, there are 30 long-history tide gauge stations along the GOM coast ( Figure 11). Long-term tide gauge records used for this study are downloaded from the NOAA COOPS (https://tidesandcurrents.noaa.gov/sltrends/sltrends.html). These measurements are monthly averages, which removes the effect of higher frequency phenomena to compute an accurate linear sea-level trend. Figure 11. Map depicting the relative sea-level (RSL) trends and vertical land movement trends derived from closely-spaced tide gauge and GNSS pairs. The RSL trends are referred to adjacent benchmarks fixed on land surface; the land movement trends are referred to GOM20. Figure 11. Map depicting the relative sea-level (RSL) trends and vertical land movement trends derived from closely-spaced tide gauge and GNSS pairs. The RSL trends are referred to adjacent benchmarks fixed on land surface; the land movement trends are referred to GOM20.
Instead of measuring absolute (eustatic) sea-level changes, long-term tide gauge measurements provide information on relative sea-level changes. This is because tide gauges measure sea levels relative to reference benchmarks fixed on the local land surface through repeat leveling surveys. Therefore, the readings of tide gauges are a combination of the absolute sea-level changes and the local land vertical movements. Absolute sea-level rise is mostly due to a combination of thermal expansion of seawater and meltwater from glaciers and ice sheets. A smaller contributor to global sea-level rise could be the decline of aquifer storage capacity as a result of excessive groundwater pumping. The pumping-compaction process squeezes water out from aquifers and leads to permanent compaction of aquifers, eventually, moving water from inland to the sea. The absolute global sea-level rise rate is critically important for evaluating the long-term risk of coastal flooding, erosion, and wetland loss at a regional scale. The stable regional reference frame and the long-history GNSS observations in the GOM region provide a fundamental infrastructure for delineating secular sea-level trends. In this section, we use closely-spaced GNSS and tide gauge measurements on Grand Isle, Louisiana (Figure 12), to describe the process of delineating absolute sea-level trends. permanent compaction of aquifers, eventually, moving water from inland to the sea. The absolute global sea-level rise rate is critically important for evaluating the long-term risk of coastal flooding, erosion, and wetland loss at a regional scale. The stable regional reference frame and the long-history GNSS observations in the GOM region provide a fundamental infrastructure for delineating secular sea-level trends. In this section, we use closely-spaced GNSS and tide gauge measurements on Grand Isle, Louisiana (Figure 12), to describe the process of delineating absolute sea-level trends. Grand Isle is a barrier island at the southern edge of Barataria Bay, one of the large interdistributary bays in the Mississippi River delta. The Grand Isle tide gauge sits atop of the barrier island, closely-spaced (250 m) with a long-term GNSS station GRIS (Figures 12a, b). The Grand Isle tide gauge was installed in 1947 and has been continuously operated for 73 years. The primary benchmark used for leveling the tide gauge measurements is a disk set on the top of a concrete sea wall adjacent to the tide gauge (https://tidesandcurrents.noaa.gov/benchmarks/8761724.html). There are a dozen of other benchmarks around the primary benchmark for assessing the long-term vertical movements of the primary benchmark. Thus, sea-level changes over time are tracked relative to a fixed station datum maintained by the benchmark network. The distance between the tide gauge and the primary benchmark is approximately 220 m; the distance between the tide gauge and the GNSS station is approximately 250 m (Figure 12c Grand Isle is a barrier island at the southern edge of Barataria Bay, one of the large inter-distributary bays in the Mississippi River delta. The Grand Isle tide gauge sits atop of the barrier island, closely-spaced (250 m) with a long-term GNSS station GRIS (Figure 12a,b). The Grand Isle tide gauge was installed in 1947 and has been continuously operated for 73 years. The primary benchmark used for leveling the tide gauge measurements is a disk set on the top of a concrete sea wall adjacent to the tide gauge (https://tidesandcurrents.noaa.gov/benchmarks/8761724.html). There are a dozen of other benchmarks around the primary benchmark for assessing the long-term vertical movements of the primary benchmark. Thus, sea-level changes over time are tracked relative to a fixed station datum maintained by the benchmark network. The distance between the tide gauge and the primary benchmark is approximately 220 m; the distance between the tide gauge and the GNSS station is approximately 250 m (Figure 12c). GRIS is operated by the Center for Geoinformatics at Louisiana State University and has accumulated 15 years (2005-2019) of continuous GNSS observations ( Figure 13). The antenna of GRIS is mounted on a two-floor concrete-brick building, 60 meters away from the primary benchmark on the concrete sea wall. Since the GNSS antenna is only 250 m away from the tide gauge site, the 15-year GNSS observations would be able to provide accurate vertical land movement trend at the tide gauge site.
Recent studies conducted by the research group at the Tulane University suggest that a significant portion of subsidence in the Mississippi River delta area occurs within the shallow portion (top 5 m) of the sediments [66][67][68]. The foundations of buildings in the soft delta area are often deeper than similar structures inland. The monuments of the tide gauges could be considerably deep in order to sustain the long-term stability of the gauges. The deeply buried foundations may exclude certain compaction of shallow sediments [69]. If the depths of the foundation of the building that GNSS is mounted on and the foundation of the tide gauge are considerably different, the GNSS antenna and tide gauge may experience different subsidence, and both could be different from the subsidence at the land surface (total compaction). Thus, the absolute sea-level changes derived from the method illustrated in Figure 13 could be biased. A further investigation of this issue will be conducted in our future study.
to sustain the long-term stability of the gauges. The deeply buried foundations may exclude certain compaction of shallow sediments [69]. If the depths of the foundation of the building that GNSS is mounted on and the foundation of the tide gauge are considerably different, the GNSS antenna and tide gauge may experience different subsidence, and both could be different from the subsidence at the land surface (total compaction). Thus, the absolute sea-level changes derived from the method illustrated in Figure 13 could be biased. A further investigation of this issue will be conducted in our future study.   The relative sea-level rise rate is 9.1 mm/year with a 95% CI of ±0.4 mm/year. That means the site has submerged to water approximately 65 cm during the past 73 years. The process for calculating the 95% CI is documented in the NOAA Technical Report [70]. Figure 13b depicts the land subsidence time-series recorded by the GNSS station GRIS and the relative sea-level rise recorded by the tide gauge during 15 years from 2005 to 2019. The land subsidence rate is 6.5 mm/year with a 95% CI of ±1 mm/year with respect to GOM20. The 95% CI is estimated by twice the velocity uncertainty obtained according to the MIDAS method [44]. Assuming the local land retained the same subsidence rate during the history of the tide gauge, the absolute sea-level trend at the tide gauge site can be obtained by subtracting out the land subsidence rate (6.5 mm/year) from the relative sea-level rise rate (9.1 mm/year). The reasonableness of the assumption will be discussed later. Figure 13c depicts the absolute sea-level rise time series with respect to GOM20. The absolute sea-level rise rate at the tide gauge site is 2.6 mm/year (with respect to GOM20). The uncertainty (95% CI) of the corrected sea-level rise rate could be approximated by the larger one between the uncertainty of the sea-level rise rate (± 0.4 mm/year) and the uncertainty of the land subsidence rate (± 1 mm/year).
The relative sea-level rise rates at 30 tide gauge stations along the GOM coast are depicted in Figure 11. The rate ranges from 2.1 mm/year (106 years: 1914-2019) at Cedar Key, Florida, to 9.7 mm/year (36 years: 1939-1974) at Eugene Island, Louisiana. Ten tide gauge stations along the Florida coast recorded relative sea-level rise rates varying between 2 to 4 mm/year. Two tide gauge stations along the Alabama coast recorded relative sea-level rise rates approximate to 4 mm/year. One tide gauge station at the Mississippi coast recorded relative sea-level rise rate of 4.6 mm/year. There are three long-history tide gauges along the coast of Louisiana. Two stations located at the frontier of the Mississippi River delta recorded significant relative sea-level rise rates: Eugene Island: 9.7 mm/year (36 years: 1939-1974) and Grand Isle: 9.1 mm/year (73 years: 1947-2019). One station located at the back of the Mississippi River delta recorded a slower sea-level rise rate: 5.4 mm/year (New Canal, 38 years: 1982-2019). There are nine long-history tide gauge stations along the coast of Texas. Five stations north of Corpus Christi recorded relative sea-level rise rates from 5 to 7 mm/year. The other four Texas stations south of Corpus Christi, as well as four gauges along the Mexico coast and one station at the coast of Cuba, recorded sea-level rise rates at 2 to 3 mm/year.
Since all tide gauges along the GOM coasts are within the same oceanic basin, the absolute sea-level change rates at different sites should be approximately the same. Thus, the variations in relative sea-level rise rates seen here ( Figure 11) primarily reflect differences of local land movements (subsidence or uplift). We selected a closely-spaced GNSS station to estimate the land vertical movement trend at each tide gauge site. The selection of GNSS stations is based on a combining consideration of the distance to the tide gauge, the observation history of the GNSS, and local tectonic movements. The selected GNSS stations for 30 tide gauge stations are depicted in Figure 11 Since there are numerous GNSS stations along the US portion of the GOM coast, a long-history GNSS station can often be found for a tide gauge station within a few kilometers. The shortest distance between a tide gauge and its corresponding GNSS is within one kilometer, such as the distance (260 m There is no considerable vertical crust deformation within the southern GOM region and there is no massive groundwater and oil gas extractions along the southern coast of GOM. Thus, the vertical ground movement trends derived from these four GNSS stations provide reliable estimates of the long-term vertical land movements at their corresponding tide gauge sites. Figure 14 depicts sea-level rise rates at 30 tide gauge sites with respect to GOM20, which are obtained according to the method illustrated in Figure 13. The average absolute sea-level rise rate with respect to GOM20 is 2.6 mm/year, with a standard deviation of approximate 0.5 mm/year (1b), which suggests a 95% spatial confidence interval of ±1 mm/year (2b). Figure 13c suggests that the temporal uncertainty (95% CI) of the absolute sea-level rise trend is also at a level of approximately ±1 mm/year. The tide gauge at Sabine Pass, Texas, exhibits the most substantial absolute sea-level rise rate (−4.86 mm/year). Sabine pass is the state boundary between Texas and Louisiana. This site is located at the Texas side of the Sabine Pass, approximately 6.4 km inland from the coastal line. Sabine pass is the connection of the Sabine Lake and GOM. The regional ocean currents at the tide gauge site may be slightly different from tide gauge stations on the coastal line. The GNSS station (TXSP) used to estimate the local land movement at Sabine Pass also had a short history of only four years (2016-2019). Accordingly, we exclude the absolute sea-level rise rates at the Sabine Pass site as an outlier in calculating the average sea-level rise rate and its corresponding standard deviation. The high sea-level rise rate (3.6 mm/year) recorded at New Canal, a shipping canal in New Orleans, Louisiana, could also be biased by regional ocean currents.
mm/year (with respect to GOM20). There is no considerable vertical crust deformation within the southern GOM region and there is no massive groundwater and oil gas extractions along the southern coast of GOM. Thus, the vertical ground movement trends derived from these four GNSS stations provide reliable estimates of the long-term vertical land movements at their corresponding tide gauge sites. Figure 14 depicts sea-level rise rates at 30 tide gauge sites with respect to GOM20, which are obtained according to the method illustrated in Figure 13. The average absolute sea-level rise rate with respect to GOM20 is 2.6 mm/year, with a standard deviation of approximate 0.5 mm/year (1б), which suggests a 95% spatial confidence interval of ±1 mm/year (2б). Figure 13c suggests that the temporal uncertainty (95% CI) of the absolute sea-level rise trend is also at a level of approximately ±1 mm/year. The tide gauge at Sabine Pass, Texas, exhibits the most substantial absolute sea-level rise rate (-4.86 mm/year). Sabine pass is the state boundary between Texas and Louisiana. This site is located at the Texas side of the Sabine Pass, approximately 6.4 km inland from the coastal line. Sabine pass is the connection of the Sabine Lake and GOM. The regional ocean currents at the tide gauge site may be slightly different from tide gauge stations on the coastal line. The GNSS station (TXSP) used to estimate the local land movement at Sabine Pass also had a short history of only four years (2016-2019). Accordingly, we exclude the absolute sea-level rise rates at the Sabine Pass site as an outlier in calculating the average sea-level rise rate and its corresponding standard deviation. The high sealevel rise rate (3.6 mm/year) recorded at New Canal, a shipping canal in New Orleans, Louisiana, could also be biased by regional ocean currents.  The method for delineating an absolute sea-level trend assumes that both relative sea-level and vertical land movements retain linear trends over time. Nevertheless, the assumption may be valid only in a limited time window. These 30 tide gauge stations have an average history of 56 years, while their corresponding GNSS stations have a much shorter history, approximately one decade. Among these 29 GNSS stations, twenty-one stations (72%) recorded vertical movement trends within ±2 mm/year; six stations (20%) recorded downward land movement trends between 2 to 3 mm/year; only two stations at the Mississippi delta coast recorded downward movement trends of approximately 6.5 mm/year (GRIS: 6.5 mm/year; LMCN: 6.4 mm/year). From the view of plate tectonics, the eastern (Florida), western (western Texas, Mexico), and southern (Mexico, Guba) portions of the GOM coast are stable over time, as demonstrated in Figure 5. Accordingly, it is reasonable to extend the GNSS derived land movement trends to the observational history of tide gauges. Since these 30 tide gauges have an average history of approximately a half-century, we suggest that the interpolation of 2.6 ± 1 mm/year (95% CI) sea-level rise rate with respect to GOM20 be limited to the past five decades (the 1970s-2010s). The 2.6 ± 1 mm/year sea-level rise rate with respect to GOM20 can be regarded as a regional adjustment to the global mean sea-level rise (GMSLR). In general, the estimates of GMSLR during the 20th-century range from 1.6 to 3.0 mm/year [71][72][73][74][75][76]. These studies also emphasized temporal variations of the tide-gauge-derived GMSLR. Holgate [71] reported global sea-level rise rates of 2.03 ± 0.35 mm/year (1904-1953), 1.45 ± 0. 34 mm/year (1954-2003), and 1.74 ± 0. 16 mm/year (1900-2000). Church and White [74] reported GMSLR of 1.7 ± 0.2 mm/year (1990-2009) and 1.9 ± 0.4 mm/year ). Hay et al. [76] reported GMSLR of 1.2 ± 0.2 mm/year , 90% CI) and 3.0 ± 0.7 mm/year (1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010). It appears that the 2.6 ± 1 mm/year seal-level rise rate along the GOM coast is slightly higher than the recently reported GMSLR. The 2.6 ± 1 mm/year is also slightly higher than recently reported average absolute sea-level rise rate along the GOM coast. Letetrel et al. [77] reported an absolute sea-level rise rate of 2.0 ± 0.4 mm/year based on long-history observations from five closely-spaced tide gauge and GNSS pairs and satellite altimetry data along the GOM coast.

GOM20 Versus NAVD88
The US National Spatial Reference System (NSRS) is currently realized by the North American Datum of 1983 (NAD83) and the North American Vertical Datum of 1988 (NAVD88). NAD83 and NAVD88 are the official horizontal and vertical datums for coastal mapping and subsidence studies in the conterminous United States (CONUS). However, both NAD83 and NAVD88 rely on physical survey marks on land that deteriorate over time. To improve NSRS, NGS will replace NAD83 and NAVD88 with the North American Terrestrial Reference Frame of 2022 (NATRF2022) and the North American-Pacific Geopotential Datum of 2022 (NAPGD2022) in 2022 (https://geodesy.noaa.gov). NAVD88 was established in 1991 by the minimum-constraint adjustment of geodetic leveling observations in North American continental covering Canada, the United States, and Mexico. In the adjustment, only the height of the primary tidal benchmark at Father Point, Rimouski, Quebec, Canada, was held fixed. The Father Point is approximately 2500 km away to the north coast of GOM. The definition of NAVD88 uses the Helmert orthometric height, which can be obtained by the geoid height at the site plus the ellipsoid height derived from GNSS observations [37]. Currently, the geoid height with respect to the NAD83 ellipsoid can be estimated by the most recent geoid model GEOID18 (https://www.ngs.noaa.gov/GEOID/GEOID18/). In practice, the geoid height at a site is regarded as a constant number over time. Thus, a trend of the orthometric height time series (with respect to NAVD88) is the same as the trend of the ellipsoidal height time series (with respect to NAD83) [37]. Figure 15a depicts the three-component displacement time series at GNSS station GRIS with respect to NAD83 and GOM20. To get the NAD83 coordinate time series, firstly, the IGS14 coordinates are transformed to IGS08 coordinates according to the method provided by IGS [12], and the IGS08 coordinates are transformed to NAD83 coordinates according to the method provided by NGS [78]. The displacement time series with respect to NAD83 indicate a movement toward the east direction with a steady velocity of 1.4 mm/year and a downward movement (subsidence) of 8.3 mm/year. The subsidence rate with respect to NAD83 is approximately 2 mm/year faster than the subsidence rate with respect to GOM20 (6.5 mm/year). GRIS is closely-spaced (250 m) with the tide gauge located on the Grand Isle, Louisiana ( Figure 12). Figure 14 suggests that the absolute sea-level rise rate at the Grand Isle, Louisiana tide gauge site is 2.6 mm/year if GOM20 is used as a reference frame. However, the absolute sea-level rise rate would be reduced to 0.8 mm/year if NAD83 (or NAVD88) is used as a reference frame. The difference of 2 mm/year is remarkable for sea-level studies. It is not a simple question as to which value is right or wrong because they specifically refer to different reference frames. A significant difference between NAVD88 and GOM20 is that GOM20 used references (points) adjacent to the Gulf coast (within 500 km to the coastal line), while NAVD88 used references are much further away from the Gulf coast.
Nevertheless, the references closer to the Gulf coast will share more common ground movements with the GOM coastal area than those references far away from the coast. Glacial isostatic adjustment (GIA) is a typical common vertical movement that can be entirely removed by the regional reference frame. As to the vertical land movements due to GIA, they are spatially uniform within the GOM region and have amplitudes less than 1 mm/year [79]. Certainly, using references adjacent to the coastal area would make more sense to appraise the long-term risks associated with coastal subsidence, sea-level rise, flooding, and erosion.
Remote Sens. 2020, 12, 350 23 of 29 Figure 15. Three-component displacement time series at GNSS station GRIS with respect to (a) GOM20 versus NAD83, and (b) with respect to GOM20 versus NA12. GRIS is located on the barrier island, Grand Isle, Louisiana (see Figure 12).

GOM20 Versus NA12
NA12 is another continental-scale reference frame that has also been applied for geological hazards and crustal motion studies in North American [34]. NA12 is designed to have no-net-rotation (NNR) with respect to the stable interior of the North American Plate. Figure 15b illustrates the threecomponent displacement time series at GNSS station GRIS with respect to GOM20 and NA12. The time series with respect to NA12 is obtained from the Nevada Geodetic Laboratory (http://geodesy.unr.edu). No outliers have been removed. The comparison indicates that the longterm ground deformation trends with respect to both reference frames are almost identical.
Nine of the 30 core reference stations used for realizing NA12 are also used as reference stations for realizing GOM20.  Table A1). One of the major differences between NA12 and GOM20 is that NA12 employs the daily 7-parameter transformation method for transforming EEF-XYZ coordinates from the global reference frame (IGS14) to the regional reference frame (NA12). In contrast, GOM20 employs the total 7-parameter transformation method for transforming EEF-XYZ coordinates from the global reference frame (IGS14) to the regional reference frame (GOM20). NA12 users need daily 7-parameters provided by the Nevada Geodesy Laboratory to convert IGS14 coordinates to NA12. In contrast, GOM20 users only need the seven parameters listed in Table 1 to convert the whole time series from IGS14 to GOM20. Furthermore, the 7-parameters used by GOM20 are unattached to the future status of these reference stations. That is to say, GOM20 provides more convenient and reliable Figure 15. Three-component displacement time series at GNSS station GRIS with respect to (a) GOM20 versus NAD83, and (b) with respect to GOM20 versus NA12. GRIS is located on the barrier island, Grand Isle, Louisiana (see Figure 12).

GOM20 Versus NA12
NA12 is another continental-scale reference frame that has also been applied for geological hazards and crustal motion studies in North American [34]. NA12 is designed to have no-net-rotation (NNR) with respect to the stable interior of the North American Plate. Figure 15b illustrates the three-component displacement time series at GNSS station GRIS with respect to GOM20 and NA12. The time series with respect to NA12 is obtained from the Nevada Geodetic Laboratory (http://geodesy.unr.edu). No outliers have been removed. The comparison indicates that the long-term ground deformation trends with respect to both reference frames are almost identical.
Nine of the 30 core reference stations used for realizing NA12 are also used as reference stations for realizing GOM20.  Table A1). One of the major differences between NA12 and GOM20 is that NA12 employs the daily 7-parameter transformation method for transforming EEF-XYZ coordinates from the global reference frame (IGS14) to the regional reference frame (NA12). In contrast, GOM20 employs the total 7-parameter transformation method for transforming EEF-XYZ coordinates from the global reference frame (IGS14) to the regional reference frame (GOM20). NA12 users need daily 7-parameters provided by the Nevada Geodesy Laboratory to convert IGS14 coordinates to NA12. In contrast, GOM20 users only need the seven parameters listed in Table 1 to convert the whole time series from IGS14 to GOM20. Furthermore, the 7-parameters used by GOM20 are unattached to the future status of these reference stations. That is to say, GOM20 provides more convenient and reliable access to users than NA12.

Broad Applications of GOM20
GOM20 provides a consistent platform to integrate positional observations over space and time from different remote sensing techniques to a unified geodetic reference and enables multidisciplinary and cross-disciplinary research. The long-term GNSS observations within the GOM region have accumulated valuable datasets for the research community. Hydrologists may use these datasets for studying ground deformation resulting from fluid withdrawal, aquifer compaction, and seasonal hydrologic and atmospheric pressure loading; geomorphologists may use these datasets for studying long-term coastal erosion and wetland loss problems along the GOM coast; oceanographic and sea-level researchers may use these datasets for calibrating long-term sea-level changes along the GOM coast; researchers in the field of remote sensing may use GNSS-derived displacement time series as "ground truth" for calibrating or validating subsidence estimations from remote sensing techniques, such as InSAR, airborne LiDAR, and UAV-photogrammetry. The long-term GNSS measurements also provide first-hand ground truth datasets for local governments to assess risks from coastal subsidence and sea-level rise and make plans for future coastal development. Hopefully, this study will further stimulate coastal hazards investigating and land-use planning discussions along the GOM coast.
GOM20 also enables an approach of using stand-alone GNSS to conduct accurate ground deformation monitoring, as well as structural health monitoring, such as the long-term stability monitoring of offshore platforms within GOM. By transforming PPP solutions to a regional reference frame, users do not need to install any ground reference stations in the field and do not need to include any reference data in their data processing. The stand-alone GNSS surveying method will substantially reduce field logistics costs and, therefore, will revolutionize the way for conducting geological hazards (landslides, subsidence, faulting) and structural health monitoring in the GOM region.

Conclusions
Current geodetic infrastructure within the GOM region makes it possible to precisely track long-term (e.g., > 3 years) ground deformation at the level of 0.5 mm/year and above using stand-alone GNSS observations. This study summarized the current GNSS stations along the GOM coast and established a stable regional reference frame. The primary product from this study is the seven parameters (Table 1) for converting the ECEF-XYZ coordinate time series from the global reference frame IGS14 to the regional reference frame GOM20. The frame stability of GOM20 is approximately 0.3 mm/year in the horizontal directions and 0.5 mm/year in the vertical direction. GOM20 can be confidently used from the 1990s to 2030 without causing positional errors larger than the precision of daily-PPP solutions. It's certainly important to be cautious when interpreting site velocities at a level of submillimeter per year. GOM20 will be incrementally improved and synchronized with future updates of IGS reference frames. According to this study, subsidence faster than 2 cm/year is ongoing in several major cities in central Mexico, and the most rapid urban subsidence is approximately 27 cm/year in Mexico City; a large portion of the Texas and Louisiana coasts is subsiding at 3 to 6.5 mm/year, and the most rapid coastal subsidence is up to 6.5 mm/year at the coastline of Mississippi River delta; the present average sea-level rise rate along the GOM coast is 2.6 ± 1mm/year (95% CI) with respect to GOM20.
Author Contributions: G.W., X.Z. and Y.B. prepared the original draft. All authors edited, reviewed, and improved the manuscript. All authors have read and agreed to the published version of the manuscript.

Acknowledgments:
The authors appreciate NGS, UNAVCO, and the Nevada Geodetic Laboratory (NGL) for sharing their GNSS products.

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