GPS-derived ground deformation ( 2005 – 2014 ) within the Gulf of Mexico region referred to a stable Gulf of Mexico reference frame

This study investigates current ground deformation derived from the GPS geodesy infrastructure in the Gulf of Mexico region. The positions and velocity vectors of 161 continuous GPS (CGPS) stations are presented with respect to a newly established local reference frame, the Stable Gulf of Mexico Reference Frame (SGOMRF). Thirteen long-term (> 5 years) CGPS are used to realize the local reference frame. The root mean square (RMS) of the velocities of the 13 SGOMRF reference stations achieves 0.2 mm yr−1 in the horizontal and 0.3 mm yr−1 in the vertical directions. GPS observations presented in this study indicate significant land subsidence in the coastal area of southeastern Louisiana, the greater Houston metropolitan area, and two cities in Mexico (Aguascalientes and Mexico City). The most rapid subsidence is recorded at the Mexico City International airport, which is up to 26.6 cm yr−1 (2008–2014). Significant spatial variation of subsidence rates is observed in both Mexico City and the Houston area. The overall subsidence rate in the Houston area is decreasing. The subsidence rate in southeastern Louisiana is relatively smaller (4.0–6.0 mm yr−1) but tends to be steady over time. This poses a potential threat to the safety of coastal infrastructure in the long-term.


Introduction
The Gulf of Mexico (GOM) region has been the heart of the U.S. energy industry because of substantial oil and gas deposits along the coast and offshore of the GOM.It is heavily populated and vulnerable to local ground deformation (faulting, subsidence, uplift) and relative sea-level rise (e.g., Day et al., 1995;Kolker et al., 2011;Thatcher et al., 2013).Land subsidence and faulting problems in the GOM region have been frequently investigated by different research groups using GPS observations (e.g., Dokka, 2011;Engelkemeir et al., 2010;Kearns et al., 2015;Khan et al., 2014;Osmanolu et al., 2011;Wang and Soler, 2013).However, it is difficult to align the results from these research groups because they used different data sets collected by different organizations during different time periods.Furthermore, they focused on localized ground deformation and applied different reference points or frames.This study aims to establish a unified local geodetic reference frame, the Stable Gulf of Mexico Reference Frame (SGOMRF), to investigate the current ground deformation within the whole GOM region during the past decade (2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014).Observations of 161 high-quality continuous GPS stations among 450 active long-term GPS stations (Fig. 1) are investigated in this study.Land subsidence and faulting in the Houston region, Mexico City, and the southeastern Louisiana region are discussed and compared.

GPS data processing
This study applies the precise point positioning (PPP) method, for solving the 24-hour average position of a GPS antenna.PPP is based on the processing of the following ionosphere-free combinations of the undifferenced code and phase observations (Zumberge et al., 1997): Published by Copernicus Publications on behalf of the European Geosciences Union.
where f 1 and f 2 are the GPS L1 and L2 frequencies; P (Li) and (Li) are the code and phase observations at the corresponding frequency; ρ is the true range; c is the speed of light; dT is the receiver clock offset; d trop is the tropospheric delays; N i is the phase ambiguity term in (Li).Therefore, the unknown parameters estimated in PPP include position coordinates, phase ambiguity, receiver clock offset and the tropospheric delays.GNSS-Inferred Positioning System and Orbit Analysis Simulation Software (GIPSY-OASIS) package (V6.3)developed at the Jet Propulsion Laboratory (JPL) is applied for calculating daily positions.The GIPSY-OASIS package provides single receiver phase ambiguity fixed PPP solutions.The single receiver phase ambiguity method uses the wide lane and phase bias estimates obtained from a global network of ground GPS stations to perform ambiguity-resolved PPP resolution (Bertiger et al., 2010).The major parameters estimated and key models applied in the PPP processing include: the VMF1 troposphere mapping model (Boehm et al., 2006), second-order ionospheric delay (Kedar, 2003), the ocean tidal loading model FES2004 (Lyard et al., 2006) calculated through the free online service operated by Onsala Space Observatory, Sweden (Free ocean tide loading provider, 2015), tropospheric gra-dient (Bar-Sever et al., 1998), zenith troposphere delay as a random walk with variance of 5 × 10 −8 km s −1 , gradient troposphere wet delay as a random walk with variance of 5 × 10 −9 km S −1 , and receiver clock as white noise with updates every measurement epoch.Station coordinates are initially provided in the loose frame of the JPL's fiducial-free GPS orbits.The coordinates are then transformed into the International GNSS Service Reference Frame of 2008 (IGS08) using the daily seven transformation parameters that are delivered with the JPL's orbit products.
The PPP processing conducts all calculations within an Earth-centered-Earth-fixed (ECEF) geocentric coordinate system (x, y, and z).In order to track land surface deformation, the ECEF geocentric coordinates are converted to cartographic (northing, easting, and ellipsoidal height) coordinates, which are referred to the GRS-80 ellipsoid.The threecomponent daily positional time series indicates the change of ground surface over time at different directions.The three components represent: north, east, and up.The up component (subsidence/uplift) measurements used in this study are obtained by differencing GPS measured ellipsoidal heights referred to the local reference frame described in the next section.Detection of regional-scale land subsidence has historically depended on surveying benchmarks periodically.This has traditionally been accomplished by differencing orthometric heights obtained from spirit leveling.A recent investigation conducted by Wang and Soler (2014) has indicated that using ellipsoid and orthometric heights would result in the same practical subsidence measurements.Accordingly, the subsidence values used in this study can be regarded as having the same "physical meaning" as the conventional subsidence measurements obtained from leveling surveys.For each time series, the outliers, defined as the days for which the uncertainty was greater than 2.0 times of the average uncertainty of the entire measurement, were removed (Firuzabadì and King, 2011;Wang, 2011).The uncertainty of each measurement was directly output by the GIPSY-OASIS program.On average, 5 % of the total samples are removed as outliers.The daily positional time series applied in this article are the "cleaned" time series.

