Houston GNSS Network for Subsidence and Faulting Monitoring: Data Analysis Methods and Products

: Harris-Galveston Subsidence District (HGSD), in collaboration with several other agencies, has been operating a dense Global Navigation Satellite System (GNSS) network for subsidence and faulting monitoring within the Greater Houston region since the early 1990s. The GNSS network is designated HoustonNet, comprising approximately 250 permanent GNSS stations as of 2021. This paper documents the methods used to produce position time series, transform coordinates from the global to regional reference frames, identify outliers and steps, analyze seasonal movements, and estimate site velocities and uncertainties. The GNSS positioning methods presented in this paper achieve 2 – 4-mm RMS accuracy for daily positions in the north – south and east – west directions and 5 – 8-mm accuracy in the vertical direction within the Greater Houston region. Five-year or longer continuous observations are able to achieve submillimeter-per-year uncertainties (95% confidence interval) for both horizontal and vertical site velocities. Two decades of GNSS observations indicate that Katy in Fort Bend County, Jersey Village in northwestern Harris County, and The Woodlands in southern Montgomery County have been the areas most affected by subsidence ( 1 – 2 cm = year) since the 2000s; the overall subsidence rate and the size of subsiding area ( > 5 mm = year) have been decreasing as a result of the groundwater regulations enforced by HGSD and other local agencies. HoustonNet data and products are released to the public through HGSD. The primary products are the daily East-North-Up (ENU) position time series and site velocities with respect to the International GNSS Service (IGS) Reference Frame 2014 (IGS14), the stable Gulf of Mexico Reference Frame (GOM20), and the stable Houston Reference Frame (Houston20). The ENU position time series with respect to Houston20 are recommended for delineating subsidence and faulting within the Greater Houston region. The ENU time series with respect to GOM20 are recommended for studying subsidence and faulting within the Gulf coastal plain and sea-level changes along the Gulf Coast. The entire HoustonNet data set is reprocessed every a few years with updated positioning software, IGS and regional reference frames, and data analysis tools. We recommend that users use the most recent release of HoustonNet data products and avoid mixing old and new positions. DOI: 10.1061/ (ASCE)SU.1943-5428.0000399. This work is made available under the terms of the Creative Commons Attribution 4.0 International license, https://creativecommons.org/licenses/by/4.0/.


Introduction
For 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 1 Professor, Dept. of Earth and Atmospheric Sciences, Univ. of Houston, Houston, TX 77204 (corresponding author).ORCID: https://orcid.org/0000-0003-3731-3839. Email: gwang@uh.edu 2 Project Chief, Harris-Galveston Subsidence District, 1660 West Bay Area Blvd., Friendswood, TX 77546.ORCID: https://orcid.org/0000-0003-0994-2566.Email: agreuter@subsidence.org 3 Dupty General Manager, Harris-Galveston Subsidence District, 1660 West Bay Area Blvd., Friendswood, TX 77546.Email: tpetersen@ subsidence.org 4 General Manager, Harris-Galveston Subsidence District, 1660 West Bay Area Blvd., Friendswood, TX 77546.Email: mturco@subsidence.orgNote.This manuscript was submitted on November 24, 2021; approved on March 8, 2022; published online on June 25, 2022.Discussion period open until November 25, 2022; separate discussions must be submitted for individual papers.This paper is part of the Journal of Surveying Engineering, © ASCE, ISSN 0733-9453.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. 2019Wang et al. , 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 andShirzaei 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 km 2 (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.

