Next Article in Journal
Automation Aspects for the Georeferencing of Photogrammetric Aerial Image Archives in Forested Scenes
Previous Article in Journal
Multilayer Perceptron Neural Networks Model for Meteosat Second Generation SEVIRI Daytime Cloud Masking
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Land Subsidence over Oilfields in the Yellow River Delta

1
Key Laboratory for Geo-Environment Monitoring of Coastal Zone of the National Administration of Surveying, Mapping and GeoInformation & Shenzhen Key Laboratory of Spatial Smart Sensing and Services, Shenzhen University, Shenzhen 518060, China
2
College of Information Engineering, Shenzhen University, Shenzhen 518060, China
3
COMET, School of Civil Engineering and Geosciences, Newcastle University, Newcastle upon Tyne NE1 7RU, UK
4
School of Geographical and Earth Sciences, University of Glasgow, Glasgow G12 8QQ, UK
5
The First Institute of Oceanography, State Oceanic Administration, Qingdao 266061, China
*
Author to whom correspondence should be addressed.
Remote Sens. 2015, 7(2), 1540-1564; https://0-doi-org.brum.beds.ac.uk/10.3390/rs70201540
Submission received: 30 September 2014 / Accepted: 27 January 2015 / Published: 2 February 2015

Abstract

:
Subsidence in river deltas is a complex process that has both natural and human causes. Increasing human activities like aquaculture and petroleum extraction are affecting the Yellow River delta, and one consequence is subsidence. The purpose of this study is to measure the surface displacements in the Yellow River delta region and to investigate the corresponding subsidence source. In this paper, the Stanford Method for Persistent Scatterers (StaMPS) package was employed to process Envisat ASAR images collected between 2007 and 2010. Consistent results between two descending tracks show subsidence with a mean rate up to 30 mm/yr in the radar line of sight direction in Gudao Town (oilfield), Gudong oilfield and Xianhe Town of the delta, each of which is within the delta, and also show that subsidence is not uniform across the delta. Field investigation shows a connection between areas of non-uniform subsidence and of petroleum extraction. In a 9 km2 area of the Gudao Oilfield, a poroelastic disk reservoir model is used to model the InSAR derived displacements. In general, good fits between InSAR observations and modeled displacements are seen. The subsidence observed in the vicinity of the oilfield is thus suggested to be caused by fluid extraction.

Graphical Abstract

1. Introduction

Nearly half a billion people live on or near deltas [1] many of which are subsiding as well as being affected by global sea-level rise [2]. Delta subsidence can be due to a range of factors and specific causes lead to spatial and temporal variability in subsidence rates [3]. For example, wetland loss in the Mississippi delta in the 1980s has been linked to hydrocarbon production in Louisiana [4] and Holocene sediment consolidation in the delta [5], amongst others [6]. In the Yellow River delta, formed as the Yellow River enters the Bohai Sea [7], subsidence is likely due to trapping of sediment in 3147 reservoirs [8] in the catchment, aquaculture [9] and hydrocarbon extraction, each of which increases the relative rate of relative sea-level rise.
Repeat pass Interferometric synthetic Aperture Radar (InSAR) is a powerful tool for subsidence mapping [10,11,12,13]. For river delta studies, InSAR has previously been used to measure subsidence due to sediment consolidation and artificial loading in several deltas, including the Pearl River [14], Nile [15] and Fraser River [16] deltas. In the Yellow River delta, Higgins et al. [9] showed subsidence rates as high as 250 mm/yr from InSAR. However, their work only focused on subsidence at aquaculture facilities in the delta.
In this work, land subsidence in oilfields over the Yellow River delta is investigated. This paper attempts to map land subsidence using InSAR time series and confirm the deformation source. We use the Small Baseline Subset (SBAS) [17] method in the Stanford method for persistent scatterers (StaMPS) [18] to investigate surface displacements from 2007 to 2010. The Yellow River delta (YRD) subsidence turns out to be a complex process involving both shallow sediment consolidation and deep oil reservoir compaction, modelled herein using a poroelastic disk reservoir applied to Gudao oilfield.

2. The Yellow River Delta Settings

After a major shift in channel course in 1855, the Yellow River re-entered the Bohai Sea (Figure 1). There have been 10 major channel adjustments from 1855 to the present [19,20], shaping a new 5000 km2 delta [21]. The current cuspate delta formed after the most recent shift of the main channel in 1976, and since 1992 has been the Yellow River Delta National Nature Reserve (YRDNNR) [22]. Extensive hydrocarbon production occurs in the Yellow River delta. Shengli oilfield, which extends from the delta into the Bohai Sea, is the fourth largest oilfield in China with annual production of 27 million tons and production facilities located within the natural reserve [23].
Figure 1. Location of the study area (black dashed rectangle). The red and blue dashed rectangles are Track 132 and Track 404 from Envisat Satellite, respectively. Sampling locations in the Yellow River delta are marked. Background is SRTM [24]. White areas are water surfaces (ocean and ponds). The location of the Yellow River is highlighted on the DEM by dikes built for flood control.
Figure 1. Location of the study area (black dashed rectangle). The red and blue dashed rectangles are Track 132 and Track 404 from Envisat Satellite, respectively. Sampling locations in the Yellow River delta are marked. Background is SRTM [24]. White areas are water surfaces (ocean and ponds). The location of the Yellow River is highlighted on the DEM by dikes built for flood control.
Remotesensing 07 01540 g001

3. Surface Displacements from InSAR

3.1. Radar Interferometry and Time Series Analysis

The differential InSAR technique uses two or more SAR complex images to generate Earth surface movements with sub-centimetre accuracy over a wide area (typically 100 ×100 km) [25]. However, the accuracy is affected by several error sources: DEM error, atmospheric effects, orbit error and other noise [26]. InSAR time series technique utilizes a series of radar interferograms to separate those error terms from line of sight (LOS) displacement time series. Those time series methods are referred to as Persistent Scattering (PS), Small Baseline Subset (SBAS), SqueeSAR or Temporarily Coherent Pixel (TCP) methods [17,18,27,28,29,30,31,32,33].
In this study, ground movements were obtained using a small baseline approach in StaMPS [17,18]. In contrast to most other time series methods, StaMPS uses phase spatial correlation rather than amplitude analysis to identify stable pixels, and does not make any prior assumption about the temporal nature of ground deformation.
Slowly decorrelating filtered phase (SDFP) pixels are selected using phase analysis, which is refined in a series of iterations. Phase analysis aims to evaluate the phase stability of a SDPF candidate. The phase analysis includes three steps: estimation of the spatially correlated part of phase values; look angle error estimation; and, statistics of the gamma values (time series coherences). Firstly, the spatially correlated deformation, atmospheric effects and orbit errors are estimated together through band-pass filtering, leaving the spatially uncorrelated phases. The strategy for step one is not to separate deformation from atmospheric effects or orbit errors. Rather, deformation, atmospheric effects and orbit errors are thought of as a single phase term. The filter consists of an adaptive phase filter and a low-pass filter with a cut-off wavelength of 800 m, which is the recommended distance for spatially correlated phase values [18]. Secondly, the spatially uncorrelated phase series and its perpendicular baselines series are used to derive the look angle error for each pixel using the least squares method. Thirdly, subtracting the spatially correlated phase and look angle error from the original phase, the residual phase is used to estimate the gamma value for each pixel which represents the noise level of pixels in the time series. SDPF candidates with higher gamma values tend to be less affected by atmospheric effects, orbit inaccuracies and look angle errors. Hence, they are selected as the SDFP pixels.
Since the interferometric phase values are modulo 2π, they are unwrapped using the three-dimensional phase unwrapping algorithm [34], which incorporates the two-dimensional SNAPHU algorithm [35]. SNAPHU is a statistical-cost, network flow algorithm. The algorithm aims to compute the most likely unwrapped solution with maximum posterior probability estimation given the observable input data [35].
For the small baseline approach, an inversion of deformation increments from the unwrapped interferograms is necessary as multiple master images are used and the phases of them are not referenced to the same time [36]. The inversion is implemented using linear equations. The inverted phase and displacement are also known as modeled phase and displacement. The long-wavelength atmospheric effects and orbit errors are estimated via a best-fit plane for the modeled phases [37]. Higher deformation consistency between adjacent descending tracks was achieved when a best-fit plane was used to account for atmospheric effects.

