Anatomy of Subsidence in Tianjin from Time Series InSAR

Groundwater is a major source of fresh water in Tianjin Municipality, China. The average rate of groundwater extraction in this area for the last 20 years fluctuates between 0.6 and 0.8 billion cubic meters per year. As a result, significant subsidence has been observed in Tianjin. In this study, C-band Envisat (Environmental Satellite) ASAR (Advanced Synthetic Aperture Radar) images and L-band ALOS (Advanced Land Observing Satellite) PALSAR (Phased Array type L-band Synthetic Aperture Radar) data were employed to recover the Earth’s surface evolution during the period between 2007 and 2009 using InSAR time series techniques. Similar subsidence patterns can be observed in the overlapping area of the ASAR and PALSAR mean velocity maps with a maximum radar line of sight rate of ~170 mm ̈ year ́1. The west subsidence is modeled for ground water volume change using Mogi source array. Geological control by major faults on the east subsidence is analyzed. Storage coefficient of the east subsidence is estimated by InSAR displacements and temporal pattern of water level changes. InSAR has proven a useful tool for subsidence monitoring and displacement interpretation associated with underground water usage.


Introduction
Subsidence caused by groundwater extraction is a worldwide problem [1][2][3].Estimated subsidence rates range up to several hundred mm per year, highlighting the importance of continuous displacement mapping.With great demand for fresh water and a limited river water resource, ground water use seems inevitable for Tianjin City (Figure 1a), which has a population of 12.93 million and is currently growing at 2.6% per year till 2010 [4].The maximum pumping rate recorded was in 1981 with the annual rate of 1.2 ˆ10 9 m 3 .The annual rates from 1981 to 2007 are in the range 0.6-0.8ˆ10 9 m 3 ¨year ´1 (Figure 1b) [5].As a consequence of this sustained high level of abstraction, an overall drop in the water table has been observed in Tianjin.Water extraction has increasingly been from deeper underground [6].This extensive use of ground water could account for the subsidence in Tianjin, which has been observed since 1959.The subsidence can be divided into three different stages: slow subsidence at 50 mm¨year ´1 from 1959 to 1977, accelerating subsidence from 1978 to 1982 at 100 mm¨year ´1 and attenuating subsidence at 20-30 mm¨year ´1 from 1983 to 1997 [7].The maximum total subsidence from 1959 to 2008 has reached 3.2 meters [5].Locally, higher rates have been recorded: in the 1980s, the maximum subsidence rate is given as 250 mm¨year ´1 in Tianjin [6].The current subsidence rate has been estimated as 8-56 mm¨year ´1 [6] and 30-40 mm¨year ´1 [5].
InSAR (Interferometric Synthetic Aperture Radar) has been used to map subsidence due to groundwater extraction in Tianjin [8][9][10] and other areas like Shanghai [11], Mekong delta [12], Ganges-Brahmaputra delta [13], Iran [14][15][16][17], Reykjanes Peninsula [18] and so on.In this study, an InSAR time series technique is also used to recover the surface displacement field over the urban area of Tianjin between 2007 and 2009.The patterns of observed subsidence are then interpreted in details and modeled.Subsidence is concentrated in particular areas of the city, and these subsidence centers are analyzed in terms of their geological controls and association with groundwater levels.The effects of groundwater volume changes on surface displacement are investigated for one subsidence center through inverse modeling of InSAR derived subsidence.
Remote Sens. 2016, 8, 266 2 of 22 for the subsidence in Tianjin, which has been observed since 1959.The subsidence can be divided into three different stages: slow subsidence at 50 mm•year −1 from 1959 to 1977, accelerating subsidence from 1978 to 1982 at 100 mm•year −1 and attenuating subsidence at 20-30 mm•year −1 from 1983 to 1997 [7].The maximum total subsidence from 1959 to 2008 has reached 3.2 meters [5].Locally, higher rates have been recorded: in the 1980s, the maximum subsidence rate is given as 250 mm•year −1 in Tianjin [6].The current subsidence rate has been estimated as 8-56 mm•year −1 [6] and 30-40 mm•year −1 [5].InSAR (Interferometric Synthetic Aperture Radar) has been used to map subsidence due to groundwater extraction in Tianjin [8][9][10] and other areas like Shanghai [11], Mekong delta [12], Ganges-Brahmaputra delta [13], Iran [14][15][16][17], Reykjanes Peninsula [18] and so on.In this study, an InSAR time series technique is also used to recover the surface displacement field over the urban area of Tianjin between 2007 and 2009.The patterns of observed subsidence are then interpreted in details and modeled.Subsidence is concentrated in particular areas of the city, and these subsidence centers are analyzed in terms of their geological controls and association with groundwater levels.The effects of groundwater volume changes on surface displacement are investigated for one subsidence center through inverse modeling of InSAR derived subsidence.

Tectonic Divisions
Tianjin Municipality is located 130 km southeast of Beijing and is on the western shore of Bohai Gulf.Tianjin lies on the North China Plain (NCP), also known as Huabei Plain.The NCP is dominated by Haihe and Cangdong fault systems (Figure 1a).The NNE-SSW oriented alignment of Cangdong fault, Tianjin fault, and the Paleogene lacuna, forms a horst-graben structure, which is transversely cut by WNW-ESE aligned faults including Baodi, Hangu, and Haihe faults [24].
In summary, Tianjin is separated by Baodi Fault into the northern Yanshan Folded Belt and the southern Bohai Bay Basin.The Bohai Bay basin contains three tectonic divisions: Jizhong depression, Cangxian uplift, and Huanghua depression.Paleogene lacuna separates Jizhong depression and Cangxian uplift.Cangdong fault separates the Cangxian uplift and the Huanghua depression [25].

Sediments
The stratigraphy in the Bohai Bay Basin is summarized in Figure 1c.The Bohai Bay Basin is a tectonically subsiding region filled with continental Tertiary and Quaternary sediments [21].The Quaternary and Tertiary sediments and aquifer system range in thickness from less than 1000 m in the Cangxian uplift [22] to 1200-9000 m in the Jizhong depression, and 900-5000 m in the Huanghua depression [20].