Stable Gulf of Mexico Reference Frame (SGMRF)
In general, a global or a continental-scale reference frame is realized with an approach of minimizing the least square residual velocities of a large number of selected reference stations.In the case of IGS08, 232 globally distributed, and well-performing GPS stations are used (Rebischung et al., 2012).The velocities at GPS sites referred to IGS08 are dominated by tectonic drift.For a regional study, a stable local reference frame is often established through Helmert transformation to exclude tectonic drift (Wang et al., 2013(Wang et al., , 2014(Wang et al., , 2015a)).It facilitates the precise physical interpretation of local ground deformation over time and space.The transformation involves seven parameters including a rotation vector, a translation vector and one scale factor.These seven transformation parameters can be estimated by comparing the positions of a group of selected reference stations referred to the new reference frame with those referred to a well-established reference frame.
In practice, at least three reference stations are needed to obtain the transformation parameters.More reference stations often result in a more reliable coordinate transformation.However, a reference station that is not locally stable will degrade the overall performance of the frame transformation.A stable site is defined as retaining zero velocities (three components) with respect to a specified reference frame.The stability (precision) of a local reference frame is therefore affected by the velocities of stable sites with respect to the reference frame.Thus, the selection of reference stations is critical for establishing a stable local reference frame.In general, there is not a fixed criterion for selecting reference stations.The selection mostly is based on the availability of long-term CGPS stations in the study area.There are over 780 CGPS stations in the GOM region (Fig. 1).As this study uses a secular frame, the linearity of daily positional time series is a critical criteria for selecting reference stations (Blewitt and Lavallée, 2002).Additionally, the geographic distribution of reference stations is also considered.The following specific criteria are initially applied for selecting reference stations: 1. Having segments of data spanning at least 5 years (installed in 2009 or earlier) with no steps (a sharp change of the mean in positional time series caused by an earthquake, equipment change, or other unknown reasons).
2. No considerable subsidence or uplifting (the linear velocity rate of vertical positional time series referred to IGS08 is less than 0.5 mm yr −1 ); 3. having less than 0.1 mm yr −1 "standard error" (σ ) of the slope (V cal , V cal ) of the geocentric coordinate time series (x, y, and z) referred to IGS08.σ is a measure of the error in the precision with which V cal has been estimated by a linear regression.A smaller σ indicates a small margin of error.Approximately 95 % of the time, the true velocity will be contained in the interval between V cal − 1.96 × σ and V cal + 1.96 × σ .
The near-coast areas could be affected by subsidence and coastal erosion problems (Simms et al., 2013;Williams et al., 1997;Yu et al., 2014).Accordingly, it is preferred that reference stations be located inland rather than within nearcoast areas.However, in order to balance the overall coverage and geometrical distribution of reference stations, one near-coast station in Florida (RMND) and one near-coast station in Mexico (TAM1) were selected as reference stations (Fig. 2).Uneven distribution of reference sites could lead to biases in frame transformation (Collilieux et al., 2010;Wang et al., 2014).Positional time series affected by steps are identified by an automated edge detection program based on the derivative of Gaussian kernel (Canny, 1986).Initially, 30 CGPS stations were selected as reference stations and the reference frame transformation were calculated.Any station that had a horizontal velocity larger than 0.5 mm yr −1 with respect to the resulted local reference frame was removed from the group of reference stations and the transformation was recalculated again.Finally, 13 CGPS stations are selected as reference stations for realizing the SGOMRF (Fig. 2).Two different approaches are often used in geodesy to transform positional time series from one reference frame to another: the daily 7-parameter Helmert transformation (Blewitt et al., 2013) and the 14-parameter similarity transformation.This study applies a 14-parameter transformation approach that has been frequently applied in the geodesy surveying community (e.g., Pearson and Snay, 2012;Wang et al., 2014).The geocentric coordinates of a station with respect to SGOMRF are calculated by the following formulas: Here, T x (t), T Y (t) and T Z (t) are translations along x, y and z axis; R X (t), R Y (t) and R Z (t) are counterclockwise rotations about three axes; s(t) is a differential scale factor between IGS08 and SGOMRF.These seven parameters at any point in time are specified relative to a reference epoch by the following linear relationships: Here, t 0 denotes a specific epoch (e.g., 2013.0).IGS08 with Eqs. ( 3) and ( 4) and 14 parameters provided by Pearson and Snay (2012) (Table 1).Both sites retain nearzero velocities (< 0.5 mm yr   2 suggest that the PPP solutions obtained in this study achieve 2 mm horizontal precision and 7 mm vertical precision, which is comparable with the overall precision of GIPSY-OASIS PPP solutions.Bertiger et al. (2010) reported that the single receiver ambiguity fixed PPP solutions achieved overall 2 mm horizontal precision and 6 mm vertical precision for the 106 worldwide IGS reference stations.It should be noticed that the 14-parameter transformation processing did not improve the precision of velocity estimates since a linear regression model was applied to the changes of the scale, three translational, and three rotational motions in time.

Horizontal ground deformation
This study investigated GPS data from 148 stations that have step-free time spans of longer than 4 years, and transform the daily positional time series to the local reference frame  (Fig. 4).The average horizontal velocity of these 148 stations is below 1 mm yr −1 , which implies that the interior of the GOM region is rigid at the level of sub-millimeter per year.Five long-term GPS stations have been moving horizontally with relatively larger velocities (> 2 mm yr −1 ).These stations are MMX1 (9.7 mm yr −1 ), FSHS (3.4 mm yr −1 ), UNIP (2.9 mm yr −1 ), TXPR (2.4 mm yr −1 ) and ROD1 (2.4 mm yr −1 ).Four of these stations are located in wellknown subsidence areas.MMX1 is located in Mexico City, Mexico.FSHS is located in Franklin, Louisiana.UNIP is located in Aguascalientes, Mexico.ROD1 is located in Houston, Texas.
In the case of groundwater withdrawal induced subsidence, a subsidence bowl can be formed by localized aquifer compaction.In such an event, it is possible that GPS stations are pulled horizontally towards the center of the subsidence bowl (e.g., Allis, 2000;Bawden et al., 2001Bawden et al., , 2012)).Depending on the position of a station relative to the subsidence center, the ratio of horizontal to vertical velocities varies.A station at the edge of the subsidence center will show a relatively large horizontal velocity.A station located closer to the center of subsidence could display large velocities in both the horizontal and vertical components.Around the center of the subsidence feature, the stations will mostly exhibit vertical movement.Figure 5 depicts the three-component positional time series at ROD1, TXCN, and TXLI referred to the local reference frame.TXLI is a stable station located in Liberty County, Texas.ROD1 shows a 2.3 mm yr −1 horizontal movement towards the northeast.TXCN shows no movement in the horizontal direction.Both ROD1 and TXCN indicate long-term subsidence.According to our previous studies on the current subsidence in the greater Houston metropolitan area (Kearns et al., 2015;Yu et al., 2014), a rapid-subsidence bowl is forming around The Woodlands area.The Woodlands is a vibrant and fast-growing business and entertainment suburban area located 43 km north of downtown Houston (Fig. 6).Groundwater is the sole water source for residential and business use in this area as of 2014.The subsidence rate in the center of the subsidence bowl is about 25 mm yr −1 (Kearns et al., 2015).ROD1 is located in the city of Spring, Texas, northern Harris County.The station is 12 km southwest to The Woodlands area (Fig. 6).The subsidence rate at ROD1 is 17 mm yr −1 derived from the whole time series from 2007 to 2014.The closely spaced contour lines represent greater spatial variation of subsidence rates in this developing subsidence bowl.The horizontal velocity vector indicates that ROD1 is affected by the differential subsidence.TXCN is located in the city of Conroe, Texas, which is 21 km north of The Woodlands.The positional time series (2008)(2009)(2010)(2011)(2012)(2013)(2014) of TXCN does not indicate any considerable horizontal movement (< 1 mm yr −1 ).However, steady land subsidence with a rate of 16 mm yr −' has been recorded at this site.The contour lines in this area are spaced far apart.The different vertical-to-horizontal velocity ratio suggests that TXCN is located at a more uniformly subsiding, therefore, flat area.whereas ROD1 is more likely located along the steep sidewall of the developing subsidence bowl.The comparison of the three-component positional time series of ROD1 with those of TXCN illustrates a good example of horizontal velocity variations around a subsidence bowl.
A similar movement pattern at two CGPS sites around a subsidence bowl is also observed in Mexico City (Fig. 7).MMX1 is located in the Mexico International Airport, eastern Mexico City.UNIP is located at the Universidad Nacional Autónoma de México, southwestern Mexico City.UNIP records a smaller rate of 2.9 mm yr −1 towards the northeast.Observations from InSAR indicated that subsidence rates in Mexico City increased eastwards towards the center of the Basin of Mexico (Chaussard et al., 2014).The horizontal movements of MMX1 and UNIP agree well with the subsidence bowl illustrated by the InSAR data: both UNIP and MMX1 are moving toward the center of the subsidence bowl.MMX1 is located much closer to the center and therefore has demonstrated higher rates of horizontal motion (Fig. 7).
Figure 8 depicts that FSHS, a permanent GPS station (2010-2014) located at Franklin, Louisiana, has been moving toward the southeast at a rate of 3.4 mm yr −1 .The antenna of FSHS is mounted on a reinforced concrete building located at Franklin High School.Subsidence at this site is consistent and could lead to more serious problem over the long term.There is no known excessive groundwater withdrawal issue in this area.It is not likely that the horizontal movement is associated with an on-going subsidence bowl.Our other GPS sites (AWES, DSTR, HOUM, GRIS, BVHS, LMCN) in southeastern Louisiana also demonstrate movements southward with rates smaller than 1 mm yr −1 .Dokka et al. (2006) observed a similar southward displacement and  proposed the detaching of the South Louisiana Allochthon which contributed to both the southward motion and subsidence in this area.Latter studies, however, suggested the major driving mechanism of subsidence in this area is the  9a).TXPR is located at Pharr, Texas (Fig. 4).The reference frame is SGOMRF.compaction of shallow strata with minor or no contribution from active faulting and deep crustal processes (Blum et al., 2008;Edrington et al., 2008;Törnqvist et al., 2008;Wolstencroft et al., 2014;Yu et al., 2012).Our data show a lower rate of southward displacement, which may indicate that the horizontal motion is mainly caused by differential compaction and minor contribution from downslope movement of listric faults.The large horizontal velocity at FSHS may be a sitespecific feature related to the obvious seasonal motion in the EW component.TXPR, a permanent GPS site (2005-2014) located at Pharr, Texas, has a horizontal southwest movement of 2.4 mm yr −1 (Fig. 8).The antenna pole is anchored on a wall of an office building owned by the Texas Department of Transportation.GPS observations also show steady subsidence (5.8 mm yr −1 ) at this site.The horizontal movement could be associated with local subsidence.Further study is needed to verify the cause of the horizontal motion.

