Ground Deformation Revealed by Sentinel-1 MSBAS-InSAR Time-Series over Karamay Oilﬁeld, China

: Fluid extraction or injection into underground reservoirs may cause ground deformation, manifested as subsidence or uplift. Excessive deformation may threaten the infrastructure of the oilﬁeld and its surroundings and may even induce earthquakes. Therefore, the monitoring of surface deformation caused by oil production activities is important for the safe production of oilﬁelds and safety assessments of surrounding infrastructure. Karamay oilﬁeld is one of the major oil and gas ﬁelds in China. In this study, we take the Karamay oilﬁeld in Xinjiang as a case study to detect surface deformation caused by subsurface ﬂuid injection. Sentinel-1A images of 32 ascending (Path 114) and 34 descending (Path 165) tracks spanning March 2017 to August 2018, were used to derive vertical and horizontal deformation over Karamay oilﬁeld using the MSBAS-InSAR method. The observed two-dimensional deformation exhibited signiﬁcant vertical and east-west deformation in this region. The maximum uplift and horizontal velocity was approximately 160 mm / year and 65 mm / year, respectively. We also modeled one of the typical deformation zones using a dislocation model in a homogenous elastic half-space.


Introduction
Oilfield fluid injection is an important means of oil extraction, which mainly uses professional equipment to inject fluid into an oil layer to improve the speed of extraction and oil recovery efficiency [1,2]. However, fluid extraction or injection cause changes in reservoir pressure, which can cause deformation of the earth surface, manifested as subsidence or uplift [3]. The deformation may damage the environment and infrastructure of the oilfield, including oil and gas pipelines, roads, railways and buildings [4]. In addition, the injection of fluid may also induce earthquakes [5]. Therefore, monitoring of oilfield surface deformation is indispensable for understanding the underground oil extraction or fluid injection process [6]. Interferometric synthetic aperture radar (InSAR) is a geospatial observation technology with advantages of large spatial coverage and high spatial resolution, having centimeter level, or higher, monitoring accuracy [7]. Due to these advantages and characteristics, InSAR technology has been successfully applied to the monitoring of volcanic motion [8,9], seismic activity [10,11], landslides [12,13], glacial drift [14,15], regional ground sink [16,17], and anthropogenic activity like tunneling or infrastructure monitoring [18][19][20]. In addition, InSAR technology has also been successfully used to monitor ground deformation over oil and geothermal fields, such as the Helliishdi geothermal field [21], the West Texas oilfield [22], and the Cushing oilfield in East Central Oklahoma [23].

Stacking Method
A stacking method was proposed by Sandwell et al. in [35]. The mathematical model is shown in Equation (1). The assumption in this method is that the atmospheric signal is random in time, and ground deformation changes linearly. By averaging the unwrapped phase sets, the atmospheric delay phase is weakened and the linear deformation rate is extracted. According to the law of error propagation, when stacking N unwrapped interferograms, the linear deformation signal is increased N times, and the atmospheric phase error is only increased √N times. When the weighted average of the unwrapped phases of the N interferograms is calculated, the atmospheric influence weakens the original value of 1/√N, thereby increasing the signal-to-noise ratio [36,37]. Consequently, the random part of the atmospheric phase is greatly weakened [38]. As a result, the stacking method represents a better choice to obtain the average rate of the ground deformation. If the selected interferograms are connected end-to-end with the master-slave image, the phase noise of the slave image in the previous interferogram and the master image of the latter interferogram can be counteracted when using Equation (1) to stack. However, we must ensure the quality of the first and last interferograms.
Ph_rat t i * ph i / (1) where ℎ_ is the average deformation rate, is the number of interferograms involved in the stack processing, is time interval between master and slave pair, and ℎ is the unwrapped phase value of the interferogram.

Stacking Method
A stacking method was proposed by Sandwell et al. in [35]. The mathematical model is shown in Equation (1). The assumption in this method is that the atmospheric signal is random in time, and ground deformation changes linearly. By averaging the unwrapped phase sets, the atmospheric delay phase is weakened and the linear deformation rate is extracted. According to the law of error propagation, when stacking N unwrapped interferograms, the linear deformation signal is increased N times, and the atmospheric phase error is only increased √ N times. When the weighted average of the unwrapped phases of the N interferograms is calculated, the atmospheric influence weakens the original value of 1/ √ N, thereby increasing the signal-to-noise ratio [36,37]. Consequently, the random part of the atmospheric phase is greatly weakened [38]. As a result, the stacking method represents a better choice to obtain the average rate of the ground deformation. If the selected interferograms are connected end-to-end with the master-slave image, the phase noise of the slave image in the previous interferogram and the master image of the latter interferogram can be counteracted when using Equation (1) to stack. However, we must ensure the quality of the first and last interferograms.
where Ph_rat is the average deformation rate, N is the number of interferograms involved in the stack processing, t i is time interval between master and slave pair, and ph i is the unwrapped phase value of the interferogram.