Hydrogeological Conditions
The ground water system in Tianjin has three main aquifers, the Quaternary, Tertiary and Lower Paleozoic-Sinian aquifers [22].Above the depth of 700 m in Tianjin, the ground water system can be sub-divided into six aquifer layers (Table 1) [26].Thickness of Aquifers I to V are in the range of 40-60, 25-60, 30-55, 30-50, and 20-30 m, respectively [6].Aquifer II is in an over consolidated late Pleistocene strata and Aquifer III is in a normally consolidated middle Pleistocene strata.In Tanggu district of Tianjin, even though aquifer II and III account for 46% and 24% of ground water exploitation, respectively, Aquifer III is the greatest deformation contribution source [7].Layers between 136 m and 300 m in depth account for more than 65% of subsidence by extensometer observation in Tanggu district of Tianjin [7].
The ground water extracted from a depth greater than 300 meters is mainly used for geothermal exchange.The Tianjin geothermal field covers approximately 11,000 km 2 [22].The horst where Tianjin is located is the most productive geothermal field and the highest measured geothermal gradients are found here.More than 200 geothermal wells, of which 10-15 are injection wells, are producing waters from 30 ˝C to boiling for room heating.Shallow geothermal wells range from 300 m in Quaternary formations to 600 m in Tertiary formations for lower temperatures.Deep wells can reach up to 2600 m in the Tertiary formation [22].

Data
In order to investigate the subsidence in Tianjin, an InSAR (Interferometric Synthetic Aperture Radar) method is used [28,29], exploring C-band Envisat-ASAR data from European Space Agency (ESA) and L-band ALOS-PALSAR data from Japan Aerospace Exploration Agency (JAXA).ASAR C band images (Table 2) are retrieved with a wavelength of 56 mm compared with PALSAR L band images (Table 3) of 236 mm.The radar line of sight look angles are 23 and 34 degrees off nadir for ASAR and PALSAR, respectively, which means a different sensitivity towards the vertical deformation.The azimuth angles are 166 and 13 degrees counter clockwise from north for ASAR and PALSAR, respectively.