Vertical ground deformation
The overall spatial variation of vertical velocities in the GOM region is much greater compared to that of horizontal velocities.There are certain stations showing extremely large downward vertical velocities in this region.Figure 4b indicates four rapid subsidence zones in the GOM region -the southeastern Louisiana area, the Houston metropolitan area, Aguascalientes, and Mexico City.The drivers of subsidence vary from place to place.Different drivers would result in different subsidence patterns -the spatial and temporal variability of subsidence rates.In this section, we discuss subsidence in southeastern Louisiana, Houston, Aguascalientes and Mexico City.

Subsidence in the southeastern Louisiana
The causes behind present subsidence in southeastern Louisiana have been controversial and heavily studied.Ram-sey and Moslow (1987) attributed 80% of the present subsidence on the coast of Louisiana to "compactional subsidence".Roberts et al. (1994) studied relationships between subsidence rates and faulting, land loss, thickness and characteristics of Holocene sediment in the Louisiana coastal area, and they concluded that sediment compaction was a primary cause of subsidence.A number of studies proposed that the present-day subsidence in the Mississippi Delta is mostly caused by the isostatic response to the delta load (e.g., Ivins et al., 2007;Jurkowski et al., 1984).However, Dokka (2006) argued the conventional opinions; using a case study conducted in the Michoud area of Orleans Parish, Louisiana, Dokka (2006) concluded that 73 % (16.9 mm yr −1 ) of subsidence during the period 1969-1971 and 50 % (7.1 mm yr −1 ) of subsidence during the period 1971-1977 was attributed to tectonism (fault movements).Dokka and his colleagues further addressed tectonic-induced subsidence in their other publications (Dixon et al., 2006;Dokka, 2011;Dokka et al., 2006).Wolstencroft et al. (2014) investigated the cause of subsidence in the Mississippi Delta through geophysical modeling and concluded that presentday basement subsidence rate due to sediment loading was less than ∼ 0.5 mm yr −1 and the glacial isostatic adjustment was likely to be the major driver of deep-seated subsidence.
Figure 9a and c illustrate the velocity vectors and positional time series at long-term GPS stations across the southeastern Louisiana area.Considerable subsidence rates are recorded at two near-coast stations: LMCN (foundation depth 36.5 m) and GRIS (foundation depth unknown) (< 10 km to the coastline).Both sites indicate steady subsidence of approximately 6 mm yr −1 over 10 years.Four inland stations FSHS (foundation depth > 5 m), AWES (foundation depth 1 m), HOUM (foundation depth > 15 m) and DSTR (foundation depth unknown) show smaller subsidence rates (2-4 mm yr −1 ).The seaward increase in the rate of subsidence may be a combined result of shallow sediment compaction and deep basement subsidence.Wolstencroft et al. (2014) demonstrated that present-day Pleistocene basement subsidence (deep subsidence) in the Mississippi Delta produced by viscoelastic deformation mechanisms increased seaward.The natural compaction of young sediments could occur in new infill and recently drained mashes (Törnqvist et al., 2008).GPS stations underlain by shallow (Holocene) sediments in this area will be subject to the ongoing compaction.Most data used in this study are from buildingbased stations.Therefore the monuments of GPS stations are building foundations which are typically 5-15 m below the land surface.In this case, the contribution from the uppermost Holocene compaction is minimum and our subsidence rate estimates should be considered minimum estimates.The foundation depth used in this study is collected from Dokka et al. (2006) and GPS log files from NGS.Unfortunately, this information is not always available in the GPS log files.The foundation depth of a GPS site directly determines its measuring target; therefore, it is of great importance to check this ).This rate is comparable to the rate (3.5 mm yr −1 ) reported by Dokka et al. (2006) and substantially smaller than the rate (5.7 mm yr −1 ) reported by Karegar et al. (2015).Note that, due to antenna changes and two large gaps (13 months and 8 months), the data prior to 2010 were not used to calculate the subsidence rate at BVHS.The choice of record length may explain the large difference with the rate reported by Karegar et al. (2015) in which the full data span is used.The reason for the smaller subsidence rate at BVHS compared to LMCN and GRIS is unclear.Morton and Bernier (2010) showed historical subsidence rates calculated from repeat leveling surveys at benchmarks along state highway LA23 between Chalmette and Venice (near the location of BVHS), state highway LA1 between Raceland and Grand Isle (near the location of GRIS), and state highway LA 56 between Houma and Cocodrie (near the location of LMCN).The 1965/1966 to 1993 average subsidence rates along LA1 and LA56 was 9.6 and 11 mm yr −1 , respectively, whereas, the subsidence rates along LA23 near the location of BVHS was much larger with greater than 25 mm yr −1 between 1964 and 1971 and decreased to about 18 mm yr −1 between 1971 to 1984.The leveling lines encompass a much larger area compared to the spot measurements by GPS stations.Therefore, the subsidence rate measured at the single spot BVHS may not be able to represent the subsidence rate in a larger area.
Despite the different subsidence rates between coastal sites and inland sites, the overall spatial variation of subsidence rates across southeastern Louisiana is relatively smaller compared to that of the Houston metropolitan area.The slight variation of subsidence rates in space and the steady subsidence in time suggest that the subsidence in southeastern Louisiana is not likely dominated by the compaction of shallow aquifers associated with groundwater pumping.Groundwater pumping induced subsidence often shows considerable spatial and/or temporal variations as illustrated by subsidence in the metropolitan area of Houston and Mexico City (discussed in the following sections).In southeast Louisiana, groundwater withdrawal is minimal because groundwater quality is affected by saltwater encroachment (Baumann et al., 2006;Meckel, 2008).The exception is for the greater New Orleans area, where groundwater is pumped from shallow upper Pleistocene aquifers.Spatial correlations between areas of large subsidence and areas with high-yield groundwater wells in the New Orleans area were reported by Dokka (2011).
Hydrocarbon production has been frequently discussed as one of the possible drivers for subsidence in southeastern Louisiana (Chang et al., 2014;Kolker et al., 2011;Meckel, 2008).The oil production in this region peaked at ∼ 446 million barrels in 1968 but has decreased consistently to less than ∼ 55 million barrels after 2005 (Kolker et al., 2011;Meckel, 2008).The data used in this study are mostly from 2005, therefore we consider the possible contribution from hydrocarbon production to the subsidence observed in our data is marginal.