3.2. Data

Twenty-four Envisat SAR images from descending Track 132 (Table A1 in Appendix A) and 13 Envisat SAR images from descending Track 404 (Table A2 in Appendix A) were used to determine surface displacements over the delta region (Figure 1). The SAR images were used to generate small baseline interferograms (Figure 2). Note that the incidence angles of tracks 132 (swath I2) and 404 (swath I1) are 4° different and so displacement differences are expected. However, it is still difficult to recover the horizontal motions since both Track 132 and Track 404 are descending tracks. In Section 3.3, line of sight displacements are given. In Section 5, the line of sight displacements in Gudao oilfield are modeled considering both horizontal and vertical movements (Equations (1) and (2)) and the look vector.
Figure 2. Small Baseline network. (a) Track 132. (b) Track 404. Circles denotes SAR images and black lines represent SAR interferograms. The spatial and temporal baseline thresholds for image pair selection are 400 m and 180 days, respectively, except when image pairs are adjusted manually to achieve a well constrained and connected network for each acquisition.
Figure 2. Small Baseline network. (a) Track 132. (b) Track 404. Circles denotes SAR images and black lines represent SAR interferograms. The spatial and temporal baseline thresholds for image pair selection are 400 m and 180 days, respectively, except when image pairs are adjusted manually to achieve a well constrained and connected network for each acquisition.
Remotesensing 07 01540 g002

3.3. InSAR Time Series Results

The processed ASAR images cover an area of about 2200 km2. However, the three areas (Areas 1, 2, & 3 in Figure 3) are not connected with each other, leaving 10 km wide gaps with only individual SDFP pixels detected. In Area 3, the YRDNNR can be seen as a group of isolated islands of SDFP pixels, some of which are facilities of Shengli oilfield located in the natural reserve. In Area 2, the aquaculture facilities produce some radar signals from the concrete edges of the ponds. The connections between SDFP pixels in aquaculture facilities are also poor. Phase jumps are seen between the three areas and even within the spatial extent of YRDNNR or aquaculture facilities after unwrapping, indicating phase unwrapping errors. Hence, Area 2 and Area 3 cannot be included in the analysis. Coherence is a limitation for InSAR application in Yellow River delta because it affects the density of stable pixels that can be detected. For phase unwrapping purposes, the phases of neighboring pixels are required to be within the same phase cycle. If the neighboring SDFP pixels are far away from each other, it is uncertain if their wrapped phases are within the same cycle.
Figure 3. SDFP pixels detected in the Yellow River delta. The SDFP pixels have three clusters: Areas 1, 2, & 3. The Gudao Town (oilfield), the Xianhe Town and the Gudong oilfield to the north of the current Yellow River channel are in Area 1. Aquaculture facilities to the south of current Yellow River channel are found in Area 2. The Yellow River Delta National Natural Reserve (YRDNNR) is in Area 3.
Figure 3. SDFP pixels detected in the Yellow River delta. The SDFP pixels have three clusters: Areas 1, 2, & 3. The Gudao Town (oilfield), the Xianhe Town and the Gudong oilfield to the north of the current Yellow River channel are in Area 1. Aquaculture facilities to the south of current Yellow River channel are found in Area 2. The Yellow River Delta National Natural Reserve (YRDNNR) is in Area 3.
Remotesensing 07 01540 g003
Figure 4. Correlation between line of sight displacement rates estimated from time series of the two tracks. Scatter plot showing all pixels common to both tracks 404 and 132. Best-fit line is shown (solid line) along with the 1:1 line (dashed). The RMS of the residuals from the best-fit line is 2.68 mm/yr.
Figure 4. Correlation between line of sight displacement rates estimated from time series of the two tracks. Scatter plot showing all pixels common to both tracks 404 and 132. Best-fit line is shown (solid line) along with the 1:1 line (dashed). The RMS of the residuals from the best-fit line is 2.68 mm/yr.
Remotesensing 07 01540 g004
Figure 5. (a) Track 132: Mean displacement rates from SBAS analysis. (b) Track 404: Displacement rates from SBAS. Positions of three 500 meter wide swaths A-A", B-B" and C-C" are marked in (a). The reference area (small) and model area (large) are shown by rectangles in (a). Background is mean amplitude of SAR images. Area M, N, and Q are the locations of Gudao Town (oilfield), Gudong Town and Xianhe Town, respectively.
Figure 5. (a) Track 132: Mean displacement rates from SBAS analysis. (b) Track 404: Displacement rates from SBAS. Positions of three 500 meter wide swaths A-A", B-B" and C-C" are marked in (a). The reference area (small) and model area (large) are shown by rectangles in (a). Background is mean amplitude of SAR images. Area M, N, and Q are the locations of Gudao Town (oilfield), Gudong Town and Xianhe Town, respectively.
Remotesensing 07 01540 g005
The mean displacement rates in Area 1 over the full time period from Tracks 132 SBAS and Track 404 SBAS cases are compared in Figure 4, and shown in Figure 5a and Figure 5b, respectively. As no independent ground truth data are available in this area, the mean values of each image are firstly used as the reference phases in the time series analysis. The area without obvious displacements on the resulting rate map was then chosen as the reference area (small rectangle in Figure 5a). Finally, the displacement time series were recalculated with the newly referenced phases.
Correlation between common SBAS pixels from Tracks 132 and 404 is 0.72 (Figure 4). Differencing the mean displacement rates between the tracks shows that 35.6% of the common pixels have absolute differences less than 1 mm/yr and 61.9% less than 2 mm/yr. The overall root mean square (RMS) difference is 2.68 mm/yr between the two tracks. The displacement differences due to different image geometry conditions from Tracks 132 and 404 are expected (Appendix B). The discrepancies seen from the two tracks can also be related to differences in temporal sampling of the deformation, atmospheric effects and orbital errors.
Figure 6. Displacement rates of the Yellow River delta estimated from small baseline interferograms for Tracks 132 and 404. Positions of three 500 meter wide swaths A-A", B-B" and C-C" are marked in Figure 5a. Note the consistency in results from the two tracks.
Figure 6. Displacement rates of the Yellow River delta estimated from small baseline interferograms for Tracks 132 and 404. Positions of three 500 meter wide swaths A-A", B-B" and C-C" are marked in Figure 5a. Note the consistency in results from the two tracks.
Remotesensing 07 01540 g006
To further assess the accuracy of the InSAR time series results, the displacements in the area of overlap between the two adjacent tracks can be compared to calculate the RMS of their differences. Since the acquisition times are different, RMS values were calculated for the differences between the displacement values measured from Track 132 and interpolated displacement values from Track 404 at the SAR acquisition times of Track 132. The RMS differences between Tracks 132 and 404 are 4.5 mm for Gudao1 and 3.8 mm for Gudong1, respectively. For these two locations, these RMS values are 6~8% of the net displacements. Pixels with time series RMS differences between the independent tracks below 3 mm, 5 mm and 10 mm, account for 11%, 60%, and 97% of the total common pixels, respectively. A 5 mm RMS is set as the variance to simulate a variance-covariance matrix (VCM) using an exponential function. The VCM is used to calculate correlated noise using Cholesky decomposition. The correlated noise is added back to the original data to provide a new noisy data set in model test below [38].
Subsidence areas can be identified from the displacement rates (Figure 5a,b). Area M is the location of Gudao Town (oilfield). Swath A-A" shows subsidence in this area with a maximum rate of about 30 mm/yr (Figure 6a). The width of the subsiding bowl cannot be determined from Swath A-A" because there is a lack of stable pixels further east. Swath B-B" also shows this area of subsidence (Figure 6b). Area N, crossed by swath C-C", is the Gudong Oilfield. Two peak areas of subsidence, from 1 km to 6 km with maximum rates of 15 mm/yr, and from 7 km to 10 km with maximum rates of 20 mm/yr, can be identified (Figure 6c).
Figure 7. (a) Time series of surface displacement at Gudao1 and Gudong1 from Tracks 132 and 404. Error bars represent the standard deviations of the displacement values for all SDFP points within a 200 m × 200 m window, centered with the given point. Displacement offset between Tracks 132 and 404 are due to different reference dates used for different tracks. (b) Displacements along Swath A-A" from Tracks 132 and 404 in three interferograms. Locations of Gudao1, Gudong1 and Swath A-A" are given in Figure 5a.
Figure 7. (a) Time series of surface displacement at Gudao1 and Gudong1 from Tracks 132 and 404. Error bars represent the standard deviations of the displacement values for all SDFP points within a 200 m × 200 m window, centered with the given point. Displacement offset between Tracks 132 and 404 are due to different reference dates used for different tracks. (b) Displacements along Swath A-A" from Tracks 132 and 404 in three interferograms. Locations of Gudao1, Gudong1 and Swath A-A" are given in Figure 5a.
Remotesensing 07 01540 g007
Both Tracks 132 and 404 show progressive displacements at Gudao1 and Gudong1 away from the satellite, indicating net subsidence of about 50 mm in three years (Figure 7a). The displacements in three different interferograms along Swath A-A" are given for the two tracks (Figure 7b). The maximum displacement of ~80 mm is reached at ~5.5 km in both tracks and is near symmetrical. Even in cases with totally symmetrical deformation, the line of sight displacement would show a greater gradient on one side due to horizontal ground motion. Figure 7b shows a greater gradient of subsidence to the west of the peak, as expected from a descending track. However, due to the lack of ascending data, the horizontal component of the displacement cannot be constrained. Alternatively, modeled horizontal displacements in Gudao oilfield reach up to 10 mm in radar LOS direction (Figure 9h in Section 5).