Methods
Single look complex image was generated using ROI_PAC [28].Interferograms were generated using Doris [30].SRTM (NASA's Shuttle Radar Topography Mission) with a 3 arc-second resolution was used for topography removal [31].The StaMPS (Stanford Method for Persistent Scatterers) software was used to process the images using Persistent Scattering (PS) methods inside the package [32,33].
In contrast to most other PS methods, StaMPS uses phase spatial correlation to identify PS pixels instead of amplitude analysis [30,34,35].The advantage of this strategy is the capability to detect persistent scattering (PS)/slowly decorrelating filtered phase (SDFP) pixels with low amplitude, which is often the case in natural terrains.The probability for a pixel to be PS/SDFP is estimated through phase analysis, which is successively refined in a series of iterations.Without any prior assumption about the temporal nature of ground deformation, StaMPS relies on the spatial correlation of deformation rather than any assumption of the temporal dependence of deformation [36].
For selected PS points, the two dimensional phases are still modulo 2π.Phase unwrapping is therefore implemented in order to derive continuous displacement fields.On each interferogram, phase unwrapping is implemented using the SNAPHU (Statistical-Cost, Network-Flow Algorithm for Phase Unwrapping) algorithm [37].SNAPHU is a statistical-cost network flow algorithm.The algorithm aims to compute the most likely unwrapped solution with maximum posterior probability estimation given the observable input data.
The topography of the NCP (North China Plain) is quite gentle because it is formed in a coastal area by sediment aggradations.Large DEM errors are not expected here.The look angle error is estimated in time series using perpendicular baselines of each scene [32].
The atmospheric effects exhibit spatially correlated pattern.Deformation pattern are also spatially correlated.However, displacement usually stays in the same area while the atmospheric phase screen (APS) pattern is variable from one image to another.For a single pixel, the deformation should be continuous in time, APS will cause unexpected signal in time series.The APS can be estimated according to its different spatial and temporal nature from displacement.The master atmospheric effect and orbital error is estimated by its presence in every interferogram.The slave atmospheric effect and orbital error is estimated by temporal high pass and spatial low pass filtering [32].As a flood plain, our research area in NCP shows gentle topography with height ranges between ´16 and +68 m.As a result, the topography dependent APS [38][39][40], which could stay in the same area and be misunderstood as displacement, are believed to be small.Long-wavelength atmospheric effects and imprecise orbit errors were reduced by a planar phase gradient.

Mean Velocities from ALOS and Envisat
The mean phase values of each image are firstly used as the reference phase in time series analysis because no ground truth data are available.The area with obvious uplift displacements on the resulting mean velocity map was then chosen as the reference area (small triangle in Figure 2).The Envisat and ALOS mean velocity rates are not uniform across the study area with similar pattern of individual subsidence zones on both maps.Two subsidence areas can be identified from the mean velocity map (Figure 2).
The west subsidence area is located between Tianjin and its northwest neighbor city Langfang.The velocity swath in the west subsidence area exhibits similar pattern between ALOS and ENIVSAT measurements (Figure 3a), which suggest that the mean velocities in the west subsidence area are reliable.
The other east subsidence area is located around the Tuanbo Reservoir to the east of Jinghai County.The newly built Jing-Hu High Speed Railway also passes through this area (Figure 2).The mean velocity swath profile along JH Railway (Figure 3b) indicates the overlap of Jing-Hu High Speed Railway and the east subsidence area e.g., between 38.8 ˝and 39.05 ˝N.

Model of West Subsidence
The western subsidence is modeled using Mogi source solution in a semi-analytical approach [44,45].The RMS differences are 9.7 mm between ASAR LOS and leveling, 8.8 mm between PALSAR and leveling, and 8.3 mm between vertical component and leveling (Figure 4a-c).Slightly improvement is seen after decomposition of horizontal movement from vertical (Figure 4c,f).The reason lies in the fact that the accuracy of Leveling data, interpolated ALOS and Envisat displacement is limited.On the one hand, the accuracy of interpolated ALOS and Envisat displacement is limited by seasonal movement.For example, ALOS displacement on 1 March 2009 is interpolated from displacements of 22 January and 9 March 2009, while Envisat displacement on 1 March 2009 is interpolated from displacements of 27 February and 3 April 2009.The water level began to drop in late March, so the interpolated Envisat displacement is greater than ALOS in March 2009.It seems that the decomposition of horizontal and vertical movement is not highly effective in the presence of seasonal displacement.On the other hand, the leveling dates available are only accurate to months, and the first day of the month is assumed as the leveling date for a campaign.However, it still can be seen that interpolated ALOS and Envisat displacements generally follow the displacement trend of benchmarks, although with less displacement (Figure 4d,e).It is likely that the vertical displacement is dominant to horizontal displacement along the leveling line, and the projection from vertical to LOS results in displacement reduction for both ascending and descending tracks.

Model of West Subsidence
The western subsidence is modeled using Mogi source solution in a semi-analytical approach [44,45].
where u z px, yq is the vertical displacement at px, yq, ∆V i is the volume change at reservoir element i, g z pr, dq is the vertical green/influence function, ν. is the Poisson's ratio, r i is horizontal distance between surface point px, yq and reservoir element px i , y i q, d is reservoir depth.The reservoir is divided into equal squares.A semi analytical inversion is performed for Mogi source array to determine the best fit distribution of source volume change for LOS annual displacements (annual rates times one year).Modeled surface displacements are the ensemble contributions from each source.Poisson's ratio of 0.25 is set for this study.Reservoir depth is fixed at about 200 meters as ground water is pumped from a depth of 100-300 m in Tianjin [7].The linear least-squares inversion procedure has been adopted to estimate the source volume change.
It is well known that InSAR observations may contain global bias due to uncertainties in reference level or other uniform signals (e.g., sediment compaction).Accordingly, a constant offset is allowed in model inversion to accommodate the global bias.
In the process of inversion, regularization on the source strength is often needed to ensure a reliable inversion with a faithful representation of the source.Directly modeled reservoir volume changes (Figure 5a) seem noisier than the regularized source (Figure 5b).This can be caused by displacement uncertainty due to noise and other error sources.Besides, the semi analytical approach allows approximation of reservoir change and efficient computations, but it does not assure the physical process that must involve a continuous contracting/inflating volume.A Laplacian regularization based process was applied to smooth the source (Figure 5b).
where G is the Green's function, L is the Laplacian smooth operator, F is the smoothing factor (weight), and their product LF is the smoothing matrix used in model inversion.S is the source volume change, O is the offset term for global bias, and D is the observed displacement.In order to balance the roughness of source volume change and model fit to displacement, the L curve method (Figure 6) is adopted to find the best smoothing factor [46].
displacement uncertainty due to noise and other error sources.Besides, the semi analytical approach allows approximation of reservoir change and efficient computations, but it does not assure the physical process that must involve a continuous contracting/inflating volume.A Laplacian regularization based process was applied to smooth the source (Figure 5b).
where is the Green's function, is the Laplacian smooth operator, is the smoothing factor (weight), and their product is the smoothing matrix used in model inversion.is the source volume change, is the offset term for global bias, and is the observed displacement.In order to balance the roughness of source volume change and model fit to displacement, the L curve method (Figure 6) is adopted to find the best smoothing factor [46].  displacement uncertainty due to noise and other error sources.Besides, the semi analytical approach allows approximation of reservoir change and efficient computations, but it does not assure the physical process that must involve a continuous contracting/inflating volume.A Laplacian regularization based process was applied to smooth the source (Figure 5b).
where is the Green's function, is the Laplacian smooth operator, is the smoothing factor (weight), and their product is the smoothing matrix used in model inversion.is the source volume change, is the offset term for global bias, and is the observed displacement.In order to balance the roughness of source volume change and model fit to displacement, the L curve method (Figure 6) is adopted to find the best smoothing factor [46].The modeled source volume change is about ´7694 to +1678 m 3 for each element (300 m square) per year (Figure 7), equivalent to water volume change of ´85,000 to +18,000 m 3 /km 2 /year, further equivalent to water storage change of ´85 to +18 mm¨year ´1 in height.Unfortunately, to the best of our knowledge, in situ measurements of fluid volume from production wells is unpublished and is not available for comparison with our modeling results.Ground water storage (GWS) change from GRACE satellite measurement is about ´17 to ´22 mm¨year ´1 in North China Plain [47][48][49], equivalent to ´17,000 to ´22,000 m 3 /km 2 /year.The GRACE measurement covers Beijing, Tianjin, Hebei and Shanxi of ~370,000 km 2 , while our model is in an area of 860 km 2 in Tianjin.Simulated recoverable groundwater storage depletion is ´30,000 m 3 /km 2 /year (equivalent to 30 mm¨year ´1 in height) from 1970 to 2008 [50], equivalent to ´2700 m 3 for each element (source) per year.Contour line of ´2700 m 3 is superimposed on source volume change (Figure 7d).Line ´2700 m 3 coincides with sharp reduction of both the water extraction volume and subsidence observed.Modeled source volume change show local water storage increase up to 1678 m 3 .InSAR observations show subsidence of ´168.5 to ´2.9 mm in model area (Figure 7a).The volume increase seems to contradict with the fact that only subsidence is observed in this area.From the perspective of semi-analytical modeling, any surface deformation is the ensemble contribution of all model sources via their influence functions.Therefore, for a surface point, the superposition of a remote contracting source with large volume decrease, and a nearby inflating source with small volume increase, may still generate a surface subsidence.From the perspective of physical process, aquifer can benefit from precipitation infiltration, irrigation return, or other vertical or horizontal recharge.For instance, rainfall in summer recharge local ground water aquifers in North China Plain by 12-29 mm seen from GRACE and 25 mm seen from ground boreholes [48].
The modeled source volume change is about −7694 to +1678 m 3 for each element (300 m square) per year (Figure 7), equivalent to water volume change of −85,000 to +18,000 m 3 /km 2 /year, further equivalent to water storage change of −85 to +18 mm•year −1 in height.Unfortunately, to the best of our knowledge, in situ measurements of fluid volume from production wells is unpublished and is not available for comparison with our modeling results.Ground water storage (GWS) change from GRACE satellite measurement is about −17 to −22 mm•year −1 in North China Plain [47-49], equivalent to −17,000 to −22,000 m 3 /km 2 /year.The GRACE measurement covers Beijing, Tianjin, Hebei and Shanxi of ~370,000 km 2 , while our model is in an area of 860 km 2 in Tianjin.Simulated recoverable groundwater storage depletion is −30,000 m 3 /km 2 /year (equivalent to 30 mm•year −1 in height) from 1970 to 2008 [50], equivalent to −2700 m 3 for each element (source) per year.Contour line of −2700 m 3 is superimposed on source volume change (Figure 7d).Line −2700 m 3 coincides with sharp reduction of both the water extraction volume and subsidence observed.Modeled source volume change show local water storage increase up to 1678 m 3 .InSAR observations show subsidence of −168.5 to −2.9 mm in model area (Figure 7a).The volume increase seems to contradict with the fact that only subsidence is observed in this area.From the perspective of semi-analytical modeling, any surface deformation is the ensemble contribution of all model sources via their influence functions.Therefore, for a surface point, the superposition of a remote contracting source with large volume decrease, and a nearby inflating source with small volume increase, may still generate a surface subsidence.From the perspective of physical process, aquifer can benefit from precipitation infiltration, irrigation return, or other vertical or horizontal recharge.For instance, rainfall in summer recharge local ground water aquifers in North China Plain by 12-29 mm seen from GRACE and 25 mm seen from ground boreholes [48].The horizontal displacement can be modeled using horizontal Green's function and the best-fit source volume change.g r pr, dq " r where g r pr, dq is the horizontal green/influence function.Modeled horizontal displacements can be substantial, with maximum annual horizontal displacements reach up to 74 mm (Figure 8).Maximum horizontal displacements are located along a semi closed oval around the modeled source center.The modeled maximum horizontal displacement is 36% of the maximum vertical displacement.Theoretical ratio of maximum horizontal to maximum vertical displacement is 38% for a single Mogi source [51].A field based aquifer test in Nevada shows that the horizontal displacements reach 8 mm and vertical displacements reach 12 mm within the first 22 days of pumping before reaching the steady state pumping [52].Thus, the horizontal displacement induced by aquifer abstraction can be significant.Horizontal displacements are believed to be a cause of ground fissures [53].Ground fissures associated with subsidence has caused severe infrastructure damage in Taiyuan Basin of China [54].Hence, it is worth measuring the horizontal displacements due to ground water depletion, including Tianjin.

( + )
where ( , ) is the horizontal green/influence function.Modeled horizontal displacements can be substantial, with maximum annual horizontal displacements reach up to 74 mm (Figure 8).Maximum horizontal displacements are located along a semi closed oval around the modeled source center.The modeled maximum horizontal displacement is 36% of the maximum vertical displacement.Theoretical ratio of maximum horizontal to maximum vertical displacement is 38% for a single Mogi source [51].A field based aquifer test in Nevada shows that the horizontal displacements reach 8 mm and vertical displacements reach 12 mm within the first 22 days of pumping before reaching the steady state pumping [52].Thus, the horizontal displacement induced by aquifer abstraction can be significant.Horizontal displacements are believed to be a cause of ground fissures [53].Ground fissures associated with subsidence has caused severe infrastructure damage in Taiyuan Basin of China [54].Hence, it is worth measuring the horizontal displacements due to ground water depletion, including Tianjin.

Impact of Reservoir Grid Size
Grid-like patterns stand out in the modeled displacement (Figure 9b).The hypotheses that these are related with reservoir grid size are examined.The 300 m reservoir is now resampled into 270 m, 225 m, and 75 m grids (Figure 9c-e).At the same time, reservoir volume changes are resampled by size to fit the surface displacement.It can be seen that surface displacements are well represented when finer reservoir grids are adopted without alternating surface grid size.Hence, the grid patterns that obscure the modeled displacement are mainly due to the reservoir grid size in this case.
It should be noted that Mogi-type modeling assumes that source radius α is much smaller than source depth ((α/d)Ñ0).A radius of 1 km is usually assumed for volcano cases [51].The maximum displacement for a finite sphere is about 1 + (α/d) 3 times of that for a point sphere [51].Accordingly, the grid-like patterns appear when the ratio of cavity radius to depth is significant (Figure 9).Hence, water storage change can be underestimated when reservoir grid size is similar to its depth.However, fine reservoir grid causes quadratic growth of source numbers, greatly increasing computational burden in model inversion.
when finer reservoir grids are adopted without alternating surface grid size.Hence, the grid patterns that obscure the modeled displacement are mainly due to the reservoir grid size in this case.
It should be noted that Mogi-type modeling assumes that source radius α is much smaller than source depth ((α/d)→0).A radius of 1 km is usually assumed for volcano cases [51].The maximum displacement for a finite sphere is about 1 + (α/d) 3 times of that for a point sphere [51].Accordingly, the grid-like patterns appear when the ratio of cavity radius to depth is significant (Figure 9).Hence, water storage change can be underestimated when reservoir grid size is similar to its depth.However, fine reservoir grid causes quadratic growth of source numbers, greatly increasing computational burden in model inversion.

Tectonic Division and Its Effect on Subsidence
There are two kinds of explanations when fault act as a subsidence barrier.First, sediment deposits separated by the fault exhibit different compressibility [55].Second, ground water hydraulic conductivity is low around the fault [56].In this section, major faults in Tianjin are examined to see if they are related with subsidence observed.The tectonic division could help address the subsidence differences seen in Tianjin.
The east subsidence, outlined by B1 in Figure 10b, is controlled by Tianjinbei fault to the west, and Dasi fault to the east.Most of the east subsidence is located in the area of Shuangyao salient.
Velocity differences (abrupt InSAR mean velocity changes) can be found between the areas to the west and to the east of lines F1, F2, and F3 (The F lines in Figure 10b are drawn by visual inspection) (Figure 10b).Line F1 coincides with a fault to the west of Baitangkou depression [57].Line F2 matches the central segment of Dasi fault [58].Dasi fault is also known as Baitangkou fault.Line F3, however, bends away from Dasi fault and is oriented in NE-SW direction.It remains unclear if a fault exists along F3.Further geological investigation is needed to confirm this.
The Tianjinbei fault is striking NNE (strike is a line representing the intersection of a fault with a horizontal plane) and dipping NW (dip is the angle of descent of a fault relative to a horizontal plane).The Tianjinnan fault is also striking NNE while dipping SE.The Tianjinbei fault and the Tianjinnan fault are 6-7 km apart.Tianjinbei fault joins Haihe fault to the north in Tianjin city center [59].A velocity strip (a long and narrow zone with mean velocities different from neighboring area) matches the location of the Tianjinnan fault (Figure 10).Dislocated stratum appears at about 47 m in depth from TJ08 seismic profile of Tianjinnan fault [59].Boreholes TS1 and TS4 near profile TJ08 show that a 9.5 m fault throw (the vertical component of the dip separation) appears at 88.3-90.3m in depth to the west of Tianjinnan fault and at 98.9-99.7 m in depth to the east of Tianjinnan fault [60].
Haihe fault separates Tianjinbei fault into north and south segments.A velocity strip matches the north segment of Tianjinbei fault (Figure 10).In the north segment, dislocated stratum appears at 132 m and 125 m in depth from TJ01 and SY02 seismic profiles of Tianjinbei fault [59].Weak velocity differences can also be observed along line F4 (Figure 10b), which can be the south segment of Tianjinbei fault.
if a fault exists along F3.Further geological investigation is needed to confirm this.
The Tianjinbei fault is striking NNE (strike is a line representing the intersection of a fault with a horizontal plane) and dipping NW (dip is the angle of descent of a fault relative to a horizontal plane).The Tianjinnan fault is also striking NNE while dipping SE.The Tianjinbei fault and the Tianjinnan fault are 6-7 km apart.Tianjinbei fault joins Haihe fault to the north in Tianjin city center [59].A velocity strip (a long and narrow zone with mean velocities different from neighboring area) matches the location of the Tianjinnan fault (Figure 10).Dislocated stratum appears at about 47 m in depth from TJ08 seismic profile of Tianjinnan fault [59].Boreholes TS1 and TS4 near profile TJ08 show that a 9.5 m fault throw (the vertical component of the dip separation) appears at 88.3-90.3m in depth to the west of Tianjinnan fault and at 98.9-99.7 m in depth to the east of Tianjinnan fault [60].
Haihe fault separates Tianjinbei fault into north and south segments.A velocity strip matches the north segment of Tianjinbei fault (Figure 10).In the north segment, dislocated stratum appears at 132 m and 125 m in depth from TJ01 and SY02 seismic profiles of Tianjinbei fault [59].Weak velocity differences can also be observed along line F4 (Figure 10b), which can be the south segment of Tianjinbei fault.Cangdong fault is striking NNE and dipping SE.Tianjin segment of Cangdong fault is cut by Baodi fault and haihe fault [62].Velocity differences are seen along lines F6 and F7 (Figure 10b).Dislocated stratum appears at 152 m in depth from TJ23 seismic profile of Cangdong fault [59].
Haihe fault is striking NW to NWW and dipping SSW [59].Distinguishable minimum depth of Haihe fault is 60-70 m [59].It crosses the Tianjin urban area where uniform velocities are observed.
The east subsidence area shows greater subsidence rates than urban areas of Tianjin on both Envisat and ALOS mean velocity maps.This is consistent with the fact that the area of ground water extraction expanded to the suburban area since 1990s and the ground water pumping volume in urban area has declined from more than 100 million in 1981 to less than 10 million cubic meters in 2007 [5].
Within a geological uplift, the east subsidence is likely the result of heavy underground water use.Geological control of the east subsidence is seen along DSF and TJBF.Subsidence is likely to be related with stratum below 125 m seen from TJBF.

Water Level and Subsidence
In the east subsidence, the time series displacement of JH-12 Well is plotted against contemporary precipitation and its water level records of Aquifer II (Figure 11).It can be seen that seasonal InSAR displacements near JH-12 generally follows the underground water levels of JH-12.North China Plain is one of the largest agricultural production area in China, accounting for about 50% of wheat and 33% of maize grain production in China [63].The boot growth stage of winter wheat is April, when water demand by irrigation reaches its peak level.The water table was the lowest from May to June each year.The irrigation return and abundant precipitation in July and August infiltrated and recharged ground water [64].The water level then rose to the highest level in March or April the following year.The water level pattern of JH-12 resembles the lateral recharge-runoff-discharge pattern in the piedmont of North China Plain [66].However, JH-12 is at least 100 km away from the piedmont plain zone besides Yanshan Mountain and Taihang Mountain, and is located in alluvial fan and flood plain [67].In central North China Plain, water level is normally low in winter and spring with strong evaporation and less precipitation [66], while JH-12 show consistently high water level in winter and spring except a tiny fluctuation of ca.0.5 m around December (Figure 11).It is likely that JH12 is laterally recharged since Tuanbo Reservoir is nearby (Figure 12b).[68].The seasonal LOS deformation, satellite incidence angle, and water level drop, yield an elastic storage coefficient of 0.0125 for JH-12 from ALOS.Although 18 mm LOS subsidence is observed between 23 May and 1 August 2008 from Envisat, no ASAR image can be found in June or July, when subsidence reaches its maximum value seen from ALOS.
In a numerical study of flow system simulation in the North China Plain, storage parameter set for Aquifer III is in the range of 0.0001-0.005with geometric mean value of 0.0012, and storage parameter set for Aquifers I and II are in the range of 0.04-0.25 with geometric mean of 0.075 [50].The elastic storage near JH-12 Well triples the maximum storage parameter set for Aquifer III, but is merely one-third of the minimum storage parameter set for Aquifers I and II [50].Hence Aquifer II and other aquifer layers are likely contribute together to the observed surface subsidence.The water level pattern of JH-12 resembles the lateral recharge-runoff-discharge pattern in the piedmont of North China Plain [66].However, JH-12 is at least 100 km away from the piedmont plain zone besides Yanshan Mountain and Taihang Mountain, and is located in alluvial fan and flood plain [67].In central North China Plain, water level is normally low in winter and spring with strong evaporation and less precipitation [66], while JH-12 show consistently high water level in winter and spring except a tiny fluctuation of ca.0.5 m around December (Figure 11).It is likely that JH12 is laterally recharged since Tuanbo Reservoir is nearby (Figure 12b).[68].The seasonal LOS deformation, satellite incidence angle, and water level drop, yield an elastic storage coefficient of 0.0125 for JH-12 from ALOS.Although 18 mm LOS subsidence is observed between 23 May and 1 August 2008 from Envisat, no ASAR image can be found in June or July, when subsidence reaches its maximum value seen from ALOS.
In a numerical study of flow system simulation in the North China Plain, storage parameter set for Aquifer III is in the range of 0.0001-0.005with geometric mean value of 0.0012, and storage parameter set for Aquifers I and II are in the range of 0.04-0.25 with geometric mean of 0.075 [50].
The elastic storage near JH-12 Well triples the maximum storage parameter set for Aquifer III, but is merely one-third of the minimum storage parameter set for Aquifers I and II [50].Hence Aquifer II and other aquifer layers are likely contribute together to the observed surface subsidence.
Subsidence rates are compared with water level of Aquifer III and Aquifer IV in 2008 [26] (Figure 12).Notably, the water levels are relative because the depth of Aquifer III and Aquifer IV are in the range of 180-400 m for the east subsidence area in their paper.The water level in Aquifer III of the east subsidence area is ´70 to ´60 m in depth, while it is >´55 m in the urban area of Tianjin in 2008.Higher subsidence rates are associated with lower water level in Aquifer III.Aquifer III is likely to be a subsidence source to the east subsidence due to similarity in map patterns between subsidence and lower water levels.Correlation between the water level of Aquifer IV and subsidence pattern are not seen.
A wavelet analysis on InSAR displacement and water level time series is implemented (Appendix).Displacement periodicities and time shift between water level and displacement time series are analyzed.12).Notably, the water levels are relative because the depth of Aquifer III and Aquifer IV are in the range of 180-400 m for the east subsidence area in their paper.The water level in Aquifer III of the east subsidence area is −70 to −60 m in depth, while it is >−55 m in the urban area of Tianjin in 2008.Higher subsidence rates are associated with lower water level in Aquifer III.Aquifer III is likely to be a subsidence source to the east subsidence due to similarity in map patterns between subsidence and lower water levels.Correlation between the water level of Aquifer IV and subsidence pattern are not seen.
A wavelet analysis on InSAR displacement and water level time series is implemented (Appendix).Displacement periodicities and time shift between water level and displacement time series are analyzed.Geological control of Tianjinnan, Tianjinbei and Dasi faults on the east subsidence can be seen from InSAR displacement rates.The east subsidence can be related with Aquifer II, seen from temporal association of subsidence and water level of JH12.The east subsidence can also be related with Aquifer III, seen from elastic storage coefficient and spatial association of subsidence and water level of Aquifer III.The west subsidence is modeled using Mogi source array.Best-fit water extraction volume is obtained through model inversion.Estimated water volume change is ´85,000 to +18,000 m 3 /km 2 /year (equivalent to ´85 to +18 mm¨year ´1 in height).High gradients of subsidence and water storage change are associated with the recoverable water volume change of ´30,000 m 3 /km 2 /year (equivalent to ´30 mm¨year ´1 in height) from other numerical simulation study.
Geological control of Tianjinnan, Tianjinbei and Dasi faults on the east subsidence can be seen from InSAR displacement rates.The east subsidence can be related with Aquifer II, seen from temporal association of subsidence and water level of JH12.The east subsidence can also be related with Aquifer III, seen from elastic storage coefficient and spatial association of subsidence and water level of Aquifer III.
This study demonstrates the capability of time series InSAR to map subsidence, model ground water volume change, identify geological control factors, and estimate storage coefficients in Tianjin.InSAR can be an effective tool to investigate ground water use conditions in Tianjin and NCP.sampling of water level time series at image acquisition time results in very coarse power spectrum.Oppositely, the non-linear component of Envisat displacement is over sampled by linear interpolation at water level time.However, the interpolation also reduces the accuracy of displacements.Envisat time series show a two-month cycle (64-day period) between April and May 2008 from CWT analysis (Figure A1b).No effective cycle can be observed for water level time series (Figure A1d).XWT analysis shows a phase rotation from about 120 ˝out of phase in the one and half month cycle (48-day period) to about 120 ˝in phase in the two month cycle (64-day period) between April and May 2008 (Figure A1e).Abrupt phase changes occur between adjacent 48-day band and 64-day band.This implies strong seasonal variations in a short time.It corresponds to the quick water level drop in April.WTC analysis shows about 90 ˝out of phase in cycles of 4-8 days and 16-32 days (Figure A1f).It means the displacement leads water level by 1-2 days or 4-8 days.Actually, this is due to lower temporal resolution of Envisat images compared with water level time series.The lowest water level is between the 4th and 5th Envisat image acquisition time, so it looks as if displacement from 4th to 5th image happens before the water level rise.In fact, the 4th image corresponds to a water drop stage while the 5th image corresponds to a water rise stage.A1b).No effective cycle can be observed for water level time series (Figure A1d).XWT analysis shows a phase rotation from about 120° out of phase in the one and half month cycle (48-day period) to about 120° in phase in the two month cycle (64-day period) between April and May 2008 (Figure A1e).Abrupt phase changes occur between adjacent 48-day band and 64-day band.This implies strong seasonal variations in a short time.It corresponds to the quick water level drop in April.WTC analysis shows about 90° out of phase in cycles of 4-8 days and 16-32 days (Figure A1f).It means the displacement leads water level by 1-2 days or 4-8 days.
Actually, this is due to lower temporal resolution of Envisat images compared with water level time series.The lowest water level is between the 4th and 5th Envisat image acquisition time, so it looks as if displacement from 4th to 5th image happens before the water level rise.In fact, the 4th image corresponds to a water drop stage while the 5th image corresponds to a water rise stage.Wavelet tools are able to find the periodicities of InSAR time series and relative phase (time) shift between InSAR and water level time series.However, the wavelet power spectrum can be misleading if the temporal resolution of time series is low.Meanwhile, the phase (time) shift from continuous wavelet analysis enhances the robustness of time delay analysis between time series.