Subsidence in the Houston area
The groundwater induced subsidence in the Houston area has been intensively investigated by researchers from the US Geological Survey (USGS) (Bawden et al., 2012;Galloway et al., 1999;Johnson et al., 2011;Kasmarek et al., 2009Kasmarek et al., , 2010Kasmarek et al., , 2012Kasmarek et al., , 2013)), National Geodetic Survey (NGS) (Zilkoski et al., 2003), and local research institutions (e.g., Engelkemeir et al., 2010;Kearns et al., 2015;Khan et al., 2014;Qu et al., 2015;Wang and Soler, 2013;Wang et al., 2015b).A recent study conducted by Yu et al. (2014) indicated that the subsidence in the Houston metropolitan area is attributed by the compaction of aquifers within about 500 m to the ground surface.Since groundwater usage changes according to the local population and land usage, the subsidence resulting from groundwater withdrawal will vary in space.The subsidence rate will also change over time in accordance to the city development and policy changes.Historically, subsidence in the Houston area had been primarily occurring in the eastern and southeastern portions.A comparison of the current subsidence and recent subsidence  mapped by the USGS (Bawden et al., 2012) indicates that the subsidence in the Houston area has been migrating to the western and northern areas since the 1990s.The overall subsidence has also been reduced significantly as a result of rigidly enforced groundwater regulation plans (Harris-Galveston Subsidence District, 2013).
Figure 9b and d illustrate the velocity vectors and positional time series at long-term GPS stations within the Houston area.The spatial variation of subsidence in the Houston area is more significant than that in southeastern Louisiana.The highest subsidence rates in the Houston area are recorded at two inland sites: TXCN and ROD1 (Table 3).These two stations are more than 100 km away from the coastline.The subsidence rate is as high as 17 mm yr −1 at TXCN.Two coastal stations TXGV and TXGA record subsidence rates of only 1.3 and 4.4 mm yr −1 , respectively.Considerable temporal variations in subsidence rates are also identified at several sites.For instance, the subsidence rate at ROD1 was 25 mm yr −1 before 2010 and has reduced to 13 mm yr −1 after 2010 due to the enforced groundwater regulation since 2010 (Fig. 5) (Harris-Galveston Subsidence District, 2013).