4. Field Investigation and Sediment Analysis

To identify possible reasons for the observed subsidence, field investigation was conducted in the Yellow River delta. The Gudao Town, Xianhe Town and Gudong oilfield to the north of the current Yellow River channel, the aquaculture facilities to the south of current channel and the Yellow River Delta National Nature Reserve (YRDNNR) in the cuspate delta, were inspected. The areas showing coherent subsidence bowls from the InSAR observations are associated with the Gudao and Gudong oilfields, in which oil extraction dates back to the late 1960s and late 1980s, respectively. Oil extraction is considered to be the major cause of the areas of coherent subsidence observed in the Yellow River delta.
Near the oil field are areas of both arable farmland and residential land. Farming, due to seasonal water use, may induce subsidence. However, the data shows linear subsidence with no seasonal pattern in Gudao oilfield. Hence, the deformation in the oilfield is related to continuous hydrocarbon extraction rather than to water use for farming.
Areas of subsidence may also be related to sediment properties so shallow sediment samples were taken from the delta (See Figure 1 for sample locations) and analysed for sediment properties (Figure C1), including particle size distribution, mineralogy, roundness and shape. The overall consistency in sediment properties (Appendix C) suggests that differences in material composition are not the cause of the non-uniform displacement signals observed from InSAR.

5. Modeling Subsidence in Gudao Oilfield

The Gudao oilfield, located in the Jiyang depression in the Bohaiwan Basin has reserves are 4 × 108 tons of oil and 47 × 108 m3 of natural gas. The tertiary sandstone reservoir is a draped anticline [39,40]. The Gudao oilfield is controlled by the Gunan fault to the south and the Gubei fault zone to the north, and is gently dipping on its east and west sides. This structure facilitates the migration of petroleum from Bonan, Gunan and Gubei hydrocarbon sources to the Gudao oilfield [39].
The cause of subsidence in Gudao oilfield is analysed here through poroelastic modeling. Both elastic and poroelastic models can be used to determine subsidence sources. Elastic models assume a continuous elastic medium in which the fluid exists within a cavity (e.g., a magma chamber) in homogenous isotropic elastic half space [41]. However, hydrocarbons in reservoir rocks exist in a more distributed way. Poroelastic models are used in these cases. The poroelastic theory independently proposed by Biot [42] and Terzaghi [43] was further developed by Nur and Byerlee [44] into an exact effective stress law for deformation of rocks with fluid. Segall [45] obtained a poroelastic analytical solution in an axisymmetric space to model the subsidence in a hydrocarbon reservoir in the Lacq gas field [46]. Similar expressions to model subsidence caused by oil extraction were also given by Geertsma [47] for other material properties.
The area of the Gudao oilfield with near axisymmetric displacement outlined in Figure 5a is modeled using a poroelastic disk reservoir model [45,46]. The vertical (uz) and radial (ur) surface displacements (z = 0) are obtained in terms of the product of zero order and first order Bessel functions J0 and J1.
u z ( r , 0 ) = α ( 1 2 v ) T Δ P μ R 0 J 1 ( k R ) J 0 ( k r ) e k d d k
u r ( r , 0 ) = + α ( 1 2 v ) T Δ P μ R 0 J 1 ( k R ) J 1 ( k r ) e k d d k
The disk radius R, thickness T, and center depth d are the geometrical parameters of the poroelastic reservoir. The reservoir’s Biot coefficient α, which is a measurement of the compressibility of the soil for a change in water pressure [42,44], Poisson’s ratio ν, shear modulus μ and the pressure change ΔP are the material parameters. The Biot coefficient is an increasing function of porosity with values approaching zero at low porosity while they approach one at high porosity. When integration is performed, k represents nuclei of strain.
Gudao oilfield is an unconsolidated sandstone reservoir with 30–33% porosity, with the depths of oil bearing Neogene Guantao formation ranging from 1120 to 1350 m, and a pressure of 12 MPa [39,48].
The Biot coefficient, Poisson’s ratio, shear modulus and the reservoir thickness of poroelastic disk reservoir model are taken from Sun et al. [49] and fixed in the model (Table 1). The center (x, y, d) of the reservoir, reservoir radius R, and the pressure change ΔP within the reservoir are modeled from InSAR observations. The modeled reservoir depth d is then compared to the published depth of 1235 m.
Table 1. Parameters used in model calculation. The shear modulus for a Poisson’s ratio of 0.26 is 298 MPa. The reservoir thickness is 230 m for 11,200–1350 depth range.
Table 1. Parameters used in model calculation. The shear modulus for a Poisson’s ratio of 0.26 is 298 MPa. The reservoir thickness is 230 m for 11,200–1350 depth range.
SymbolQuantityValue
αBiot coefficient0.85
νPoisson’s ratio0.26
μShear modulus298 MPa
TReservoir thickness230 m
The inversion is implemented using a simulated annealing algorithm to find the optimal parameters within a space limited by lower and upper boundaries by minimizing the misfit between model predictions and data [50]. Simulated noise was added to the observed displacement field and the inversion was repeated. The 5 mm RMS estimated from InSAR time series displacement is used to construct a full variance-covariance matrix (VCM). The VCM matrix is decomposed to simulate spatially correlated noise [38]. A total of 100 trials were used to generate 100 groups of best fit model parameters, from which the joint sensitivity of the estimated parameters was assessed.
When the InSAR displacement is inverted directly, the estimated depth is 1722 m, which is 487 m deeper than the published depth of 1235 m. This overestimation could be caused by an unknown global bias due to an inaccurate reference level in the InSAR data [51]. In this case, the global bias could be incorporated into the model by co-estimating an offset between the model and the InSAR displacement [52] (case A). Alternatively, the mean phase of each image is again used as the reference phase for InSAR (case B). The two methods actually produce similar results (Table 2).
Table 2. Best Fit Poroelastic disk reservoir model parameters and uncertainty. A: offset between InSAR and model estimated. B: mean value as InSAR reference.
Table 2. Best Fit Poroelastic disk reservoir model parameters and uncertainty. A: offset between InSAR and model estimated. B: mean value as InSAR reference.
Quantity (unit)Disk Centre East (km)Disk Centre North (km)Disk Centre Depth (km)Disk Reservoir Radius (km)Pressure Change (MPa)Offset (mm)Misfit (mm)
SymbolxydRΔPCRMSE
Offset estimated1.56±0.052.48±0.051.08±0.090.83±0.110.80±0.2120.0±0.48.1±0.2
Mean value reference1.57±0.052.41±0.041.14±0.090.91±0.090.73±0.13N/A8.3±0.2
After consideration of the global bias, the estimated depths are of 1080 and 1140 m in categories A and B, 100–160 m shallower than the published depth of 1235 m. These differences are close to the uncertainty level of 90 m (Table 2). The modeled pressure changes are 0.8 and 0.73 MPa in Categories A and B, respectively, but pressure change data are unavailable for comparison.
Figure 8. Matrix plot for 100 groups of best fit poroelastic disk reservoir parameters for InSAR derived displacement with mean phase value as reference and with added noise in Gudao Town (oilfield). A well confined parameter shows symmetric distribution in histogram or dots plot. A wide scatter means the parameter is unstable.
Figure 8. Matrix plot for 100 groups of best fit poroelastic disk reservoir parameters for InSAR derived displacement with mean phase value as reference and with added noise in Gudao Town (oilfield). A well confined parameter shows symmetric distribution in histogram or dots plot. A wide scatter means the parameter is unstable.
Remotesensing 07 01540 g008
The parameters in Category B (mean value as InSAR reference) are plotted against each other to examine potential trade-offs (Figure 8). The model parameters are well-constrained with symmetrical distributions meaning that there are no major trade-offs, except as follows. Note that there is a linear trade-off between the pressure decline and the depth of the reservoir when the reservoir gets deeper. Generally, a deeper reservoir needs greater pressure change to produce the same amount of surface displacement as a shallow one. Another trade-off is observed between the pressure decline and the radius of the reservoir when the lateral dimension of the reservoir gets greater. Accordingly, a greater reservoir requires less pressure change for the same amount of surface displacement than a small reservoir. Actually, the volume (equivalent to the radius here) and pressure change are treated together in an elastic Mogi type model and their product is defined as the source strength [51].
The model underestimates the displacement in the SW corner of the mapped area (Figure 9). This area of displacement is separated from the central area of subsidence and may be related to a separate source. The central area of near symmetrical displacement of 3 × 3 km with deformation up to 67 mm is relatively well modeled. Along the profile PP", we can group observations into three clusters from SW to NE (Figure 9d). The left cluster (0–1 km along the profile) shows 10 mm subsidence. The model systematically overestimates subsidence in this region (positive residuals). The central cluster (1–2.5 km) shows clear subsidence and exhibits a displacement gradient. The model generally underestimates subsidence leading residuals up to 10 mm. The right cluster (2.5–4 km) has only a few coherent pixels and these show less subsidence than the central cluster. Overall, the three clusters suggest a bowl shape subsidence field, consistent with the model, although with steeper sides than the model predicts. Profile SS" shows a higher displacement gradient and the model systematically overestimates displacement away from the center of deformation. This misfit is related to the non-axisymmetric pattern of deformation and could be due to heterogeneity of the pressure distribution in the reservoir that could be solved by considering non-uniform pressure gradients. This requires that the actual pressure change distributions within the reservoir to be known, which is not the case here.
Figure 9. (a) InSAR derived line of sight cumulative displacement field in Gudao Town between Feb 7, 2007 and Jan 21, 2010 using the mean value as InSAR reference. (b) Best approximated displacement field using a disk reservoir model. (c) Residuals between InSAR and modeled displacements. (d) Profile P-P” through the data (triangle), the model (line), and the residuals (crosses). (e) Same for the profile S-S”. (f) Modeled vertical and horizontal displacement seen from radar LOS direction. (g) Modeled vertical displacement seen from radar LOS direction. (h) Modeled horizontal displacement seen from radar LOS direction. (i) Profile P-P” through both vertical and horizontal (black), only vertical (red), and only horizontal (blue) displacements. (j) Same for the profile S-S”.
Figure 9. (a) InSAR derived line of sight cumulative displacement field in Gudao Town between Feb 7, 2007 and Jan 21, 2010 using the mean value as InSAR reference. (b) Best approximated displacement field using a disk reservoir model. (c) Residuals between InSAR and modeled displacements. (d) Profile P-P” through the data (triangle), the model (line), and the residuals (crosses). (e) Same for the profile S-S”. (f) Modeled vertical and horizontal displacement seen from radar LOS direction. (g) Modeled vertical displacement seen from radar LOS direction. (h) Modeled horizontal displacement seen from radar LOS direction. (i) Profile P-P” through both vertical and horizontal (black), only vertical (red), and only horizontal (blue) displacements. (j) Same for the profile S-S”.
Remotesensing 07 01540 g009aRemotesensing 07 01540 g009b