Figure 1 .
Figure 1.(a) The tectonic framework of the Northern China Plain.Boundary of North China Plain is in thick red line.Piedmont plain, alluvial fan and flood plain, and coastal plain in North China Plain are bounded by thick yellow and blue lines.Major faults and subsidiary faults in North China Plain are shown in thin black and red solid lines, respectively [19].Palaeogene lacuna (missing line of Paleogene system) is in thin dotted line, while further details of the fault system in Tianjin are shown by thin yellow solid lines [20].Locations of major cities in this area are shown.Envisat ASAR frame, ALOS PALSAR frame, and political boundary of Tianjin are outlined with dashed lines; (b) Annual groundwater pumping volume and rainfall in Tianjin City from 1980 to 2007 (after [5]); (c) Stratigraphy in Bohai Bay Basin (after [21-23]).The upper Tertiary formations can be divided into Neogene Minghuazhen (Nm) and Guantao (Ng) groups.The lower Tertiary formations can be divided into Eogene Dongying (Ed), Shahejie (Es), and Kongdian (Ek) formations.

Figure 1 .
Figure 1.(a) The tectonic framework of the Northern China Plain.Boundary of North China Plain is in thick red line.Piedmont plain, alluvial fan and flood plain, and coastal plain in North China Plain are bounded by thick yellow and blue lines.Major faults and subsidiary faults in North China Plain are shown in thin black and red solid lines, respectively [19].Palaeogene lacuna (missing line of Paleogene system) is in thin dotted line, while further details of the fault system in Tianjin are shown by thin yellow solid lines [20].Locations of major cities in this area are shown.Envisat ASAR frame, ALOS PALSAR frame, and political boundary of Tianjin are outlined with dashed lines; (b) Annual groundwater pumping volume and rainfall in Tianjin City from 1980 to 2007 (after [5]); (c) Stratigraphy in Bohai Bay Basin (after [21-23]).The upper Tertiary formations can be divided into Neogene Minghuazhen (Nm) and Guantao (Ng) groups.The lower Tertiary formations can be divided into Eogene Dongying (Ed), Shahejie (Es), and Kongdian (Ek) formations.