HoustonNet
As 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 GNSS
HGSD 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 GNSS
The 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.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 Flow
The 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).
Producing Daily Solutions with Respect to the IGS Reference Frame GNSS 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.0IGS14 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 sevenparameter 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. (2010Bertiger et al. ( , 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. (2013Wang et al. ( , 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 Outliers
The 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)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020)(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): 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 μ AE 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 AE3σ (or AE4.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 Transformation
A 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:  to IGS14 and Houston20 or GOM20; T 0 x , T 0 y , and T 0 z are constant parameters indicating rates (one-time derivatives) of three translational shifts along x, y, z-coordinate axes; and R 0 x , R 0 y , and R 0 z = 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 t 0 , T 0 x , T 0 y , T 0 z , R 0 x , R 0 y , and R 0 z are listed in Table 1.To 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 10 6 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: where λ 0 and φ 0 = initial geodetic longitude and latitude of site corresponding to initial XYZ coordinate ðX 0 ; Y 0 ; Z 0 Þ, 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 AE 0.3 mm=year with respect to Houston20 and 3.2 AE 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 AE 0.3 mm=year with respect to Houston20 and 5.8 AE 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 Hous-ton20 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 AE 0.2 mm= year, which is unrealistically large.The horizontal velocity of 1.1 AE 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 Pole
The 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 R 0 x , R 0 y , and R 0 z used to realize the regional reference frame (Table 1) Table 1.Seven parameters for realizing Houston20 and GOM20 Parameter a IGS14 to Houston20 IGS14 to GOM20 a These seven parameters are used to transform ECEF-XYZ coordinates from IGS14 to Houston20 and GOM20 according to Eq. ( 2).Counterclockwise rotations around the x-, y-, and z-axes are positive.
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 (R 0 x , R 0 y , R 0 z ) 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 (R 0 x , R 0 y , and R 0 z ) 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 welldetermined 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 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 Steps
A 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 Houston-Net 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 Modeling
Seasonal 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) 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) where ; the amplitude of the semiannual signal can be estimated as . The peak-to-trough amplitude of the seasonal motions can be estimated as 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 AE 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) (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 (2015MRHK ( -2020) ) and 12 mm=year at UHCR (2015)(2016)(2017)(2018)(2019)(2020).The peakto-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 sitespecific 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 Uncertainties
In 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 Applications
The 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 Subsidence
As 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).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.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 deepseated 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 Activities
Fig. 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 (< AE 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.

Discussion and Concluding Remarks
The geodesy research community continuously has been improving the current models and developing new methods for high-accuracy high-precision GNSS positioning.Current positioning methods employed in HoustonNet routine processing are able to achieve 2-4 mm RMS accuracy in the horizontal directions and 5-8 mm RMS accuracy in the vertical direction.Five-year or longer continuous observations can achieve submillimeter-per-year uncertainties (95% CI) for site velocity estimates in both horizontal and vertical directions.The entire HoustonNet data are reprocessed every 2 or 3 years to integrate the state-of-the-art achievements in global and regional GNSS geodetic infrastructures.As a result, the positional time series from the most recent processing may have higher precision and fewer outliers than those published a few years ago, which would result in more-reliable site velocity estimates.Accordingly, we strongly encourage data users to use the most recent data products rather than merging the new products with the previous products.GOM20 and Houston20 are the key components of the regional geodetic infrastructure, which will be improved incrementally and synchronized with the updates of the IGS reference frame.The next significant change to the HoustonNet data processing will be the update of regional reference frames.ITRF is expecting to replace ITRF14 with ITRF20 in a few years.Subsequently, IGS will update IGS14 to IGS20.In turn, the Stable Houston Reference Frame and the Stable Gulf of Mexico Reference Frame will be updated to align with IGS20.In the future, it is likely that full GNSS processing (GPS + GLONASS + Galileo + BeiDou) will be adopted in HoustonNet processing.The new GNSS signals will complement the GPS-only signals and help identify certain systemic errors and biases associated with GPS.Preparations for a transition to full GNSS analysis are already underway.
This paper summarizes the ongoing subsidence, faulting, and seasonal motions within the Greater Houston region.The amplitudes of seasonal ground deformation in subsiding areas are remarkably larger than the amplitudes in nonsubsiding areas, which is evident in GNSS-derived displacement time series and ready for further analysis and interpretation.HoustonNet and its data products have the potential for broad applications in urban geohazards and regional sea-level studies.The data analysis methods introduced in this paper are applicable to other regional GNSS networks.

Fig. 1 .
Fig. 1.Map showing the geographic locations of HoustonNet and reference GNSS stations used for realizing the Stable Houston Reference Frame (Houston20) and the Stable Gulf of Mexico Reference Fame (GOM20).The horizontal site velocity vectors (2005-2020) refer to the International GNSS Service (IGS) Reference Frame 2014 (IGS14).

Fig. 2 .Fig. 4 .Fig. 3 .
Fig. 2. Locations of HoustonNet GNSS stations.Squares represent 120 Harris-Galveston Subsidence District (HGSD) stations.Triangles represent 65 continuous stations operated by the University of Houston (UH).Circles represent approximately 65 continuous stations operated by other agencies.The color of each station indicates its observational history.Texas coast aquifer system outcrops in the northwest corner of the Greater Houston region (Casarez 2020).The color areas are regulated by the Harris-Galveston Subsidence District (HGSD Areas 1, 2, 3), the Fort Bend Subsidence District (FBSD Areas A and B), the Brazoria County Groundwater Conservation District (BCGCD), and the Lone Star Groundwater Conservation District (LSGCD).

where
Fig. 6.Plots showing the methods for identifying and removing outliers: (a) three steps for identifying and removing outliers from the geocentric coordinate time series (y-component); and (b) removed outliers shown in the ENU time series.

Fig. 7 .
Fig. 7. Displacement time series compared with the Stable Houston Reference Frame (Houston20) and the Stable Gulf of Mexico Reference Frame (GOM20) at three permanent GNSS stations located inside and outside the Greater Houston region: (a) TXLM (Houston); (b) LMCN (New Orleans); and (c) ZMA1 (Miami).Locations of LMCN and ZMA1 are marked in Fig. 1.The location of TXLM is marked in Fig. 2.

Fig. 8 .
Fig. 8. Plots illustrating the potential steps detected by the automated change point detection employed in the routine data analysis: (a) JGS2 is a continuous GNSS station located in Dayton, Liberty County; and (b) P043 is a PAM GNSS station situated on the southwestern end of Galveston Island.The ENU displacement time series are aligned with Houston20.Locations of JGS2 and P043 are marked in Fig. 2.

Fig. 10 .Fig. 11 .
Fig. 10.Seasonal vertical ground movements at (a) nonsubsiding sites; and (b) subsiding sites.The time series (solid circles) are the residuals after removing the linear and nonlinear trends [Eq.(8)].Locations of these sites are marked in Fig. 2.
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.

Fig. 14 .
Fig. 14.Map depicting horizontal velocity vectors at 220 GNSS sites with respect to Houston20 derived from decadal observations (2011-2020) within the Greater Houston region.Uncertainties (95% CI) of the horizontal velocities are within AE0.5 mm=year.The rectangle indicates the Long Point Fault GNSS array.The lines represent active faults mapped by USGS(Shah and Lanning-Rush 2005).The solid polygons represent salt domes mapped by USGS(Huffman et al. 2004).

Fig. 15 .
Fig. 15.Horizontal displacement time series and trajectories with respect to Houston20 at three sites adjacent to active faults: (a) P038; (b) P075; and (c) P090.Locations of the GNSS sites and faults are marked in Fig. 14.

Table S1 . General Information of HoustonNet GNSS (as of October 2021)
The vertical displacement occurred during the entire observational time span of the station.** HGSD: Harris-Galveston Subsidence District; FBSD: Fort Bend Subsidence District; LSGCD: Lone Star Groundwater Conservation District; BCGCD: Brazoria County Groundwater Conservation District; UH: University of Houston; SmartNet: SmartNet North America; USCG: U.S. Coastal Guard; COH: City of Houston; USF: University of South Florida; LCRA: Lower Colorado River Authority; NOAA: National Oceanic and Atmospheric Administration; RODS: RODS Surveying Inc.; TXAM: Texas A&M University; TxDOT: Texas Department of Transportation; FAA: Federal Aviation Administration.