Ground Deformation Revealed by Sentinel-1 MSBAS-InSAR Time-Series over Karamay Oilfield, China

Abstract

Fluid extraction or injection into underground reservoirs may cause ground deformation, manifested as subsidence or uplift. Excessive deformation may threaten the infrastructure of the oilfield 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 oilfields and safety assessments of surrounding infrastructure. Karamay oilfield is one of the major oil and gas fields in China. In this study, we take the Karamay oilfield in Xinjiang as a case study to detect surface deformation caused by subsurface fluid 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 oilfield using the MSBAS-InSAR method. The observed two-dimensional deformation exhibited significant 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.

1. 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].

Karamay oilfield, as one of the major oil and gas fields in China, started using underground water injection to improve oil recovery rate in 1985 [1], and horizontal well drilling technology is used to further improve single well production [24]. The Karamay conglomerate reservoir was discovered in 1955. In 1960, oil production was 1.64 million tons, followed by steady-state growth. In 2002, oil production exceeded 10 million tons [25]. Several studies have examined ground deformation characteristics caused by oil and gas exploitation before 2010 in this area, by using Envisat ASAR and ALOS PALSAR images [6,26]. These studies revealed that there was significant deformation in this area, with maximum accumulated uplift of 200 mm between 2006 and 2011. However, there have been few studies on surface deformation characteristics of this area in recent years. Therefore, monitoring recent surface deformations over the Karamay oilfield area is of great significance for guiding the continued exploitation for oil and gas. Liquid injection may cause local stress changes, which can lead to earthquakes [21,23,27]. Are the low-magnitude earthquakes within and around the oilfield related to oil and gas exploration in this area? In this paper, surface deformation characteristics of the Karamay oilfield during the period March 2017 to August 2018 were monitored with 32 ascending and 34 descending images from Sentinel-1 A/B. The two-dimensional deformation field of the Karamay region was obtained following application of the multi-dimensional small baseline subset technique (MSBAS) [28,29,30]. Detailed surface deformation processes were analyzed and discussed. This paper is organized as follows. First, background information on the oilfield is introduced in Section 2. Then, the stacking and MSBAS methods and SAR image processing are presented in Section 3. The InSAR results are presented in Section 4. Analysis and discussion of the results are implemented in Section 5. Finally, some conclusions are drawn in Section 6.

2. Study Area

Karamay region in which the oilfield is located covers an area of 9500 ${\mathrm{km}}^2$ and consists of four administrative districts, including Karamay, Baijiantan, Urho and Dushanzi. The altitude is high to the northwest and low to the southeast. The central and eastern part of field are flat, with an average elevation of 265 m (Figure 1). Karamay oilfield is located in the central and eastern parts of the region [26]. The climate is temperate continental, with an average annual precipitation of 110 mm and an annual evaporation of 2690 mm. Karamay oilfield is a typical oil and gas accumulation zone in the piedmont thrust zone of the foreland basin. Its reservoir type is a structural lithologic monoclinic reservoir with fault occlusion [31], and the depth range is 300 to 3000 m [32]. As a typical glutenite reservoir, Karamay oilfield has obvious geological features. First, the reservoir is a sandstone body deposited by alluvial fans, and the rock types are diverse, including conglomerate, sandstone, siltstone, and mudstone. Second, the sedimentary environment is mainly the interbedded deposition of sand, gravel and mudstone. The pores formed by the conglomerate skeleton are partially or completely filled with gravel or clay. Conglomerate reservoirs have low porosity and permeability and are highly heterogeneous [33,34].

3. Methods and Data

3.1. 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 $\sqrt{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=\left({{\displaystyle \sum }}_{i=1}^N{t}_i*p{h}_i\right)/{\displaystyle \sum }_i^N{t}_i^2$$

(1)

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 $p{h}_i$ is the unwrapped phase value of the interferogram.

3.2. 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:

$$\left(\begin{array}{c}-\cos \alpha \sin \theta \Delta t\cos \theta \Delta t\\ \lambda I\end{array}\right)\left(\begin{array}{c}{V}_E\\ V_U\end{array}\right)=\left(\begin{array}{c}D\\ 0\end{array}\right)$$

(2)

where $\alpha$ is the radar azimuth angle, $-$14.2$°$ (ascending) and $-$165.2$°$ (descending) in our case, and $\theta$ 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 $\lambda$ is a regularization parameter, $I$ is an identity matrix. The regularization parameter $\lambda$ is used to regularize the solution. The larger the $\lambda$ 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 $\left(t_1-t_0\right)/\left(t_1-t_{-1}\right)$; similarly the last ascending interferogram is compensated by a factor of $\left(t_5-t_4\right)/\left(t_6-t_4\right)$. 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.

3.3. 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].

4. Results

4.1. 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.

4.2. 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 (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 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.

5. Analysis and Discussion

5.1. 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 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.

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.

5.2. 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.

5.2.1. 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. Optimized parameters estimated with the Mogi model.

Parameter	Optimization (B1)	Optimization (C1)
X (km)	1.61 ± 0.02	0.57 ± 0.01
Y (km)	1.60 ± 0.03	0.67 ± 0.02
Depth (m)	580 ± 30	370 ± 20
Volume change ($\Delta V$, m3)	1.97$\times {10}^5$ ± 1.1 $\times$ ${10}^4$	3.87$\times {10}^4$ ± 5.9 $\times {10}^3$

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 (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.

5.2.2. 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. Optimized parameters estimated with the Sill model.

Parameter	Optimization (B1)	Optimization (C1)
X (km)	1.68 ± 0.02	0.61 ± 0.01
Y (km)	1.63 ± 0.02	0.67 ± 0.01
Depth (m)	606 ± 30	422 ± 30
Strike (°)	270° ± 3°	189° ± 2°
Length (km)	1.00 ± 0.07	0.29 ± 0.03
Width (km)	0.23 ± 0.03	0.26 ± 0.01
Opening (m)	1.26 ± 0.2	0.39 ± 0.05

. 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.

Comparing the results from the two models (Table 1 and Table 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⁵ and 2.15 × 10⁵ m³, respectively) were very similar to each other. Similar to the C1 area, the volume changes (3.2 × 10⁴ to 4.4 × 10⁴ m³) from Mogi model were also close to the volume changes (2.2 × 10⁴ to 3.5 × 10⁴ m³) from Sill model. Unfortunately, we cannot compare our modeled volumes with the actual volume mined, because it is not available in the published data.

5.3. 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 Figure 1 and Figure 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.

6. Conclusions

In this study, we monitored deformation induced by processes associated with oil production at the Karamay oilfield by analyzing Sentinel images acquired from March 2017 to August 2018. We have reached the following main conclusions:

1. The surface deformation caused by oil mining was still significant between 2017 and 2018. Four uplift zones and one settling zone were found in the oilfield. The maximum lifting deformation rate exceeded 160 mm/year, and the maximum subsidence was approximately 32 mm/year. Time series deformation results show that enhanced injection and production activities may cause nonlinear surface deformation.

2. Two-dimensional (vertical and east–west) time series deformation maps were derived by the MSBAS. The results show that the east–west deformation induced by oil mining cannot be ignored in the Karamay oilfield.

3. There was no direct relationship between the occurrence of earthquakes near the oilfield and the activities associated with oil production in Karamay oilfield from 2017 to 2018.

4. If combined with liquid injection data and oil extraction data, time series InSAR results could play a significant role in evaluating future uplift hazards and alerting oil companies of the possible risks.

Unfortunately, the liquid injection data and oil extraction data are not available at this time because of their proprietary. It should be noted that these data should be collected for future research.