Figure 2 .
Figure 2. LOS (line of sight) rates from ALOS (17 January 2007 to 9 September 2009) and Envisat (4 January 2008 to 12 June 2009) satellites.Negative values mean moving away from the satellite (e.g., subsidence).The straight line "TJLF (Tianjin Langfang)" is located between Tianjin City and Langfang City.The poly line "HSR (High Speed Railway)" is the path of Jing-Hu high speed railway.The displacement rates along the two swathes are shown in Figure 3. BM2010 and BM2011, along with other benchmarks, are shown in white dots.Leveling results for all benchmarks are compared with InSAR displacements in Figure 4. Reference points for InSAR and leveling are shown in black triangle and dot respectively.Model area is outlined by rectangle.Counties and districts in Tianjin are outlined by thin dashed line.Geothermal fields with gradient above 3.5 (°C/100 m) are outlined by thick dashed line.

Figure 3 .
Figure 3. (a) LOS rates against longitude in TJLF swath with buffer distance of 500 meters to both sides of the line; and (b) LOS rates against latitude in a swath along the Jing-Hu High Speed Railway with buffer distance of 500 meters to the railway at both sides.

Figure 2 .
Figure 2. LOS (line of sight) rates from ALOS (17 January 2007 to 9 September 2009) and Envisat (4 January 2008 to 12 June 2009) satellites.Negative values mean moving away from the satellite (e.g., subsidence).The straight line "TJLF (Tianjin Langfang)" is located between Tianjin City and Langfang City.The poly line "HSR (High Speed Railway)" is the path of Jing-Hu high speed railway.The displacement rates along the two swathes are shown in Figure 3. BM2010 and BM2011, along with other benchmarks, are shown in white dots.Leveling results for all benchmarks are compared with InSAR displacements in Figure 4. Reference points for InSAR and leveling are shown in black triangle and dot respectively.Model area is outlined by rectangle.Counties and districts in Tianjin are outlined by thin dashed line.Geothermal fields with gradient above 3.5 ( ˝C/100 m) are outlined by thick dashed line.