6. Discussion

6.1. Subsidence Patterns in Gudao Town (Oilfield), Gudong Oilfield and Xianhe Town

InSAR was used to monitor the surface movements in Gudao Town (oilfield), Gudong Town and Xianhe Town of the Yellow River delta. Subsidence bowls were observed in Gudao and Gudong oilfields with linear subsidence rates over three years. The mean rates reached about 30 mm/yr in Gudao oilfield with a bowl shaped rate map where rates gradually decrease towards the margins. The data do not suggest any other cause for surface subsidence than oil extraction. In the Gudao oilfield, using a poroelastic disk reservoir model, a deformation offset of 20 mm is well reproduced.
Subsidence in Gudong oilfield does not have a regular shape, and can be thought of as two separate subsidence centers with mean rates of about 15 mm/yr (northern) and about 20 mm/yr (southern). These subsidence areas are prone to flooding and are difficult to drain resulting in significant hazards.
In Xianhe Town, the majority of pixels show subsidence of 8~12 mm/yr (within the spatial gradient of D-InSAR and above the noise threshold for D-InSAR time series), and that do not have the bowl shape that characterizes areas of fluid extraction [53]. It suggests that other factors like sediment consolidation and loading may also contribute to the total subsidence in this part of the Yellow River delta. These subsidence rates will decrease through time if they are caused by loading [16].

6.2. Subsidence Model in Gudao Oilfield from Surface Displacement

The poroelastic disk reservoir model predicts depths for the fluid extraction that are close to the actual oilfield depths. Other parameters derived from the model inversions need to be examined against field data. In the Lacq gas field, it took about ~30 years to have a 60 MPa pressure change and 60 mm maximum subsidence with a linear maximum rate of 1 mm/MPa in the center field and pressure change of 2 MPa/yr [46]. The three years of subsidence due to oil extraction in Gudao oilfield resulted in 62 mm of displacement, on the same order as the subsidence in Lacq. However, the estimated pressure change of 0.73–0.8 MPa in Gudao oilfield is much smaller than in Lacq. The shear modulus used for the Gudao and Lacq reservoirs are 0.3 GPa (high porosity sandstone reservoir) and 23 GPa (carbonate reservoir), respectively. So, it is suggested that less pressure change is needed in the sandstone reservoir of Gudao oilfield than in the Lacq reservoir to produce the same subsidence.
Improvements can be made by improving the geometry assumed in the poroelastic model. Such an approach is usually implemented using finite element numerical models (FEM). Several studies using finite element modelling for oilfield subsidence have been published [54,55,56,57,58]. Implementation of FEM requires a static 3D model of the reservoir and of the surrounding region, describing all its geological, lithological, stratigraphical and petrophysical aspects (e.g., the shapes of the layers and the trend of the faults, initial porosity and permeability), and a dynamic model regarding the characteristics of the fluids, the rock and the well system [58]. This approach requires extensive data that remains unknown or is commercially confidential for the Yellow River delta.

6.3. Limitations for Applying InSAR in the Yellow River Delta

It should be noted that the SDFP pixels detected in the Yellow River delta region mainly occur at the locations of roads, structures, oil fields and towns. Most of the natural and farmed vegetated areas give no stable radar back scattering, leaving empty zones in the InSAR deformation maps. Hence, the subsidence rates in areas with no or sparsely distributed SDFP points remain unknown. Accordingly, subsidence at low rates due to sediment compaction in these areas may be missed. Nevertheless greater coherence can be achieved using a longer wavelength satellite (e.g., ALOS data archive from 2006 to 2011 and its current follow-on mission ALOS-2). Alternatively, using shorter revisit intervals can also improve the coherence. The newly launched Sentinel-1 satellite with a 12 day repeat cycle could enhance the performance of C band data in the Yellow River delta. Hence the application of ALOS1/2 and Sentinel-1 data in Yellow River delta will be of interest to detect displacement occurring over a wider area than is reported in this paper.

