Investigating the Ground Deformation and Source Model of the Yangbajing Geothermal Field in Tibet , China with the WLS InSAR Technique

Ground deformation contains important information that can be exploited to look into the dynamics of a geothermal system. In recent years, InSAR has manifested its strong power in the monitoring of ground deformation. In this paper, a multi-temporal InSAR algorithm, WLS InSAR, is employed to monitor and characterize the Yangbajing geothermal field in Tibet, China, using 51 ENVISAT/ASAR images acquired from two overlapping descending tracks. The results reveal that the WLS InSAR algorithm can suppress the adverse effects of seasonal oscillations, associated with the freezing-thawing cycle of the permafrost in the Qinghai-Tibet Plateau. Deformations of up to 2 cm/yr resulting from the exploitation of the geothermal resource have been detected in the southern part of the Yangbajing field between 2006 and 2010. A source model inversion of the subsurface geothermal fluids was carried out based on the elastic half-space theory using the accumulated deformations. It was found that most geothermal fluid loss has occurred in the southern part of the shallow reservoir as the pore space beneath the northern part of field was recharged by the ascending flow from the deep layers of the reservoir through well-developed faults in the region.


Introduction
The Yangbajing geothermal field, located 94 km northwest from Lhasa City, the capital of the Tibet Autonomous Region, is characterized by the highest reservoir temperature in China.Since its initial exploration in 1976, the Yangbajing geothermal field has played an important role for the electricity power supply in Lhasa [1].It is widely acknowledged that the production from the geothermal reservoir has caused considerable surface deformation in the Yangbajing basin, which could cause environmental and structural hazards, such as damage to the power plants [2].However, the surface deformation provides useful information about the volume, pressure and temperature variations within the geothermal system, which makes it possible to study its dynamics [3].The relationship between ground deformation and geothermal dynamics is straightforward, i.e., stress change caused by injection or production of fluid leads to surface deformation.Therefore, it is essential to employ geodetic techniques for monitoring the surface deformation in order to minimize the associated damage and infer the dynamics of the geothermal field and the injection/extraction processes.
Recent decades have witnessed a significant development of geodetic tools for ground deformation measurements.In addition to ground-based measurements (e.g., leveling and GPS), active geothermal fields can now be monitored using InSAR technology and satellite data.The advantages of InSAR, e.g., large-scale, high spatial resolution, low cost, all-day and all-weather imaging capability, have been widely proven in monitoring ground deformation caused by natural disasters [4][5][6][7] and anthropological activities [8][9][10].In addition, many advanced InSAR algorithms have been developed since the beginning of this century [11][12][13][14][15].By analyzing a time series of SAR images acquired along the same track, we find that these advanced InSAR algorithms enhance the accuracy and reliability of the derived average deformation and enable us to measure the evolution of surface deformation over a long period of time.Many geothermal fields have also been successfully observed by InSAR.For instance, the deformation and seismicity in the Coso geothermal field in eastern California have been observed by processing 10 ERS-1/2 images acquired between 1993 and 1998 [16].The ground subsidence at the Tauhara and Ohaaki geothermal fields in New Zealand have been analyzed with 12 ALOS PALSAR images between 2007 and 2009 [17].The deformation and its temporal variation at the Geysers geothermal field in California has been resolved by using two distinct sets of InSAR data (i.e., ERS-1/2 and TerraSAR-X), which agreed well with the coupled numerical modeling results [18].
In this paper, we present a study of the Yangbajing geothermal field during the time period 2006-2010 by analyzing a set of ESA ENVISAT C-band SAR data collected from two overlapping descending tracks.The Yangbajing geothermal field has already been monitored by the SBAS algorithm with C-band SAR data and the PS algorithm with X-band SAR data, respectively [19,20].However, the seasonal oscillations of the permafrost, which have been detected in the Qinghai-Tibet Plateau due to the freezing-thawing cycle of the permafrost [21], were not previously considered in the retrieval of the displacements caused by geothermal exploitation.In this study, a WLS InSAR algorithm [22], a modified version of the SBAS algorithm, is used to distinguish the geothermal exploitation-related displacements and seasonal oscillations.Furthermore, the resolved deformation is used to characterize the source model, which can describe the behavior of the Yangbajing geothermal system.The source model parameters are searched by using a nonlinear optimization method (i.e., genetic algorithm), followed by a linear least square adjustment.