Multidimensional Small Baseline Subset Method
The combination of ascending and descending track images not only can improve the time density of surface observations, but can also generate a two-dimensional deformation field of the target points [39,40]. In this paper, the multidimensional small baseline subset (MSBAS) method was used to process the interferometric displacement results of the ascending and descending images [28,30]. The precondition for solving two-dimensional time series deformation is that the SAR images from different orbits have the same imaging time and time span; these conditions are not met for SAR images from different satellites or orbits. However, studies from Samsonov and co-workers have shown that SAR satellites with different orbits can solve the east-west and vertical time series well, if the time span has good consistency [28,30]. They have also achieved good results in monitoring the surface deformation possibly caused by CO2 storage with MSBAS [41]. Compared with the one-dimensional deformation model, the unknown number of the solution is multiplied in the two-dimensional deformation model, which leads to the need for more observations to solve the rank defect of the equation. Here we use the regularization procedure to add the identity matrix to alleviate the ill-posed problem of the equation and to bring the equation to the full rank condition. The two-dimensional mathematical model was calculated with the MSBAS method: − cos α sin θ∆t cos θ∆t λI where α is the radar azimuth angle, −14.2 • (ascending) and −165.2 • (descending) in our case, and θ is the average incidence angle, 33.5 • (ascending) and 33.7 • (descending) in our case. ∆t is the time interval between the two images of the interferogram and λ is a regularization parameter, I is an identity matrix. The regularization parameter λ is used to regularize the solution. The larger the λ value, the smoother the solution. The L-curve method can be used to select the optimal value [42]. V E and V U are the average deformation rates in the vertical and east-west direction, D is the observed (unwrapped) interferometric displacement. As shown in Figure 2, the first image from the descending path is acquired earlier than that of the ascending image, and we compensate the first descending interferogram by a factor of (t 1 − t 0 )/(t 1 − t −1 ); similarly the last ascending interferogram is compensated by a factor of (t 5 − t 4 )/(t 6 − t 4 ). After applying boundary correction, the first and last image acquisition dates are assumed to be at t 0 and t 5 for the ascending and descending images [29]. Then, we can obtain more detailed ground surface deformation caused by fluid extraction or injection by MSBAS method.

Multidimensional Small Baseline Subset Method
The combination of ascending and descending track images not only can improve the time density of surface observations, but can also generate a two-dimensional deformation field of the target points [39,40]. In this paper, the multidimensional small baseline subset (MSBAS) method was used to process the interferometric displacement results of the ascending and descending images [28,30]. The precondition for solving two-dimensional time series deformation is that the SAR images from different orbits have the same imaging time and time span; these conditions are not met for SAR images from different satellites or orbits. However, studies from Samsonov and co-workers have shown that SAR satellites with different orbits can solve the east-west and vertical time series well, if the time span has good consistency [28,30]. They have also achieved good results in monitoring the surface deformation possibly caused by CO2 storage with MSBAS [41]. Compared with the onedimensional deformation model, the unknown number of the solution is multiplied in the twodimensional deformation model, which leads to the need for more observations to solve the rank defect of the equation. Here we use the regularization procedure to add the identity matrix to alleviate the ill-posed problem of the equation and to bring the equation to the full rank condition. The two-dimensional mathematical model was calculated with the MSBAS method: where α is the radar azimuth angle, 14.2° (ascending) and 165.2° (descending) in our case, and θ is the average incidence angle, 33.5° (ascending) and 33.7° (descending) in our case. ∆t is the time interval between the two images of the interferogram and λ is a regularization parameter, I is an identity matrix. The regularization parameter λ is used to regularize the solution. The larger the λ value, the smoother the solution. The L-curve method can be used to select the optimal value [42]. V and V are the average deformation rates in the vertical and east-west direction, D is the observed (unwrapped) interferometric displacement. As shown in Figure 2, the first image from the descending path is acquired earlier than that of the ascending image, and we compensate the first descending interferogram by a factor of t t / t t ; similarly the last ascending interferogram is compensated by a factor of t t / t t . After applying boundary correction, the first and last image acquisition dates are assumed to be at t and t for the ascending and descending images [29]. Then, we can obtain more detailed ground surface deformation caused by fluid extraction or injection by MSBAS method. acquisitions at time t-1 to t6 are marked with small black circles and blue circles, respectively. The horizontal solid lines between two dots represent interferograms, which are marked with I1 to I6. ∆t1 to ∆t5 represent the time intervals between acquired adjacent images (modified from Samsonov [29]).

