IntroductionFor almost 1 century, the Greater Houston region has been impacted by ground deformation associated with subsidence and creeping faults. Subsidence and faulting in urban areas have caused widespread damages to residential, commercial, industrial buildings, and public infrastructure. Damages associated with subsidence and faulting often go unnoticed until a broad area is affected and substantial indirect risks have been induced. Prior to the 1990s, subsidence within the Greater Houston region was measured using leveling surveys and borehole extensometers. GPS technology was implemented into land surveying in the Houston region by the joint efforts of the National Geodetic Survey (NGS) and Harris-Galveston Subsidence District (HGSD) in the late 1980s. At that time, the full GPS constellation (minimum 24 satellites) had not been completed. Nevertheless, Houston is one of the first places that employed GPS technology for urban geological hazards monitoring. GPS gradually has replaced conventional leveling surveying since the 1990s and has become the primary tool for subsidence monitoring in the Greater Houston region.Global Navigation Satellite System (GNSS) has become the standard generic term for satellite navigation systems, including GPS, GLONASS, Galileo, BeiDou, and other regional satellite-navigation systems. The majority of HoustonNet stations record only GPS signals during their entire operational history. A few stations started to record GLONASS, Galileo, and BeiDou signals recently. As of 2021, HoustonNet uses GPS-only signals, specifically L1 and L2 code and phase observations (L1, C1, L2, and P2), to calculate daily static positions. The methodology developed through this study applies to the observations from other satellite systems. Accordingly, we use the umbrella term GNSS throughout this paper.Land subsidence in the Greater Houston region is caused primarily by excessive groundwater withdrawals within shallow aquifers, which consist of interbedded clays, silts, sands, and gravels (e.g., Gabrysch 1982; Kearns et al. 2015). The major aquifers belong to the Gulf Coast Aquifer system (Fig. 1), from the land surface downward, including the Chicot aquifer, the Evangeline aquifer, the Burkeville confining unit, the Jasper aquifer, and the Catahoula sandstone (Baker 1979). The Gulf Coast Aquifer system outcrops in the northwest corner of the Greater Houston region (Casarez 2020) (Fig. 2). Accurate monitoring of ground deformation over a long period and a large area is vital to managing ongoing urban geological hazards and is pertinent to plans for future urban development. Long-term GNSS observations also provide ground deformation truth and calibration data for other remote sensing techniques and subsidence modeling.The Texas Legislature created HGSD in 1975, Fort Bend Subsidence District (FBSD) in 1989, Lone Star Groundwater Conservation District (LSGCD) in 2001, and Brazoria County Groundwater Conservation District (BCGCD) in 2005. HGSD and FBSD were created to regulate groundwater withdrawals to prevent land subsidence in areas within their respective jurisdictions; LSGCD and BCGCD were created to manage groundwater resources within their jurisdictions (Fig. 2). As of 2021, HGSD is divided into three regulatory areas—Area 1, Area 2, and Area 3—and FBSD is divided into two regulatory areas—Area A and Area B. In collaboration with the University of Houston (UH), FBSD, LSGCD, BCGCD, and other local institutes, HGSD has integrated approximately 250 permanent GPS stations within the Greater Houston region into its routine subsidence monitoring (Fig. 2). The combination of all publicly available permanent GNSS stations within Greater Houston is called HoustonNet. The Geodetic Laboratory at UH has been responsible for the HoustonNet data process and analysis since 2016.A modern regional geodetic infrastructure comprises three fundamental components: a network of permanent GNSS stations (hardware), sophisticated software packages for raw data processing (software), and rigorous regional reference frames (firmware) (e.g., Wang et al. 2019, 2020). Firmware refers to the regional reference frame, also known as soft-hardware or embedded software. Establishing a rigorous regional reference frame requires long-term (e.g., >7 years) accumulation of continuous GNSS observations from many permanent GNSS stations in the region. Hardware, software, and firmware work together to build a sophisticated regional geodetic infrastructure for high-accuracy and high-precision GNSS monitoring. HoustonNet data and its products have been used widely by researchers from subsidence districts, the USGS, universities, and private consulting companies for subsidence (e.g., Engelkemeir et al. 2010; Bawden et al. 2012; Qu et al. 2015; Yu and Wang 2016; Thornhill and Keester 2020), faulting (e.g., Khan et al. 2014; Liu et al. 2019; Qu et al. 2019), coastal flooding (e.g., Miller and Shirzaei 2019, 2021), and sea-level change studies (e.g., Liu et al. 2020; Zhou et al. 2021). Data analysis methods have been mentioned sporadically in previous publications. Because of the frequent updates of software packages and reference frames, the details of HoustonNet’s data processing and analysis methods are adjusted and improved continuously. The results from previous processing could be slightly different from the results of new processing. Stringent users may need to understand the details of data analysis methods. This paper summarizes the current GNSS geodetic infrastructure and documents the details of the GNSS data analysis methodology.The Greater Houston region described in this paper covers an area of approximately 26,000 km2 (160×160 km) comprising nine counties: Harris, Fort Bend, Brazoria, Galveston, Montgomery, Waller, Liberty, Chambers, and Austin (Fig. 2). Greater Houston has been among the fastest-growing metropolitan areas in the US since the 2010s. As of 2020, the population of the Greater Houston region is over 7 million, according to the 2020 census estimates (US Census Bureau 2020). Major cities located in this area include Houston, Galveston, Texas City, Dickinson, La Marque, League City, La Porte, Baytown, Pasadena, Spring, The Woodlands, Conroe, Cypress, Jersey Village, Katy, Rosenberg, Richmond, Sugar Land, Missouri City, and Pearland.HoustonNetAs of 2021, HoustonNet comprises approximately 250 permanent GNSS stations (Fig. 2), including 230 active stations and 20 decommissioned stations (no data since 2018). This network is expanding continuously. These GNSS stations are operated by HGSD, FBSD, UH, and several other local agencies. Approximately 25 stations have data spanning more than 20 years, 80 stations have data spanning between 10 and 20 years, and 100 stations have data spanning between 5 and 10 years (Fig. 3). The basic information of these stations is listed in Table S1.HGSD GNSSHGSD started to install permanent GNSS stations for land subsidence monitoring in the early 1990s. The early permanent GNSS stations, known as Port-A-Measure (PAM) stations (Zilkoski et al. 2003), were designed for periodic surveys rather than continuous surveys, to overcome the high cost of GNSS units (receiver plus antenna) at that time. Several PAM stations have been switched to continuous operation, such as P026, P034, P043, P049, P079, P080, P081, and P096. The first group of PAM stations (P001, P002, P003, and P004) were installed in 1994. PAM stations were designed for a campaign-style long-term monitoring solution. The antenna pole is a 5-cm-diameter steel pipe anchored at the bottom of a borehole 10.5 m below ground surface. The top 6 m of the borehole is lined with a PVC (polyvinyl chloride) pipe. Thus, ground deformation resulting from the shrink–swell behavior of shallow clay-rich soils within the top 6 m is excluded from PAM measurements. The antenna height above the land surface is approximately 2.5 m [Fig. 4(a)]. The PAM GNSS network has expanded to 105 permanent stations as of 2021 (Greuter and Petersen 2021). HGSD and FBSD personnel routinely mount GNSS antennas on permanent antenna poles to monitor land subsidence. On average, GNSS data were collected continuously for 1 week every month prior to 2006 and 1 week every 2 months since 2006. The raw data of HGSD stations are archived at HGSD and are open to the public.University of Houston GNSSThe Geodetic Laboratory at UH continuously has been building a permanent GNSS network within the Greater Houston region since 2013 (Wang et al. 2015b). As of 2021, the UH GNSS network comprises approximately 65 permanent stations (Fig. 2 and Table S1). The primary purpose of the UH GNSS is to provide a platform for studying multiple urban natural hazards, such as subsidence, faulting, salt dome uplift, coastal erosion, flooding, hurricanes (monitoring water vapor intensity), and urban heat island effects. Most UH GNSS stations are installed on one- or two-story buildings on the campuses of public schools and colleges [Fig. 4(b)]. Raw GNSS data from these stations are available to the public through the data archiving facility at UNAVCO.Other GNSSThe Texas DOT (TxDOT) operates a real-time GNSS network with 194 permanent GNSS stations (as of 2021) in Texas (TxDOT 2021). TxDOT has been using GNSS technology to support its surveying activities since the early 1980s. TxDOT stations are installed in the free-field [Fig. 4(c)] or on office buildings [Fig. 4(d)]. As of 2021, HoustonNet has integrated 20 TxDOT stations with an observational history of over 2 years into its routine data analysis.SmartNet operates a massive network of over 4,500 permanent GNSS stations worldwide for providing real-time kinematic (RTK) corrections for high-accuracy navigation, mapping, and monitoring applications (SmartNet North America 2021). SmartNet operates over 1,400 permanent GNSS stations in North America. As of 2021, HoustonNet has integrated 25 SmartNet stations with an observational history of over 3 years into its routine data analysis. Several Smartnet stations are operated jointly by SmartNet and UH. Raw data are archived at SmartNet and UNAVCO.The City of Houston, US Coast Guard, Federal Aviation Administration (FAA), National Oceanic and Atmospheric Administration (NOAA), RODS Surveying, and several other agencies operate a few permanent continuous GNSS stations in the Greater Houston region (Table S1).According to Yang et al.’s (2016) investigation of several pairs of closely spaced ground-based and building-based GNSS stations, there is no considerable difference between building-based (one- to two-floor buildings) and ground-based GNSS positions regarding the GNSS-derived site velocities and their uncertainties as long as the buildings are stable. Accordingly, we did not distinguish building-based and ground-based (also called free-field) GNSS data in our data analysis.Data and Data FlowThe GNSS receivers in the field record multifrequency pseudoranges and phase circles, and signal-to-noise (SNR) measurements for each satellite being tracked. HGSD’s PAM data are collected at 1 sample/30 s. The majority of continuous GNSS data of HoustonNet are collected at a rate of 1 sample/15 s. GNSS data generally are stored in receivers as daily files that span one GPS-time day in a vendor-specific proprietary format. HGSD and UH stations are Trimble (Sunnyvale, California) antennas and receivers, mostly Trimble NetR9 receivers and Geodetic Zephyr II antennas. In general, TxDOT stations are equipped with Trimble antennas and receivers, and SmartNet (Norcross, Georgia) stations are equipped with Leica (Norcross, Georgia) antennas and receivers. Raw data from the majority of HGSD stations are downloaded manually and archived at HGSD. UH GNSS stations are downloaded automatically and archived at UH and UNAVCO. Most TxDOT stations in the Greater Houston region participate in the NGS Continuously Operating Reference Station (CORS) program, and raw data are archived at NGS. TxDOT also archives raw data at its FTP site. Most other GNSS stations are archived at NGS, UNAVCO, or their own archiving facilities. In addition, a few HoustonNet stations stream 1-Hz data in real-time for high-accuracy land surveying and other real-time applications. Real-time GNSS data and their applications are not addressed in this paper.Fig. 5 illustrates the HoustonNet data flow schema. The data flow consists of the following primary components: (1) collection and transfer of raw GPS data (e.g., *.T02) from field instruments to its data center, (2) generation of receiver independent exchange (RINEX) data, including file translation, quality checks, archiving, distribution, and metadata management at UH, (3) generation of Precise Point Positioning (PPP) daily solutions referred to the current International GNSS Service (IGS) Reference Frame, and (4) transforming the daily solutions to the regional reference frames (Houston20 and GOM20) and generating final data products, including position time series, site velocities and their uncertainties, and contour maps. The raw data and RINEX files are archived at UH. TEQC software is used for conducting data quality-check and generating RINEX (version 2.11) files (Estey and Meertens 1999). HoustonNet data products currently are distributed to the public via the Web services at HGSD (2021).Analysis MethodsThe major steps for GNSS data (daily RINEX files) postprocessing include (1) producing Earth-Centered-Earth-Fixed (ECEF) Cartesian coordinates (XYZ) from RINEX files, (2) removing outliers and correcting steps, (3) transforming the ECEF-XYZ coordinates from the global to regional reference frames, (4) converting ECEF-XYZ coordinates to East-North-Up (ENU) coordinates, and (5) calculating site velocities and their uncertainties.Producing Daily Solutions with Respect to the IGS Reference FrameGNSS data processing algorithms generally implement two approaches to achieving high-precision GNSS positions: relative and absolute positioning (e.g., Herring et al. 2016; Wang et al. 2017). A relative positioning method uses simultaneous observations from two or more GNSS units; at least one of these antennas is fixed at a known location with respect to a specific reference frame. The position of a rover station can be determined relative to the fixed station by applying a carrier-phase double difference method. The relative method also is called the differential method. In contrast to the differential method, the absolute method involves only one GPS station to determine its coordinates with respect to a global reference frame. The PPP method is a typical absolute positioning method. The methodology and algorithms of PPP were described by Zumberge et al. (1997). The PPP method requires only a single GPS station at the user’s end. Users do not need to install any reference stations in the field, nor do they need to include any data from reference stations during post processing. PPP techniques, combined with stable regional reference frames, have attracted broad interest in geological hazards monitoring (e.g., Wang 2013; Murray and Svarc 2017; Zhao et al. 2020) and structural health monitoring (SHM) (e.g., Bao et al. 2018) because of the operational simplicity and the consistency of the positioning accuracy over time and space.The differential method was used in the early study of land subsidence in the Houston region by HGSD (e.g., Zilkoski et al. 2003) and the geodesy research community (e.g., Wang and Soler 2013). Before 2015, HGSD’s GNSS data were processed primarily by the software package Program for Adjustment of GPS Ephemerides (PAGES) developed by NGS (Schenewerk and Hilla 1999). PAGES also is the core GNSS data processing engine of the Online Positioning User Service (OPUS), which solves a static position by averaging three separate differential solutions processed using three selected reference stations (e.g., Wang and Soler 2012). The disadvantage of the differential method is that the stability of the reference stations needs to be justified independently and precisely, which often is tricky in practice (e.g., Guo et al. 2019). Subsidence at rover sites will be underestimated if the reference station also is subsiding. It has become more and more difficult to find stable reference stations since the 1990s because of the spread of subsidence in the Greater Houston region. Thus, it is challenging to obtain accurate and reliable subsidence measurements using the differential method within the Greater Houston region. Another disadvantage of the differential method is that the accuracy of the relative positioning depends on the length of the baseline, i.e., the distance between reference and rovers. Thus, relative positioning could result in inconsistent accuracy of positions within the Greater Houston region (Soler and Wang 2016).HGSD changed its data processing strategy from using the relative positioning to using the absolute positioning method in 2016. Jet Propulsion Laboratory’s (JPL) GipsyX/RTGx software package, previously GIPSY/OASIS, has been employed in HoustonNet data processing since 2016. As of October of 2021, HoustonNet uses the single-receiver phase-ambiguity-fixed PPP in GipsyX (v.1.7) to generate ECEF-XYZ coordinates with respect to the International GNSS Service Reference Frame 2014 (IGS14), the IGS realization of the International Terrestrial Reference Frame 2014 (ITRF14) (Altamimi et al. 2016; Rebischung et al. 2016). IGS replaced IGS14 with IGb14 on May 17, 2020 (Rebischung 2020), the realization of IGb14 used 5-year longer data sets than the realization of IGS14. As of 2021, JPL’s products are aligned with IGb14. However, most geodesy publications continue to use the term IGS14 rather than IGb14. For this reason, we use the term IGS14 throughout this paper.The GipsyX processing uses JPL’s repro3.0 IGS14 final no-net rotation (NNR) products to solve an initial position for RINEX files. The NNR products contain files aligned to the NNR reference frame, including orbital state estimates, transmitter clock estimates, Earth rotation parameter (ERP) estimates, spacecraft attitude information, and wide-lane phase biases (WLPB) information. A seven-parameter Helmert transformation is performed using the x-file provided in the JPL products to transform the NNR position into the current IGS reference frame. The x-file contains the seven Helmert parameters (three small rotations, three translations, and one scale parameter of the day) needed to transform the NNR frame to the IGS reference frame. To set up the GipsyX-PPP processing, HoustonNet uses a parameter tree very similar to the default tree included in the GipsyX-1.7 release. To account for ocean tide loading effects in daily positions, HoustonNet uses the ocean loading displacements coefficients calculated by the ocean tide loading Web service operated by Chalmers University of Technology, Sweden, using the Finite Element Solution ocean tide model FES2004 (Bos and Scherneck 2021). To account for ionospheric effects, HoustonNet uses Ionosphere Exchange (IONEX) files (Schaer et al. 1998) for second-order ionospheric corrections. To correct for tropospheric effects, HoustonNet uses the GPT2w-based nominal troposphere mapping functions (Böhm et al. 2015). To correct for antenna changes, HoustonNet uses calibrations made available by the IGS14 ANTEX calibration table (Schmid et al. 2007).The PPP processing results in 24-h average positions (daily solutions) from the daily RINEX files. Each daily file from each station is processed independently. Thus, the coordinate results are insulated from potential data problems at other days or other sites. The accuracy of GNSS static positioning has improved significantly during the last 3 decades with the continuous improvement of global and regional geodetic infrastructures. Bertiger et al. (2010, 2020) reported that the PPP resolutions achieve daily RMS accuracy of approximately 2 mm in the horizontal direction and 6–7 mm in the vertical direction. In general, the RMS accuracy of PPP solutions in the Greater Houston region is slightly worse than the global average because of the humid climate and the fluctuations of groundwater levels. According to Wang et al. (2013, 2017) and this study, the RMS accuracy of daily PPP solutions in the Houston region is about 2–4 mm in the horizontal direction and 5–8 mm in the vertical direction.For GNSS stations for which the completed RINEX data sets are not archived at HoustonNet, the IGS14-XYZ solutions provided by the Nevada Geodetic Laboratory (NGL) at the University of Nevada, Reno are used (Blewitt et al. 2018). NGL also employs the GipsyX software for their routine GNSS data processing. HoustonNet periodically compares the daily solutions with NGL’s daily solutions. The IGS14 XYZ time series produced by UH and NGL are the same.Removing OutliersThe PPP solutions are defined in the ECEF-XYZ coordinate system referring to a global reference frame, currently IGS14. Although the accuracy and precision (repeatability) of GNSS solutions have been improving over time, certain large variations (anomalies) often are superimposed onto the daily coordinate time series. These obvious anomalies, often referred to as outliers, do not reflect actual physical motions of the antenna. Removing outliers is necessary before further analysis; otherwise, one may obtain unrealistic site velocities and significant uncertainties. The effect of outliers on site-velocity estimates could be significant for short-period data sets. In our analysis, the outliers were identified and removed in the IGS14-XYZ time series before conducting reference frame transformation and calculating the ENU position time series. The causes of outliers are various. We found that the primary reason is the short period of observations in the field. RINEX files comprising a few hours of observations often present as outliers in the 24-h position time series. Power outages and equipment maintenance in the field, flooding, and extreme weather conditions often are the causes of outliers in the daily position time series. Weather fronts and heavy rains are typical weather conditions that cause outliers in GNSS time series (e.g., Wang 2013).There are numerous outlier detection algorithms for GNSS time series in literature (e.g., Gökalp et al. 2008; Ordoñez et al. 2011; Wang 2011). In general, most methods work well for time series following a linear trend. However, most GNSS stations in the Greater Houston region are affected by nonlinear subsidence. Furthermore, HoustonNet comprises both campaign and continuous data. The campaign data sets are much noisier than the continuous data sets. In practice, it is challenging to apply a unified outlier removing method to both continuous and campaign data sets.The static PPP provides the x-, y-, and z-coordinates and their uncertainties. The corresponding uncertainties of 24-h XYZ coordinates usually are at a level of a few millimeters. Occasionally, the uncertainty can be tremendous, up to several centimeters, even a few decimeters, which indicates that the carrier-phase ambiguities were not fixed successfully to their correct integer number. Therefore, an initial screening is conducted to remove those epochs with substantial uncertainties in the x-, y-, and z-components. The threshold of the uncertainty is set to 0.2 m. If the uncertainty in one component is larger than 0.2 m, this epoch is removed from all three components. For RINEX files from continuous GNSS stations, the initial screening rarely identifies outliers. For HGSD’s campaign data, this initial screening can identify certain outliers occasionally. In this section, we use the data set from P026 (2002–2021) to illustrate the method for identifying and removing outliers. P026 is a PAM GNSS station operated as a campaign station before 2017 and as a continuous station since 2017. The location of P026 is marked in Fig. 2. P026 is an extreme case in which the time series comprises mixed campaign and continuous data with a sharp curve associated with the severe drought from 2011 to 2014. The initial screening removed 32 outliers from the original 2,491 measurements.Fig. 6(a) illustrates the steps for removing outliers at P026 after the initial screening. The first step is to smooth the XYZ time series with locally weighted scatterplot smoothing (LOWESS) to obtain the trend of the time series. LOWESS uses a nonparametric technique to fit a smooth curve through points in a time series (Cleveland 1981). A smoothing parameter (F) ranging from 0 to 1 is designed to determine the proportion of observations for local regression. The value of F controls the smoothness of the curve. A larger value of F results in a smoother curve, and a smaller F captures more short-period signals. The value of F was set as 0.1 in our analysis to catch sharp curves in time series. A larger F may result in too many outliers in a curved time series. However, for time series that follow a linear trend, the selection of F rarely affects the number of identified outliers. The LOWESS algorithm employs an iterative process to smooth the time series. The number of iterations was set to 2 in our processing. The second step is to obtain the residual time series by removing the trend time series from the original time series. The third step is to calculate the median of absolute deviations (MAD) of the residual time series and identify and remove outliers. MAD is a robust measure of the variability of a time series which is more resilient to outliers than the standard deviation. The criteria for identifying outliers are discussed subsequently. When an outlier is identified from one component (x, y, or z), the coordinates of this day are removed from all three components, although the coordinates in other components may be normal. This explains why more outliers were removed from the y-component in Step 3 than the outliers identified in Step 2.Assuming that the residuals after removing the trend roughly fit a Gaussian probability density function, the standard deviation (σ) of the residual time series can be estimated robustly by the following equation (Wilcox 2013, p. 75): (1) According to the three-sigma rule for a Gaussian distribution, also called the 68–95–99.7 rule, only about 0.3% of observations fall beyond the first three standard deviations μ±3σ, where μ denotes the mean of the residual time series. Accordingly, in our routine data processing, we set the threshold for identifying outliers to 4.5 times the MAD, which is approximately 3σ. Of course, GNSS time series are not stationary, and the residuals do not exactly follow a Gaussian distribution. In practice, the 3σ threshold often results in more outliers than 0.3% of the observations.Fig. 6(b) illustrates the final displacement time series and identified outliers at P026 in the north–south (NS), east–west (EW), and up–down (UD) directions. Forty-seven measurements among 2,459 total measurements (approximately 2%) were identified as outliers. A cursory judgment suggests that the ±3σ (or ±4.5 MAD) threshold does a reasonable job of cleaning the GNSS time series. In geophysics, it is well known that one person’s noise may be another person’s useful signal. Accordingly, we removed only the exceptional outliers in our routine processing. For HoustonNet data products, approximately 1%–3% of the total samples are removed as outliers in the routine data process. HoustonNet also provides whole data (XYZ and ENU time series) without removing any outliers. Users may consider removing more or fewer outliers for their specific applications.Reference Frame TransformationA global geodetic reference frame is realized with an approach of minimizing the overall movements of a group of selected reference stations distributed worldwide. In practice, site velocities with respect to a global reference frame are difficult to interpret visually from a regional- or local-scale geophysical perspective. In the Greater Houston region, the changes of the ENU positions with respect to IGS14 are dominated by the long-term drift and rotation of the North American plate (Fig. 1). Stable sites within the Greater Houston region generally retain approximately 14 mm/year horizontal movement in the southwest direction and 2 mm/year downward movement with respect to the global reference frame. Localized and temporal ground deformation, such as subsidence and fault creeping, could be obscured or biased by the regional motions. A stable regional- or local-scale reference frame is needed to exclude those common ground movements and highlight localized ground deformation.In routine processing, HoustonNet transforms the IGS14-XYZ coordinates to the Stable Houston Reference Frame and the Stable Gulf of Mexico (GOM) Reference Frame. The initial version of the Stable Houston Reference Frame was realized in 2013 by 10 stations located outside of the Greater Houston region (Wang et al. 2013). The initial reference frame was tied to the International GNSS Service Reference Frame 2008 (IGS08). The Stable Houston Reference Frame was updated several times with more reference stations and longer observational time spans (Wang et al. 2015b; Kearns et al. 2019). The most recent update was conducted in 2020 using 25 continuous reference stations with a minimum observational history of 8 years (Agudelo et al. 2020). The updated Stable Houston Reference Frame is designated Houston20. The Stable Gulf of Mexico Reference Frame initially was established in 2016 (Yu and Wang 2017) and updated in 2020, designated GOM20 (Wang et al. 2020). GOM20 was realized using 55 continuous GNSS stations located on the stable portion of the Gulf Coastal Plain. The average observational period of these reference stations was 13.5 years. Locations of reference stations used for realizing Houston20 and GOM20 are depicted in Fig. 1. Both reference frames currently are tied to IGS14 and will be improved incrementally and synchronized with the update of the IGS reference frame.The detailed methods for realizing a stable regional reference frame and the criteria for selecting reference stations were addressed in previous publications (e.g., Kearns et al. 2019; Wang et al. 2020). The ECEF-XYZ coordinates with respect to IGS14 are transformed to Houston20 or GOM20 according to the following equations: (2) [X(t)Y(t)Z(t)]RegionalRF=[TX′·(t−t0)TY′·(t−t0)TZ′·(t−t0)]+[1RZ′·(t−t0)−RY′·(t−t0)−RZ′·(t−t0)1RX′·(t−t0)RY′·(t−t0)−RX′·(t−t0)1]×[X(t)Y(t)Z(t)]IGS14where X(t)IGS14, Y(t)IGS14, and Z(t)IGS14 = geocentric XYZ coordinates (at epoch t) of a site with respect to IGS14; X(t)R, Y(t)R, and Z(t)R = XYZ coordinates of site with respect to Houston20 or GOM20 at epoch t; X(t)IGS14, Y(t)IGS14, and Z(t)IGS14 = IGS14 positions obtained from the PPP solutions (m); and t0 = epoch that aligns coordinates with respect to the two reference frames, where t0 is fixed at 2016.0 (year) for Houston20 and at 2015.0 for GOM20, and a site retains identical XYZ coordinates at epoch t0 with respect to IGS14 and Houston20 or GOM20; Tx′, Ty′, and Tz′ are constant parameters indicating rates (one-time derivatives) of three translational shifts along x, y, z-coordinate axes; and Rx′, Ry′, and Rz′ = rates of three rotations between two reference frames around the x, y, z-coordinate axes. Counterclockwise rotations are regarded as positive. The frame stability of Houston20 and GOM20 is at a submillimeter-per-year level in all three directions (Agudelo et al. 2020; Wang et al. 2020). These seven parameters t0, Tx′, Ty′, Tz′, Rx′, Ry′, and Rz′ are listed in Table 1.Table 1. Seven parameters for realizing Houston20 and GOM20Table 1. Seven parameters for realizing Houston20 and GOM20ParameteraIGS14 to Houston20IGS14 to GOM20t0 (years)2,016.02,015.0Tx′ (m/year)1.4040400×10−0027.1281610×10−004Ty′ (m/year)9.6139040×10−0045.6136741×10−004Tz′ (m/year)7.2404862×10−0032.9287337×10−003Rx′ (rad/year)−9.8590126×10−010−4.0941604×10−010Ry′ (rad/year)−1.7311089×10−009−3.1975595×10−009Rz′ (rad/year)1.3205311×10−009−2.3610546×10−010To study ground deformation at the Earth’s surface, the geocentric XYZ coordinates are converted to a geodetic orthogonal curvilinear coordinate system (longitude, latitude, and ellipsoid height) referenced to the GRS80 ellipsoid. The conversion from XYZ coordinates to geodetic longitude (λ) is straightforward. However, the conversions for the latitude (φ) and ellipsoid height (h) involve complicated calculations. The computation must be done iteratively, and it is sensitive to high accuracy (small errors) because the magnitudes of radius of the ellipsoid and an ellipsoid height (h) are about 106 apart. We used the algorithm introduced by Heiskanen and Moritz (1967, Chapter 5.3) to calculate latitudes and ellipsoid heights.The longitude and latitude coordinates then are projected onto a two-dimensional (2D) local horizontal plane for tracking surface ground deformation in the north–south and east–west directions at each site. The change of ellipsoid heights over time is used to depict the land surface displacement in the up–down direction. The following equations are used to obtain the site-specific topocentric ENU displacement time series: (3) [E(t)N(t)U(t)]=[−sinλ0cosλ00−cosλ0·sinφ0−sinλ0·sinφ0cosφ0cosλ0·cosφ0sinλ0·cosφ0sinφ0]×[X(t)−X0Y(t)−Y0Z(t)−Z0]where λ0 and φ0 = initial geodetic longitude and latitude of site corresponding to initial XYZ coordinate (X0,Y0,Z0), which is the position at the first epoch of the time series.The detection of land subsidence historically has depended on periodical leveling surveys of benchmarks. This traditionally has been accomplished by differencing orthometric heights obtained from spirit leveling. According to the investigation of Wang and Soler (2014), for practical use, the vertical displacements derived from ellipsoid heights are the same as those derived from orthometric heights. In HoustonNet data products, negative vertical displacements indicate subsidence, and positive vertical displacements indicate uplift. Many GNSS sites in southeastern Houston had minor land uplift over a decadal period after Chicot and Evangeline groundwater levels recovered to their preconsolidation heads (Kearns et al. 2015).Fig. 7 compares the ENU series (2005–2020) with Houston20 and GOM20. TXLM (Houston), LMCN (New Orleans), and ZMA1 (Miami) are three long-history GNSS stations located inside and outside the Greater Houston region (Fig. 1). TXLM is located in La Marque, Galveston County, south of Houston (Fig. 2). The average subsidence rate at this site during the last 16 years was 2.9±0.3 mm/year with respect to Houston20 and 3.2±0.3 mm/year with respect to GOM20. The difference is 0.3 mm/year, which is at the level of the 95% confidence interval (CI) of the velocity estimates. The difference of horizontal site velocities with respect to Houston20 and GOM20 is below the 95% CI (0.2 mm/year). We checked many other stations within the Greater Houston region. The difference of site velocities with respect to Houston20 and GOM20 is below 1 mm/year within the Greater Houston region.LMCN is a long-period GNSS located on top of a Louisiana University Marine Consortium building in Cocodrie, Louisiana, about 150 km southwest of New Orleans and about 500 km from Houston (Fig. 1). GNSS-derived subsidence at LMCN has been applied frequently to delineate long-term coastal subsidence and absolute sea-level rise rates along the Mississippi delta coast (e.g., Kuchar et al. 2018; Keogh and Törnqvist 2019; Wang et al. 2020). This site has a subsidence rate of 4.2±0.3 mm/year with respect to Houston20 and 5.8±0.3 mm/year with respect to GOM20. The difference of 1.6 mm/year is significant for sea-level studies. Using the coastal subsidence rate with respect to Houston20 will exaggerate the absolute-sea-level rise rate estimate at this site. ZMA1 is located in Miami, which is about 1,500 km from Houston (Fig. 1). ZMA1 has a near-zero (<0.5 mm/year) site velocity with respect to GOM20 in all three directions, which is consistent with the geological understanding that the Florida peninsula is a stable portion of the Gulf Coastal Plain. The vertical velocity at ZMA1 with respect to Houston20 is 3.2±0.2 mm/year, which is unrealistically large. The horizontal velocity of 1.1±0.2 mm/year (NS) with respect to Houston20 also is remarkable. These unphysically large velocities with respect to Houston20 are caused by the rotations associated with the coordinate transformation when sites are far from the area covered by the reference stations. The effect of rotations is investigated in the next section.Locating Euler PoleThe movement of a rigid body on the surface of the Earth with respect to a global reference frame can be described as a rotation around a fixed pole that passes through the center of the Earth. This pole of rotation is known as an Euler pole. The coordinate transformation from a global reference frame to a regional reference frame comprises both translations and rotations [Eq. (2)]. The translations often are at a level of a few millimeters per year (Table 1). Therefore, the coordinate transformation from a global reference frame to a regional reference frame can be approached by applying a rotation against the Euler-pole rotation. Accordingly, a regional reference frame also can be described as an Euler-pole reference frame. The rotation is called anti-Euler-pole rotation. The location of the Euler pole location and the rotation rate can be estimated according to the three rotation rates Rx′, Ry′, and Rz′ used to realize the regional reference frame (Table 1) (4) (5) φp=tan−1(Rz′Rx′2+Ry′2)(6) where λp (longitude) and φp (latitude) = location of Euler pole (rad); and ω = rotation rate (rad/year). A positive latitude estimate (φp) indicates north of the equator, and a negative latitude estimate indicates south of the equator. A positive longitude estimate (λp) indicates east of the prime meridian, and a negative longitude estimate indicates west of the prime meridian. Eq. (4) results in a longitude (λp) between −π/2 and π/2. However, the real longitude of the Euler pole could be −π+λp if λp is positive, and could be π+λp if λp is negative. Fortunately, it is not difficult to determine the correct longitude (λp, −π+λp, or π+λp) based on the geographic locations of reference stations and their rotation directions. In practice, a counterclockwise rotation is represented by a positive ω, and a clockwise rotation is represented by a negative ω.The three rotational parameters (Rx′, Ry′, Rz′) of GOM20 result in an Euler pole (with respect to IGS14) at 97.3°W and 4.2°S with a rotation rate of 0. 185°/million years around the counterclockwise direction. The Euler pole location and rotation rate roughly are comparable with the Euler poles of the North American Plate published by previous researchers. For example, Blewitt et al. (2013) estimated that the plate motion aligned with the North America Reference Frame 2012 (NA12) (with respect to IGS08) is a counterclockwise rotation of 0.177°/million years around the Euler pole at 88.6°W and 9.9°S; Ding et al. (2019) estimated the North American plate motion (with respect to IGS08) to be a counterclockwise rotation of 0.194°/million years around an Euler pole at 85.5°W and 5.0°S; Kreemer et al. (2018) estimated the North American plate motion (with respect to IGS08) to be a counterclockwise rotation of 0.201°/million years around an Euler pole at 86.0°W and 2.3°S. The coverage area of GOM20 is much smaller than the coverage of these continental-scale reference frames. The three rotational parameters (Rx′, Ry′, and Rz′) of Houston20 result in an Euler pole (with respect to IGS14) at 119.7°W and 33.5°N with a rotation rate of 0.137°/million years, which is far from the GOM20 Euler pole. The method for realizing a stable regional reference frame is to find the six best parameters (Table 1) to minimize the overall motions among selected reference stations with respect to the new reference frame. Those parameters do not necessarily reflect the physical motions (shifts and rotations) of the regional crustal block in a precise way. We found that slightly different site velocities can result in considerably different Euler poles. It also is true that very different Euler poles can result in quite similar velocities in certain areas. Only large tectonic plates have well-determined Euler poles. Accordingly, the estimated Euler pole vector (location and rotation) associated with a regional reference frame, particularly a local-scale reference frame, should be interpreted with caution.The horizontal ground movement velocity (v) at a site resulting from the rotation around the Euler pole can be estimated by (7) where ω = rotation rate of Euler pole vector (°/year); R = radius of the Earth (m); and θ = distance from site to Euler pole along the surface of the Earth, measured in terms of angular length (rad). The site velocity increases with the distance to the Euler pole (i.e., θ changes from 0 to π/2).For stable sites within the area covered by the reference network, the ground movements associated with the anti-Euler-pole rotation would offset the ground movement with respect to IGS14. Thus, stable sites would retain near-zero (e.g., <1 mm/year) site velocities with respect to the regional reference frame. However, ground movements produced by the anti-Euler-pole rotation could be significantly large or small if a site is too far from or too close to the Euler pole [Eq. (7)]. As a result, a physically stable site far from the area covered by the network of reference stations may not retain a near-zero site velocity with respect to the regional reference frame. In this case, the reference frame transformation makes no sense in practice. Accordingly, we suggest that the application scope of a regional reference frame should be limited to the area covered by the network of reference stations. We recommend that users targeting ground deformation within the Greater Houston region should use Houston20; users targeting regional studies, such as delineating coastal erosion, creeping of growth faults, subsidence, and sea-level changes along the GOM coast, should use GOM20. HoustonNet data products provide the ENU time series with respect to three reference frames: IGS14, GOM20, and Houston20. Users are obliged to select the correct time series according to their specific research purposes.Detecting StepsA well-known issue associated with the GNSS-derived displacement time series is the presence of instant position changes, known as steps or breaks (e.g., Williams 2008; Gazeaux et al. 2013; Griffiths and Ray 2016; Heflin et al. 2020). For the displacement time series in the Greater Houston region, most steps are associated with equipment changes (e.g., antenna, receiver, cable, firmware), severe drought, flooding, and temporal multipaths. These steps need to be identified and treated carefully before performing linear regression and seasonal analysis. Depending on their locations in the time series, undetected, and therefore uncorrected, steps may have a detrimental effect on linear regression. Although HoustonNet tries to maintain an accurate field log for equipment changes and field maintenance, there still are certain undocumented antenna and receiver changes. In practice, it is not easy to perfectly reoccupy a site in the field, even using the same antenna. Coordinate shifts up to a few centimeters can occur for campaign GNSS surveys involving different personnel, antennas and receivers, and antenna mounts. This explains why the displacement time series from PAM stations are considerably noisier than those from continuous GNSS stations (Figs. 6 and 8). Severe drought and flooding events also could cause step-like ground deformation in the vertical direction (e.g., Zhou et al. 2021).We used a change point detection (CPD) program to detect potential steps in the displacement time of each component of each station. The automated CPD is written in Python and available to the public through GitHub (Wang 2021). A change point is an abrupt change in the distributional properties of data. Over the years, numerous CPD algorithms have been proposed to identify abrupt changes in time series (e.g., Guo et al. 2019; Militino et al. 2020). In general, CPD methods comprise three elements: a cost function, a search method, and a constraint on the number of changes to detect. There are two types of CPD algorithms, offline and online (or real-time). The offline algorithms use the whole data to find the change points. In our routine process, we used the pruned exact linear time (PELT) algorithm (Jackson et al. 2005; Killick et al. 2012; Truong et al. 2020), which detects changepoints by minimizing a cost function over possible numbers and locations of change points.Fig. 8(a) depicts the CPD-detected steps in the ENU time series of JGS2. The steps are represented by vertical lines. JGS2 is a permanent GNSS station operated by SmartNet. It is located in Dayton, Liberty County (Fig. 2). This station was installed in October 2010 with a Leica GRX1200 GG Pro receiver and a Leica S10 antenna [Fig. 4(e)]. There was a receiver and antenna change on February 24, 2017. The receiver was changed to a Leica GR30 receiver, and the antenna was changed to a Leica AR10 antenna. The receiver and antenna change on February 24, 2017, was considered in the GipsyX processing with corresponding receiver and antenna parameters. The CPD program found a change point on this day in the vertical direction. The height of the step is 5 mm, which is the typical error associated with the reoccupation of GNSS antennas. The automated program detected another step on November 18, 2015, in the vertical direction. This step also can be identified visually, but it does not coincide with any known cause. The CPD detected two steps in the NS component in January 2016 and January 2019. The causes of these steps are unknown.Fig. 8(b) depicts the CPD-detected steps in the ENU time series of P043, a PAM station located at the southwestern end of Galveston Island (Fig. 2). P043 was a campaign station with data collection of 1 week/month on average before 2017. This station was switched to a continuous station in 2017. There were four change points in the vertical component, on December 5, 2010, July 28, 2013, March 20, 2015, and October 1, 2018. There were no receiver or antenna changes on these four days. The last change point is a minor step (<2 mm) that barely can be verified visually. The first three steps were related to the prolonged drought in southern Texas from fall 2010 to spring 2015 (Zhou et al. 2021). The GNSS-derived site velocity has returned to the predrought level (i.e., before 2011) since 2016.Numerous CPD techniques have been developed in the fields of data mining, statistics, digital signal process, and computer science. In general, most algorithms can identify precisely obvious steps (e.g., >1 cm), but often miss minor steps or result in many false steps. A visual check always is necessary to verify real steps, particularly minor steps (e.g., <5 mm). Nevertheless, an automated CPD program helps to identify minor steps and determine precisely the exact initial epochs (days) of these steps. The CPD program (Wang 2021) does have a function to remove those identified steps. However, HoustonNet does not activate the automated step removal in its routine data process. Obvious steps related to equipment changes are corrected; steps related to real ground movements or unknown reasons are only marked. Stringent users may pay attention to those steps in their own analysis.Seasonal ModelingSeasonal oscillations, particularly in the vertical direction, have been observed widely from GNSS-derived displacement time series (e.g., Dong et al. 2002). In practice, seasonal signals could affect site velocity estimates considerably, which in turn could affect site stability assessments, particularly for observations spanning less than 3 years (e.g., Wang et al. 2017; Wang 2022). A seasonal model often is established to model the seasonal motions associated with seasonal changes of temperature, terrestrial hydrosphere, atmospheric pressure, soil moisture, groundwater pumping, and modeling errors. In our data analysis, each ENU time series was decomposed into four components (Fig. 9) (8) D(i)=L(i)+NL(i)+S(i)+R(i)where L(i) = linear trend obtained by applying a linear regression over the whole time series; NL(i) = nonlinear trend obtained by applying LOWESS filtering to the de-linear-trended time series; S(i) = seasonal motion; and R(i) = residuals after removing the linear trend, the nonlinear trend, and the seasonal motion. The seasonal motion in each direction can be modeled by a combination of annual and semiannual sinusoids with constant amplitudes and phases (e.g., Bao et al. 2021) (9) S(i)=c1cos(2π×ti)+d1sin(2π×ti)+c2cos(4π×ti)+d2sin(4π×ti)where ti = time at epoch i [decimal years (e.g., 2,013.6433)]; and c1, d1, c2, and d2 are coefficients determining amplitudes of annual and semiannual signals. The amplitude of the annual signal can be estimated as p1=c12+d12 ; the amplitude of the semiannual signal can be estimated as p2=c22+d22. The peak-to-trough amplitude of the seasonal motions can be estimated as P=2×p12+p22.Long-term GNSS observations provided by HoustonNet have accumulated rich data sets for studying the seasonal motions in the Houston region. In general, the amplitudes of seasonal motions (up and down) at stable sites within the Greater Houston region are smaller than the amplitudes of seasonal motions in the other regions, such as south Alaska (Wang et al. 2015a) and north China (Wang et al. 2018). Fig. 10(a) depicts the seasonal motions at four sites (TXLI, UHCL, TXLM, TXGA) in the nonsubsiding areas (<2 mm/year). Locations of these four sites are marked in Fig. 2. TXLI is located in Liberty County, one of the most rural Texas counties, with no heavy groundwater pumping during the last century, and in which groundwater levels have been fairly stable. The GNSS-derived site velocity was 0.6±0.5 mm/year from 2015 to 2020. The peak-to-trough amplitude of the seasonal subsidence and uplift is 2 mm, much smaller than the RMS (5.5 mm) of the residuals. UHCL is located on the UH Clear Lake campus in Pasadena, Harris County. TXLM is located in the City of La Marque, Galveston County. The Clear Lake and La Marque areas experienced rapid subsidence during the 1960s and 1970s because of abundant groundwater withdrawals (e.g., Gabrysch 1982). Subsidence slowed beginning in the 1980s and finally has ceased (<2 mm/year) since the 2000s as a result of the groundwater regulations enforced by HGSD. The peak-to-trough amplitudes of seasonal subsidence and uplift at these two sites also are minor (<4 mm). TXGA is located on Galveston Island. According to Zhou et al. (2021), the ongoing subsidence of approximately 2.0 mm/year [2015–2020 (GOM20)] is dominated by the natural compaction of deep aquifers (Evangeline and Jasper). The peak-to-trough amplitude of the seasonal subsidence and uplift at this site is 1.6 mm, which is almost invisible from the time series. The warm climate throughout the year and the thick unconsolidated sediments may result in minor seasonal vertical ground movements at these sites, at which the seasonal fluctuations of groundwater levels are also minor.Fig. 10(b) depicts seasonal motions at four sites (ROD1, CFHS, MRHK, UHCR) in rapidly subsiding (>10 mm/year) areas. The locations of these four sites are marked in Fig. 2. ROD1 is located in Spring, northern Harris County. The average subsidence rate at this site was 11 mm/year from 2006 to 2020 and 6 mm/year from 2015 to 2020. CFHS is located in Cypress, northwestern Harris County. The subsidence rate at this site was 15 mm/year from 2015 to 2020. MRHK and UHCR are located in Katy, Fort Bend County. The subsidence rate was 17 mm/year at MRHK (2015–2020) and 12 mm/year at UHCR (2015–2020). The peak-to-trough amplitudes of seasonal motions at these four sites are outstanding, approximately 1 cm at ROD1, CFHS and MRHK, and 2 cm at UHCR. In addition, we checked more stations in both nonsubsiding and subsiding areas. The seasonal motions are site-specific and dominated by the seasonal fluctuation of groundwater levels in the Chicot and Evangeline aquifers. Accordingly, we do not intend to develop a unified model to characterize the seasonal ground motions within the Greater Houston region. Instead, a seasonal model [Eq. (9)] was determined for each component of each site for calculating the uncertainty of the linear trend according to the method documented by Wang (2022).Calculating Site Velocities and Their UncertaintiesIn the field of geological hazards monitoring, users mostly rely on GNSS-derived site velocities for site stability assessment and risk analysis, rather than on the positions or displacements. The uncertainties of site velocities have become a critical parameter for geological hazards monitoring. In general, accurate positions do not guarantee reliable velocity estimates. The accuracy (uncertainty) of site velocities does not depend entirely on the hardware (antennas and receivers) used for collecting data, but largely relies on the length of the time span and the rigor of the regional reference frame.The conventional approach for estimating a site velocity from the GNSS-derived displacement time series is to apply a linear regression for the whole time series. GNSS-derived vertical displacement time series in the Greater Houston region are affected considerably by long-term groundwater level changes. As a result, long-term displacement time series often have obvious nonlinear features (Fig. 11). In practice, a linear regression is a robust tool for assessing the magnitude of land subsidence over time. Therefore, linear trends derived from the ENU time series have become the key products from HoustonNet. A linear trend and its 95% CI are calculated for the entire time series after removing obvious outliers and correcting steps. The detailed method for calculating a linear trend and its 95% CI from GNSS-derived daily displacement time series was documented by Wang (2022).Fig. 11 depicts the linear regression within 5-year windows at P001 and ROD1. The locations of P001 and ROD1 are marked in Fig. 2. P001 was the first PAM GNSS station installed in Jersey Village, northwestern Harris County. ROD1 is a CORS located in Spring, northern Harris County. Both sites have pronounced nonlinear features of subsidence over time. The subsidence time series within a 5-year time window can be modeled rationally by a linear regression line. The 95% CI of the subsidence rate derived from a 5-year time series varies from about 0.5 to over 1.0 mm/year. Campaign data result in larger site-velocity uncertainties than do continuous data.Major Products and ApplicationsThe major products provided by HoustonNet comprise the ENU time series and site velocities with respect to Houston20, GOM20, and IGS14. The displacement time series with respect to Houston20 provides a fundamental data set for delineating spatial and temporal variations of subsidence and faulting within the Greater Houston region.Tracking Urban SubsidenceAs the population has grown outward from Houston to the north and the northwest since the 1990s, groundwater pumping, and in turn land subsidence, has followed the same pattern of urban expansion (e.g., Coplin and Galloway 1999; Turco and Petrov 2015). GNSS-derived subsidence has provided first-hand data sets for HGSD to assess present groundwater regulations and plan for future groundwater regulation (Petersen et al. 2020; Greuter et al. 2021). Fig. 12 depicts subsidence contour lines during the last two decadal periods, from 2001 to 2010 and from 2011 to 2020. Site velocities at 85 sites were used to create the contour lines for the period from 2001 to 2010, and site velocities at 220 sites were used to create the subsidence contour lines for the period from 2011 to 2020.The subsidence contour lines indicate that subsidence in downtown Houston and Galveston County (HGSD Areas 1, 2) has ceased since the 2000s. Katy in Fort Bend County, Jersey Village in northwestern Harris County (HGSD Area 3), and The Woodlands in southern Montgomery County have been the most affected areas by subsidence (>1 cm/year) since the 2000s. There is a clear trend that the overall size of the subsiding area is shrinking, and the overall subsidence rate decreased from the 2000s to the 2010s. This is the result of the groundwater regulations enforced by HGSD, FBSD, LSGCD, and other local groundwater management agencies.Compared with Harris, Fort Bend, and Galveston Counties, there are few GNSS stations in other counties. There were only three stations [TXCN, P013, and P012 (Fig. 12)] in Montgomery County from 2001 to 2010. Therefore, the 2001–2010 subsidence contour lines in Montgomery County (Fig. 12) were deduced from very limited measurements and should be interpreted with caution. Fig. 13 depicts the subsidence time series at TXCN (2005–2021), P013 (2000–2021), and P012 (2000–2021). P012 is adjacent to the southeast border of Montgomery County, within the Kingwood area, the largest master-planned residential community in northern Harris County and southern Montgomery County. P012 has had a steady subsidence rate of approximately 5 mm/year since 2007. TXCN is located in Conroe, central Montgomery County. P013 is located in The Woodlands, southern Montgomery County. TXCN and P013 are approximately 17 km apart and have had very similar subsidence trends over time (2005–2021). The subsidence rate was about 14–16 mm/year (2005–2015) before 2015, and it has decreased to about 5–6 mm/year (2016–2021) since 2016. The reduction of subsidence rates in 2016 was coincident with the availability of surface water for public use in Montgomery County since 2015 (Wang et al. 2021). As a part of the groundwater usage reduction plan enforced by LSGCD (LSGCD 2013), treated surface water from Lake Conroe has been transported to urban areas, such as the City of Conroe and The Woodlands, in the central and southern county since 2015.The GNSS antennas at five sites—Addicks (ADKS), −549 m, below land surface; Lake Houston (LKHU), −591 m below land surface; Northeast (NETP), −661 m below land surface), Clear Lake (TXEX), −936 m below land surface; and Katy (UHKD), −904 m below land surface)—are mounted on the top of the inner poles of deep extensometer boreholes (Yu et al. 2014). The locations of these five stations are marked in Fig. 12. A photo of the deep-seated GNSS site (LKHU) is shown in Fig. 2(f). These deep-seated GNSS antennas were installed to provide stable references for calculating site deformation at other GNSS sites. The antennas continuously record the vertical movements at the bottom of each borehole, not the movements at the land surface (e.g., Yu et al. 2014; Wang et al. 2014). Accordingly, the GNSS-derived subsidence rates at these five sites should not be used for creating land surface subsidence contour maps.Monitoring Fault ActivitiesFig. 14 depicts GNSS-derived horizontal velocity vectors within the Greater Houston region during the last decade (2011–2020), and shows active faulting traces in the Houston region mapped by the USGS using high-resolution airborne light detection and ranging (lidar) data collected in October 2001 (Shah and Lanning-Rush 2005). These faults were formed millions of years ago during the formation of GOM, and belong to a class of geologic structures known as growth faults cutting the pre-Miocene-age sediments (Holocene-, Pleistocene-, and Pliocene-age sediments) along the GOM coast. These faults have been documented at depths from 1,000 to 4,000 m based on extensive investigations using geophysical well logs and deep seismic surveys (Verbeek et al. 1979).The velocity vectors indicate that the majority of GNSS sites are horizontally stable (<1 mm/year) with respect to Houston20. Only a few places experienced localized horizontal movements larger than 2 mm/year during the last decade. There are no spatially coherent movement trends except in the Long Point Fault area as recorded by the Long Point GNSS array, which comprises 13 permanent stations distributed along the 2 sides of the main trace of the Long Point Fault. These sites have a coherent horizontal movement toward the northwest with an average speed of approximately 2–3 mm/year (2014–2021). Further detailed investigation indicated no considerable difference between stations at two sides of the fault trace with regard to the direction and magnitude of site velocities. According to Liu et al. (2019), the coherent movement toward the northwest is associated with the developing of a subsidence bowl in the Jersey Village area since the 1990s, rather than with deep-seated or tectonic-controlled fault movements. The ongoing subsidence within the subsidence bowl is approximately 1–2 cm/year [Figs. 11(a) and 12].Fig. 15 illustrates the GNSS-derived displacement time series at three sites (P038, P075, and P090) that show remarkable horizontal velocities. These sites are adjacent to known faulting traces (Fig. 14). The ongoing horizontal site velocity is approximately 10 mm/year (2012–2021) at P038 toward the northeast, 4 mm/year (2012–2021) at P075 toward the northwest, and 3 mm/year (2016–2021) at P090 toward the northeast. These sites have been stable (<±3 mm/year) in the vertical direction since the 2000s (Fig. 12). The horizontal displacement time series at each site reveals steady linear movements, suggesting that the movements may be associated with faulting activities. Nevertheless, the movement at a single location does not indicate movements along the whole fault. In fact, several stations adjacent to these sites do not show any considerable horizontal movements. Long-term and dense GNSS observations are needed to understand the relationships between horizontal and vertical ground movements, between faulting activities and groundwater withdrawals, between faulting and subsidence, between ground deformation and flooding, and between ground deformation and structure damages.References Agudelo, G., G. Wang, Y. Liu, Y. Bao, and M. J. Turco. 2020. “GPS geodetic infrastructure for subsidence and fault monitoring in Houston, Texas, USA.” Proc. Int. Assoc. Hydrol. Sci. 382 (Apr): 11–18. https://doi.org/10.5194/piahs-382-11-2020. Altamimi, Z., P. Rebischung, L. Métivier, and X. Collilieux. 2016. “ITRF2014: A new release of the International Terrestrial Reference Frame modeling nonlinear station motions.” J. Geophys. Res.: Solid Earth 121 (8): 6109–6131. https://doi.org/10.1002/2016JB013098. Baker, E. T. J. 1979. Stratigraphic and hydrogeologic framework of part of the coastal plain of Texas. Texas Department of Water Resources Rep. No. 236. Austin, TX: Texas Dept. of Water Resources. Bao, Y., W. Guo, G. Wang, W. Gan, M. Zhang, and S. Shen. 2018. “Millimeter-accuracy structural deformation monitoring using stand-alone GPS: Case study in Beijing, China.” J. Surv. Eng. 144 (1): 05017007. https://doi.org/10.1061/(ASCE)SU.1943-5428.0000242. Bao, Y., X. Yu, G. Wang, H. Zhou, X. Ding, G. Xiao, S. Shen, R. Zhao, and W. Gan. 2021. “SChina20: A stable geodetic reference frame for ground movement and structural deformation monitoring in South China.” J. Surv. Eng. 147 (3): 04021006. https://doi.org/10.1061/(ASCE)SU.1943-5428.0000352. Bawden, G., M. Johnson, M. Kasmarek, J. Brandt, and C. S. Middleton. 2012. Investigation of land subsidence in the Houston-Galveston region of Texas by using the global positioning system and interferometric synthetic aperture radar, 1993–2000. USGS Scientific Investigation Rep. No. 2012-5211. Washington, DC: USGS. Bertiger, W., S. D. Desai, B. Haines, N. Harvey, A. W. Moore, S. Owen, and J. P. Weiss. 2010. “Single receiver phase ambiguity resolution with GPS data.” J. Geod. 84 (5): 327–337. https://doi.org/10.1007/s00190-010-0371-9. Blewitt, G., C. Kreemer, W. Hammond, and J. M. Goldfarb. 2013. “Terrestrial reference frame NA12 for crustal deformation studies in North America.” J. Geodyn. 72 (Dec): 11–24. https://doi.org/10.1016/j.jog.2013.08.004. Böhm, J., G. Möller, M. Schindelegger, G. Pain, and R. Weber. 2015. “Development of an improved empirical model for slant delays in the troposphere (GPT2w).” GPS Solutions 19 (3): 433–441. https://doi.org/10.1007/s10291-014-0403-7. Casarez, I. R. 2020. Aquifer extents in the coastal lowlands aquifer system regional groundwater availability study area in Texas, Louisiana, Mississippi, Alabama, and Florida. Washington, DC: USGS. https://doi.org/10.5066/P9BH2KG2. Cleveland, W. S. 1981. “LOWESS: A program for smoothing scatterplots by robust locally weighted regression.” Am. Stat. 35 (1): 54. https://doi.org/10.2307/2683591. Coplin, L., and D. Galloway. 1999. “Houston Galveston, Texas: Managing coastal subsidence.” Vol. 1182 of Land subsidence in the United States, 35–48. Washington, DC: USGS. Ding, K., J. T. Freymueller, P. He, Q. Wang, and C. Xu. 2019. “Glacial isostatic adjustment, intraplate strain, and relative sea level changes in Eastern United States.” J. Geophys. Res.: Solid Earth 124 (6): 6056–6071. https://doi.org/10.1029/2018JB017060. Dong, D., P. Fang, Y. Bock, M. K. Cheng, and S. Miyazaki. 2002. “Anatomy of apparent seasonal variations from GPS-derived site position time series.” J. Geophys. Res. 107 (B4): ETG 9-1–ETG 9-16. https://doi.org/10.1029/2001JB000573. Gabrysch, R. K. 1982. Ground-water withdrawals and land-surface subsidence in the Houston-Galveston region, Texas, 1906-80. USGS Open-File Rep. No. 82-571. Washington, DC: USGS. Gazeaux, J., et al. 2013. “Detecting offsets in GPS time series: First results from the detection of offsets in GPS experiment.” J. Geophys. Res.: Solid Earth 118 (5): 2397–2407. https://doi.org/10.1002/jgrb.50152. Gökalp, E., O. Güngör, and Y. Boz. 2008. “Evaluation of different outlier detection methods for GPS networks.” Sensors 8 (11): 7344–7358. https://doi.org/10.3390/s8117344. Greuter, A., M. J. Turco, C. M. Petersen, and G. Wang. 2021. “Impacts of groundwater withdrawal regulation on subsidence in Harris and Galveston counties, Texas, 1978-2020.” GeoGulf Trans. 71: 109–118. Griffiths, J., and J. Ray. 2016. “Impacts of GNSS position offsets on global frame stability.” Geophys. J. Int. 204 (1): 480–487. https://doi.org/10.1093/gji/ggv455. Guo, W., G. Wang, Y. Bao, P. Li, M. Zhang, Q. Gong, R. Li, Y. Gao, R. Zhao, and S. Shen. 2019. “Detection and monitoring of tunneling-induced riverbed deformation using GPS and BeiDou: A case study.” Appl. Sci. 9 (13): 2759. https://doi.org/10.3390/app9132759. Heflin, M., A. Donnellan, J. Parker, G. Lyzenga, A. Moore, L. G. Ludwig, J. Rundle, J. Wang, and M. Pierce. 2020. “Automated estimation and tools to extract positions, velocities, breaks, and seasonal terms from daily GNSS measurements: Illuminating nonlinear Salton Trough deformation.” Earth Space Sci. 7 (7): e2019EA000644. https://doi.org/10.1029/2019EA000644. Heiskanen, W. A., and H. Moritz. 1967. Physical geodesy, 364. London: W. H. Freeman. Herring, T. A., T. I. Melbourne, M. H. Murray, M. A. Floyd, W. M. Szeliga, R. W. King, D. A. Phillips, C. M. Puskas, M. Santillan, and L. Wang. 2016. “Plate boundary observatory and related networks: GPS data analysis methods and geodetic products.” Rev. Geophys. 54 (4): 759–808. https://doi.org/10.1002/2016RG000529. HGSD (Harris-Galveston Subsidence District). 2021. “Subsidence rates in Harris County, Texas and surrounding areas.” Accessed April 1, 2022. https://hgsubsidence.org. Huffman, A. C., Jr., S. A. Kinney, L. R. H. Biewick, H. R. Mitchell, and G. L. Gunther. 2004. “Salt diapirs in the Gulf Coast.” In Gulf Coast geology (GCG) online—Miocene of Southern Louisiana, Version 1. Washington, DC: USGS. Jackson, B., J. D. Sargle, D. Barnes, S. Arabhi, A. Alt, P. Gioumousis, E. Gwin, P. Sangtrakulcharoen, L. Tan, and T. T. Tsai. 2005. “An algorithm for optimal partitioning of data on an interval.” IEEE Signal Process Lett. 12 (2): 105–108. https://doi.org/10.1109/LSP.2001.838216. Kearns, T. J., G. Wang, Y. Bao, J. Jiang, and D. Lee. 2015. “Current land subsidence and groundwater level changes in the Houston metropolitan area (2005–2012).” J. Surv. Eng. 141 (4): 05015002. https://doi.org/10.1061/(ASCE)SU.1943-5428.0000147. Kearns, T. J., G. Wang, M. Turco, J. Welch, V. Tsibanos, and H. Liu. 2019. “Houston16: A stable geodetic reference frame for subsidence and faulting study in the Houston metropolitan area, Texas, US.” Geod. Geodyn. 10 (5): 382–393. https://doi.org/10.1016/j.geog.2018.05.005. Keogh, M. E., and T. E. Törnqvist. 2019. “Measuring rates of present-day relative sea-level rise in low-elevation coastal zones: A critical evaluation.” Ocean Sci. 15 (1): 61–73. https://doi.org/10.5194/os-15-61-2019. Khan, S. D., Z. Huang, and A. Karacy. 2014. “Study of ground subsidence in northwest Harris County using GPS, LiDAR, and InSAR techniques.” Nat. Hazard. 73 (3): 1143–1173. https://doi.org/10.1007/s11069-014-1067-x. Kreemer, C., W. C. Hammond, and G. Blewitt. 2018. “A robust estimation of the 3-D intraplate deformation of the North American plate from GPS.” J. Geophys. Res.: Solid Earth 123 (5): 4388–4412. https://doi.org/10.1029/2017JB015257. Kuchar, J., G. Milne, M. Wolstencroft, R. Love, L. Tarasov, and M. Hijma. 2018. “The influence of sediment isostatic adjustment on sea level change and land motion along the US Gulf Coast.” J. Geophys. Res.: Solid Earth 123 (1): 780–796. https://doi.org/10.1002/2017JB014695. Liu, Y., J. Li, J. Fasullo, and D. L. Galloway. 2020. “Land subsidence contributions to relative sea level rise at tide gauge Galveston Pier 21, Texas.” Sci. Rep. 10 (1): 1–11. https://doi.org/10.1038/s41598-020-74696-4. Liu, Y., X. Sun, G. Wang, M. J. Turco, G. Agudelo, Y. Bao, R. Zhao, and S. Shen. 2019. “Current activity of the long point fault in Houston, Texas constrained by continuous GPS measurements (2013–2018).” Remote Sens. 11 (10): 1213. https://doi.org/10.3390/rs11101213. Militino, A. F., M. Moradi, and M. D. Ugarte. 2020. “On the performances of trend and change-point detection methods for remote sensing data.” Remote Sens. 12 (6): 1008. https://doi.org/10.3390/rs12061008. Miller, M. M., and M. Shirzaei. 2019. “Land subsidence in Houston correlated with flooding from Hurricane Harvey.” Remote Sens. Environ. 225 (May): 368–378. https://doi.org/10.1016/j.rse.2019.03.022. Miller, M. M., and M. Shirzaei. 2021. “Assessment of future flood hazards for southeastern Texas: Synthesizing subsidence, sea-level rise, and storm surge scenarios.” Geophys. Res. Lett. 48 (8): e2021GL092544. https://doi.org/10.1029/2021GL092544. Murray, J. R., and J. Svarc. 2017. “Global positioning system data collection, processing, and analysis conducted by the US Geological Survey earthquake hazards program.” Seismol. Res. Lett. 88 (3): 916–925. https://doi.org/10.1785/0220160204. Petersen, C., M. J. Turco, A. Vinson, J. A. Turco, A. Petrov, and M. Evans. 2020. “Groundwater regulation and the development of alternative source waters to prevent subsidence, Houston region, Texas, USA.” Proc. Int. Assoc. Hydrol. Sci. 382 (Apr): 797–801. https://doi.org/10.5194/piahs-382-797-2020. Qu, F., Z. Lu, J.-W. Kim, and W. Zheng. 2019. “Identify and monitor growth faulting using InSAR over northern greater Houston, Texas, USA.” Remote Sens. 11 (12): 1498. https://doi.org/10.3390/rs11121498. Qu, F., Z. Lu, Q. Zhang, G. W. Bawden, J.-W. Kim, C. Zhao, and W. Qu. 2015. “Mapping ground deformation over Houston–Galveston, Texas using multi-temporal InSAR.” Remote Sens. Environ. 169 (Nov): 290–306. https://doi.org/10.1016/j.rse.2015.08.027. Schaer, S., W. Gurtner, and J. Feltens. 1998. “IONEX: The IONosphere map exchange format version 1.” In Proc., 1998 IGS Analysis Centers Workshop, ESOC. Darmstadt, Germany: European Space Operations Centre. http://ftp.aiub.unibe.ch/ionex/draft/ionex11.pdf. Schmid, R., P. Steigenberger, G. Gendt, M. Ge, and M. Rothacher. 2007. “Generation of a consistent absolute phase-center correction model for GPS receiver and satellite antennas.” J. Geod. 81 (12): 781–798. https://doi.org/10.1007/s00190-007-0148-y. Shah, S., and J. Lanning-Rush. 2005. Principal faults in the Houston, Texas, metropolitan area. USGS Science Investigation Map 2874. Washington, DC: USGS. Thornhill, M. R., and M. R. Keester. 2020. “Subsidence investigations—Phase 1: Assessment of past and current investigations, prepared for lone star groundwater conservation district.” Accessed April 1, 2022. https://www.lonestargcd.org/subsidence. Turco, M. J., and A. Petrov. 2015. “Effects of groundwater regulation on aquifer-system compaction and subsidence in the Houston-Galveston Region, Texas, USA.” Proc. Int. Assoc. Hydrol. Sci. 372 (Nov): 511–514. https://doi.org/10.5194/piahs-372-511-2015. Verbeek, E. R., K. W. Ratzlaff, and U. S. Clanton. 1979. Faults in parts of north-central and western Houston metropolitan area, Texas. US Geological Survey Miscellaneous Field Studies Map 1136. Washington, DC: USGS. Wang, G. 2011. “GPS landslide monitoring: Single base vs. network solutions—A case study based on the Puerto Rico and Virgin Islands permanent GPS network.” J. Geodetic Sci. 1 (3): 191–203. https://doi.org/10.2478/v10156-010-0022-3. Wang, G. 2013. “Millimeter-accuracy GPS landslide monitoring using precise point positioning with single receiver phase ambiguity (PPP-SRPA) resolution: A case study in Puerto Rico.” J. Geodetic. Sci. 3 (1): 22–31. https://doi.org/10.2478/jogs-2013-0001. Wang, G., Y. Bao, Y. Cuddus, X. Jia, J. Serna, and Q. Jing. 2015a. “A methodology to derive precise landslide displacements from GPS observations in tectonically active and cold regions: A case study in Alaska.” Nat. Hazard. 77 (3): 1939–1961. https://doi.org/10.1007/s11069-015-1684-z. Wang, G., Y. Bao, W. Gan, J. Geng, G. Xiao, and J. S. Shen. 2018. “NChina16, a stable geodetic reference frame for geological hazard studies in North China.” J. Geodyn. 115 (Apr): 10–22. https://doi.org/10.1016/j.jog.2018.01.003. Wang, G., H. Liu, G. S. Mattioli, M. M. Miller, K. Feaux, and J. Braun. 2019. “CARIB18: A stable geodetic reference frame for geological hazard monitoring in the Caribbean region.” Remote Sens. 11 (6): 680. https://doi.org/10.3390/rs11060680. Wang, G., J. Welch, T. J. Kearns, L. Yang, and J. Serna Jr. 2015b. “Introduction to GPS geodetic infrastructure for land subsidence monitoring in Houston, Texas, USA.” Proc. Int. Assoc. Hydrol. Sci. 372 (Nov): 297–303. https://doi.org/10.5194/piahs-372-297-2015. Wang, G., J. Yu, T. J. Kearns, and J. Ortega. 2014. “Assessing the accuracy of long-term subsidence derived from borehole extensometer data using GPS observations: Case study in Houston, Texas.” J. Surv. Eng. 140 (3): 05014001. https://doi.org/10.1061/(ASCE)SU.1943-5428.0000133. Wang, G., J. Yu, J. Ortega, G. Saenz, T. Burrough, and R. Neill. 2013. “A stable reference frame for ground deformation study in the Houston metropolitan area, Texas.” J. Geodetic Sci. 3 (3): 188–202. https://doi.org/10.2478/jogs-2013-0021. Wang, G., X. Zhou, K. Wang, X. Ke, Y. Zhang, R. Zhao, and Y. Bao. 2020. “GOM20, a stable geodetic reference frame for subsidence, faulting, and sea-level rise studies along the coast of the Gulf of Mexico.” Remote Sens. 12 (3): 350. https://doi.org/10.3390/rs12030350. Wang, K., G. Wang, B. Cornelison, H. Liu, and Y. Bao. 2021. “Land subsidence and aquifer compaction in Montgomery County, Texas, US: 2000–2020.” Geoenviron. Disasters 8 (1): 28. https://doi.org/10.1186/s40677-021-00199-7. Yang, L., G. Wang, Y. Bao, T. J. Kearns, and J. Yu. 2016. “Comparisons of ground-based and building-based CORS: A case study in the region of Puerto Rico and Virgin Islands.” J. Surv. Eng. 142 (3): 05015006. https://doi.org/10.1061/(ASCE)SU.1943-5428.0000155. Yu, J., and G. Wang. 2016. “GPS derived ground motions (2005–2014) within the Gulf of Mexico region referred to a stable Gulf of Mexico reference frame.” Nat. Hazards Earth Syst. Sci. 16 (7): 1583–1602. https://doi.org/10.5194/nhess-16-1583-2016. Yu, J., G. Wang, T. J. Kearns, and L. Yang. 2014. “Is there deep-seated subsidence in the Houston-Galveston area?” J. Geophys. 2014 (Jan): 1–11. https://doi.org/10.1155/2014/942834. Zhao, R., G. Wang, X. Yu, X. Sun, Y. Bao, G. Xiao, W. Gan, and S. Shen. 2020. “Rapid land subsidence in Tianjin, China derived from continuous GPS observations (2010–2019).” Proc. Int. Assoc. Hydrol. Sci. 382 (Apr): 241–247. https://doi.org/10.5194/piahs-382-241-2020. Zhou, X., G. Wang, K. Wang, H. Liu, H. Lyu, and M. J. Turco. 2021. “Rates of natural subsidence and submergence along the Texas Coast derived from GPS and tide gauge measurements (1904–2020).” J. Surv. Eng. 147 (4): 04021020. https://doi.org/10.1061/(ASCE)SU.1943-5428.0000371. Zilkoski, D. B., L. W. Hall, G. J. Mitchell, V. Kammula, A. Singh, W. M. Chrismer, and R. J. Neighbors. 2003. The Harris Galveston coastal subsidence district/national geodetic survey automated global positioning system subsidence monitoring project. USGS Open File Rep. No. 03-308. Washington, DC: USGS. Zumberge, J. F., M. B. Heflin, D. C. Jefferson, M. M. Watkins, and F. H. Webb. 1997. “Precise point positioning for the efficient and robust analysis of GPS data from large networks.” J. Geophys. Res.: Solid Earth 102 (B3): 5005–5017. https://doi.org/10.1029/96JB03860.