Salt tectonics
The Gulf Coast area is one of the world's largest salt dome regions.Over 500 salt domes have been discovered onshore and under the sea floor of the GOM (Beckman and Williamson, 1990).Long-term accumulation of the salt movements could exert an impact on surface morphological features and cause fault growth.The vertical velocity vectors illustrated in this study do not show any considerable vertical movement that can be resulted from salt dome uplift.Jackson and Seni (1983) showed that the maximum net growth rate of diapirs in East Texas is 150 to 230 m per million years, which equals ∼ 0.2 mm yr −1 .The ground deformation rate at this level is below the limit that can be identified with the current GPS geodesy infrastructure in this region.

Subsidence in central Mexico
The CGPS sites in two central Mexico cities -Aguascalientes and Mexico City show extremely rapid subsidence.The causality between groundwater extraction and land subsidence in Mexico City was first investigated in the 1930s and then the 1940s (Carrillo, 1947;Cuevas, 1936).Groundwater accounts for nearly half of Mexico City's water usage (Sosa-Rodriguez, 2010).As of 2011, shallow aquifers in this region had been seriously overexploited (Engel et al., 2011).The ground water level has been declining at aver-Table 3. Horizontal (V h ) and vertical velocities (V v ) of GPS stations plotted in Fig. 9.