Figure 2 .
Figure 2. LOS (line of sight) rates from ALOS (17 January 2007 to 9 September 2009) and Envisat (4 January 2008 to 12 June 2009) satellites.Negative values mean moving away from the satellite (e.g., subsidence).The straight line "TJLF (Tianjin Langfang)" is located between Tianjin City and Langfang City.The poly line "HSR (High Speed Railway)" is the path of Jing-Hu high speed railway.The displacement rates along the two swathes are shown in Figure 3. BM2010 and BM2011, along with other benchmarks, are shown in white dots.Leveling results for all benchmarks are compared with InSAR displacements in Figure 4. Reference points for InSAR and leveling are shown in black triangle and dot respectively.Model area is outlined by rectangle.Counties and districts in Tianjin are outlined by thin dashed line.Geothermal fields with gradient above 3.5 (°C/100 m) are outlined by thick dashed line.

Figure 3 .
Figure 3. (a) LOS rates against longitude in TJLF swath with buffer distance of 500 meters to both sides of the line; and (b) LOS rates against latitude in a swath along the Jing-Hu High Speed Railway with buffer distance of 500 meters to the railway at both sides.

Figure 3 .
Figure 3. (a) LOS rates against longitude in TJLF swath with buffer distance of 500 meters to both sides of the line; and (b) LOS rates against latitude in a swath along the Jing-Hu High Speed Railway with buffer distance of 500 meters to the railway at both sides.