7. Conclusions

Two tracks of Envisat ASAR images are employed to investigate the subsidence in the Yellow River delta. Using InSAR time series techniques, three subsidence locations (Gudao, Gudong and Xianhe) in the Yellow River delta are detected. Cross validation between the two independent tracks show that the displacement rates and time series are reliable. Maximum subsidence rates of 30 mm/yr, 20 mm/yr and 12 mm/yr are observed for Gudao, Gudong and Xianhe, respectively.
Shallow sedimentary samples from the Yellow River delta are analysed for sediment properties, including particle size distribution, mineral composition, roundness and shape. Consistent sediment properties reveal that shallow delta materials are not the cause of InSAR observed displacement differences.
The primary cause of land subsidence in two locations (Gudao and Gudong) is suggested to be the oil extraction. Oil extraction reduces the oil reserve in the reservoir, decreasing the pressure within the reservoir. Reduced pore pressure increases the effective stress and cause the reservoir to shrink. Reservoir compaction will in turn cause subsidence at the Earth surface. Poroelastic models well represent this kind of deformation mechanism. The disk reservoir is in good agreement with the geometry of Gudao oilfield and its analytical solutions to the Earth surface displacement are available. Hence, a poroelastic disk reservoir is used to model subsidence due to oil extraction in Gudao oilfield.
The best fit poroelastic disk reservoir parameters for Gudao oilfield are estimated from InSAR observed displacement through model inversion. The accuracy of model parameters is examined when the inversion was repeated. Similitude between InSAR observed subsidence and the best fit poroelastic model confirms the dominance of oil exploitation in total subsidence of Gudao oilfield.
Oil extraction induced pressure decline of 0.73–0.8 MPa is estimated in the sandstone reservoir of Gudao oilfield. Unfortunately, actual pressure decline is unavailable for comparison. Oil extraction induced maximum subsidence of 62 mm is estimated in Gudao from the best fit model parameters. More uniform signals like the 20 mm deformation offset estimated in Gudao from the model, and subsidence in Xianhe (8~12 mm/yr) are related to shallow sediment compaction or loading.
Subsidence mapping in Yellow River delta are limited by the number of SDFPs due to spatial and temporal decorrelations. Radar data with longer wavelength and shorter revisit time will exhibit better coherence and raise density of SDPFs across the delta. Hence, longer wavelength L band ALOS 1/2 images and shorter revisit time C band Sentinel-1 images are of interest to map subsidence in a wider area of the Yellow River delta. Advanced subsidence models for the Yellow River delta can be implemented using numerical methods if currently confidential data become available.

Acknowledgments

This work is supported by a joint University of Glasgow/China Scholarship Council (GU/CSC) scholarship program to PL. The ENVISAT images were supplied through the ESA-MOST Dragon 2 Cooperation Program (ID: 5343). Field work was funded by a grant from the Carnegie Trust for the Universities of Scotland to ZL. QL is supported by grants from Shenzhen Scientific Research and Development Funding Program (No. ZDSY20121019111146499, No. JSGG20121026111056204), and Shenzhen Dedicated Funding of Strategic Emerging Industry Development Program (No. JCYJ20121019111128765). We thank JPL/Caltech for the use of ROI PAC, TU-Delft for DORIS and Andy Hooper for StaMPS in our data processing and analysis. Figure 1 were prepared using the public domain Generic Mapping Tools [59]. We are grateful to Callum Dowds, Peter Chung and Kenny Roberts for carrying out particle size and SEM analyses, and Andrew Singleton for helpful discussions. Constructive comments from L. Zhang, E. Chaussard, and other anonymous reviewers significantly helped to improve the manuscript.

Author Contributions

Peng Liu carried out the experiments and drafted the manuscript. Qingquan Li, Zhenhong Li, and Trevor Hoey provided useful suggestions and improved the manuscript. Yanxiong Liu hosted the field trip in Yellow River delta and facilitated sediment sampling. Chisheng Wang contributed to the model inversion.

Conflicts of Interest

The authors declare no conflict of interest

Appendix A

Table A1. SAR images used from Track 132 (Reference date: 26-Jul-2007).
Table A1. SAR images used from Track 132 (Reference date: 26-Jul-2007).
Scene No.Acquisition DateBaseline Perpendicular (m)Scene No.Acquisition DateBaseline Perpendicular (m)
101-Feb-2007−3721301-Jan-2009158
208-Mar-20074811405-Feb-2009−292
312-Apr-2007 −1851512-Mar-2009520
417-May-2007−1361616-Apr-2009−208
521-Jun-2007 71 1721-May-2009−24
626-Jul-200701825-Jun-2009342
708-Nov-2007 2021930-Jul-2009−6
817-Jan-2008−412003-Sep-2009 312
921-Feb-2008−3242108-Oct-2009−99
1027-Mar-20083472212-Nov-2009289
1123-Oct-20082442317-Dec-2009−308
1227-Nov-2008−2502421-Jan-2010267
Table A2. SAR images used from Track 404 (Reference date: 05-September-2009).
Table A2. SAR images used from Track 404 (Reference date: 05-September-2009).
Scene No.Acquisition DateBaseline Perpendicular (m)Scene No.Acquisition DateBaseline Perpendicular (m)
116-Jan-2007592707-Oct-200875
227-Mar-2007535816-Dec-2008−93
305-Jun-2007 125931-Mar-2009 742
401-Jan-2008−2751005-May-20090
505-Feb-20082001109-Jun-2009336
615-Apr-20083351227-Oct-200976

Appendix B

Supposing that the actual and two line of sight displacements are d, d1 and d2 respectively, the LOS displacements differences (Δd) between two tracks, and the ratio (rΔθ) of LOS displacement differences to the actual displacement are:
Δ d = d 1 d 2 = d c o s ( θ i θ d ) d c o s ( θ i θ d + Δ θ )
r Δ θ = Δ d / d = c o s ( θ i θ d ) c o s ( θ i θ d + Δ θ )
where θ i is radar incidence angle for one track, while Δθ is the radar incidence angle difference between two tracks, and θd is the deformation vector angle off nadir (Figure B1a). In this case, θ i is in the range (15.0, 22.9) for Track 132, Δθ is ~4° for Track 132 and 404, while θ d is in the range (−180°, 180°) for any possible deformations. Hence the range for r Δ θ is (−0.07, 0.07) (Figure B1b). The displacement differences due to different image geometry conditions from Tracks 132 and 404 are less than 7% of the actual displacement.
Figure B1. The ratio (rΔθ) of LOS displacement differences to the actual displacement. (a) incidence angle geometries. (b). rΔθ in terms of θd and θi in this case. (c) rΔθ in profile HH". (d). rΔθ in profile VV".
Figure B1. The ratio (rΔθ) of LOS displacement differences to the actual displacement. (a) incidence angle geometries. (b). rΔθ in terms of θd and θi in this case. (c) rΔθ in profile HH". (d). rΔθ in profile VV".
Remotesensing 07 01540 g010

Appendix C