Houston
Southeastern Louisiana Reference frame: SGOMRF Reference frame: SGOMRF  , 1995). Figure 10 illustrates the threecomponent positional time series at three CGPS sites: INEG, MMX1, and UNIP.INEG is located in the city of Aguascalientes, Mexico.This station shows a steady subsidence rate of 25.7 mm yr −1 .MMX1 records steady land subsidence at a rate of 266.3 mm yr −1 , which could be the most rapid subsidence rate ever recorded by a CGPS station.The subsidence rate at the center of the subsidence bowl could be even larger.In fact, a subsidence rate as high as 370 mm yr −1 (1996-2005) was derived from InSAR studies in Mexico City (Cabral-Cano et al., 2008).UNIP records a small subsidence rate of 2.7 mm yr −1 .The distance between MMX1 and UNIP is only about 17.8 km.The subsidence rate at MMX1 is about 100 times higher than that at UNIP.It demonstrates the significant spatial variation of the subsidence rate in Mexico City.

Conclusions
This study utilizes the current GPS geodesy infrastructure in the GOM region to investigate the ground deformation associated with subsidence and faulting.A sophisticated regional GPS geodesy infrastructure should include three components: individual GPS stations, a stable local reference frame, and sophisticated positioning software packages.Currently, a unified "local reference frame" does not exist in the GOM region.This study established a stable local reference frame (SGOMRF) to fill the gap.In the first release of the SGOMRF, the 14 Helmert transformation parameters for converting coordinates from the IGS08 to SGOMRF are provided (Table 1).The SGOMRF will be incrementally improved and periodically updated to synchronize with the updates of the IGS reference frame.The potential applica- tions of the local reference frame include providing a consistent framework for precisely monitoring coastal hazards over space and time, studying long-term coastal erosion and wetland loss, studying sea level rise, and comparing research results from different research groups.The stable reference frame will also be useful for long-term health monitoring of dams, sea walls, high-rise buildings, long-span bridges and for planning and designing coastal restoration projects.GPS observations show significant land subsidence in the coastal area of southeastern Louisiana, the Houston metropolitan area, Aguascalientes, and Mexico City.Significant spatial variations of subsidence rates due to differences in groundwater withdraw and clay layer thickness are observed in the Houston area and Mexico City.The decrease of the subsidence rate over time is also observed at GPS stations located in Houston.The GPS sites in the southeastern Louisiana area show steady subsidence and a general south-ward horizontal movement toward the GOM.This may suggest a deep tectonic process associated with faulting.Subsidence resulting from faulting would be difficult to stop through human efforts.As a result, the smaller but steady subsidence (4-6 mm yr −1 ) would cause considerable damage to the coastal protection infrastructure (e.g., sea walls, levees, flood walls, storm surge barriers) in the long term.
GPS observations presented in this study do not show any considerable ground movements that could be associated with salt tectonics or faults.The magnitude of salt dome uplift and faulting may be below the level that can be identified by the current GPS geodesy infrastructure with a time span less than a decade.A denser CGPS network and a longer period of data accumulation are crucial for a more comprehensive study of local ground deformation within the GOM region.
Appendix A: The GPS velocities with respect to SGOMRF