Data and Processing
SAR images covering the study area were acquired by C-band (5.6 cm wavelength and 5.4 GHz) radar on board the European Space Agency (ESA) Sentinel-1 A/B satellites. Sentinel A and Sentinel B are the latest generation of C-band satellites launched by ESA in 2014 and 2016, respectively; the coherence between images is increased compared with previous ESA satellites, because of a reduced revisit time (12 days for Sentinel A, and 6 days for the combination of Sentinel A and B) [43]. The image in wide swath mode provides a wide coverage of about 250 km with a slant range resolution Simplified schematics of two-dimensional time series. Ascending and descending SAR acquisitions at time t −1 to t 6 are marked with small black circles and blue circles, respectively. The horizontal solid lines between two dots represent interferograms, which are marked with I 1 to I 6 . ∆t 1 to ∆t 5 represent the time intervals between acquired adjacent images (modified from Samsonov [29]).

Data and Processing
SAR images covering the study area were acquired by C-band (5.6 cm wavelength and 5.4 GHz) radar on board the European Space Agency (ESA) Sentinel-1 A/B satellites. Sentinel A and Sentinel B are the latest generation of C-band satellites launched by ESA in 2014 and 2016, respectively; the coherence between images is increased compared with previous ESA satellites, because of a reduced revisit time (12 days for Sentinel A, and 6 days for the combination of Sentinel A and B) [43]. The image in wide swath mode provides a wide coverage of about 250 km with a slant range resolution of 5 m and an azimuth resolution of 20 m. In this study, we used GAMMA software to process 32 ascending ( Figure 3a) and 34 descending (Figure 3b) datasets, gathered between March 2017 and August 2018 [44]. There are two kinds of orbital data files provided by ESA, including correctional orbital parameters with accuracy within 10 cm and precise orbital ephemeris parameters with accuracy within 5 cm [45]. The latter were used to improve the accuracy of the satellite position and decrease the perpendicular baseline error. We chose spatial perpendicular baselines smaller than 100 m and temporal baselines smaller than 100 days as the thresholds of the interferograms. Interferograms with good coherence and small atmospheric errors were selected for further processing. As a result, we obtained 82 ascending ( Figure 3a) and 131 descending (Figure 3b) interferograms. A digital elevation model (DEM) with 1 arc-second (about 30 m resolution) from the Shuttle Radar Topography Mission (SRTM) (http://srtm.csi.cgiar.org) was used to remove topographic effects from the interferograms. Adaptive filtering was used to filter the interferograms [46], and the minimum cost flow (MCF) algorithm was used to unwrap them [47].
Remote Sens. 2019, 11, x FOR PEER REVIEW 5 of 18 of 5 m and an azimuth resolution of 20 m. In this study, we used GAMMA software to process 32 ascending ( Figure 3a) and 34 descending (Figure 3b) datasets, gathered between March 2017 and August 2018 [44]. There are two kinds of orbital data files provided by ESA, including correctional orbital parameters with accuracy within 10 cm and precise orbital ephemeris parameters with accuracy within 5 cm [45]. The latter were used to improve the accuracy of the satellite position and decrease the perpendicular baseline error. We chose spatial perpendicular baselines smaller than 100 m and temporal baselines smaller than 100 days as the thresholds of the interferograms. Interferograms with good coherence and small atmospheric errors were selected for further processing. As a result, we obtained 82 ascending ( Figure 3a) and 131 descending (Figure 3b) interferograms. A digital elevation model (DEM) with 1 arc-second (about 30 m resolution) from the Shuttle Radar Topography Mission (SRTM) (http://srtm.csi.cgiar.org) was used to remove topographic effects from the interferograms. Adaptive filtering was used to filter the interferograms [46], and the minimum cost flow (MCF) algorithm was used to unwrap them [47].

Deformation Velocity
The annual average deformation rate ( Figure 4) in the line of sight (LOS) over the Karamay oilfield area was obtained using the stacking method (Equation (1)). The deformation rates from ascending and descending orbits exhibited a good consistency in magnitude. There were four obvious uplift areas (marked as A, B, C and D in Figure 4) and one subsidence center. The maximum rate of subsidence was about 32 mm/year. The maximum rate of uplift, exceeding 100 mm/year, appeared in area A.

Deformation Velocity
The annual average deformation rate (Figure 4) in the line of sight (LOS) over the Karamay oilfield area was obtained using the stacking method (Equation (1)). The deformation rates from ascending and descending orbits exhibited a good consistency in magnitude. There were four obvious uplift areas (marked as A, B, C and D in Figure 4) and one subsidence center. The maximum rate of subsidence was about 32 mm/year. The maximum rate of uplift, exceeding 100 mm/year, appeared in area A. Compared with previous research results [6,25], the uplift deformation rate changed from 30 mm/year between 2003 and 2010 to 150 mm/year from 2017 to 2018. The region with the largest deformation rate occurred in Zone B from 2003 to 2010, while the region with the largest deformation rate occurred in Zone A from 2017 to 2018. We suspect that these phenomenon may be caused by the transfer of main oil exploration area from Zone B to Zone A, which needs to be confirmed by the oilfield mining company. Compared with previous research results [6,25], the uplift deformation rate changed from 30 mm/year between 2003 and 2010 to 150 mm/year from 2017 to 2018. The region with the largest deformation rate occurred in Zone B from 2003 to 2010, while the region with the largest deformation rate occurred in Zone A from 2017 to 2018. We suspect that these phenomenon may be caused by the transfer of main oil exploration area from Zone B to Zone A, which needs to be confirmed by the oilfield mining company.

Two-Dimensional Time Series Deformation
In order to derive the vertical and horizontal components of deformation of the ground surface, we applied the MSBAS method (Equation (2)) to obtain the two-dimensional time series deformation field (east-west and vertical) of the study area [40,48]. The vertical and the east-west deformation rates are presented in Figure 5a,b, respectively. Both vertical and east-west displacements occurred in the four deformation zones and are marked in Figure 5  In order to derive the vertical and horizontal components of deformation of the ground surface, we applied the MSBAS method (Equation (2)) to obtain the two-dimensional time series deformation field (east-west and vertical) of the study area [40,48]. The vertical and the east-west deformation rates are presented in Figures 5a and 5b, respectively. Both vertical and east-west displacements occurred in the four deformation zones and are marked in Figure 5 (A, B, C and D).  The maximum uplift rate was up to 160 mm/year, and the maximum eastward movement rate was 65 mm/year in area A ( Figure 5). Area B was also a typical uplift area with horizontal motion. The magnitude of the uplift was about 110 mm/year, and the east flank of the uplifted center moved eastward at a rate of 30 mm/year, and the west side of the uplifted center moved westward at a rate of 15 mm/year. The uplift rate at area C was approximately 80 mm/year and there was also a small area with subsidence rate of 16 mm/year. Area D exhibited both uplift and obvious subsidence; uplift was at a rate of approximately 110 mm/year and the subsidence rate was approximately 25 mm/year. Although there were varying degrees of east-west and vertical deformation in the four study areas, the vertical displacement was dominant.
In order to study the deformation process of the uplift zone, area B1 ( Figure 5) was selected for examination of its two-dimensional time series deformation in the vertical ( Figure 6) and east-west directions (Figure 7). The maximum uplift rate was up to 160 mm/year, and the maximum eastward movement rate was 65 mm/year in area A ( Figure 5). Area B was also a typical uplift area with horizontal motion. The magnitude of the uplift was about 110 mm/year, and the east flank of the uplifted center moved eastward at a rate of 30 mm/year, and the west side of the uplifted center moved westward at a rate of 15 mm/year. The uplift rate at area C was approximately 80 mm/year and there was also a small area with subsidence rate of 16 mm/year. Area D exhibited both uplift and obvious subsidence; uplift was at a rate of approximately 110 mm/year and the subsidence rate was approximately 25 mm/year. Although there were varying degrees of east-west and vertical deformation in the four study areas, the vertical displacement was dominant.
In order to study the deformation process of the uplift zone, area B1 ( Figure 5) was selected for examination of its two-dimensional time series deformation in the vertical ( Figure 6) and east-west directions (Figure 7).  The vertical deformation rate of the area increased slowly with time before March 2018, but significantly accelerated thereafter ( Figure 6). The magnitude of the uplift reached approximately 160 mm by 15 August 2018. Similarly, the east-west deformation rate was also slow before March 2018 and the horizontal motion intensified after March 2018 (Figure 7). The amounts of westward displacement on the western flank of the uplift reached 30 mm, and eastward movement on the other side was closed to 60 mm. Several previous studies have shown that fluid injection, a commonly used method in oil and gas production, can lead to ground uplift [6,26]. This may be related to enhanced water injection and production activities. The increase of pore pressure with fluid injection leads to a decrease of the effective stress of the reservoir, which leads to uplift of the ground surface [27]. Therefore, we speculated that the deformation in this area was likely to be related to liquid injection activities during oil and gas production. Rapid ground deformation may pose a threat to underground pipeline facilities and surface constructions [49]. Therefore, oil and gas production managers should pay attention to it. The vertical deformation rate of the area increased slowly with time before March 2018, but significantly accelerated thereafter ( Figure 6). The magnitude of the uplift reached approximately 160 mm by 15 August 2018. Similarly, the east-west deformation rate was also slow before March 2018 and the horizontal motion intensified after March 2018 (Figure 7). The amounts of westward displacement on the western flank of the uplift reached 30 mm, and eastward movement on the other side was closed to 60 mm. Several previous studies have shown that fluid injection, a commonly used method in oil and gas production, can lead to ground uplift [6,26]. This may be related to enhanced water injection and production activities. The increase of pore pressure with fluid injection leads to a decrease of the effective stress of the reservoir, which leads to uplift of the ground surface [27]. Therefore, we speculated that the deformation in this area was likely to be related to liquid injection activities during oil and gas production. Rapid ground deformation may pose a threat to underground pipeline facilities and surface constructions [49]. Therefore, oil and gas production managers should pay attention to it.

Deformation Analysis
GPS measurements at station XJKL (triangle in Figure 5) were obtained from (ftp://ftp.cgps.ac.cn/products/position/gamit/XJKL.ios_gamit_detrend.neu) (Global Navigation Satellite System Data Products of China Earthquake Administration). They were used to evaluate the reliability of the InSAR results (Figure 8). This comparison shows that the InSAR and GPS observations are in good agreement. Keeping in mind that the GPS station XJKL is a station from China Mainland Tectonic Environment Monitoring Network Project. The horizontal displacement of XJKL refers to the International GNSS Service (IGS) station information of neighboring countries or regions [50]. Therefore, the horizontal displacement on the XJKL station represents the tectonic  Figure 5) were obtained from (ftp://ftp.cgps.ac. cn/products/position/gamit/XJKL.ios_gamit_detrend.neu) (Global Navigation Satellite System Data Products of China Earthquake Administration). They were used to evaluate the reliability of the InSAR results ( Figure 8). This comparison shows that the InSAR and GPS observations are in good agreement. Keeping in mind that the GPS station XJKL is a station from China Mainland Tectonic Environment Monitoring Network Project. The horizontal displacement of XJKL refers to the International GNSS Service (IGS) station information of neighboring countries or regions [50]. Therefore, the horizontal displacement on the XJKL station represents the tectonic movement over a large area. However, the horizontal displacement from InSAR is very local. This may result in the horizontal displacements from InSAR being less than that in XJKL. movement over a large area. However, the horizontal displacement from InSAR is very local. This may result in the horizontal displacements from InSAR being less than that in XJKL. In order to study the time-varying characteristics of surface deformation, time series InSAR results at six points were extracted ( Figure 5). In order to avoid the influence of inconsistent precision at individual pixels, the average value of 5 × 5 pixels, i.e., 150 m × 150 m around each point was used to plot the time series results (Figure 9).  In order to study the time-varying characteristics of surface deformation, time series InSAR results at six points were extracted ( Figure 5). In order to avoid the influence of inconsistent precision at individual pixels, the average value of 5 × 5 pixels, i.e., 150 m × 150 m around each point was used to plot the time series results (Figure 9). movement over a large area. However, the horizontal displacement from InSAR is very local. This may result in the horizontal displacements from InSAR being less than that in XJKL. In order to study the time-varying characteristics of surface deformation, time series InSAR results at six points were extracted ( Figure 5). In order to avoid the influence of inconsistent precision at individual pixels, the average value of 5 × 5 pixels, i.e., 150 m × 150 m around each point was used to plot the time series results (Figure 9).  Point P3 (in Figure 5) located at Jinlong Town is away from the obvious deformation zone, and the other selected points are near the center of uplifting of each deformation zone. The horizontal and vertical displacements at point P3 (Figure 9c) oscillated slightly around 0 with time, indicating that the results we obtained were reliable. The uplift center P1 (Figure 9a), located 6 km southwest of the railway station, showed close to a linear rate of uplift with time. By August 2018, the relative uplift was about 240 mm. The maximum deformation rate reached 160 mm/year from 2017 to 2018. This result is different from a previously reported maximum deformation rate of 50 mm/year from 2006 to 2010 [6,26]. According to the report, the oil reservoirs in the Karamay oilfield area generally entered a high moisture cut period after 2015. In order to ensure oil production, increasing water injection was taken after 2015 [51]. Therefore, continuous fluid injection activity may be a factor that caused continued uplift. Unlike P1, P2 in area B (Figure 5), located to the west of Jinlong Town, exhibited a slow uplift before March 2018, reaching approximately 30 mm (Figure 9b). However, the magnitude of the uplift increased to 130 mm from March 2018 to August 2018, which may have been related to enhanced injection and production activities in this area.

GPS measurements at station XJKL (triangle in
We selected images from Figure 5a,b at regular time intervals to make time series profiles of uplift deformation and east-west displacement. Deformation time series along profile i-j in Figure 5 highlight the spatiotemporal progress of uplift and east-west movement ( Figure 10) in the vertical (a) and east-west (b) directions. Similar deformation modes were apparent points P4 and P5 on the west side of the Baijiantan (Figure 9d,e). The displacement of P4 and P5 increased with time, and the deformation was relatively small after March 2018. Uplift in the early stage and a tendency to subsidence in the later stage was apparent at point P6, located to the east of Baijiantan (Figure 9f). This phenomenon may be due to the injection of liquid in the early stage, leading to the uplift. The oil exploitation caused the later subsidence. Deformation in area D occurred near the national highway G217 ( Figure 5) and has potential to threaten road stability and the nearby oil/gas pipelines. Point P3 (in Figure 5) located at Jinlong Town is away from the obvious deformation zone, and the other selected points are near the center of uplifting of each deformation zone. The horizontal and vertical displacements at point P3 (Figure 9c) oscillated slightly around 0 with time, indicating that the results we obtained were reliable. The uplift center P1 (Figure 9a), located 6 km southwest of the railway station, showed close to a linear rate of uplift with time. By August 2018, the relative uplift was about 240 mm. The maximum deformation rate reached 160 mm/year from 2017 to 2018. This result is different from a previously reported maximum deformation rate of 50 mm/year from 2006 to 2010 [6,26]. According to the report, the oil reservoirs in the Karamay oilfield area generally entered a high moisture cut period after 2015. In order to ensure oil production, increasing water injection was taken after 2015 [51]. Therefore, continuous fluid injection activity may be a factor that caused continued uplift. Unlike P1, P2 in area B (Figure 5), located to the west of Jinlong Town, exhibited a slow uplift before March 2018, reaching approximately 30 mm (Figure 9b). However, the magnitude of the uplift increased to 130 mm from March 2018 to August 2018, which may have been related to enhanced injection and production activities in this area.
We selected images from Figures 5a and 5b at regular time intervals to make time series profiles of uplift deformation and east-west displacement. Deformation time series along profile i-j in Figure  5 highlight the spatiotemporal progress of uplift and east-west movement ( Figure 10) in the vertical (a) and east-west (b) directions. Similar deformation modes were apparent points P4 and P5 on the west side of the Baijiantan (Figure 9d, 9e). The displacement of P4 and P5 increased with time, and the deformation was relatively small after March 2018. Uplift in the early stage and a tendency to subsidence in the later stage was apparent at point P6, located to the east of Baijiantan (Figure 9f). This phenomenon may be due to the injection of liquid in the early stage, leading to the uplift. The oil exploitation caused the later subsidence. Deformation in area D occurred near the national highway G217 ( Figure 5) and has potential to threaten road stability and the nearby oil/gas pipelines.

Inversion Modeling
Geophysical inversion in our study is based on monitoring data and applying a corresponding physical model to explore model parameters and reservoir characteristics. Because the ground uplift was induced by reservoir expansion, it provides valuable information allowing us to evaluate the dynamic changes in reservoir properties [52]. Using the deformation information acquired by InSAR and the corresponding geophysical model, reservoir parameter values can be inverted to better understand the underground fluid injection and surface dynamics [3]. Compared to other existing models, only physical models from Mogi [53] and Sill [54] are commonly used to explain ground surface deformation caused by changes in underground filling. Therefore, both models were selected to perform parameter inversion for uplift areas B1 and C1 (see Figure 5), in order to explain surface uplift behavior in this study.

Mogi Modeling
The Mogi model is a simple and effective inversion model of surface deformation arising from pressure changes of underground reservoirs [53]. The simulated parameters include point source center position (X, Y), reservoir depth and volume change. The vertical cumulative deformation maps ( Figure 7) were used to invert the reservoir parameters. The root mean square error (RMSE) between the observed and modeled values was used as inversion criterion. Table 1 shows the best fitting parameters obtained by the Mogi model. Volume changes in areas B1 and C1 are 1.97 × 10 5 m 3 and 3.87 × 10 4 m 3 , respectively. Most of the residuals were smaller than 10 mm (Figure 11c,f). However, a relatively large residual occurred on the long axis of the elliptical deformation for B1. This observation indicates that the point source Mogi model is suitable for circular surface deformation, but when the deformation of the surface approximates an ellipse, it does not perform well anymore. For area C1, the deformation characteristics are close to circular, so the Mogi model performs adequately. Geophysical inversion in our study is based on monitoring data and applying a corresponding physical model to explore model parameters and reservoir characteristics. Because the ground uplift was induced by reservoir expansion, it provides valuable information allowing us to evaluate the dynamic changes in reservoir properties [52]. Using the deformation information acquired by InSAR and the corresponding geophysical model, reservoir parameter values can be inverted to better understand the underground fluid injection and surface dynamics [3]. Compared to other existing models, only physical models from Mogi [53] and Sill [54] are commonly used to explain ground surface deformation caused by changes in underground filling. Therefore, both models were selected to perform parameter inversion for uplift areas B1 and C1 (see Figure 5), in order to explain surface uplift behavior in this study.

Mogi Modeling
The Mogi model is a simple and effective inversion model of surface deformation arising from pressure changes of underground reservoirs [53]. The simulated parameters include point source center position (X, Y), reservoir depth and volume change. The vertical cumulative deformation maps ( Figure 7) were used to invert the reservoir parameters. The root mean square error (RMSE) between the observed and modeled values was used as inversion criterion. Table 1 shows the best fitting parameters obtained by the Mogi model. Volume changes in areas B1 and C1 are 1.97 10 m and 3.87 10 m , respectively. Most of the residuals were smaller than 10 mm (Figures 11c,f). However, a relatively large residual occurred on the long axis of the elliptical deformation for B1. This observation indicates that the point source Mogi model is suitable for circular surface deformation, but when the deformation of the surface approximates an ellipse, it does not perform well anymore. For area C1, the deformation characteristics are close to circular, so the Mogi model performs adequately.

Sill Modeling
We also modeled oilfield deformation to explain the subsurface processes and explain the integral uplift pattern of the ground with the Sill model [54]. The Sill model, embedded in an elastic half-space, is defined by seven parameters: Location (X, Y), depth, strike, length, width, and opening. The RMSE between the observed and modeled values was used as inversion criterion. The optimal parameters of the Sill simulation are shown in Table 2. The calculated reservoir depth in areas B1 and C1 were 606 m and 422 m, respectively. One source with a rectangular sill of size 996 × 230 m, below the deformation area B1, was directed toward 270 • from the north (anticlockwise). Another source with a rectangular sill 289 × 256 m, below area C1, was directed 189 • from the north (anticlockwise). The residuals (Figure 12c,f) were mostly concentrated between −15 mm and 15 mm. Although we used a uniform elastic half-space model [54], the geology of the study area has some heterogeneities, which may affect the inversion results to some extent. We also modeled oilfield deformation to explain the subsurface processes and explain the integral uplift pattern of the ground with the Sill model [54]. The Sill model, embedded in an elastic half-space, is defined by seven parameters: Location (X, Y), depth, strike, length, width, and opening. The RMSE between the observed and modeled values was used as inversion criterion. The optimal parameters of the Sill simulation are shown in Table 2. The calculated reservoir depth in areas B1 and C1 were 606 m and 422 m, respectively. One source with a rectangular sill of size 996 × 230 m, below the deformation area B1, was directed toward 270° from the north (anticlockwise). Another source with a rectangular sill 289 × 256 m, below area C1, was directed 189° from the north (anticlockwise). The residuals (Figures 12c,f) were mostly concentrated between 15 mm and 15 mm. Although we used a uniform elastic half-space model [54], the geology of the study area has some heterogeneities, which may affect the inversion results to some extent.  Comparing the results from the two models (Tables 1 and 2), the depths were close and were consistent with the depth of the oil well in this area [55]. In area B1, the maximum volume changes from the Mogi model and the minimum volume change from the Sill model (2.08×10 5 and 2.15×10 5 m 3 , respectively) were very similar to each other. Similar to the C1 area, the volume changes (3.2 ×10 4 to 4.4×10 4 m 3 ) from Mogi model were also close to the volume changes (2.2×10 4 to 3.5×10 4 m 3 ) from Sill model. Unfortunately, we cannot compare our modeled volumes with the actual volume mined, because it is not available in the published data.

Deformation and Seismicity
Prior to liquid injection, the stress state is determined by local tectonic stress, which may cause new geomechanical processes after liquid injection [27]. Increased pore pressure due to fluid injection reduces the effective stress, which causes the reservoir to expand and act on the surface to cause ground uplift. This mechanism is easier to understand with the Mohr-Coulomb circles (Figure 13), Comparing the results from the two models (Tables 1 and 2), the depths were close and were consistent with the depth of the oil well in this area [55]. In area B1, the maximum volume changes from the Mogi model and the minimum volume change from the Sill model (2.08 × 10 5 and 2.15 × 10 5 m 3 , respectively) were very similar to each other. Similar to the C1 area, the volume changes (3.2 × 10 4 to 4.4 × 10 4 m 3 ) from Mogi model were also close to the volume changes (2.2 × 10 4 to 3.5 × 10 4 m 3 ) from Sill model. Unfortunately, we cannot compare our modeled volumes with the actual volume mined, because it is not available in the published data.

Deformation and Seismicity
Prior to liquid injection, the stress state is determined by local tectonic stress, which may cause new geomechanical processes after liquid injection [27]. Increased pore pressure due to fluid injection reduces the effective stress, which causes the reservoir to expand and act on the surface to cause ground uplift. This mechanism is easier to understand with the Mohr-Coulomb circles (Figure 13), which is modified from Teatini [26]. According to Terzaghi's principle, when a fluid is injected into the tiny pores of the rock, the pore pressure will increase, resulting in a decrease in the effective stress of the rock. The Mohr circle shifts to the left as effective stress decreases. There are maybe two outcomes when the Mohr circle continues to move left: (1) A shear failure may happen if the circle is tangent to the envelope; and (2) a tensile damage may occur when the circle intersects the envelope [27]. The new stress is closer to the failure envelope than before the injection because of the injection of fluid, which may induce seismicity [56].
We used seismic data of magnitude 1.0 or higher within 70 km of the oilfield from March 2017 to August 2018 from the China Earthquake Data Center (http://data.earthquake.cn). These data are plotted in Figures 1 and 14. The locations of seismicity are mainly distributed near faults, especially near the Karamay fault ( Figure 1). The earthquake epicenters mainly occur near the Baibai fault, West Mahu fault, and the Baiquankou fault. In Figure 14, the depth of the inversion reservoir is about 600 m (red line), and as mentioned earlier, the depth of the reservoir is from 300 (black line) to 3000 meters (blue line). However, the depth of the earthquakes are greater than 3000 m. So based on the location and depth of earthquakes in study area, we suggest that the occurrence of the earthquakes in this region is not related to the activities associated with oil production. which is modified from Teatini [26]. According to Terzaghi's principle, when a fluid is injected into the tiny pores of the rock, the pore pressure will increase, resulting in a decrease in the effective stress of the rock. The Mohr circle shifts to the left as effective stress decreases. There are maybe two outcomes when the Mohr circle continues to move left: (1) A shear failure may happen if the circle is tangent to the envelope; and (2) a tensile damage may occur when the circle intersects the envelope [27]. The new stress is closer to the failure envelope than before the injection because of the injection of fluid, which may induce seismicity [56]. We used seismic data of magnitude 1.0 or higher within 70 km of the oilfield from March 2017 to August 2018 from the China Earthquake Data Center (http://data.earthquake.cn). These data are plotted in Figures 1 and 14. The locations of seismicity are mainly distributed near faults, especially near the Karamay fault ( Figure 1). The earthquake epicenters mainly occur near the Baibai fault, West Mahu fault, and the Baiquankou fault. In Figure 14, the depth of the inversion reservoir is about 600 m (red line), and as mentioned earlier, the depth of the reservoir is from 300 (black line) to 3000 meters (blue line). However, the depth of the earthquakes are greater than 3000 m. So based on the location and depth of earthquakes in study area, we suggest that the occurrence of the earthquakes in this region is not related to the activities associated with oil production. Figure 13. Mohr-Coulomb circles. τ1 an τ2 are the current and maximum allowable shear stress, respectively. σ and σ are the maximum and minimum principal stress, respectively. The Mohr-Coulomb circles move to the left when pore pressure (p0) increases (modified from Teatini [27]). Figure 13. Mohr-Coulomb circles. τ 1 an τ 2 are the current and maximum allowable shear stress, respectively. σ 1 and σ 2 are the maximum and minimum principal stress, respectively. The Mohr-Coulomb circles move to the left when pore pressure (p 0 ) increases (modified from Teatini [27]).