Twenty-one sediment samples were collected from water channels, pond banks, reed marsh, cotton fields and soil blocks left by excavators in the Yellow River delta in September 2010 (Figure 1). Grain size was determined using a Coulter LS 230 laser diffraction particle size analyser [60], each sample being run three times. At least two sub-samples from each sample were analysed with further tests carried out for samples with unusual particle size distribution patterns. Scanning Electron Microscopy (SEM) imaging was used to analyse the mineralogy, shape, and roundness of the samples.
Samples P15, P33 and P34 include occasional coarse grains ranging up to 600 µm, while 97.8–99.9% of the sediment lies in the range from 0.04 to 200 µm (Figure C1).
The mean particle size of 18 samples lies within the silt classification (5.6–46.1µm) with modal sizes from 4 to 60 µm. The mean particle size for sediment samples collected over the delta is 25.9 µm (Table C1). The particle size and mineralogy supports the view that the provenance of sediments in the Yellow River delta is the Loess Plateau [8].
Figure C1. (a) SEM image showing a quartz grain in sample P24. This grain exhibits the dominant angular, tabular form of quartz grains from the delta. Some clay particles are seen on the central grain and on the silt-sized particle behind. The growth on the top right particles can be salt crystallisation. (b) Differential volume particle size for all samples. Differential volume particle size distribution is the fractional volume of particles of given size, and is equivalent to a distribution by mass if particle density is equal for all sizes [61].
Figure C1. (a) SEM image showing a quartz grain in sample P24. This grain exhibits the dominant angular, tabular form of quartz grains from the delta. Some clay particles are seen on the central grain and on the silt-sized particle behind. The growth on the top right particles can be salt crystallisation. (b) Differential volume particle size for all samples. Differential volume particle size distribution is the fractional volume of particles of given size, and is equivalent to a distribution by mass if particle density is equal for all sizes [61].
Remotesensing 07 01540 g011
The mineralogy of the delta sediments was assessed from SEM images of six samples (Table C2) that were quantified using ImageJ [62], Quartz (Figure C1a) is found to be the dominant mineral (60~65% abundance) which agrees well with the mineralogy of loess seen by Nemecz et al. [63]. Calcite, feldspar and biotitie mica are also found with abundances of 5–8%, 5–10% and 2–5%, respectively. Clay minerals account for 20–25%, mainly smectite and illite.
About 10% of the particles are high sphericity and the other 90% particles have low sphericity, e.g., the tabular quartz particle in Figure C1 [64]. Within the low sphericity class, the sub-angular, angular and very angular particles account for 16%, 38% and 36% of the particles, respectively (Table C3).
All of the delta samples show similar physical properties to the Loess Plateau deposits, including the dominant quartz composition of 60~65% (Figure C1), dominance of silt size particles (4–60 μm) and their angular, tabular shape. This suggests that (i) the Loess Plateau is the source of the delta sediments, and (ii) the particles are not significantly rounded during transport from the source to the delta, although the increase in quartz abundance does indicate some increase in sediment maturity. The dominance of angular particles suggests limited abrasion during transport, and Russell [65] noted that fine particles round less by abrasion than do larger ones. In addition, some shape-sorting occurs during transport. The dominance of low-sphericity grains here suggests that suspension is the dominant transport mode for the delta sediments [65,66,67].
Table C1. Particle size statistics for shallow samples from the Yellow River delta. Sorting, skewness and kurtosis were calculated using the formulae given by Folk and Ward [68]. Note that phi units are defined as p h i = l o g 2 D ie log to the base 2 of the size D, which is in mm. Sorting is the standard deviation of particle size. The skewness characterizes the asymmetry of a distribution: skewness > 0 signifies a distribution with a tail towards larger sizes; skewness < 0 has a tail towards smaller sizes [69]. The kurtosis is a measure of grains around the median [70].
Table C1. Particle size statistics for shallow samples from the Yellow River delta. Sorting, skewness and kurtosis were calculated using the formulae given by Folk and Ward [68]. Note that phi units are defined as p h i = l o g 2 D ie log to the base 2 of the size D, which is in mm. Sorting is the standard deviation of particle size. The skewness characterizes the asymmetry of a distribution: skewness > 0 signifies a distribution with a tail towards larger sizes; skewness < 0 has a tail towards smaller sizes [69]. The kurtosis is a measure of grains around the median [70].
SampleMean Particle Size (µm/microns)Sorting (phi units)SkewnessKurtosis
P1125.61.39−0.383.08
P11B28.91.66−0.454.09
P1240.41.45−0.413.19
P1345.01.36−0.422.98
P145.642.03−0.107.53
P1533.81.41−0.363.01
YRD_BANK_P23.751.45−0.103.87
YRD_BANK_P2B35.71.51−0.463.39
P219.211.85−0.335.71
P2217.51.59−0.374.12
P2346.11.19−0.412.08
P2418.81.73−0.315.13
P255.032.08−0.067.69
P266.992.07−0.187.78
P2710.62.10−0.298.18
P31A66.31.33−0.492.39
P31B72.41.04−0.381.59
P326.041.92−0.236.64
P3324.11.80−0.415.25
P3432.31.43−0.433.02
P359.611.83−0.295.88
Table C2. Mineralogy of six samples: P11, P15, P24, P27, P31A, and P34.
Table C2. Mineralogy of six samples: P11, P15, P24, P27, P31A, and P34.
MineralQuartzCalciteK-feldspar and plagioclase feldsparBiotite micaClay minerals
Abundance60–65%5–8%5–10%2–5%20–25%
Table C3. Roundness of 384 particles from four samples: P11, P24, P31A, and P34. The roundness of a particle depends on the sharpness of the edges and corners rather than the shape [64].
Table C3. Roundness of 384 particles from four samples: P11, P24, P31A, and P34. The roundness of a particle depends on the sharpness of the edges and corners rather than the shape [64].
High sphericityLow sphericity
Very angular0137
Angular24145
Sub angular1460
Sub rounded13
Rounded00
Well rounded00