Study Area
The Yangbajing basin (Figure 1a), bordered by Pangduo Mountain to the southeast and Nyenchen Tonglha Mountain to the northwest, is one of the most developed and largest fault basins in the Qinghai-Tibet Plateau, with a length of ~90 km and an area of ~450 km 2 [23].The terrain of the basin is characterized by steep slopes in the southern area of the Zangbo River and gentle relief in the northwestern area, where several creeks have developed [2].The China-Nepal Highway and Zangbo River pass through the basin.There has been intensive tectonic activity in the basin.For instance, 30 earthquakes with a magnitude ě4.75 have occurred from 1921-1991 in the Yangbajing geothermal field [24].As no evidence shows any occurrence of Quaternary volcanic eruptions in the basin, the Yangbajing geothermal field belongs to the non-volcanic high-temperature fields in the Qinghai-Tibet Plateau [24].The basin is covered by different types of frozen soil, including seasonally frozen ground and permafrost [23].
The Yangbajing geothermal system includes two reservoirs at distinct depths [1,24,25].The shallow reservoir, divided into northern and southern parts by the China-Nepal Highway, covers an area of approximately 14.6 km 2 at a depth of 180 ~280 m (Figure 1b).The temperature in the shallow reservoir ranges from 150 ˝C to 160 ˝C, and the northern part is generally a little hotter than the southern part.As the primary electricity generation source, the shallow reservoir yields thermal water at a rate of 1200 tons/day to the Yangbajing geothermal power plant [24].The deep reservoir is a typical bedrock-fractured geothermal reservoir and covers an area of approximately 3.8 km 2 in the northern part of the field.The deep geothermal fluid is stored in fracture zones and tectonic fissures in the rocks, and migrates along the faults in the region.Downhole temperature measurements show that the deep reservoir is characterized by higher temperatures and has at least two layers [2].The upper layer, below a depth of at least 950 m, has a temperature of 250~275 ˝C and a thickness of more than 350 m.The lower layer is found below a depth of 1850 m and has a temperature higher than 300 ˝C.Although the shallow and deep reservoirs are distinguished depthwise by different temperatures, they belong to the same geothermal system.It has been confirmed that the occurrence of the Yangbajing geothermal field is due to a magmatic heat source [26].The recharge of the Yangbajing geothermal fluids occurs because of the snowmelt water from the Nyenchen Tonglha Mountains at altitudes of 4400~5800 m [1].The primary discharge paths were the hot springs before the power plant was built.Since the Yangbajing geothermal field started to produce electricity in 1977, intense exploitation dominates the discharge of the geothermal fluids [25].

of 14
The upper layer, below a depth of at least 950 m, has a temperature of 250~275 °C and a thickness of more than 350 m.The lower layer is found below a depth of 1850 m and has a temperature higher than 300 °C.Although the shallow and deep reservoirs are distinguished depthwise by different temperatures, they belong to the same geothermal system.It has been confirmed that the occurrence of the Yangbajing geothermal field is due to a magmatic heat source [26].The recharge of the Yangbajing geothermal fluids occurs because of the snowmelt water from the Nyenchen Tonglha Mountains at altitudes of 4400~5800 m [1].The primary discharge paths were the hot springs before the power plant was built.Since the Yangbajing geothermal field started to produce electricity in 1977, intense exploitation dominates the discharge of the geothermal fluids [25].

Datasets and Processing
A total of 51 SAR images acquired from two ENVISAT descending tracks are analyzed in the investigation of the Yangbajing geothermal field.The datasets span a time interval of four and a half years from April, 2006-September, 2010.The footprints of the two sets of SAR images can be seen in Figure 1a.Clearly, more than half of the frames overlap, providing an opportunity to conduct a quantitative comparison.The datasets include 29 scenes from Track 176 and 22 scenes from Track

Datasets and Processing
A total of 51 SAR images acquired from two ENVISAT descending tracks are analyzed in the investigation of the Yangbajing geothermal field.The datasets span a time interval of four and a half years from April 2006-September 2010.The footprints of the two sets of SAR images can be seen in Figure 1a.Clearly, more than half of the frames overlap, providing an opportunity to conduct a quantitative comparison.The datasets include 29 scenes from Track 176 and 22 scenes from Track 405, allowing us to generate 121 and 68 interferograms, respectively, after applying thresholds for the spatial-temporal baselines (Figure 2).The selected thresholds are a temporal baseline smaller than 700 days (approximately 2 years) and a perpendicular baseline shorter than 150 m.Precise orbital data provided by ESA are employed to refine the orbit vectors of the ENVISAT SAR images.
Remote Sens. 2016, 8, 191 4 of 14 405, allowing us to generate 121 and 68 interferograms, respectively, after applying thresholds for the spatial-temporal baselines (Figure 2).The selected thresholds are a temporal baseline smaller than 700 days (approximately 2 years) and a perpendicular baseline shorter than 150 m.Precise orbital data provided by ESA are employed to refine the orbit vectors of the ENVISAT SAR images.The standard two-pass differential InSAR approach is used to process each small baseline interferogram.The 90-m SRTM DEM data were used to remove the topographic contribution in the interferograms.A multi-look processing (i.e., 2 looks in range and 10 looks in azimuth directions) is adopted to suppress the decorrelation noise.In addition, the decorrelation noise in the generated differential interferograms is further minimized by using an improved Goldstein filter [27].After a reference pixel is selected, the differential interferograms are then unwrapped by a minimum-cost flow algorithm in order to retrieve the real phase at the coherent pixels.