Figure 1 .
Figure 1.Map showing current CGPS stations in the GOM region.Blue triangles represent current CGPS installed in and after 2009.Red circles represent current CGPS installed before 2009.Grey squares represent decommissioned CGPS that have data spanned for more than 5 years.Data are available through NGS and UNAVCO.

Figure 2 .
Figure 2. Map showing the locations and velocity vectors with 95 % confidence ellipses of the 13 reference stations used to define SGOMRF.Black vectors are referred to NAD83; blue vectors are referred to IGS08; red vectors are referred to SGOMRF.

Figure 3
Figure 3 depicts the three-component positional time series of two stations (OKAN and SG05) referred to SGOMRF, North American Datum of 1983 (NAD 83) (2011), and IGS08.The NAD 83 is the horizontal control datum for the US, Canada, Mexico, and Central America(Schwarz, 1989;Snay and Soler, 2000).It is widely used as a North American plate-fixed reference frame in the practice of surveying.The continental-scale reference frame has been updated for several times.The most recent realization is referred as NAD83(2011) at epoch 2010.0.The positional coordinates referred to NAD83(2011) are transformed from −1 ) with respect to the local reference frame SGOMRF.The three-component velocities derived from the 10-year(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014) continuous observations at OKAN are −2.3 mm yr −1 (north), −13.4 mm yr −1 (east), and 0.2 mm yr −1 (up) with respect to IGS08 and 0.6 mm yr −1 (north), 1.6 mm yr −1 (east), and −0.2 mm yr −1 (up) with respect to NAD83.The velocities at SG05 are generally the same (Fig.3).The horizontal velocity (1.5 cm yr −1 ) referred to IGS08 can be explained by the tectonic movement of the GOM region with respect to the global reference frame.However the minor horizontal movements (∼ 2 mm yr −1 ) with respect to the supposedly continent-fixed reference frame (NAD83) can be misleading.The same minor horizontal movements with respect to NAD83(2011) can also be observed at the 13 reference stations, even though the nearzero velocities with respect to SGOMRF indicate they are stable (Fig.2).The coherent horizontal movement referred to the NAD83(2011) could be incorrectly interpreted as the result of local active faulting.In fact, the 2 mm yr −1 horizontal velocity does not represent any local ground movement.It indicates the low precision of the NAD83 reference frame within the GOM region.The instability of the continentalscale reference frame will overlook or bias minor local horizontal ground deformation signals.The root-mean-square (RMS) value of the linear velocities of these 13 reference stations are 0.15 mm yr −1 for the east component, 0.19 mm yr −1 for the north component, and 0.25 mm yr −1 for vertical component with respect to the local

