Identify and Monitor Growth Faulting Using InSAR over Northern Greater Houston , Texas , USA

Growth faults are widely distributed in the Greater Houston (GH) region of Texas, USA, and the existence of faulting could interrupt groundwater flow and aggravate local deformation. Faulting-induced property damages have become more pronounced over the last few years, necessitating further investigation of these faults. Interferometric synthetic aperture radar (InSAR) has been proved to be an effective way for mapping deformations along and/or across fault traces. However, extracting short-wavelength small-amplitude creep signal (about 10–20 mm/yr) from long time span interferograms is extremely difficult, especially in agricultural or vegetated areas. This study aims to position, map and monitor the rate, extent, and temporal evolution of faulting over GH at the highest spatial density using Multi-temporal InSAR (MTI) technique. The MTI method, which maximizes usable signal and correlation, has the ability to identify and monitor faulting and provide accurate and detailed depiction of active faults. Two neighboring L-band Advanced Land Observing (ALOS) tracks (2007–2011) are utilized in this research. Numerous areas of sharp phase discontinuities have been discerned from MTI-derived velocity map. InSAR measurements allow us to position both previously known faults traces as well as nucleation of new fractures not previously revealed by other ground/space techniques. Faulting damages and surface scarps were evident at most InSAR-mapped fault locations through our site investigations. The newly discovered fault activation appears to be related to excessive groundwater exploitation from the Jasper aquifer in Montgomery County. The continuous mining of groundwater from the Jasper aquifer formed new water-level decline cones over Montgomery County, corroborating the intensity of new fractures. Finally, we elaborate the localized fault activities and evaluate the characteristics of faulting (locking depth and slip rate) through modeling MTI-derived deformation maps. The SW–NE-oriented faults pertain to normal faulting with an average slip rate of 7–13 mm/yr at a shallow locking depth of less than 4 km. Identifying and characterizing active faults through MTI and deformation modeling can provide insights into faulting, its causal mechanism and potential damages to infrastructure over the GH.