WLS InSAR Estimation
The phase at a pixel in a differential interferogram, including the contributions from the surface deformation, topographic residuals, orbital errors, atmospheric artifacts, and remaining phase noises, can be written as: where ϕ is the phase in the differential interferogram, ∆d denotes the phase component due to the LOS surface deformation between the master and slave acquisitions, ∆ is phase component due to the topographic residuals, ϕ and ϕ are the phase components due to the inaccurate orbit and atmospheric delays, respectively (which usually behave as phase ramp throughout the whole image and can be removed by the polynomial fitting method [28,29]), ϕ is the phase component due to the remained noises, λ is the wavelength, is the perpendicular baseline, is the slant range between the sensor and the target, and θ is the incidence angle.
According to Equation (1), it is obviously difficult to resolve the deformation of interest from other contributions (e.g., the topographic residuals), if only a single differential interferogram were used.However, when many differential interferograms along the same track are available, a time series analysis can be conducted to estimate the LP temporal component of the surface deformation.This is how most current advanced InSAR algorithms work.Although the production-induced deformation can generally be characterized by gradual accumulation, the seasonal variations require consideration in this study due to the freezing-thawing cycle of the permafrost in the Qinghai-Tibet Plateau.Therefore, the LP deformation component can be modeled by the following superimposed function [22]: The standard two-pass differential InSAR approach is used to process each small baseline interferogram.The 90-m SRTM DEM data were used to remove the topographic contribution in the interferograms.A multi-look processing (i.e., 2 looks in range and 10 looks in azimuth directions) is adopted to suppress the decorrelation noise.In addition, the decorrelation noise in the generated differential interferograms is further minimized by using an improved Goldstein filter [27].After a reference pixel is selected, the differential interferograms are then unwrapped by a minimum-cost flow algorithm in order to retrieve the real phase at the coherent pixels.

WLS InSAR Estimation
The phase at a pixel in a differential interferogram, including the contributions from the surface deformation, topographic residuals, orbital errors, atmospheric artifacts, and remaining phase noises, can be written as: where φ int is the phase in the differential interferogram, ∆d los denotes the phase component due to the LOS surface deformation between the master and slave acquisitions, ∆ε is phase component due to the topographic residuals, φ orbit and φ atmo are the phase components due to the inaccurate orbit and atmospheric delays, respectively (which usually behave as phase ramp throughout the whole image and can be removed by the polynomial fitting method [28,29]), φ noise is the phase component due to the remained noises, λ is the wavelength, B K is the perpendicular baseline, R is the slant range between the sensor and the target, and θ is the incidence angle.
According to Equation (1), it is obviously difficult to resolve the deformation of interest from other contributions (e.g., the topographic residuals), if only a single differential interferogram were used.However, when many differential interferograms along the same track are available, a time series analysis can be conducted to estimate the LP temporal component of the surface deformation.This is how most current advanced InSAR algorithms work.Although the production-induced deformation can generally be characterized by gradual accumulation, the seasonal variations require consideration in this study due to the freezing-thawing cycle of the permafrost in the Qinghai-Tibet Plateau.Therefore, the LP deformation component can be modeled by the following superimposed function [22]: The first term on the right side of this equation is the linear deformation component, the second and third terms reflect the seasonal deformation component, ∆ t is the time interval between the master and slave acquisitions, T is the period of the seasonal deformation (one year in this study), and α 1 , α 2 and α 3 are the amplitudes of the base functions requiring estimation.Assuming there are m differential interferograms available along a track, a linear system with four unknowns and m observations can be constructed as: where Φ " " φ int,1 , φ int,2 , ¨¨¨, φ int,m ‰ T is the observation matrix, X " rα 1 , α 2 , α 3 , ∆εs T is the unknown matrix, and B is the design matrix This linear system can be resolved using the WLS adjustment when m ě 4: where Σ Φ reflect the stochastic properties of the observation matrix, composed by the variances of the differential interferograms.For each pixel, the variance estimate is based on the probability density function of the differential phase, which depends on the interferometric coherence and multi-look number [22].
The phase residuals, contributed by the HP temporal phase components and the unmodeled atmospheric phase components (i.e., turbulent atmosphere), can be separated from the phase measurements after removing the estimated LP deformation components and topographic residuals.A least-squares inversion is then applied to retrieve the time series of the phase residuals, corresponding to the SAR images, assuming its mean is zero.Subsequently, spatial median and temporal triangle filtering are carried out to separate the turbulent atmospheric phase components from the time series of the phase residuals.Thus the estimated turbulent atmospheric phase components and topographic residuals are subtracted from the differential interferograms.Finally, the refined differential interferograms are processed using the WLS adjustment again, in order to estimate the time series of surface displacements.