Figure 3 .
Figure 3. Comparisons of time series at two CGPS sites, OKAN and SG05, referred to three reference frames: IGS08 (black), NAD83 (blue), and SGOMRF (red).OKAN is located at Antlers, Oklahoma.SG05 is located at Melbourne, Florida.

Figure 5 .
Figure 5. Displacement time series of two rapidly subsiding CGPS sites (ROD1 and TXCN) in northern Houston.The displacement time series of a stable site (TXLI) are plotted for comparative purposes.Locations of these three stations are plotted in Fig. 6.The reference frame is SGOMRF.

Figure 6 .
Figure 6.Locations of three CGPS stations and cities within the greater Houston area.Black and blue velocity vectors represent horizontal and vertical velocities with respect to SGOMRF, respectively.Red lines are growth faults and yellow dots are salt domes(Garrity and Soller, 2009).Contours lines are subsidence rate from 2005 to 2012(Kearns et al., 2015).

Figure 7 .
Figure 7. Horizontal (red) and vertical (green) velocity vectors at MMX1 and UNIP with respect to the local reference fame (SGOMRF).The color patterns represent the average subsidence rate derived from InSAR analysis (Chaussard et al., 2014).The three-component positional time series of UNIP and MMXI are plotted in Fig. 10.

Figure 8 .
Figure8.Three-component displacement time series from two CGPS sites with considerable horizontal movements.FSHS is located at Franklin, Louisiana (Fig.9a).TXPR is located at Pharr, Texas (Fig.4).The reference frame is SGOMRF.

Figure 9 .
Figure 9. Top maps show vertical and horizontal velocity vectors in (a) southeastern Louisiana and (b) Houston.Bottom plots show vertical positional time series of subsiding stations in (c) southeastern Louisiana (> 1.5 mm yr −1 ) and (d) Houston-Galveston (> 4.5 mm yr −1 ).The reference frame is SGOMRF.

Figure 10 .
Figure 10.Three-component displacement time series of three CGPS stations in central Mexico.The reference frame is SGOMRF.The locations of UNIP and MMX1 are marked in Fig. 7 and the location of INEG is marked in Fig. 4.

Table 1 .
Fourteen parameters for reference frame transformations from IGS08 to SGOMRF and from IGS08 to NAD83(2011).

Table 2 .
Geodetic coordinates (longitude, latitude, ellipsoidal height) of the 13 reference sites with respect to SGOMRF and their RMS accuracy (repeatability).

Table A1 .
Locations and velocities of 148 GPS stations within the Gulf of Mexico region with respect to SGOMRF.

Table A2 .
Locations and velocities of 13 reference stations used to establish SGOMRF.