4. 2 .
InSAR and Leveling Comparison Leveling data are available in May 2007, March 2008, September 2008 and March 2009 from 13 benchmarks along a 20 km section of the Jing-Hu High Speed railway in this area.InSAR results are compared with the leveling survey results.Comparison between InSAR and leveling is not straightforward due to different reference systems.Leveling reference point for Tianjin is the Liqizhuang bedrock point located in Tianjin institution of surveying and mapping.Liqizhuang bedrock point is a deep stratum benchmark, with its base 1088 m underground in Cambrian formations [41].The InSAR results are referenced to the mean value of phases first.Then, the areas with the greatest uplift signals are chosen as the InSAR reference area.Comparison between leveling and InSAR displacements can be made through a double difference in time and space [42].The spatial reference for InSAR and leveling comparison is leveling point BM2010, which exhibited the least subsidence among all leveling points.The temporal reference for InSAR and leveling comparison is set at the first leveling campaign in May 2007.Firstly, LOS InSAR displacements are linearly interpolated at the leveling times.Secondly, Interpolated LOS InSAR results are double differenced using the spatial and temporal reference.Finally, LOS Envisat and ALOS results are decomposed into east-west and vertical components by ignoring north-south horizontal movement for a polar orbit SAR satellite [43].To be consistent with InSAR results, the leveling results are re-referenced to BM2010 at May 2007 as well.It should be noted that first Envisat acquisition is later than first leveling campaign, so InSAR and leveling results are temporally re-referenced to March 2008 when necessary.Remote Sens. 2016, 8, 266 8 of 22

Figure 4 .
Figure 4. (a) Correlation between ALOS LOS displacements and four leveling campaigns; (b) correlation between Envisat LOS displacements and three leveling campaigns; (c) correlation between decomposed vertical displacements and three leveling campaigns; (d) interpolated ALOS displacements and leveling displacement; (e) interpolated Envisat displacements and leveling displacement; and (f) decomposed vertical displacement and leveling displacement.

Figure 4 .
Figure 4. (a) Correlation between ALOS LOS displacements and four leveling campaigns; (b) correlation between Envisat LOS displacements and three leveling campaigns; (c) correlation between decomposed vertical displacements and three leveling campaigns; (d) interpolated ALOS displacements and leveling displacement; (e) interpolated Envisat displacements and leveling displacement; and (f) decomposed vertical displacement and leveling displacement.

Figure 6 .
Figure 6.L curve to determine source smoothness.The horizontal axis denotes the modeled source smoothness logp||LS|| 2 q, where L and S are defined in Equation (4).The vertical axis denotes source fit to InSAR observations logp||GS ´D|| 2 q , where G and D are defined in Equation (4).The curve is formed by 13 different Laplacian regularized sources S.Each S corresponds to the best-fit result of each inversion.Smoothing strength for the source is controlled by LF in model inversion.There are 13 smoothing factors F range between 10 ´4 to 10 ´1 with exponential step of 0.25.Note that horizontal axis LS is calculated without multiplying smoothing factor F. The square centered with an asterisk represents that a value of 10 ´3 is employed as the smoothing factor for final result.