Source Model Inversion
Once the surface deformation due to the dynamics of the geothermal system is accurately measured, it should be inverted to retrieve the geothermal parameters at depth.In this study, we assume a volumetric source component for the Yangbajing geothermal system based on the elastic half-space theory [3].This means that the surface deformation at the geothermal field is caused by the movement of fluids and gases generated by the geothermal activity beneath the field.It should be noted that although not all tectonic stresses can be described, the volumetric source component is an acceptable assumption for a geothermal system [30].Therefore, when the Earth deforms elastically, a linear relationship between the surface displacement and the subsurface volume change can be constructed as [31]: where d i pxq is the ith component of surface displacement at pixel x (i " 1, 2, 3 correspond to the up, east, and north components, respectively), ∆v pyq is the fractional volume change of block y with a source volume V, and G i px, yq is the Green's function of the observation pixel x and source block y using a homogeneous elastic half-space assumption [32], the Green's function can be expressed as follows: where S " b px 1 ´y1 q 2 `px 2 ´y2 q 2 `px 3 ´y3 q 2 is the distance between the observation pixel x and the source block y, and ν is Poisson's ratio (0.25 in this study).
For more details about the presented forward modeling of the volumetric source component, see [33].However, it should be noted that other sources, such as finite rectangular source [34], nucleus-of-strain formulation [35] and its developers in geophysics [36] and reservoir engineering [37], can also be used to construct the relationship between the surface displacement and the subsurface dynamics.
Since the InSAR-derived LOS measurement is related to the up, east, and north deformation components according to the radar imaging geometry [38,39], Equation ( 6) can be rewritten as: where u 1 , u 2 and u 3 are the components of unit vector, projecting from the up, east, and north directions onto the LOS direction, respectively: where θ is the radar incidence angle and α is the azimuth angle (clockwise from north).Since a number of displacement measurements are provided by InSAR, the set of Equation ( 6) can be written into a vector-matrix form: D " GV, (10) where D is the observation vector containing the InSAR LOS displacement measurements at the surface points, and V is the unknown vector containing the fractional volume changes within a grid of subsurface blocks.Apparently, this system can be solved by the linearization approach (i.e., least-squares adjustment) if the design matrix G is accurately determined and the number of the observations is larger than that of the surface blocks.The latter condition can be easily satisfied, since we can divide the subsurface volume into a coarse gird (several times the InSAR pixel size) and a few layers (fewer than three layers).However, the former condition depends on the a priori geometrical information about the source volume, i.e., the depth and thickness of the geothermal system.We can employ a nonlinear optimization method (e.g., genetic algorithm) to estimate the best-fitted geometrical parameters in case they are not accurately provided.In order to balance the estimation accuracy and the computation burden of the nonlinear inversion we assume that the fractional volume changes are uniform at all blocks, so that only three unknowns (depth, thickness, and uniform volume change) require estimation.In addition, a quadtree down-sampling algorithm is applied to reduce the large number of InSAR point observations, while retaining sufficient information for inversion.
As the geometric parameters of the geothermal system are either provided in advance or estimated by the genetic algorithm, the fractional volume change at each block can be estimated by applying the least-squares adjustment to Equation (10).Note that the original sampling InSAR point observations are used in the linear inversion in order to make use of as many observations as possible.However, although the linear system is over-determined, the solution of the fractional volume changes could be unstable, due to the observation errors and the errors of design matrix.A roughness penalty should be included in the system to stabilize the least-squares inversion [3,40]: where Π is the penalty matrix that estimates the spatial derivative of the fractional volume change, and W r is a weighting factor controlling the relative importance of the roughness penalty term with respect to the data misfit term, which can be solved by the method of trial and error in a number of inversions [3].Therefore, continuity can be preserved for the estimated fractional volume changes, which ensures that the geothermal system represents a connected network rather than a collection of isolated blocks.The roughness penalty is however not necessary if the a priori knowledge about the reservoir is provided [41,42].
Besides the presented two-step scheme, the inversion can also be conducted by some integrated schemes (e.g., the Ensemble-based approaches [42,43]), which can avoid the possible inappropriate results yielded by the nonlinear inversion.The two-step scheme is preferred in this study since its applicability has been widely demonstrated in many geophysical inversion problems related to such seismic events.In addition, it is quite simple and efficient to conduct the two-step scheme.

Results
The linear LOS surface displacement velocity maps are first estimated by the WLS InSAR algorithm for Tracks 176 and 405, respectively, as shown in Figure 3a,b.Note that the exposed rock at the foot of the Nyenchen Tonglha Mountain is selected as the reference area, where the deformation can be neglected (see black square in Figure 4a).In addition, the areas with low coherence (<0.35) have been masked out to avoid unwrapping errors induced by large decorrelation noise.Figure 3 shows that most areas in the Yangbajing basin were rather stable, with LOS displacement rate below 0.5 cm/yr.However, the linear displacement velocity estimated by the data from Track 176 (Figure 3a) shows that the Yangbajing geothermal field (black box in Figure 1a), has experienced obvious local surface deformation.This ground subsidence is also detected by the data from Track 405, although only the right half of the geothermal field is covered (Figure 3b).As mentioned above, a linear ground deformation is expected to be caused by the withdrawal of subsurface fluid from the Yangbajing geothermal system.Figure 4a shows the superimposed displacement velocity fields estimated by both sets of SAR data, where the overlapping regions are difficult to distinguish due to their good resemblance.An enlarged view of the Yangbajing geothermal field is shown in Figure 5a.Most of the ground subsidence occurs in the southern part of the Yangbajing geothermal field, with a maximum LOS velocity of ~2 cm/yr.The northern part experiences gentle, smaller-scale ground subsidence, with LOS velocities generally smaller than 1 cm/yr.It is also observed that the ground subsidence increases from the margins to the center of the geothermal field.However, the displacements near the hydrothermal explode hollow and the Zangbo River are unavailable due to decorrelation of the InSAR observations.This is somewhat expected since the moist and changeable features of the ground surface at the discharge area of the geothermal field is quite unfavorable for InSAR observations.It is possible that even larger ground subsidence have occurred in that area.From the LOS surface displacement evolution estimated by the WLS InSAR algorithm for Tracks 176 and 405, we found that the whole Yangbajing basin experiences obvious seasonal oscillations, i.e., ground uplift in the cold season and subsidence in the warm season.This is to be expected, since this area is covered by the permafrost, deforming with the freezing-thawing cycle of the active layer [23].However, the surface displacements induced by geothermal production are somewhat obscured due to the seasonal oscillations.In order to present the details of the local deformation evolution, the deformations at two selected points marked A and B in Figure 5a, are shown in Figure 6a,b, respectively.It is observed that the displacement evolutions estimated from Track 176 and 405 agree quite well with each other.Point A shows seasonal variations in the margins of the geothermal field, where the cumulative displacement velocities are smaller than 2 mm/yr.Point B at the center of the geothermal field is characterized by both seasonal variation and cumulative displacement.The cumulative displacement velocity at point B, about ´1.5 cm/yr, reveals that the Yangbajing geothermal field experiences obvious ground subsidence.Profile CC' across the field, whose location is marked with a dashed line in Figure 5a, is also shown in Figure 6c, exhibiting the progress of the deformation along the profile.In order to better represent the cumulative deformation associated with the production of the Yangbajing geothermal system, four deformation evolutions acquired on nearly the same dates (but different years), i.The estimated LOS linear surface displacements at the Yangbajing geothermal field are further down-sampled by the quadtree algorithm to provide less-dense observations for the nonlinear inversion of the geometrical information about the source volume.In the quadtree down-sampling, the LOS displacement field is first divided by the blocks with a dimension size of 64 × 64.These initial blocks are then divided into four equal-sized square blocks if the difference between the maximum and the minimum displacements of the block elements is less than 2 mm, until the dimension size of block reaches 4 × 4. The result is shown in Figure 5b, achieving 844 data points, 2% of the original 43,684 displacement observations.Using the down-sampled InSAR displacement observations, the best fitted geometrical parameters are inverted by using a genetic algorithm with 51 iterations.These estimates are 189.1 and 81.0 m for the reservoir depth and thickness, respectively.Subsequently, the shallow reservoir is divided into a rectangular grid of size 100 × 128, i.e., two times coarser than that of the InSAR pixel grid.This means that the grid spacing of the shallow reservoir is ~90 m.The fractional volume change at each block is then determined by using a least-squares adjustment with a roughness penalty, as shown in Figure 7.The total value of the fractional volume change is estimated at −21.5, which is comparable to that of the uniform volume change (i.e., −24.3).Since no a priori information about the reservoir can be used to constrain the distribution of volume change, the roughness penalty term W of 10 is used in the inversion, on the basis of an examination of the data misfit.In other words, we pick W such that the difference between the InSAR-measured and the simulated displacement velocity fields is insignificant.In addition, we adopt an iterative strategy to ensure the robustness of the inversion.In the initial processing, the variety range of uniform volume change is relatively large to cover most situations.In the iterative processing, the variety range of uniform volume change is set to smaller by using the average value of the fractional volume change estimated from the last iteration as the reference.The iterative processing will be terminated when the standard deviation of the difference between the two consecutively estimated fractional volume changes is smaller than a pre-determined threshold (0.002 in this study).The estimated LOS linear surface displacements at the Yangbajing geothermal field are further down-sampled by the quadtree algorithm to provide less-dense observations for the nonlinear inversion of the geometrical information about the source volume.In the quadtree down-sampling, the LOS displacement field is first divided by the blocks with a dimension size of 64 ˆ64.These initial blocks are then divided into four equal-sized square blocks if the difference between the maximum and the minimum displacements of the block elements is less than 2 mm, until the dimension size of block reaches 4 ˆ4.The result is shown in Figure 5b, achieving 844 data points, 2% of the original 43,684 displacement observations.Using the down-sampled InSAR displacement observations, the best fitted geometrical parameters are inverted by using a genetic algorithm with 51 iterations.These estimates are 189.1 and 81.0 m for the reservoir depth and thickness, respectively.Subsequently, the shallow reservoir is divided into a rectangular grid of size 100 ˆ128, i.e., two times coarser than that of the InSAR pixel grid.This means that the grid spacing of the shallow reservoir is ~90 m.The fractional volume change at each block is then determined by using a least-squares adjustment with a roughness penalty, as shown in Figure 7.The total value of the fractional volume change is estimated at ´21.5, which is comparable to that of the uniform volume change (i.e., ´24.3).Since no a priori information about the reservoir can be used to constrain the distribution of volume change, the roughness penalty term W r of 10 is used in the inversion, on the basis of an examination of the data misfit.In other words, we pick W r such that the difference between the InSAR-measured and the simulated displacement velocity fields is insignificant.In addition, we adopt an iterative strategy to ensure the robustness of the inversion.In the initial processing, the variety range of uniform volume change is relatively large to cover most situations.In the iterative processing, the variety range of uniform volume change is set to smaller by using the average value of the fractional volume change estimated from the last iteration as the reference.The iterative processing will be terminated when the standard deviation of the difference between the two consecutively estimated fractional volume changes is smaller than a pre-determined threshold (0.002 in this study).The linear displacement velocity field simulated by the inverted source model is shown in Figure 5c and then compared with the InSAR-estimated one.The two displacement velocity fields are rather similar.Figure 5d shows the differences between the InSAR-estimated and the simulated displacement velocities.Most of the differences range from −1-1 mm, demonstrating that the inverted source model for the shallow reservoir can well explain the linear surface displacements measured in the Yangbajing geothermal field.Despite the overall agreement, some relatively large residuals The linear displacement velocity field simulated by the inverted source model is shown in Figure 5c and then compared with the InSAR-estimated one.The two displacement velocity fields are rather similar.Figure 5d shows the differences between the InSAR-estimated and the simulated displacement velocities.Most of the differences range from ´1 to 1 mm, de monstrating that the inverted source model for the shallow reservoir can well explain the linear surface displacements measured in the Yangbajing geothermal field.Despite the overall agreement, some relatively large residuals (up to 4 mm) can be found around the decorrelated areas.They are characterized by a grid-shape, mostly due to the InSAR observation noises and errors associated with the sparse grid used for the geothermal model inversion.
As expected, obvious volume changes are found under the ground where subsidence occurs (Figure 7).The total volume change of the Yangbajing geothermal system was estimated at ´1.41 ˆ10 7 m 3 /yr, indicating that a significant volume of geothermal fluids had been extracted from the shallow reservoir during the investigated period.Most of the source volume reduction occurred in the southern part of the field, between the Nepal Highway and the Zangbo River.

Discussion
Field surveys have been carried out in the Yangbajing field since 1983 to monitor the ground subsidence associated with the geothermal exploitation.The observed mean annual velocity reached ~3.5 cm/yr at the southern part of the field [2].Field work was suspended in 1994 due to shortage of financial support, so there are no in situ data to assess the accuracy of the InSAR-measured surface deformation.However, the decrease of the mean annual velocity during the investigated period is reasonable, since the main geothermal exploitation gradually transferred from the southern to the northern field [44].The observed deformation velocities are comparable to other measurements with the InSAR technique [19,20].Furthermore, we performed quantitative comparisons in the overlapping area between the displacement fields observed using SAR data from Tracks 176 and 405.As shown in Figure 4b, the differences basically range from ´2 to 2 mm/yr, with a standard deviation (STD) of 0.1 mm/yr and a mean of 0.4 mm/yr.This good agreement indicates that the WLS InSAR algorithm is suitable to estimate the linear surface displacements in the geothermal field.The residuals can be explained by the differences between the incidence angles of the two sets of SAR images.Compared to conventional survey methods, InSAR requires much less human and financial support, and therefore can be used as a routine tool in the investigation of the geothermal field.
Also, a 4-km portion of the China-Nepal Highway passes through the subsidence region.Although the highway and railway constructed in the permafrost area have been designed to prevent hazards related to the freezing-thawing cycle of the active layer, the cumulative inhomogeneous deformation induced by the geothermal production represents a more serious threat to road infrastructure.Therefore, the InSAR-derived displacements should be taken into account in any further exploitation plan for the Yangbajing geothermal field.
Although it has been reported that the Yangbajing geothermal system consists of shallow and deep reservoirs [25], in this study we only consider the shallow layer for the source model inversion.This simplification is realistic for three reasons.First, most thermal water used for electricity generation was extracted from the shallow reservoir during the study period [24].Second, the deep geothermal fluid is stored in the fracture zones and tectonic fissures at least 950 beneath the surface m, where the space and pressure are fairly steady.Therefore, it has quite limited effect on the ground surface deformation, which depends on the combination of pressure reduction and the effective elastic modulus of the reservoir.Thirdly, the deep reservoir is expected to be located in the northern part of the Yangbajing field and to generate the shallow reservoir in the Quaternary pores through a side flow recharge of hot fluids [2].Therefore, monitoring the dynamics of the shallow reservoir reflects in a certain degree the dynamics of the deep reservoir.
The production wells are distributed uniformly in the Yangbajing field, 24 in the northern part and 31 in the southern part (see dots in Figure 7).The localization of the volume reduction in the southern part indicates that the deep reservoir is located in the northern field and recharges the shallow reservoir there.This result supports the conceptual reservoir models for the Yangbajing geothermal system proposed in [1,2].Red dashed lines in Figure 7 denote the observed and inferred faults in the Yangbajing field.We found that rhombic block structures are crossly formed along the northeast and the northwest striking faults.The longest northeast striking fault is located near the Nepal Highway, serving as a boundary between the northern and southern geothermal systems.The mature faults in the northern part provide ideal conditions for the ascending flows from the deep reservoir.However, only one of the northwest striking faults extends to the southern field.This could explain why the fluid loss in the southern part of the shallow reservoir cannot be immediately replenished by the flow from the deep reservoir.

Conclusions
Monitoring the surface deformation induced by the geothermal production is essential to prevent possible geo-hazards, as well as to infer the dynamics of the subsurface flow fluids.In this study we employ an advanced InSAR algorithm, termed WLS InSAR, to investigate the Yangbajing geothermal field in the Qinghai-Tibet Plateau.We used a set of ENVISAT ASAR data acquired from two adjacent descending tracks from 2006 to 2010.By introducing a seasonal deformation model, the accumulated displacements caused by the geothermal production are separated from the seasonal oscillations associated with the freezing-thawing cycle of the permafrost.The elastic half-space theory is then applied to characterize the source model of the Yangbajing geothermal system from the WLS InSAR displacement measurements.Two main conclusions can be summarized from this study: (1) The WLS InSAR algorithm can be applied to detect and monitor ground deformation occurring in geothermal fields.Although in situ observations are unavailable during the study period, the pattern and amplitude of the InSAR-measured displacements are reasonable according to the field survey between 1983 and 1994.A cross-validation has been conducted between the two adjacent tracks' InSAR measurements, indicating an accuracy of millimeter or even sub-millimeter level.It is observed that the southern part of the Yangbajing field experiences considerable linear deformation along the LOS direction, of up to ~2 cm/yr, which behaves as ground subsidence and radial contraction.In addition, almost 8 cm peak-to-peak seasonal oscillations have been inferred from the deformation evolution results.
(2) InSAR-measured displacements are very useful for the inversion of the source model of the subsurface geothermal fluids.The geometrical parameters and the fractional volume changes of the Yanbajing geothermal system have been determined by a genetic algorithm and least squares adjustment, respectively.The results reveal that geothermal fluids of 1.41 ˆ10 7 m 3 had been extracted from the shallow reservoir every year in the period 2006-2010.The fluid loss in the southern part of the shallow reservoir dominates the source volume reduction.In the northern part, the discharged fluid has been refilled by ascending flow from the deep reservoir.The reason for these assumptions is that the deep reservoir is located in the northern field and has several developed faults serving as natural pipelines for the thermal fluids.
Since field surveying is always a time-and labor-consuming work, InSAR should be viewed as a routine tool to monitor and characterize the Yangbajing geothermal field in the future.In particular, the newer TerraSAR-X, COSMO-SkyMed and Sentinel-1 SAR data have finer resolution and shorter revisit period than the ENVISAT ASAR data, and the use of them will significantly enhance the spatial and temporal resolutions of the measured surface displacements.Although the volume change in the shallow reservoir is proven sufficient to fit the InSAR-measured displacements, the interactive behavior of the fluid flows between the shallow and deep reservoirs should be considered in order to investigate the interior dynamics of a geothermal system, especially when more surface displacement observations and subsurface geometrical information are provided.

Figure 1 .
Figure 1.(a) Optical image of the Yangbajing Basin.Dashed boxes outline the footprints of the ENVISAT ASAR descending images.The black rectangle indicates the study area.The inset provides the location of the study area in the Qinghai-Tibet Plateau; (b) Optical image of the study area, i.e., the Yangbajing geothermal field.

Figure 1 .
Figure 1.(a) Optical image of the Yangbajing Basin.Dashed boxes outline the footprints of the ENVISAT ASAR descending images.The black rectangle indicates the study area.The inset provides the location of the study area in the Qinghai-Tibet Plateau; (b) Optical image of the study area, i.e., the Yangbajing geothermal field.

Figure 2 .
Figure 2. Spatial-temporal baselines of the SAR data: (a) Track 176; and (b) Track 405.Circles and lines represent the SAR images and small-baseline interferograms.

Figure 2 .
Figure 2. Spatial-temporal baselines of the SAR data: (a) Track 176; and (b) Track 405.Circles and lines represent the SAR images and small-baseline interferograms.

Figure 3 .
Figure 3. Linear LOS displacement velocity maps for the Yangbajing basin, which are superimposed on the SRTM-derived shaded relief map.(a) Track 176; (b) Track 405.Black squares represent the reference area.

Figure 4 .
Figure 4. (a) Superimposed LOS linear displacement velocity map.The black rectangle marks the study area (Yangbajing geothermal field); (b) Histogram of the differences between the displacement results estimated from the SAR data of Tracks 176 and 405.

Figure 5 .
Figure 5. Maps of measured and simulated displacement velocities at the Yangbajing geothermal field, as superimposed on the optical map.(a) Linear displacement velocity map from InSAR; (b) Down-sampled linear displacement velocity; (c) Simulated linear displacement velocity; (d) Differences between the InSAR-estimated and the simulated linear displacement velocity.

Figure 3 . 14 Figure 3 .
Figure 3. Linear LOS displacement velocity maps for the Yangbajing basin, which are superimposed on the SRTM-derived shaded relief map.(a) Track 176; (b) Track 405.Black squares represent the reference area.

Figure 4 .
Figure 4. (a) Superimposed LOS linear displacement velocity map.The black rectangle marks the study area (Yangbajing geothermal field); (b) Histogram of the differences between the displacement results estimated from the SAR data of Tracks 176 and 405.

Figure 5 .
Figure 5. Maps of measured and simulated displacement velocities at the Yangbajing geothermal field, as superimposed on the optical map.(a) Linear displacement velocity map from InSAR; (b) Down-sampled linear displacement velocity; (c) Simulated linear displacement velocity; (d) Differences between the InSAR-estimated and the simulated linear displacement velocity.

Figure 4 .
Figure 4. (a) Superimposed LOS linear displacement velocity map.The black rectangle marks the study area (Yangbajing geothermal field); (b) Histogram of the differences between the displacement results estimated from the SAR data of Tracks 176 and 405.
e., 15 April 2007, 4 May 2008, 19 April 2009, and 9 May 2010, are selected to minimize the effects of the seasonal oscillations.The reference time of the four deformation evolutions is 30 April 2006.As displayed in Figure 6c a subsidence funnel was clearly formed in the Yangbajing field that kept growing with time.

Figure 4 .
Figure 4. (a) Superimposed LOS linear displacement velocity map.The black rectangle marks the study area (Yangbajing geothermal field); (b) Histogram of the differences between the displacement results estimated from the SAR data of Tracks 176 and 405.

Figure 5 .
Figure 5. Maps of measured and simulated displacement velocities at the Yangbajing geothermal field, as superimposed on the optical map.(a) Linear displacement velocity map from InSAR; (b) Down-sampled linear displacement velocity; (c) Simulated linear displacement velocity; (d) Differences between the InSAR-estimated and the simulated linear displacement velocity.

Figure 5 .
Figure 5. Maps of measured and simulated displacement velocities at the Yangbajing geothermal field, as superimposed on the optical map.(a) Linear displacement velocity map from InSAR; (b) Down-sampled linear displacement velocity; (c) Simulated linear displacement velocity; (d) Differences between the InSAR-estimated and the simulated linear displacement velocity.Triangles mark selected points A and B (see text).The dashed line shows the location of profile CC'.Black and white solid lines mark the China-Nepal Highway and Zangbo River, respectively.

Figure 6 .
Figure 6.Point time series of deformation and deformation profiles.(a,b) Evolutions of LOS surface deformation at points A and B marked in Figure 5a.Circles and triangles represent the results of

Figure 6 .
Figure 6.Point time series of deformation and deformation profiles.(a,b) Evolutions of LOS surface deformation at points A and B marked in Figure 5a.Circles and triangles represent the results of Track 176 and 405, respectively; (c) LOS surface deformation (with respect to 30 April 2006) along profile CC'.Deformations acquired on similar dates in four consecutive years are presented, in order to minimize the effects of seasonal variations.
Remote Sens. 2016, 8, 191 10 of 14 Track 176 and 405, respectively; (c) LOS surface deformation (with respect to 30 April 2006) along profile CC'.Deformations acquired on similar dates in four consecutive years are presented, in order to minimize the effects of seasonal variations.

Figure 7 .
Figure 7. Fractional volume changes of the Yangbajing geothermal system inverted from the InSAR displacement observations.Black and white solid lines represent the China-Nepal Highway and Zangbo River, respectively.Red dashed lines and black dots mark fault locations and the locations of production wells [2].

Figure 7 .
Figure 7. Fractional volume changes of the Yangbajing geothermal system inverted from the InSAR displacement observations.Black and white solid lines represent the China-Nepal Highway and Zangbo River, respectively.Red dashed lines and black dots mark fault locations and the locations of production wells [2].