Introduction
Greater Houston, Texas, USA (hereafter GH), comprising Montgomery, Grimes, Waller, Harris, Liberty, Chambers, Austin, Fort Bend, Brazoria and Galveston Counties, located on the Gulf coastal plain, is laced by numerous growth faults [1][2][3] (Figure 1).The surficial processes of the deep-seated ancient faults could be depicted as a cut-and-fill scarp, which have impacted and damaged buildings and infrastructures [4,5].Hundreds of paved roads and homes in the Houston area are being offset by faults and require frequent maintenance.Airports, factories, commercial establishments, and railroads have also been affected by fault movement, and it costs several million dollars to repair such structures a year [6].There are more subtle or hidden effects of faulting which have never been recognized and many breaks may not be related to faulting because the existence of faults is unknown [6,7]. million dollars to repair such structures a year [6].There are more subtle or hidden effects of faulting which have never been recognized and many breaks may not be related to faulting because the existence of faults is unknown [6,7].[8], where colors represent generalized lithology, and black letters (Ql-Lissie Formation; Pow-Willis Formation; Qb-Beaumont Formation; Qbs-Beaumont Formation (sand); Qbc-Beaumont Formation (clay); Qal-Alluvium; Qt-Fluviatile terrace deposits) and grey lines show distribution of the major geologic units [1].Faults (marked by black lines) [9,10], hydrocarbon wells (small grey dots [11]), GPS benchmarks (green triangles), salt dome positions (green polygons) and county boundaries (white dashed lines, county names are labeled by purple letters) are superimposed.The blue and black rectangles show coverage of two neighboring Advanced Land Observing (ALOS) data.The Texas State district is shown as an inset map on the upper-right corner.
While the lithology of GH is primarily unconsolidated, normal faulting does exist along with the tectonic evolution of the Gulf of Mexico basin and associated sediments of the Gulf Coast aquifer since Middle Jurassic [1].Numerous growth faults occur parallel to the United States Gulf Coast, dipping from northwest to southeast [1].It has been suggested that faulting in the GH area is mainly related to the development of the Gulf of Mexico, its regional faults, sedimentary history, salt movement and fluid (oil, gas, and groundwater) extraction [12][13][14].Movement of faults in the Gulf Coast was first described in 1926 following local land-surface subsidence around the Goose Creek oil field about 20 miles east of Houston [15].However, it was not until the 1970s that the public began to recognize their importance as geologic hazards, and numerous papers devoted to various aspects of faulting in the GH area raised a dramatic and continuing discussion about the fault movement mechanisms (e.g., [13,16,17]).Beginning in the late 1970s and early 1980s, the U.S. Geological Survey (USGS) launched an extensive fault study of the GH area [7,13,[17][18][19].Over 300 active faults with an aggregate length of up to 500 km were discovered over the HG region [5,17].These faults generally move as aseismic fault creep, with a typical creep rate ranging from 4 to 27 mm/yr [20][21][22][23][24][25][26].[8], where colors represent generalized lithology, and black letters (Ql-Lissie Formation; Pow-Willis Formation; Qb-Beaumont Formation; Qbs-Beaumont Formation (sand); Qbc-Beaumont Formation (clay); Qal-Alluvium; Qt-Fluviatile terrace deposits) and grey lines show distribution of the major geologic units [1].Faults (marked by black lines) [9,10], hydrocarbon wells (small grey dots [11]), GPS benchmarks (green triangles), salt dome positions (green polygons) and county boundaries (white dashed lines, county names are labeled by purple letters) are superimposed.The blue and black rectangles show coverage of two neighboring Advanced Land Observing (ALOS) data.The Texas State district is shown as an inset map on the upper-right corner.
While the lithology of GH is primarily unconsolidated, normal faulting does exist along with the tectonic evolution of the Gulf of Mexico basin and associated sediments of the Gulf Coast aquifer since Middle Jurassic [1].Numerous growth faults occur parallel to the United States Gulf Coast, dipping from northwest to southeast [1].It has been suggested that faulting in the GH area is mainly related to the development of the Gulf of Mexico, its regional faults, sedimentary history, salt movement and fluid (oil, gas, and groundwater) extraction [12][13][14].Movement of faults in the Gulf Coast was first described in 1926 following local land-surface subsidence around the Goose Creek oil field about 20 miles east of Houston [15].However, it was not until the 1970s that the public began to recognize their importance as geologic hazards, and numerous papers devoted to various aspects of faulting in the GH area raised a dramatic and continuing discussion about the fault movement mechanisms (e.g., [13,16,17]).Beginning in the late 1970s and early 1980s, the U.S. Geological Survey (USGS) launched an extensive fault study of the GH area [7,13,[17][18][19].Over 300 active faults with an aggregate length of up to 500 km were discovered over the HG region [5,17].These faults generally move as aseismic fault creep, with a typical creep rate ranging from 4 to 27 mm/yr [20][21][22][23][24][25][26].
Early detection of active faults can help mitigate impacts to buildings and infrastructure.The active faults of GH have been delineated by several researchers (e.g., [10,13,18,19,27,28]). Fault detection can be accomplished through various means.The aerial photographs interpretation followed by a ground inspection provides an efficient way to identify the specific fault-affected areas in a broad scale (e.g., [29,30]).Aerial photos are often helpful for identifying the points of contact between the drier upthrown side and the wetter downthrown side of a fault because of the changes in vegetation.However, the applications of aerial photographs are limited over areas of rapid expansion, significant topographic relief, and/or dense vegetation or tree cover.Thus, ground surveys are the principal means to identify and map surface fault fractures over much of north GH, including northeastern Harris County and Montgomery County (e.g., [31,32]).Ground investigations locate faults on the ground through surface clues such as topographic scarps, sharp changes in vegetative communities, offset stream meanders and some other ground representations, but most faults over GH may not show phenomenon on the ground surface and/or the processes are too subtle to be detected under the field investigations.Drilling in hydrocarbon exploration activity could provide crucial evidences, such as the location and orientation, for the deeper faults, and perceive the movement of surface faults over the well location.Various ground-based geophysical surveys, such as geophysical logs [33], seismic reflection [27,28,34], resistivity [24,25,29,[35][36][37], conductivity [35,36], magnetic imaging [35,36], ground-penetrating radar (GPR) [24,27,28,35], and gravity [28,[34][35][36][37] have all been used to identify fault zones with varying degrees of success over small areas.Among the above techniques, GPR has been used for faults investigation in the GH region to a certain extent, but their applications were limited over some areas with high soil moisture [34].None of the geophysical methods have been entirely satisfactory to apply to all the faults and most of them are expensive and time consuming.Airborne LiDAR (Light Detection and Ranging) could produce high-resolution DEMs with typically 10-30 cm vertical precision depending on the land cover and terrain conditions.Two observed DEMs can then be differentiated to generate a time-lapse picture and obtain elevation changes across the faults, i.e., DoD (DEM of difference), which is a kind of efficient methods in identifying fault scarps/ruptures.Researchers from the University of Houston have published an updated faults map of the GH region with a number of new discovered faults using LiDAR [27,28].Unfortunately, none of the methods above indicate the historical and/or recent activity of faults over time.Site reconnaissance using the permanent constructed global positioning systems (GPS) and leveling sites across fault traces could precisely monitor local fault motions and is quit applicable in revealing the fault segment movements under the presuppose that well awareness of fault locations (e.g., [38][39][40][41][42][43][44]).
Extensive fault study has been conducted and more than 350 known growth faults were discovered in the Harris County metropolitan area during the later last century, and three main fault systems in Houston region are the Long Point Faults system, the Addicks Faults system, and the Hockley-Conroe Faults system, including the Long Point, Hockley, Addicks, Tomball, Willow Creek, and Eureka Faults [19,[23][24][25]36].Most of the mapped faults over Houston region are predominantly listric normal faults, while the main structures in northern Greater Houston, mainly north Harris County and Montgomery County, are normal faulting and dip-slip down mostly to the coast with about 60 • -75 • dip angle, suggesting the rate of horizontal movement component is less than one-fourth of that from vertical motion [3,6,17,23,34,37,45].The study of faults slowed beginning in the 1990s, but many faults, especially the extremely active ones over north GH, have not been completely mapped out or researched.However, they are still active and represent potential geohazards to their surrounding environment, as is evident by reports from local residents in Montgomery County claiming property damage due to faulting during the last decade (e.g., [46,47]).Multiple attempts to acquire seismic surveys in this area were unsuccessful due to the heavy vegetation [34].The mechanism of such geo-hazards (i.e., active faults) is still unclear, with many relatively small faults or fissures still unrecognized over GH.Though Growth faults in the GH area are aseismic, and their motions do not give rise to substantial seismic events, the potential risks of active faults are still high, with potential to cost millions of dollars in property and infrastructure [3].In some cases, the location of specific fault segments remained unknown until they caused extensive damage to local neighborhoods.Accurately locating active faults remains crucial for protecting people and infrastructures from severe damage.Project developers and government planners can avoid higher risk areas or accommodate potential ground shifts in their construction plans with prior knowledge of the location of faults.A more detailed study of fault segments is necessary to determine the cause of increased fault activity and therefore reduce the potential hazard of the active faults.
Interferometric Synthetic Aperture Radar (InSAR) is an established and reliable technique to study surface displacements.InSAR's ability to retrieve historic deformation, identify small sized surface anomalies, comprehensively reveal the spatial dimension of ground displacement and detect deformation borders in a high spatial resolution of meters and a centimeter-to millimeter-measurement accuracy are some of the advantages of this technique.Discontinuity in deformation gradient across fault lines offers the opportunity to map faults in InSAR deformation maps, because two sides of faults generally move past each other and present different deformation rate (e.g., [48][49][50][51][52][53]). InSAR velocity fields have been not only used to describe the geometry of the fault rupture during earthquake, but also have been employed in evaluating the long-wavelength deformation associated with interseismic strain accumulation.Though InSAR has proven successful in monitoring many large-scale long-wavelength seismic faults (e.g., [49][50][51]54,55]), few researches focused on the aseismic fault creep signal, and most of which are still large-scale and can be related to some seismic activities (e.g., [56][57][58][59]).Almost no research, to our knowledge, was mainly focus on identifying the active aseismic creep signal, which has potential hazard but have not been mapped out.Challenges remain particularly in the plant-covered lands or coastal regions, such as extracting a short wavelength, small scale creep signal of about 10 mm/yr, from long time period InSAR interferograms.Several scholars investigated the subsidence of GH area utilizing InSAR technique, but there is no particular InSAR research focus on the continuous faulting anomaly over this region [10,26,53,60,61].The authors of [26,60,61] examined land subsidence of Houston (mainly in Harris County) during the period of 1990s using C-band ERS-1/2 interferometry respectively and comprehensive analysis the subsidence bowl centering at Jersey Village by integrating InSAR observation with extensometer and GPS measurements.The displacement anomalies along Long Point Fault was identified and discussed in both [26] and [61] using InSAR technology.The authors of [10] published PS-InSAR (Persistent Scatterer InSAR) deformation over a 55*5 km 2 rectangle area in northwest Harris County using 25 ERS-1/2 images.These works located subsidence of Houston from 1992 to 2002 through either analyzing individual interferometric phase signals or PS-InSAR in a localized area, containing amount of noise due to the phase decorrelation, atmosphere, and other phase errors.The authors of [53] explored ground deformation features over Houston-Galveston region during the period of 1993-2011 utilizing C-band ERS-1/2, Envisat, and L-band ALOS data by MTI (Multi-temporal InSAR) method.They observed 5 mm/yr to 40 mm/yr differential deformation rate along numerous fault locations over the study area.
Growth faulting has influenced a wide variety of geological conditions in GH, and is known to have impacted and damaged buildings, highways, wells, and pipelines [2,4].Locating and characterizing active faults is crucial for protecting people and infrastructure from severe damage.The objective of this work is to carry out an integrated study of the active faults over northern GH region, including establishing position and monitoring the distribution, velocity, and temporal development of faulting utilizing MTI technique.Furthermore, we use the analysis to identify the key driving mechanisms.To do so, we utilized L-band ALOS datasets, which can reach to the ground partially penetrating through vegetation to obtain ground surface information and reveal the faulting creep over northern GH region.First, we estimated the long-term deformation rate from 2007 to 2011 by MTI to characterize the spatial distribution of active faults.Second, we mapped fault fractures by phase jumps/discontinuities identified from our InSAR average deformation map.Our InSAR mapped fault traces were validated through LiDAR, geophysical survey and field investigated observations.Third, we derived the fault slip models from two independent InSAR displacement measurements to study the characteristics of faulting.Finally, we expound how the faulting over GH is in connection with regional faults, sedimentary history of the Gulf of Mexico, salt movement and fluid extraction.

Geologic Background and Hydrologic Setting
GH lies largely in the smooth and low-lying Gulf coastal plain with elevation gradually rising inland from about 0 to 110 m.The east and north GH are densely forested, while the west and south parts are mainly covered with prairie grassland, and the coastal regions are covered in prairie and sand [1].The surficial geology of GH is complex, largely constituted by poorly-cemented sands, clay shales and unconsolidated clays, into few kilometers deep.Figure 1 shows the geological map of GH [8].GH is underlain by sediment deposited along Gulf of Mexico during the Cenozoic era; large amounts of unconsolidated sands and sandy clays were deposited throughout the Pleistocene and Holocene periods [1].Both the Lissie Formation and the Beaumont Formation are part of the Houston lithology, and the Willis Formation underlies uncomfortably the Lissie Formation stratigraphically.The older Pleistocene Willis Formation mainly contains clays with less volume of sands and silts.The Pleistocene Lissie formation is primarily composed of sands with fewer silts and clays, while the Beaumont formation is comprised of finer clays with silt [8,62].The contacts regions among formations are usually with the characteristic of low cohesion and will be easily to turn into normal faults, and many faults of the GH region lie at the contact of these two formations, according to the geologic map (Figure 1).For example, the Hockley fault within the Hockley-Conroe fault system lies along the contact of the Willis (clay-dominated) upper lithologic unit at the upthrown block, with the Lissie (sand-dominated) upper unit on the downthrown side of the fault [3,23,34].
GH areas receive all or part of their water from the Gulf Coast Aquifer system, which spreads across the Gulf of Mexico coastline from the U.S.-Mexico border to the Texas-Louisiana border [63], and extensive amount groundwater of about 1.3 billion cubic meters are exploited from this aquifer per year [1,63].The Gulf Coast Aquifer according to their facies and hydraulic characteristics has been divided into five hydrological units: the Chicot aquifer (top layer), the Evangeline aquifer, the Burkeville confining layer, the Jasper aquifer, and the Catahoula aquifer (at the bottom), dipping from the northwest to southeast [62].Three primary water-bearing layers in this aquifer system (Chicot, Evangeline, and Jasper) are made up of laterally discontinuous deposit sediments of gravel, sand, silt, and clay.The uppermost Chicot aquifer is comprised of Pleistocene-and Holocene-age sediments, and the Evangeline aquifer is composed of Miocene-and Pliocene-age sediments.The Jasper aquifer and the Burkeville confining layer consist of Miocene-age sediments [62,64].Since 1890s, groundwater withdrawals from both Chicot and Evangeline aquifers have been the major water-source of for development in Harris, Galveston, Fort Bend, Montgomery, and Brazoria Counties, while the Jasper aquifer is also source of water in Montgomery County and in north Harris County [65].Groundwater withdrawal from the deeper Jasper aquifer has been increasing since 2000 as urban growth spreads northward [66].A majority of wells serving as freshwater resources have been drilled to the depth of 300~700 m in the underground aquifers, where declining of pumping-well water tables may results in compacting of the fine-grained clay layers within the aquifer [67,68].

InSAR Datasets and Processing
InSAR technique has the capability of measuring land surface displacement with a centimeterto millimeter-level at a spatial resolution of a few meters to tens of meters over a broad scale.The conventional individual interferograms analysis approach has been widely and successfully applied in the researches of land subsidence investigation, however the reliability of conventional InSAR method was limited by the temporal and/or spatial decorrelation between repeated acquisition of images as well as the orbital and atmospheric phase artifacts (e.g., [48,[69][70][71][72]).MTI methods, which currently include both Persistent Scatterer InSAR (PSInSAR) (e.g., [71]) and Small BAseline Subset (SBAS) algorithm (e.g., [73]) techniques, overcome the above mentioned limitations through utilizing series of radar images.PSInSAR approach generates interferograms between the selected master image, on account of favorable Doppler, spatial and temporal baselines, and all the other images, which discerns coherent persistent scatterers from the incoherent subscriptions, thereby getting physically meaningful observations.PSInSAR method has been successfully applied in city regions where the resolution elements are principally dominated by strong reflecting manufactured constructions and [71,74].However, this technique is restricted to phase pixels that are constant over time with adequate coherence, which generally results in sparse PS in non-urban regions.The SBAS processing forms interferograms between SAR acquisitions with short temporal intervals and spatial baselines to minimize the decorrelation, and detects pixels decorrelate little over short time intervals after filtering, which have been termed as slow varying filtered phase pixels (SFP) [75].The PS pixels and SFP pixels form two partitions, but overlapping sets of pixels.By combining both PS and SBAS InSAR techniques maximizes the spatial density (both PS and SFP pixels) of usable signal.Conventional InSAR processing and their applications over the GH face two main problems: decorrelation (over non-urban areas) and tropospheric artifacts.Few to no studies are conducted for the nonurban areas, which occupy more than 70% of the GH, with substantial amounts of faults, water/oil/gas wells, and distributed salt domes.
To study faulting of GH, we used the MTI approach proposed by [75] incorporating both PS and SBAS InSAR to maximize interferometric phase correlation.Details of this method can be found in [74][75][76].The GH region is covered by 23 images from two adjacent ascending L-band ALOS PALSAR tracks (wavelength is 23.6 cm): 175 (12 scenes, date range: 20070926-20110104) and 176 (11 scenes, date range: 20070713-20110121).The ascending radar satellite travels with a heading of −10.3 • and incidence angle of about 38 • at the center of the track.The topographic phase was removed from the interferograms using a simulated external DEM from ~30 m SRTM [77].In order to avoid independent interferogram clusters, consequently ensure the continuity of time-series analysis, we included some interferograms with larger spatial and/or temporal baselines and constructed 68 interferograms totally (Figure 2) using Doris software [78].The two datasets were then processed using StaMPS MTI method respectively.We utilized Toolbox for Reducing Atmospheric InSAR Noise (TRAIN) package to produce near real-time water-vapor measurements for select SAR images [79].We then estimated the atmospheric artifact and baseline error by high-pass filtering of the unwrapped phase data in time, followed by low-pass filtering in space [74].

InSAR-Derived Annual Deformation Rate
We processed SAR images from ALOS PALSAR datasets utilizing the MTI method described in Section 3 and generated two annual line-of-sight (LOS) deformation maps (track 175 and 176) of GH area during 2007 and 2011.InSAR measurements from two neighboring ALOS tracks were then averaged and mosaicked to show the overall distribution of deformation.An oblique corner shape deformation pattern with its corner at Katy, TX can be detected, including two large-scale subsidence zones (Figure 3a).The main broad-scale subsidence zone covers an elongated northeastward area of roughly 80 km by 30 km, which spans from Katy (Harris County), across Jersey Village (Harris County), to The Woodlands (southern Montgomery County) and Spring (northern Harris County).This subsidence cone centered at Spring and The Woodlands with a maximum subsidence rates of about 30 mm/yr (Figure 3a).The other broad-scale deformation cone is elongated northwestward, extending from Arcola to near Katy (Fort Bend County).A relatively lower subsidence rate of about 15 mm/yr is observed over this feature during the period of 2007-2011.Meanwhile, the southeastern Harris County presents ~20 mm/yr uplifts during the same 2007-2011 timespan owing to the effective groundwater regulation implementation.According to [53], the center of maximum subsidence was migrating northeastward, transferring from Jersey Village (~53 mm/yr) in 1990s to Spring and The Woodlands after 2000, where a relatively stable rate of subsidence of about 30 mm/yr was observed from the 1990s to 2011 [53].Figure 3a shows coincident subsiding distributions and magnitudes between the observations from two ALOS tracks.As the same spatial-temporal cover of the two datasets, we compared the vertical deformation rates in the area of overlap from two ALOS displacement maps.The InSAR derived deformations along LOS direction of SAR, we then projected InSAR measurements to vertical direction referring to the corresponding incidence angles, by assuming the deformation is fully vertical.Figure 3b shows the difference map of two adjacent tracks while an inset on the upper-right corner of Figure 3b shows the corresponding histogram.The standard deviation of the whole common area is approximately 5.8 mm (Figure 3).
Figure 4 shows profiles of the LOS annual displacement (a,b) and cumulative displacement (c) along the cross-fault lines of P#P#' marked on Figure 3a.Sharp displacement gradient changes or discontinuous are evident across a majority of faults, and two sides of a fault normally have different surface displacements (Figure 4).Fault activities during 2007-2011 are not as significant as that before 2000 at the Long Point fault system and Addicks fault system, and [53] concluded that the rate of creep is about 10-15 mm/yr.The faults in the northern Greater-Houston were still active, with a creep rate of 10-20 mm/yr at Hockley fault system and about 10 mm/yr at the Conroe fault system (Figure 4). Figure 3a shows coincident subsiding distributions and magnitudes between the observations from two ALOS tracks.As the same spatial-temporal cover of the two datasets, we compared the vertical deformation rates in the area of overlap from two ALOS displacement maps.The InSAR derived deformations along LOS direction of SAR, we then projected InSAR measurements to vertical direction referring to the corresponding incidence angles, by assuming the deformation is fully vertical.Figure 3b shows the difference map of two adjacent tracks while an inset on the upper-right corner of Figure 3b shows the corresponding histogram.The standard deviation of the whole common area is approximately 5.8 mm (Figure 3).
Figure 4 shows profiles of the LOS annual displacement (a,b) and cumulative displacement (c) along the cross-fault lines of P#P#' marked on Figure 3a.Sharp displacement gradient changes or discontinuous are evident across a majority of faults, and two sides of a fault normally have different surface displacements (Figure 4).Fault activities during 2007-2011 are not as significant as that before 2000 at the Long Point fault system and Addicks fault system, and [53] concluded that the rate of creep is about 10-15 mm/yr.The faults in the northern Greater-Houston were still active, with a creep rate of 10-20 mm/yr at Hockley fault system and about 10 mm/yr at the Conroe fault system (Figure 4).

Fault Traces Identification by InSAR Deformation
All InSAR observations from ALOS, Envisat and ERS indicate that two blocks of a fault always present different movements, which may weaken or modify the local displacement trend, resulting steep subsidence gradients across the mapped faults positions (Figures 3 and 4) [53].The measured fault-related gradients are products of continuous differential erosion and/or deposition across fault traces but are hardly noticed on sites.However, reports from local residents in The Woodlands (Montgomery County) area claimed damages on their property due to faulting that has become more evident during the last decade (e.g., [46,47]).
From the long term velocity maps we can see numerous sharp phase discontinuity along faults traces over the GH region, especially in the northern GH (i.e., north Harris County and south Montgomery County), seen as sharp color contrasts (pointed by white arrows in Figure 5a,b).These areas are covered by both neighboring ALOS tracks (i.e., overlap of two rectangles in Figure 1).Both independent tracks show remarkably consistent velocity steps and gradients (Figure 5a,b).Based on the inference that steep phase gradients and/or discontinuities should be caused by the ongoing differential movement of two sides of faults, fault lines can be depicted along the steep phase gradients and/or discontinuities as fault traces.At least three broad-scale active faulting zones can be imaged over northern GH region as outlined by the purple dashed rectangles marked on Figure 5c: the Hockley fault System, the Big Barn fault System and the Conroe Fault System (from south to north) (Figure 5), all within the larger Hockley-Conroe Fault System that runs northeast from the town of Hockley to the city of Conroe.The Hockley Fault System was identified by [13] as a continuously active growth fault with a dip angle of 60-70 degrees toward the Gulf Coast [23,36].As one of the most active faults in GH, the Hockley fault System is more than 40 km long, extending from the Hockley salt dome, across Highway 290 and terminating northeast of Hufsmith, Texas [19,[23][24][25]36].At least 5 segments of the Hockley fault System could be mapped by InSAR (Figure 5c), including the Hockley fault that has been studied by many geophysical engineers (e.g., [23][24][25]28,34,36].The Big Barn Fault System in Montgomery County, Texas is located approximately 20 miles north of Houston, Texas, extending from the very end of the Hockley Fault

Fault Traces Identification by InSAR Deformation
All InSAR observations from ALOS, Envisat and ERS indicate that two blocks of a fault always present different movements, which may weaken or modify the local displacement trend, resulting steep subsidence gradients across the mapped faults positions (Figures 3 and 4) [53].The measured fault-related gradients are products of continuous differential erosion and/or deposition across fault traces but are hardly noticed on sites.However, reports from local residents in The Woodlands (Montgomery County) area claimed damages on their property due to faulting that has become more evident during the last decade (e.g., [46,47]).
From the long term velocity maps we can see numerous sharp phase discontinuity along faults traces over the GH region, especially in the northern GH (i.e., north Harris County and south Montgomery County), seen as sharp color contrasts (pointed by white arrows in Figure 5a,b).These areas are covered by both neighboring ALOS tracks (i.e., overlap of two rectangles in Figure 1).Both independent tracks show remarkably consistent velocity steps and gradients (Figure 5a,b).Based on the inference that steep phase gradients and/or discontinuities should be caused by the ongoing differential movement of two sides of faults, fault lines can be depicted along the steep phase gradients and/or discontinuities as fault traces.At least three broad-scale active faulting zones can be imaged over northern GH region as outlined by the purple dashed rectangles marked on Figure 5c: the Hockley fault System, the Big Barn fault System and the Conroe Fault System (from south to north) (Figure 5), all within the larger Hockley-Conroe Fault System that runs northeast from the town of Hockley to the city of Conroe.The Hockley Fault System was identified by [13] as a continuously active growth fault with a dip angle of 60-70 degrees toward the Gulf Coast [23,36].As one of the most active faults in GH, the Hockley fault System is more than 40 km long, extending from the Hockley salt dome, across Highway 290 and terminating northeast of Hufsmith, Texas [19,[23][24][25]36].At least 5 segments of the Hockley fault System could be mapped by InSAR (Figure 5c), including the Hockley fault that has been studied by many geophysical engineers (e.g., [23][24][25]28,34,36].The Big Barn Fault System in Montgomery County, Texas is located approximately 20 miles north of Houston, Texas, extending from the very end of the Hockley Fault to its terminus at Interstate 45.At least 5 segments of the Big Barn Fault System could be imaged by InSAR (Figure 5c), which have not been formally named by USGS in their database but are mentioned and/or studied by several researchers [31,37,80].Extensive deformation and damages to nearby businesses and residences within the vicinity of the fault were reported by local residents within the last decade, but no fault activities were reported in this area before.The full Conroe Fault System has not been reported or studied by other researchers, and two faults with a length of 34.3 km and 11.8 km were identified by InSAR maps for the first time, running northeast from the northeast of Magnolia to the city of Conroe.In addition, two segments of the InSAR-mapped faults (within the top purple rectangle in Figure 5c) were also shown by [28], however, without names and discussions by them (Figure 5d).While the Conroe Fault System has not been formally named by USGS in their database, [80] referred one segment (the top dashed pink line on Figure 5d) as the Conroe fault.The agreement between InSAR-mapped faults traces with those from [80], which indicates that the InSAR-discovered fractures constitute parts of the Big Barn Fault Systems and Conroe Fault Systems.So the discovered faults are referred as Conroe Fault System in this study.to its terminus at Interstate 45.At least 5 segments of the Big Barn Fault System could be imaged by InSAR (Figure 5c), which have not been formally named by USGS in their database but are mentioned and/or studied by several researchers [31,37,80].Extensive deformation and damages to nearby businesses and residences within the vicinity of the fault were reported by local residents within the last decade, but no fault activities were reported in this area before.The full Conroe Fault System has not been reported or studied by other researchers, and two faults with a length of 34.3 km and 11.8 km were identified by InSAR maps for the first time, running northeast from the northeast of Magnolia to the city of Conroe.In addition, two segments of the InSAR-mapped faults (within the top purple rectangle in Figure 5c) were also shown by [28], however, without names and discussions by them (Figure 5d).While the Conroe Fault System has not been formally named by USGS in their database, [80] referred one segment (the top dashed pink line on Figure 5d) as the Conroe fault.The agreement between InSAR-mapped faults traces with those from [80], which indicates that the InSAR-discovered fractures constitute parts of the Big Barn Fault Systems and Conroe Fault Systems.So the discovered faults are referred as Conroe Fault System in this study.

Validation of InSAR-Mapped Faults
Our InSAR deformation maps have successfully imaged significant LOS displacement anomalies across the known fault location (Figure 5).The discontinuity lines run just on the LiDAR-mapped fault traces at some points, indicating our results are basically consistent with the faults mapped by LiDAR (white lines in Figure 5d), but with the added benefits of low-cost and very wide coverage.The Hockley Fault was extensively investigated with a variety of geophysical methods to illuminate the subsurface character of this fault: seismic [28,34,81], magnetic [36], conductivity [36], gravity [28,34,36], GPR [24,25,28] and resistivity imaging [24,25,36].All methods imaged fault-like signatures, such as a steep slope over the fault location from the conductivity data, different magnetic values associated with the fault sides, and fault-related gravity anomalies, across the Hockley Fault (at the intersection with Highway 290), which all confirm the position of the

Validation of InSAR-Mapped Faults
Our InSAR deformation maps have successfully imaged significant LOS displacement anomalies across the known fault location (Figure 5).The discontinuity lines run just on the LiDAR-mapped fault traces at some points, indicating our results are basically consistent with the faults mapped by LiDAR (white lines in Figure 5d), but with the added benefits of low-cost and very wide coverage.The Hockley Fault was extensively investigated with a variety of geophysical methods to illuminate the subsurface character of this fault: seismic [28,34,81], magnetic [36], conductivity [36], gravity [28,34,36], GPR [24,25,28] and resistivity imaging [24,25,36].All methods imaged fault-like signatures, such as a steep slope over the fault location from the conductivity data, different magnetic values associated with the fault sides, and fault-related gravity anomalies, across the Hockley Fault (at the intersection with Highway 290), which all confirm the position of the Hockley Fault trace.Figure 6a displays a fault anomaly whose location is consistent with the Hockley Fault mapped by previous researchers at precisely the same location [36].
The Big Barn fault was first mentioned in a field investigation guide book by the Houston Geological Society [80].Fugro Consultants, Inc. (FCI) [31] and Tolunay-Wong Engineers, Inc (TWEI) [32] acknowledged the presence of faulting and mapped where the faulting occurred through their field investigations over the Big Barn Fault System.Figure 6b shows enlarged InSAR-derived deformation, where supposed fault traces, LiDAR-mapped faults and faults imaged by FCI and TWEI are superimposed.The dots show some surface scarps published by FCI and TWEI through their geological survey for the purpose of designing network of water pipelines.All published faults' information is consistent with the InSAR-mapped fault traces at either Big Barn fault or Egypt fault (Figure 6b).The authors of [37] conducted an integrated geophysical investigation of the Big Barn fault, including the use of electrical resistivity and gravity techniques, to delineate the fault and define the geology of the upper rock units.Six of the seven field sites in their study were superimposed onto the enlarged InSAR deformation map, and sharp contrasts in InSAR deformation map from the upthrown side (smaller deformation) of the fault to the downthrown side (larger deformation) are evident at all the field sites (Figure 6c,d).Thus we confirmed that the discontinuities in the MTI-derived InSAR images do correspond to the surface ruptures.All geophysical techniques showed the different lithology features of two sides of the fault, i.e., Lissie Sand in higher accumulation on the faults' downthrown blocks and the Willis Clay in higher accumulation on the upthrown sides.The authors of [37] regarded these truncation features as normal faulting.
Hockley Fault trace.Figure 6a displays a fault anomaly whose location is consistent with the Hockley Fault mapped by previous researchers at precisely the same location [36].
The Big Barn fault was first mentioned in a field investigation guide book by the Houston Geological Society [80].Fugro Consultants, Inc. (FCI) [31] and Tolunay-Wong Engineers, Inc (TWEI) [32] acknowledged the presence of faulting and mapped where the faulting occurred through their field investigations over the Big Barn Fault System.Figure 6b shows enlarged InSAR-derived deformation, where supposed fault traces, LiDAR-mapped faults and faults imaged by FCI and TWEI are superimposed.The dots show some surface scarps published by FCI and TWEI through their geological survey for the purpose of designing network of water pipelines.All published faults' information is consistent with the InSAR-mapped fault traces at either Big Barn fault or Egypt fault (Figure 6b).The authors of [37] conducted an integrated geophysical investigation of the Big Barn fault, including the use of electrical resistivity and gravity techniques, to delineate the fault and define the geology of the upper rock units.Six of the seven field sites in their study were superimposed onto the enlarged InSAR deformation map, and sharp contrasts in InSAR deformation map from the upthrown side (smaller deformation) of the fault to the downthrown side (larger deformation) are evident at all the field sites (Figure 6c,d).Thus we confirmed that the discontinuities in the MTI-derived InSAR images do correspond to the surface ruptures.All geophysical techniques showed the different lithology features of two sides of the fault, i.e., Lissie Sand in higher accumulation on the faults' downthrown blocks and the Willis Clay in higher accumulation on the upthrown sides.The authors of [37] regarded these truncation features as normal faulting.[37] (the other site is beyond the scope of (b)).White lines represent the faults mapped by LiDAR [28], while newly revealed fractures by our MTI processing are in black lines.(e-g) show the photos taken in the field investigation, whose locations are marked as colored stars in (b) and in Figure 5c (green star for (e), light blue star for (f) and blue star for (g)).[37] (the other site is beyond the scope of (b)).White lines represent the faults mapped by LiDAR [28], while newly revealed fractures by our MTI processing are in black lines.(e-g) show the photos taken in the field investigation, whose locations are marked as colored stars in (b) and in Figure 5c (green star for (e), light blue star for (f) and blue star for (g)).
The Conroe Fault was mentioned in a field trip guide book published by the Houston Geological Society, in which [80] delineated the Conroe Fault according to a short surface rupture of about one kilometer.They indicated that there was no movement of the Conroe Fault during 1985 and 1987, while the rate of movement was 18 mm/yr for the year 1987.Damage to a swimming pool at the Conroe Aquatics Center (green star on Figure 5d), which is located on the InSAR-mapped fault line, was reported early in 2018 [82], indicating the activation of Conroe Fault.
The correlations between surface ruptures and subsurface faults activities have been confirmed at some faults positions over Houston region [6].A field survey was conducted in the study area, and visible ground cracks were evident at a considerable portion of the locations, in other areas the field expression was subtle and the presence of a fault was difficult to confirm. Figure 6e-g displays photos taken at places of the InSAR-observed phase discontinuities through our filed trip, whose locations have been marked as colorful stars on Figures 5c and 6b.The photograph (Figure 6e) taken at Hufsmith fault, viewed from the southwest shows that the downthrown side drops down relative to the upthrown one.The photographs (f) and (g), viewed from northeast and southwest, respectively, also indicate that the hanging wall moves downward compare to the upthrown one on the southeast, exhibiting the normal fault properties.The field investigation proves the existence of surface changes occurring along the Hockley and Big Barn faults as identified in our InSAR deformation map, where the LiDAR data were not applicable.InSAR-derived deformation measurements at 6 locations were compared with continuous GPS measurements.InSAR could only measure a projection vector of the three dimensional displacements on the LOS direction (i.e., deformation away from or towards the satellites), thus we projected the 3-dimensional GPS measurements into the corresponding InSAR LOS direction according to their local incidence angles.We selected PS points situating within the scope of the 100 m (PAM13, PAM17, rod1 from track 175) or 200 m (PAM11, PAM18, PAM48 from track 176) of each GPS stations and averaged the InSAR measurements around each benchmarks to compare their LOS displacements (Figure 7).The time-series deformations from GPS and InSAR measurements agree well in both the magnitude and tendency with an average RMSE of about 9 mm between the two measurements.However, InSAR measured relative displacement and different spatial datum was chosen for the InSAR and GPS data processing, so we compared the relative displacement between GPS stations, locating at the hanging wall and footwall separately, to evaluate the intensity of fault activity.The GPS observations from stations PAM 11 and PAM 18 located, respectively, on the upthrown and downthrown block of the Hockley Fault System indicate a deferential vertical displacement of −13.8 mm/yr [28], while InSAR presents a deferential vertical displacement of −13.0 mm/yr between stations PAM 11 and PAM 18. DoD generated by two versions of airborne LiDAR measurement from 2001 and 2008 (with vertical accuracy of 11.6 cm) indicated the slip rate of Hockley Fault System was −10.9 mm/yr [28].

Rate, Extent, and Temporal Evolution of Growth Faulting
The InSAR time-series deformation results also indicate the continuous differential movement of Hockley and Conroe Fault System (Figure 4c).The most active fault, the Big Barn Fault, presents a high deformation gradient of about 80 mm during the ALOS period (3.5 years) from 2007 to 2011.The cumulative deformation gradients of other faults are around 40-80 mm.The influence area of each fault is around 0.3-1 km (Figure 4).The mapped faults generally parallel the Gulf Coast line, dipping from northwest to southeast, showing consistency with the scarps of the surrounding terrain and the regional faults pattern [1].

Modeling Faulting Parameters
The surface faulting over GH shows neither reverse-slip nor strike-slip motion, which are firmly dip-slip normal faults [3].The InSAR deformation displays typical fault anomaly with steep slope over the fault location (Figure 5).Many acknowledged surface faults in the GH could be traced to certain depths (about hundreds to thousands of meters), which are part of the Gulf Coast geologic structures rather than purely a surficial phenomenon [17].Geodetic data inversion provides an independent way to explore the subsurface faulting characters, and in this paper we attempted to generate a fault slip model using two independent InSAR vertical displacement measurements.The continuous mining of underground aquifer at Montgomery County formed new water-level decline cones, reflecting the nucleation of new fractures.Meanwhile the large scale subsidence field caused by groundwater withdrawal makes it difficult to separate fault-related deformation for the purpose of studying fault mechanism.So, over the Conroe Fault lines we could only select a 2-3 km buffer zone that is far from the water-level decline cones to minimize the effects of the large-scale subsidence caused by groundwater withdrawal on fault parameter modeling.We utilized the linear regression model in [53] and groundwater level changes contours at the Jasper aquifer for the period 2000-2011 from USGS to establish an approximate deformation trend surface induced by groundwater withdrawal to reduce the influence of induced subsidence on the fault slip model.
To improve calculative efficiency, we subsampled InSAR measurements with spacing of 50 m, and a two-segment fault model was discretized into 216 and 144 patches for F1 and F2 of a uniform size (500 m × 500 m).We used a finite dislocation model in the elastic half-space [83], hypothesizing a Poisson ratio of 0.25, to compute Green's function for each patch and utilized an iterative constrained least-squares algorithm, Steepest Descent Method (SDM) [84][85][86], to estimate the distributed slip model using InSAR measurements.The homogeneous elastic half-space model may not the best model to represent mechanisms of the observed fault creeping, but it could be a simple and efficient way to delineate the features of underground fault activities through the surface displacement.A smoothing parameter was employed to constrain the distribution of slipping to avoid unrealistic fluctuations of the inversion result.The optimal smoothing factor of 0.08 was resolved by the analysis of the tradeoff curve between the squared data misfit and the squared slip roughness [87].The slip rate angles of F1 and F2 were allowed to vary within a reasonable variation range of −93° to −87° (by assuming purely normal component) in order to avoid excessive increase of

Modeling Faulting Parameters
The surface faulting over GH shows neither reverse-slip nor strike-slip motion, which are firmly dip-slip normal faults [3].The InSAR deformation displays typical fault anomaly with steep slope over the fault location (Figure 5).Many acknowledged surface faults in the GH could be traced to certain depths (about hundreds to thousands of meters), which are part of the Gulf Coast geologic structures rather than purely a surficial phenomenon [17].Geodetic data inversion provides an independent way to explore the subsurface faulting characters, and in this paper we attempted to generate a fault slip model using two independent InSAR vertical displacement measurements.The continuous mining of underground aquifer at Montgomery County formed new water-level decline cones, reflecting the nucleation of new fractures.Meanwhile the large scale subsidence field caused by groundwater withdrawal makes it difficult to separate fault-related deformation for the purpose of studying fault mechanism.So, over the Conroe Fault lines we could only select a 2-3 km buffer zone that is far from the water-level decline cones to minimize the effects of the large-scale subsidence caused by groundwater withdrawal on fault parameter modeling.We utilized the linear regression model in [53] and groundwater level changes contours at the Jasper aquifer for the period 2000-2011 from USGS to establish an approximate deformation trend surface induced by groundwater withdrawal to reduce the influence of induced subsidence on the fault slip model.
To improve calculative efficiency, we subsampled InSAR measurements with spacing of 50 m, and a two-segment fault model was discretized into 216 and 144 patches for F1 and F2 of a uniform size (500 m × 500 m).We used a finite dislocation model in the elastic half-space [83], hypothesizing a Poisson ratio of 0.25, to compute Green's function for each patch and utilized an iterative constrained least-squares algorithm, Steepest Descent Method (SDM) [84][85][86], to estimate the distributed slip model using InSAR measurements.The homogeneous elastic half-space model may not the best model to represent mechanisms of the observed fault creeping, but it could be a simple and efficient way to delineate the features of underground fault activities through the surface displacement.A smoothing parameter was employed to constrain the distribution of slipping to avoid unrealistic fluctuations of the inversion result.The optimal smoothing factor of 0.08 was resolved by the analysis of the tradeoff curve between the squared data misfit and the squared slip roughness [87].The slip rate angles of F1 and F2 were allowed to vary within a reasonable variation range of −93 • to −87 • (by assuming purely normal component) in order to avoid excessive increase of free parameters during the inversion.The strike directions of the two faults were directly derived from the surface rupture, LiDAR-mapped faults and InSAR-derived faults traces.We adopt a standard linear least-squares optimization method to find the best-fit slip distribution parameters that match the surface InSAR observations.The model fits the observed interferogram reasonably well with a RMS misfit of 5 mm.
We compared results derived by the estimated slip model with InSAR observations in Figure 8, and both independent InSAR tracks, i.e., track 176 (Figure 8a-c) and 175 (Figure 8d-f), show consistent fault parameters results.The overall picture of land surface deformation was well reconstructed through the best-fit slip distribution model (Figure 8b,e).The RMS of the residuals is ∼5 mm, which drops into the range of observations uncertainties.The result reveals slip rates of F1 and F2 are heterogeneous along the fault traces, which range from negligible dip motion to a peak value of 27 mm/yr with average slip rates of 13 mm/yr for F1 and 7 mm/yr for F2, respectively.The average locking depths are about 2 km (up to 4 km, F1) and 0.75 km (up to 2 km, F2) along the fault parallel direction indicating that the Conroe Fault is active within a very shallow locking depth (Figure 8g-h).While no published research exactly proves the modeled slip rates, there are good aspect of demonstrations on other segments of the Conroe Fault and surrounding faults (system).The authors of [3] concluded that the fault movements of GH were intermittent throughout any given year with an average rate of 12.7 mm/yr from 1966 to present, and indicated a slip rate of 18 mm/yr at the very east section of the Conroe Fault (top pink line on Figure 5d) in 1987 [80].The slip rate of the nearby Hockley Fault System was suggested to be 10.9 mm/yr (by differentiating the 2001 and 2008 LiDAR measurements) and −13.8 mm/yr (GPS observations) [28].GPR data indicated the nearby Hockley Fault System faulting to at least 300 m deep [28], and the nearby Big Barn Faults were thought to be active above depths of 1200 to 1500 m below ground surface through the study of well logs [3].The borehole and well logs, as well as the seismic reflection data delineated a majority of faults presented at 1000-4000 m below the surface in the GH region [17,88].
Remote Sens. 2019, 11, x FOR PEER REVIEW 14 of 23 free parameters during the inversion.The strike directions of the two faults were directly derived from the surface rupture, LiDAR-mapped faults and InSAR-derived faults traces.We adopt a standard linear least-squares optimization method to find the best-fit slip distribution parameters that match the surface InSAR observations.The model fits the observed interferogram reasonably well with a RMS misfit of 5 mm.We compared results derived by the estimated slip model with InSAR observations in Figure 8, and both independent InSAR tracks, i.e., track 176 (Figure 8a-c) and 175 (Figure 8d-f), show consistent fault parameters results.The overall picture of land surface deformation was well reconstructed through the best-fit slip distribution model (Figure 8b,e).The RMS of the residuals is ∼5 mm, which drops into the range of observations uncertainties.The result reveals slip rates of F1 and F2 are heterogeneous along the fault traces, which range from negligible dip motion to a peak value of 27 mm/yr with average slip rates of 13 mm/yr for F1 and 7 mm/yr for F2, respectively.The average locking depths are about 2 km (up to 4 km, F1) and 0.75 km (up to 2 km, F2) along the fault parallel direction indicating that the Conroe Fault is active within a very shallow locking depth (Figure 8g-h).While no published research exactly proves the modeled slip rates, there are good aspect of demonstrations on other segments of the Conroe Fault and surrounding faults (system).The authors of [3] concluded that the fault movements of GH were intermittent throughout any given year with an average rate of 12.7 mm/yr from 1966 to present, and indicated a slip rate of 18 mm/yr at the very east section of the Conroe Fault (top pink line on Figure 5d) in 1987 [80].The slip rate of the nearby Hockley Fault System was suggested to be 10.9 mm/yr (by differentiating the 2001 and 2008 LiDAR measurements) and −13.8 mm/yr (GPS observations) [28].GPR data indicated the nearby Hockley Fault System faulting to at least 300 m deep [28], and the nearby Big Barn Faults were thought to be active above depths of 1200 to 1500 m below ground surface through the study of well logs [3].The borehole and well logs, as well as the seismic reflection data delineated a majority of faults presented at 1000-4000 m below the surface in the GH region [17,88].

Discussion About the Causes of Fault Movement over GH
The MTI-derived deformation maps identify several SW-NE-orientation fault fractures between Hockley community and Conroe city, Texas, named Hockley-Conroe Fault System.Faulting in the GH area is mainly related to the development of the Gulf of Mexico, its regional faults, sedimentary history, salt movement and fluid extraction (oil and gas, and groundwater) [12,14].
Most of the geomorphic steps have been demonstrated as the surface expressions of known subsurface faults in the Gulf Coast, and the geological forces of nature play certain roles in actuating detected fault activities over northern GH [89].Origins of Gulf Coast faults can be attributed to the accumulation of thick sedimentary deposits and sediment masses slid basinward on top of mobile salt (and/or shale) substrata [90].The land surface in the GH area is almost entirely constituted of unconsolidated clays, clay shales, and poorly cemented sands of a young geological age, i.e., tens of thousands of years.Increased rates of sediment loading since Pleistocene time have apparently reactivated older fault systems over the study area [1].For instance, the topographic maps surveyed by USGS showed visible surface scarps induced by prehistoric fault movements along Eureka Heights fault and Long Point fault [91], which were regarded as natural process-triggered faults, given that the movement occurred before large-scale anthropogenic activity, such as fluid extraction, had affected the GH region.Large-scale forces, such as load-induced crustal warping, may also affect regional faults, and consequently play significant roles in fault movement [92].Three regional fault zones pass through the northern GH area: Yegua fault zone, Hitchcock fault zone, and the Wilcox Fault Zone, which runs across Montgomery County [3].
The regional faults might transfer stresses to neighboring region where are experiencing stresses, creating new fractures areas from some distance apart and stimulating movement of existing faults sections.The best proof of that is most of the mapped faults of GH tend to trend parallel to the coast, generally have greater recognized length and show a regional fault downthrown toward the coast [1,4,12,26,93].For example, the east-most section of the Conroe Fault has been attributed to the extensive, deep regional fault system [80].
The activation of these faults on the topographic surface may also have resulted from natural geologic processes such as salt movement [12,14,94,95].Thick-bed of evaporate minerals (primarily of salt) at certain depth make inroads into overlying rock layers, forming a large structure of active diapirism which is referred as salt dome [96].Salt tectonics not only plays important roles in petroleum geology, but also affects the regional groundwater flow.Large proportion of salt domes occurred at the U.S. Gulf Coast, with cross section diameters of 1-10 km and depths extending as far down as 6.5 km [19,97,98].The salt was deposited in Jurassic age that slowly move upwards, passing through the younger sediments, and reaching the earth surface in some cases [88,99].Salt domes are irrefutably among of the most important elements for inducing faulting, because almost 80% of the known faults in GH region exist over diapir of salt domes, with many of them showing a radial pattern [5,9,13,15,27].While some do display radial patterns around the position of salt domes, these faults show the general northeast-southwest trending orientations, which roughly parallel to the coast line of Gulf of Mexico [1,4,15,19,26,93].With regards to our study area, the Hockley Fault System locates between the Hockley salt dome at the west end and Tomball salt dome at the east end, while the Big Barn fault is situated at the southwest flank of the Conroe salt dome [80], suggesting the potential relationship between the fault activities and the rising salt domes.
The concept of faulting triggered by fluid extraction (such as oil, gas, and groundwater) induced subsidence has gained favor, and the faulting of GH was also linked to the fluid extraction in some literature since the 1970s (e.g., [6,16,29,100,101]).The state of stress on/around faults might be altered along with the pore pressure declines within the underground reservoirs that caused by groundwater and/or hydrocarbon extraction [102].There was always a connection between petroleum and faulting.Movement of faults was first described in 1926 as a result of oil extraction induced subsidence at the Goose Creek oil field on Galveston Bay [15].In northern GH, the east-most Conroe Fault was not only correlated to the deeper regional faults but also related to the hydrocarbon exploration activities at the Grand Lake-Risher Field, located at western Conroe city [80].No other faults in this region were reported to be directly linked to the hydrocarbon production activity.In addition, oil and gas exploration wells are sparsely distributed within our InSAR-mapped active faults region.
Substantial decline in the potentiometric surface could possibly be another contributing factor for acceleration of the faults slip rate over northern GH after 2000.Large scale groundwater removal can lower piezometric level, causing changes in the vertical effective stress in the vicinity of faults, thereby contributing to the acceleration of pre-existing faults motion [19].Since the 1890s, groundwater withdrawals from the Chicot and Evangeline aquifers have been the primary water source for development in Brazoria, Montgomery, Harris, Fort Bend, and Galveston Counties, while Jasper aquifer is also source of water for the whole Montgomery County and northern part of Harris County [65].The rates of piezometric surface decline in Chicot and Evangeline aquifers after Harris and Fort Bend Counties implemented groundwater withdrawal regulations in the 1970s.The groundwater levels were monitored by 13extensometers over GH, and all sites showed general rise of water level for the past 50 years [103].Accordingly, reports from USGS showed the consequent aquifer compaction rates have decreased at most extensometer sites, and even mildly recovered at some sites since 1990, addressing a delayed effect of 15 years on aquifer compaction associated with groundwater withdrawal over Houston region [43,53,65].Differential movements were observed across numerous mapped faults over Houston region as one fault side usually moves with respect to the other side, indicating the faulting activities correlate strongly with the InSAR-derived displacements in the GH area (Figure 4).Faults in northwestern Harris County were active during 1993-2000, with a rate (differential movement) of about 15 mm/yr at the Hockley fault system, 15-25 mm/yr at the Long Point fault system, and 40 mm/ year at the Addicks fault system [53].Little to no deferential movement was observed for the other Hockley-Conroe fault system in the 1990s.The observed fault activities were in according with the severe land subsidence cones over western GH in the 1990s [53].
Groundwater withdrawal from the Jasper aquifer has been increasing since 2000 as urban growth spreads northward, especially the extensive urban development in Montgomery County, and the total groundwater withdrawal continuously increased from 7.84 MG/d in 1976 to 43.6 MG/d in 2000, and then to 64.2 MG/d By 2010 [103].As a result, groundwater levels rose by 60 m in both the top two aquifers (i.e., Evangeline and Chicot aquifers), while groundwater levels declined by 60 m in the Jasper aquifer system from 1976 to 2015, and formed new water-level decline/subsidence cones over northern GH.The rate of faults slip dropped to about 5-20 mm/yr along the Addicks fault system and Long Point fault system from 2005 to 2011 on account of the uplifted groundwater level in Harris County and the subsidence cones migrated northward.Faulting receded or ceased in the water-level recovered areas with stable or uplift deformation, but sustain unabated over the ongoing subsidence areas with continued groundwater level declines, i.e., Hockley fault system.Reports from local residents in The Woodlands area claimed property damage due to faulting that has become more evident during the last decade (e.g., [31,47]).Faulting activities were in connection with the spatial distribution and density of water-level decline and ground subsidence.Figure 9 compares the groundwater level changes of the Jasper aquifer for the period 2000-2011 with the InSAR-derived deformation rate during 2007-2011.Subsidence near The Woodlands, where groundwater level declined most in a decade and Big Barn Fault System lies, is highly correlated with the spatial pattern of groundwater level declines (Figure 9).Not merely does accelerate activity of faulting in the time of excessive pumping of groundwater, but also the differential movements along faults in the region with maximum decline of piezometric surface are more significant than that within areas with lesser piezometric decline.
There is still no uniform consensus of origin(s) for surface faulting over northern GH.In general, regional faulting, salt tectonics, as well as subsidence caused by underground fluids withdrawal are suggested activation triggers for these faults, and all three mechanisms are often interactive with each other to an extent.However, according to the InSAR-discovered ground surface faulting over northern GH, the high spatial-temporal correlation between withdrawal of underground fluids and the accelerated faulting activities is beyond dispute.Sharp displacement gradients along fault traces decrease the surface tension of soil and facilitate the relative movement between two blocks of faults.Meanwhile, faults interrupt flow of groundwater, limit the horizontal expansion of localized subsidence zones, and accelerate depression and fault activities.Furthermore, the inverted fault locking depth described in Section 5 is slightly correlated with the groundwater withdrawal operations.Our study seems to validate that subsidence and related shallow subsurface fault activities in northern GH relates to mining of aquifers.There is still no uniform consensus of origin(s) for surface faulting over northern GH.In general, regional faulting, salt tectonics, as well as subsidence caused by underground fluids withdrawal are suggested activation triggers for these faults, and all three mechanisms are often interactive with each other to an extent.However, according to the InSAR-discovered ground surface faulting over northern GH, the high spatial-temporal correlation between withdrawal of underground fluids and the accelerated faulting activities is beyond dispute.Sharp displacement gradients along fault traces decrease the surface tension of soil and facilitate the relative movement between two blocks of faults.Meanwhile, faults interrupt flow of groundwater, limit the horizontal expansion of localized subsidence zones, and accelerate depression and fault activities.Furthermore, the inverted fault locking depth described in Section 5 is slightly correlated with the groundwater withdrawal operations.Our study seems to validate that subsidence and related shallow subsurface fault activities in northern GH relates to mining of aquifers.

Conclusions
The MTI method combining aspects of both PSInSAR and SBAS shows the extraordinary ability of identifying and monitoring faulting activity though detailed analysis of displacement maps.We augmented the MTI technique with GPS data and deformation modeling to delineate the extent and study the mechanisms of the faulting over northern GH.First, the InSAR results enabled us to identify and map fault fractures according to the steep phase gradients and/or discontinuities caused by differential movement of two blocks of faults from InSAR velocity maps.Three broad-scale active faulting zones located within the Hockley-Conroe Fault System (runs northeast from the town of Hockley to the city of Conroe) were imaged in detail: The Hockley fault System, the Big Barn fault System and the Conroe Fault System.Second, InSAR-mapped fault positions were then compared to

Conclusions
The MTI method combining aspects of both PSInSAR and SBAS shows the extraordinary ability of identifying and monitoring faulting activity though detailed analysis of displacement maps.We augmented the MTI technique with GPS data and deformation modeling to delineate the extent and study the mechanisms of the faulting over northern GH.First, the InSAR results enabled us to identify and map fault fractures according to the steep phase gradients and/or discontinuities caused by differential movement of two blocks of faults from InSAR velocity maps.Three broad-scale active faulting zones located within the Hockley-Conroe Fault System (runs northeast from the town of Hockley to the city of Conroe) were imaged in detail: The Hockley fault System, the Big Barn fault System and the Conroe Fault System.Second, InSAR-mapped fault positions were then compared to the published fault maps [28,80], geophysical surveys [36,37], surface scarps by geological surveys [31,32], and field investigation photos; all available evidence supports the reliability of the approach.Third, the accuracy of our InSAR time-series subsidence was reported to be approximately 9 mm/yr by comparing with GPS observations.Fourth, we characterize the faulting activity over North GH by estimating the fault parameters (locking depth and slip rate) through modeling of InSAR-derived deformation maps.The SW-NE-oriented faults on ground surface pertain to normal-slip movement with an average slip rate of 7-13 mm/yr at locking depth of < 4 km.Finally, we reviewed and assessed vast quantities of information on the features and origins of faulting of GH, their obvious affiliation with growth

Figure 1 .
Figure 1.Geological map of GH[8], where colors represent generalized lithology, and black letters

Figure 1 .
Figure 1.Geological map of GH[8], where colors represent generalized lithology, and black letters (Ql-Lissie Formation; Pow-Willis Formation; Qb-Beaumont Formation; Qbs-Beaumont Formation (sand); Qbc-Beaumont Formation (clay); Qal-Alluvium; Qt-Fluviatile terrace deposits) and grey lines show distribution of the major geologic units[1].Faults (marked by black lines)[9,10], hydrocarbon wells (small grey dots[11]), GPS benchmarks (green triangles), salt dome positions (green polygons) and county boundaries (white dashed lines, county names are labeled by purple letters) are superimposed.The blue and black rectangles show coverage of two neighboring Advanced Land Observing (ALOS) data.The Texas State district is shown as an inset map on the upper-right corner.

Figure 2 .
Figure 2. SAR interferograms spatial-temporal baseline distributions: (a) shows PS-InSAR (single master) interferograms of both two tracks: blue segments indicate 11 interferograms generated by ALOS path 175, with 20071227 as the master image; black segments represent 10 interferograms generated by ALOS path 176, with 20100605 as the master image; (b,c) SBAS (multiple master) interferograms with temporal and perpendicular baseline thresholds of 350 days and 1200 m, respectively.

Figure 2 .
Figure 2. SAR interferograms spatial-temporal baseline distributions: (a) shows PS-InSAR (single master) interferograms of both two tracks: blue segments indicate 11 interferograms generated by ALOS path 175, with 20071227 as the master image; black segments represent 10 interferograms generated by ALOS path 176, with 20100605 as the master image; (b,c) SBAS (multiple master) interferograms with temporal and perpendicular baseline thresholds of 350 days and 1200 m, respectively.

Figure 3 .
Figure 3. (a) Annual line-of-sight (LOS) deformation map derived from ALOS datasets, which shows a mosaic of two ALOS tracks 175 and 176.The negative values represent ground surface movement away from the satellite.The pink box in (a) indicates our main study area used for the discussion of InSAR results in Figure 5; the dark red lines are profiles across the active faults (labeled as P#P#').(b)Vertical displacement difference derived from path 175 and 176 by assuming the observed InSAR deformation is vertical; the corresponding statistical histogram of (b) is an inset on the upper-right corner.

Figure 3 .
Figure 3. (a) Annual line-of-sight (LOS) deformation map derived from ALOS datasets, which shows a mosaic of two ALOS tracks 175 and 176.The negative values represent ground surface movement away from the satellite.The pink box in (a) indicates our main study area used for the discussion of InSAR results in Figure 5; the dark red lines are profiles across the active faults (labeled as P#P#').(b)Vertical displacement difference derived from path 175 and 176 by assuming the observed InSAR deformation is vertical; the corresponding statistical histogram of (b) is an inset on the upper-right corner.

Figure 4 .
Figure 4. Annual (a,b) and time-series (c) subsidence result from interferometric synthetic aperture radar (InSAR) along three profiles whose positions are shown on Figure 3 as P#P#'.Dark yellow lines show the corresponding surface height.The rough sketch of some integrated faults' geometry is shown at the bottom of (b).Green polygon in (a) shows approximate location of Tomball salt dome, whose underground depth is not clear.Dash lines represent fault locations.

Figure 4 .
Figure 4. Annual (a,b) and time-series (c) subsidence result from interferometric synthetic aperture radar (InSAR) along three profiles whose positions are shown on Figure 3 as P#P#'.Dark yellow lines show the corresponding surface height.The rough sketch of some integrated faults' geometry is shown at the bottom of (b).Green polygon in (a) shows approximate location of Tomball salt dome, whose underground depth is not clear.Dash lines represent fault locations.

Figure 5 .
Figure 5. Enlarged deformation maps over northwestern Houston, whose location is outlined by pink dashed rectangles in Figure 3, from two adjacent ALOS-1 paths: (a) Path 176; (b) Path 175; (c,d) Both Path 175 and Path 176.The white arrows in (a) and (b) represent zones of high deformation velocity gradient (differential movement).White lines represent the faults mapped by LiDAR [10] and pink dashed lines show the faults published by Norman and Elsbury (1991), while newly revealed fractures by our Multi-temporal InSAR (MTI) processing are in black lines.The purple boxes in (c) indicate the three discovered faults systems, while the middle one (the same as the bottom purple box in (d)) is also used for InSAR results validation in Figure 6b.The top purple box in (d) shows an area that will be used for fault activity model analysis in Section 5.The colored stars in (c) and (d) mark positions for field survey and validation in Figure 6a,c-g.

Figure 5 .
Figure 5. Enlarged deformation maps over northwestern Houston, whose location is outlined by pink dashed rectangles in Figure 3, from two adjacent ALOS-1 paths: (a) Path 176; (b) Path 175; (c,d) Both Path 175 and Path 176.The white arrows in (a) and (b) represent zones of high deformation velocity gradient (differential movement).White lines represent the faults mapped by LiDAR [10] and pink dashed lines show the faults published by Norman and Elsbury (1991), while newly revealed fractures by our Multi-temporal InSAR (MTI) processing are in black lines.The purple boxes in (c) indicate the three discovered faults systems, while the middle one (the same as the bottom purple box in (d)) is also used for InSAR results validation in Figure 6b.The top purple box in (d) shows an area that will be used for fault activity model analysis in Section 5.The colored stars in (c) and (d) mark positions for field survey and validation in Figure 6a,c-g.

Figure 6 .
Figure 6.(a) Enlarged deformation map of Hockley fault, whose location is shown as a pink star in Figure 5d.(b) Enlarged deformation map around Big Barn Fault System, whose location is outlined by the large purple dashed rectangle in Figure 5d.Scarp positions are from FCI and TWEI, and the faults they mapped are based on the scarp locations in the field.(c,d) show the enlarged deformation maps of Part 1 and Part 2 (labeled on (b)) of Big Barn faults, respectively, and the pink lines and numbers display six of the seven field sites of geophysical survey conducted in[37] (the other site is beyond the scope of (b)).White lines represent the faults mapped by LiDAR[28], while newly

Figure 6 .
Figure 6.(a) Enlarged deformation map of Hockley fault, whose location is shown as a pink star in Figure 5d.(b) Enlarged deformation map around Big Barn Fault System, whose location is outlined by the large purple dashed rectangle in Figure 5d.Scarp positions are from FCI and TWEI, and the faults they mapped are based on the scarp locations in the field.(c,d) show the enlarged deformation maps of Part 1 and Part 2 (labeled on (b)) of Big Barn faults, respectively, and the pink lines and numbers display six of the seven field sites of geophysical survey conducted in[37] (the other site is beyond the scope of (b)).White lines represent the faults mapped by LiDAR[28], while newly revealed fractures by our MTI processing are in black lines.(e-g) show the photos taken in the field investigation, whose locations are marked as colored stars in (b) and in Figure5c(green star for (e), light blue star for (f) and blue star for (g)).

Figure 7 .
Figure 7.Comparison between InSAR-mapped time-series ground surface displacement and GPS observations at 6 stations.

Figure 7 .
Figure 7.Comparison between InSAR-mapped time-series ground surface displacement and GPS observations at 6 stations.

Figure 8 .
Figure 8. Observed (a and d), modeled (b and e), and residual (c and f) average vertical deformation maps of the two independent tracks (i.e., 176 and 175) of InSAR data.The black line shows the location of the fault.(g and h) show slip rate distribution maps of fault F1 and F2, respectively.

Figure 8 .
Figure 8. Observed (a and d), modeled (b and e), and residual (c and f) average vertical deformation maps of the two independent tracks (i.e., 176 and 175) of InSAR data.The black line shows the location of the fault.(g and h) show slip rate distribution maps of fault F1 and F2, respectively.

23 Figure 9 .
Figure 9. InSAR-observed deformation rate from 2007 to 2011 and the contoured groundwater elevation change at Jasper aquifer for the period 2000-2011 [104].Black triangles show location of groundwater wells in Jasper aquifer that were used to derive the water level change counters.

Figure 9 .
Figure 9. InSAR-observed deformation rate from 2007 to 2011 and the contoured groundwater elevation change at Jasper aquifer for the period 2000-2011 [104].Black triangles show location of groundwater wells in Jasper aquifer that were used to derive the water level change counters.