References and Notes

  1. Wang, J.; Gao, W.; Xu, S.; Yu, L. Evaluation of the combined risk of sea level rise, land subsidence, and storm surges on the coastal areas of Shanghai, China. Climatic Change 2012, 115, 1–22. [Google Scholar] [CrossRef]
  2. Syvitski, J.P.M.; Kettner, A.J.; Overeem, I.; Hutton, E.W.H.; Hannon, M.T.; Brakenridge, G.R.; Day, J.; Vorosmarty, C.; Saito, Y.; Giosan, L.; et al. Sinking deltas due to human activities. Nature Geosci. 2009, 2, 681–686. [Google Scholar] [CrossRef]
  3. Chaussard, E.; Amelung, F.; Abidin, H.; Hong, S.-H. Sinking cities in Indonesia: ALOS PALSAR detects rapid subsidence due to groundwater and gas extraction. Remote Sensing Environ. 2013, 128, 150–161. [Google Scholar] [CrossRef]
  4. Morton, R.; Bernier, J.; Barras, J. Evidence of regional subsidence and associated interior wetland loss induced by hydrocarbon production, gulf coast region, USA. Environ. Geol. 2006, 50, 261–274. [Google Scholar] [CrossRef]
  5. Tornqvist, T.E.; Wallace, D.J.; Storms, J.E.A.; Wallinga, J.; van Dam, R.L.; Blaauw, M.; Derksen, M.S.; Klerks, C.J.W.; Meijneken, C.; Snijders, E.M.A. Mississippi Delta subsidence primarily caused by compaction of holocene strata. Nature Geosci. 2008, 1, 173–176. [Google Scholar] [CrossRef]
  6. Boesch, D.F.; Josselyn, M.N.; Mehta, A.J.; Morris, J.T.; Nuttle, W.K.; Simenstad, C.A.; Swift, D.J.P. Scientific assessment of coastal wetland loss, restoration and management in Louisiana. J. Coastal Res. 1994, Special Issue 20, i–v, 1–104. [Google Scholar]
  7. Cui, B.-L.; Li, X.-Y. Coastline change of the yellow river estuary and its response to the sediment and runoff (1976–2005). Geomorphology 2011, 127, 32–40. [Google Scholar] [CrossRef]
  8. Jia, X.; Wang, H. Mineral compositions and sources of the riverbed sediment in the desert channel of yellow river. Environ. Monit. Assess. 2011, 173, 969–983. [Google Scholar] [CrossRef] [PubMed]
  9. Higgins, S.; Overeem, I.; Tanaka, A.; Syvitski, J.P.M. Land subsidence at aquaculture facilities in the Yellow River Delta, China. Geophys. Res. Lett. 2013, 40, 3898–3902. [Google Scholar] [CrossRef]
  10. Chaussard, E.; Wdowinski, S.; Cabral-Cano, E.; Amelung, F. Land subsidence in Central Mexico detected by ALOS InSAR time-series. Remote Sens. Environ. 2014, 140, 94–106. [Google Scholar] [CrossRef]
  11. Samsonov, S.; d'Oreye, N.; Smets, B. Ground deformation associated with post-mining activity at the French-German border revealed by novel insar time series method. Int. J. Appl. Earth Obs. Geoinf. 2013, 23, 142–154. [Google Scholar] [CrossRef]
  12. Yang, C.-S.; Zhang, Q.; Zhao, C.-Y.; Wang, Q.-L.; Ji, L.-Y. Monitoring land subsidence and fault deformation using the small baseline subset InSAR technique: A case study in the Datong Basin, China. J. Geodynam. 2014, 75, 34–40. [Google Scholar] [CrossRef]
  13. Plattner, C.; Wdowinski, S.; Dixon, T.H.; Biggs, J. Surface subsidence induced by the Crandall Canyon Mine (Utah) collapse: InSAR observations and elasto-plastic modelling. Geophys. J. Int. 2010, 183, 1089–1096. [Google Scholar] [CrossRef]
  14. Wang, H.; Wright, T.J.; Yu, Y.; Lin, H.; Jiang, L.; Li, C.; Qiu, G. InSAR reveals coastal subsidence in the Pearl River Delta, China. Geophys. J. Int. 2012, 191, 1119–1128. [Google Scholar] [CrossRef]
  15. Becker, R.H.; Sultan, M. Land subsidence in the Nile Delta: Inferences from radar interferometry. Holocene 2009, 19, 949–954. [Google Scholar] [CrossRef]
  16. Mazzotti, S.; Lambert, A.; Van der Kooij, M.; Mainville, A. Impact of anthropogenic subsidence on relative sea-level rise in the Fraser River Delta. Geology 2009, 37, 771–774. [Google Scholar] [CrossRef]
  17. Hooper, A. A multi-temporal InSAR method incorporating both persistent scatterer and small baseline approaches. Geophys. Res. Lett. 2008, 35. [Google Scholar] [CrossRef]
  18. Hooper, A.; Segall, P.; Zebker, H. Persistent scatterer interferometric synthetic aperture radar for crustal deformation analysis. J. Geophys. Res. 2007, 112. [Google Scholar] [CrossRef]
  19. Shi, C.; Zhang, D. Processes and mechanisms of dynamic channel adjustment to delta progradation: The case of the mouth channel of the Yellow River, China. Earth Surf. Process. Landf. 2003, 28, 609–624. [Google Scholar] [CrossRef]
  20. Chu, Z.X.; Sun, X.G.; Zhai, S.K.; Xu, K.H. Changing pattern of accretion/erosion of the modern Yellow River (Huanghe) Subaerial Delta, China: Based on remote sensing images. Mar. Geol. 2006, 227, 13–30. [Google Scholar] [CrossRef]
  21. Wang, H.; Bi, N.; Saito, Y.; Wang, Y.; Sun, X.; Zhang, J.; Yang, Z. Recent changes in sediment delivery by the Huanghe (Yellow River) to the sea: Causes and environmental implications in its estuary. J. Hydrol. 2010, 391, 302–313. [Google Scholar] [CrossRef]
  22. Fang, H.; Xu, J. Land cover and vegetation change in the Yellow River Delta nature reserve analyzed with Landsat Thematic Mapper data. Geocarto Int. 2000, 15. [Google Scholar] [CrossRef]
  23. Chen, K.; Yuan, J.; Yan, C. China Shandong Yellow River Delta national nature reserve. Available online: http://www.ramsar.org/pdf/lib/hbk4-07cs05.pdf (accessed on 30 September 2014).
  24. Farr, T.G.; Rosen, P.A.; Caro, E.; Crippen, R.; Duren, R.; Hensley, S.; Kobrick, M.; Paller, M.; Rodriguez, E.; Roth, L.; et al. The shuttle radar topography mission. Rev. Geophys. 2007, 45. [Google Scholar] [CrossRef]
  25. Massonnet, D.; Feigl, K.L. Radar interferometry and its application to changes in the earthʼs surface. Rev. Geophys. 1998, 36, 441–500. [Google Scholar] [CrossRef]
  26. Hanssen, R.F. Radar Interferometry, Data Interpretation and Error Analysis; Kluwer Academic Publishers: New York, NY, USA, 2001; Volume 2, p. 298. [Google Scholar]
  27. Ferretti, A.; Prati, C.; Rocca, F. Permanent scatterers in SAR interferometry. IEEE Trans. Geosci. Remote Sens. 2001, 39, 8–20. [Google Scholar] [CrossRef]
  28. Berardino, P.; Fornaro, G.; Lanari, R.; Sansosti, E. A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms. IEEE Trans. Geosci. Remote Sens. 2002, 40, 2375–2383. [Google Scholar] [CrossRef]
  29. Hooper, A.; Bekaert, D.; Spaans, K.; Arıkan, M. Recent advances in SAR interferometry time series analysis for measuring crustal deformation. Tectonophysics 2012, 514–517, 1–13. [Google Scholar] [CrossRef]
  30. Ferretti, A.; Fumagalli, A.; Novali, F.; Prati, C.; Rocca, F.; Rucci, A. A new algorithm for processing interferometric data-stacks: SQUEESAR. IEEE Trans. Geosci. Remote Sens. 2011, 49, 3460–3470. [Google Scholar] [CrossRef]
  31. Li, Z.; Fielding, E.J.; Cross, P. Integration of InSAR time-series analysis and water-vapor correction for mapping postseismic motion after the 2003 Bam (Iran) earthquake. IEEE Trans. Geosci. Remote Sens. 2009, 47, 3220–3230. [Google Scholar] [CrossRef]
  32. Zhang, L.; Ding, X.; Lu, Z. Modeling PSInSAR time series without phase unwrapping. IEEE Trans. Geosci. Remote Sens. 2011, 49, 547–556. [Google Scholar] [CrossRef]
  33. Zhang, L.; Ding, X.; Lu, Z.; Jung, H.-S.; Hu, J.; Feng, G. A novel multitemporal InSAR model for joint estimation of deformation rates and orbital errors. IEEE Trans. Geosci. Remote Sens. 2014, 52, 3529–3540. [Google Scholar] [CrossRef]
  34. Hooper, A.; Zebker, H. Phase unwrapping in three dimensions with application to InSAR time series. J. Opt. Soc. Am. A 2007, 24, 2737–2747. [Google Scholar] [CrossRef]
  35. Chen, C.W.; Zebker, H.A. Two-dimensional phase unwrapping with use of statistical models for cost functions in nonlinear optimization. J. Opt. Soc. Am. A 2001, 18, 338–351. [Google Scholar] [CrossRef]
  36. Schmidt, D.A.; Bürgmann, R. Time-dependent land uplift and subsidence in the Santa Clara Valley, California, from a large interferometric synthetic aperture radar data set. J. Geophys. Res. 2003, 108, 2416. [Google Scholar] [CrossRef]
  37. Ofeigsson, B.G.; Hooper, A.; Sigmundsson, F.; Sturkell, E.; Grapenthin, R. Deep magma storage at Hekla Volcano, Iceland, revealed by InSAR time series analysis. J. Geophys. Res. 2011, 116, B05401. [Google Scholar]
  38. Parsons, B.; Wright, T.; Rowe, P.; Andrews, J.; Jackson, J.; Walker, R.; Khatib, M.; Talebian, M.; Bergman, E.; Engdahl, E.R. The 1994 Sefidabeh (Eastern Iran) earthquakes revisited: New evidence from satellite radar interferometry and carbonate dating about the growth of an active fold above a blind thrust fault. Geophys. J. Int. 2006, 164, 202–217. [Google Scholar] [CrossRef]
  39. Lu, X.; Shu, Q.; Zeng, X.; Liu, T. Geology of Gudao Oilfield; (in Chinese). Petroleum Industry Press: Beijing, China, 2005. [Google Scholar]
  40. Zhang, S.; Wang, Y.; Shi, D.; Xu, H.; Pang, X.; Li, M. Fault-fracture mesh petroleum plays in the jiyang superdepression of the Bohai Bay Basin, Eastern China. Mar. Pet. Geol. 2004, 21, 651–668. [Google Scholar] [CrossRef]
  41. Davis, P.M. Surface deformation due to inflation of an arbitrarily oriented triaxial ellipsoidal cavity in an elastic half-space, with reference to Kilauea Volcano, Hawaii. J. Geophys. Res. 1986, 91, 7429–7438. [Google Scholar] [CrossRef]
  42. Biot, M.A. General theory of three dimensional consolidation. J. Appl. Phys. 1941, 12, 155–164. [Google Scholar] [CrossRef]
  43. Terzaghi, K. Theoretical Soil Mechanic; Wiley: New York, NY, USA, 1943. [Google Scholar]
  44. Nur, A.; Byerlee, J.D. An exact effective stress law for elastic deformation of rock with fluids. J. Geophys. Res. 1971, 76, 6414–6419. [Google Scholar] [CrossRef]
  45. Segall, P. Induced stresses due to fluid extraction from axisymmetric reservoirs. Pure Appl. Geophys. 1992, 139, 535–560. [Google Scholar] [CrossRef]
  46. Segall, P.; Grasso, J.-R.; Mossop, A. Poroelastic stressing and induced seismicity near the Lacq Gas field, Southwestern France. J. Geophys. Res. 1994, 99, 15423–15438. [Google Scholar] [CrossRef]
  47. Geertsma, J. Land subsidence above compacting oil and gas reservoirs. J. Pet. Technol. 1973, 25, 734–744. [Google Scholar] [CrossRef]
  48. Ren, H.-Y.; Zhang, X.-J.; Song, Z.-Y.; Rupert, W.; Gao, G.-J.; Guo, S.-X.; Zhao, L.-P. Comparison of microbial community compositions of injection and production well samples in a long-term water-flooded petroleum reservoir. PLoS ONE 2011, 6, e23258. [Google Scholar] [CrossRef] [PubMed]
  49. Sun, F.; Xue, S.; Ge, H.; Song, L. Numerical simulation of stress-induced wellbore damage in unconsolidated sandstone reservoirs. (In Chinese)Petrol. Drill. Tech. 2008, 36, 47–50. [Google Scholar]
  50. Kirkpatrick, S. Optimization by simulated annealing: Quantitative studies. J. Stat. Phys. 1984, 34, 975–986. [Google Scholar] [CrossRef]
  51. Biggs, J.; Anthony, E.Y.; Ebinger, C.J. Multiple inflation and deflation events at Kenyan Volcanoes, East African Rift. Geology 2009, 37, 979–982. [Google Scholar] [CrossRef]
  52. Muntendam-Bos, A.; Kroon, I.; Fokker, P. Time-dependent inversion of surface subsidence due to dynamic reservoir compaction. Math. Geosci. 2008, 40, 159–177. [Google Scholar] [CrossRef]
  53. Fialko, Y.; Simons, M. Deformation and seismicity in the Coso geothermal area, Inyo County, California: Observations and modeling using satellite radar interferometry. J. Geophys. Res. 2000, 105, 21781–21793. [Google Scholar] [CrossRef]
  54. Lewis, R.W.; Sukirman, Y. Finite element modelling for simulating the surface subsidence above a compacting hydrocarbon reservoir. Int. J. Numer. Anal. Methods Geomech. 1994, 18, 619–639. [Google Scholar] [CrossRef]
  55. Lewis, R.W.; Sukirman, Y. Finite element modelling of three-phase flow in deforming saturated oil reservoirs. Int. J. Numer. Anal. Methods Geomech. 1993, 17, 577–598. [Google Scholar] [CrossRef]
  56. Kosloff, D.; Scott, R.F.; Scranton, J. Finite element simulation of Wilmington oil field subsidence: II. Nonlinear modelling. Tectonophysics 1980, 70, 159–183. [Google Scholar] [CrossRef]
  57. Kosloff, D.; Scott, R.F.; Scranton, J. Finite element simulation of Wilmington oil field subsidence: I. Linear modelling. Tectonophysics 1980, 65, 339–368. [Google Scholar] [CrossRef]
  58. Capasso, G.; Mantica, S. Numerical simulation of compaction and subsidence using abaqus. In Proceedings of 2006 ABAQUS Usersʼ Conference, Cambridge, MA, USA, 23 May 2006.
  59. Wessel, P.; Smith, W.H.F. New, improved version of generic mapping tools released. Eos Trans. AGU 1998, 79, 579. [Google Scholar] [CrossRef]
  60. Blott, S.J.; Croft, D.J.; Pye, K.; Saye, S.E.; Wilson, H.E. Particle size analysis by laser diffraction. Geological Society, London, Special Publications 2004, 232, 63–73. [Google Scholar] [CrossRef]
  61. Hackley, V.A.; Lum, L.-S.; Gintautas, V.; Ferraris, C.F. Particle Size Analysis by Laser Diffraction Spectrometry: Application to Cementitious Powders; NISTIR 7097; NIST: Gaithersburg, MD, USA, 2004. [Google Scholar]
  62. Ferreira, T.; Rasband, W. Imagej; US National Institutes of Health: Bethesda, MD, USA, 2006. [Google Scholar]
  63. Nemecz, E.; Pécsi, M.; Hartyáni, Z.; Horváth, T. The origin of the silt size quartz grains and minerals in loess. Quaternary Int. 2000, 68–71, 199–208. [Google Scholar] [CrossRef]
  64. Powers, M.C. A new roundness scale for sedimentary particles. J. Sediment. Res. 1953, 23, 117–119. [Google Scholar] [CrossRef]
  65. Russell, R.D. Effects of transportation on sedimentary particles; Part 1. Transportation. In Recent Marine Sediments; Trask, P., Ed.; SEPM (Society for Sedimentary Geology): Tulsa, OK, USA, 1955; Volume 4, pp. 32–47. [Google Scholar]
  66. Wang, S.J.; Hassan, M.A.; Xie, X.P. Relationship between suspended sediment load, channel geometry and land area increment in the Yellow River Delta. CATENA 2006, 65, 302–314. [Google Scholar] [CrossRef]
  67. Qiao, S.; Shi, X.; Zhu, A.; Liu, Y.; Bi, N.; Fang, X.; Yang, G. Distribution and transport of suspended sediments off the Yellow River (Huanghe) mouth and the nearby Bohai Sea. Estuar. Coast. Shelf Sci. 2010, 86, 337–344. [Google Scholar] [CrossRef]
  68. Folk, R.L.; Ward, W.C. Brazos river bar [texas]; A study in the significance of grain size parameters. J. Sediment. Res. 1957, 27, 3–26. [Google Scholar] [CrossRef]
  69. Press, W.H.; Teukolsky, S.A.; Vetterling, W.T.; Flannery, B.P. Numerical Recipes in Fortran 77: The Art of Scientific Computing; Cambridge University Press: Cambridge, UK, 1997; Volume 1. [Google Scholar]
  70. Blott, S.J.; Pye, K. Gradistat: A grain size distribution and statistics package for the analysis of unconsolidated sediments. Earth Surf. Process. Landf. 2001, 26, 1237–1248. [Google Scholar] [CrossRef]

Share and Cite

MDPI and ACS Style

Liu, P.; Li, Q.; Li, Z.; Hoey, T.; Liu, Y.; Wang, C. Land Subsidence over Oilfields in the Yellow River Delta. Remote Sens. 2015, 7, 1540-1564. https://0-doi-org.brum.beds.ac.uk/10.3390/rs70201540

AMA Style

Liu P, Li Q, Li Z, Hoey T, Liu Y, Wang C. Land Subsidence over Oilfields in the Yellow River Delta. Remote Sensing. 2015; 7(2):1540-1564. https://0-doi-org.brum.beds.ac.uk/10.3390/rs70201540

Chicago/Turabian Style

Liu, Peng, Qingquan Li, Zhenhong Li, Trevor Hoey, Yanxiong Liu, and Chisheng Wang. 2015. "Land Subsidence over Oilfields in the Yellow River Delta" Remote Sensing 7, no. 2: 1540-1564. https://0-doi-org.brum.beds.ac.uk/10.3390/rs70201540

Article Metrics

Back to TopTop