Figure 7 .
Figure 7. (a) ALOS LOS annual displacements (annual rates times one year); (b) modeled ALOS LOS annual displacement; (c) residuals between observed and modeled displacements for ALOS; (d) Envisat LOS annual displacements; (e) modeled Envisat LOS annual displacement; (f) residuals between observed and modeled displacements for Envisat; and (g) modeled Source array volume changes from ALOS displacement.Dashed line is a contour of ´2700 m 3 .

Figure 8 .
Figure 8. Modeled annual: (a) horizontal displacements; (b) horizontal displacements in E-W direction; and (c) horizontal displacements in N-S direction.

Figure 11 .
Figure 11.Envisat and ALOS LOS displacement at JH-12 Well (the triangle in Figure 12b) in east subsidence area.Precipitation data are available from China Meteorological Administration, specifically meteorology station No. 54527 in Tianjin [65].Water level of JH-12 Well in Aquifer II is available to the end of 2008 [27].
Water level of JH-12 declines sharply from −62.57m to −69.03 m between 23 March and 26 April 2008.Consequent subsidence between 21 April and 22 July in 2008 reaches about 63 mm in ALOS LOS direction.The short term aquifer compaction should be elastic as the water level recovers to a nearly identical level next march.Hence we evaluate the elastic storage coefficient S * ke near JH-12 by relating aquifer compaction and effective stress change.Aquifer compaction and stress change are reflected by ALOS displacement and hydraulic head change respectively

Figure 11 .
Figure 11.Envisat and ALOS LOS displacement at JH-12 Well (the triangle in Figure 12b) in east subsidence area.Precipitation data are available from China Meteorological Administration, specifically meteorology station No. 54527 in Tianjin [65].Water level of JH-12 Well in Aquifer II is available to the end of 2008 [27].
Water level of JH-12 declines sharply from ´62.57m to ´69.03 m between 23 March and 26 April 2008.Consequent subsidence between 21 April and 22 July in 2008 reaches about 63 mm in ALOS LOS direction.The short term aquifer compaction should be elastic as the water level recovers to a nearly identical level next march.Hence we evaluate the elastic storage coefficient S ke near JH-12 by relating aquifer compaction and effective stress change.Aquifer compaction and stress change are reflected by ALOS displacement and hydraulic head change respectively

Figure 12 .
Figure 12.Water levels in 2008 superimposed on LOS subsidence rates from Envisat (4 January 2008 to 12 June 2009): (a) Aquifer III; and (b) Aquifer IV.Location of JH-12 Well is marked in triangle.Tuanbo Reservoir is labeled.
Subsidence up to 170 mm•year −1 is observed in the Tianjin Municipality from Envisat/ASAR and ALOS/PALSAR images collected between 2007 and 2009.The subsidence rates are similar across the overlap area between Envisat and ALOS results.Subsidence along the JH High Speed Railway has also been verified by four leveling campaigns between 2007 and 2009 with a RMS difference of 8-9 mm between InSAR and leveling.The west subsidence is modeled using Mogi source array.Best-fit water extraction volume is obtained through model inversion.Estimated water volume change is −85,000 to +18,000 m 3 /km 2 /year (equivalent to −85 to +18 mm•year −1 in height).High gradients of subsidence and water storage change are associated with the recoverable water volume change of −30,000 m 3 /km 2 /year (equivalent to −30 mm•year −1 in height) from other numerical simulation study.

Figure 12 .
Figure 12.Water levels in 2008 superimposed on LOS subsidence rates from Envisat (4 January 2008 to 12 June 2009): (a) Aquifer III; and (b) Aquifer IV.Location of JH-12 Well is marked in triangle.Tuanbo Reservoir is labeled.
Subsidence up to 170 mm¨year ´1 is observed in the Tianjin Municipality from Envisat/ASAR and ALOS/PALSAR images collected between 2007 and 2009.The subsidence rates are similar across the overlap area between Envisat and ALOS results.Subsidence along the JH High Speed Railway has also been verified by four leveling campaigns between 2007 and 2009 with a RMS difference of 8-9 mm between InSAR and leveling.

Figure A1 .
Figure A1.(a) Detrend of Envisat time series at JH-12 Well; (b) continuous wavelet power of the Envisat time series; (c) water level of JH-12 Well; (d) continuous wavelet power of water level time series; (e) cross wavelet transform of Envisat and water level time series at JH-12 Well; and (f) wavelet coherence of Envisat and water level time series for JH-12.Images (g-l) are the same for ALOS and water level time series at JH-12 Well.The thick contour is the 95% confidence level.The cone of influence (COI) is shown in light shadow.

Figure A1 .
Figure A1.(a) Detrend of Envisat time series at JH-12 Well; (b) continuous wavelet power of the Envisat time series; (c) water level of JH-12 Well; (d) continuous wavelet power of water level time series; (e) cross wavelet transform of Envisat and water level time series at JH-12 Well; and (f) wavelet coherence of Envisat and water level time series for JH-12.Images (g-l) are the same for ALOS and water level time series at JH-12 Well.The thick contour is the 95% confidence level.The cone of influence (COI) is shown in light shadow.Common interval for ALOS and water level time series is from 17 January 2007 to 7 December 2008.CWT analysis of ALOS time series shows a four-month cycle (128 days) between the 8th (6 June 2008) and the 10th (6 September 2008) ALOS image.CWT analysis for water level time series also shows a 128-day cycle, though at a lower confidence level.XWT analysis show 90 ˝in phase in the 80-day band and 90 ˝out of phase in the 128-day band, from the 6th (20 January 2008) to the 9th (22 July 2008) ALOS image.This means strong seasonality occurs between the two images.WTC analysis shows about 90 ˝in phase in cycles of 4-8 days near the 7th ALOS image (21 April 2008).This means water level leads displacement by 1-2 days, which is beyond the 46-day repetition frequency of ALOS acquisition.WTC analysis also shows 90 ˝in phase in cycles of 64 days between 7th (21 April 2008) and 10th (6 September 2008) ALOS image, although with lower confidence level.It means water level leads displacement by 16 days.The displacement began to stabilize at the 10th (6 September 2008) ALOS image, while the water level's recent rise stops at 8 August 2008, which is about 29 days ahead of the

Table 1 .
Bottom depth of aquifers in Tianjin.Baodi fault and Plain of Tianjin are depicted in Figure1a.Tanggu district is shown in Figure2.

Table 2 .
Envisat/ASAR archive with acquisition time and length of perpendicular baseline relative to the reference image on 10 October 2008.

Table 3 .
ALOS/PALSAR archive with acquisition time and length of perpendicular baseline relative to the reference image on 7 December 2008.
[26]idence rates are compared with water level of Aquifer III and Aquifer IV in 2008[26](Figure