Land Subsidence over Oilfields in the Yellow River Delta

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, OPEN ACCESS Remote Sens. 2015, 7 1541 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.


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.

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 km 2 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].[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.

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.

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.

InSAR Time Series Results
The processed ASAR images cover an area of about 2200 km 2 .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.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 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.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).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).

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 nonuniform displacement signals observed from InSAR.

Modeling Subsidence in Gudao Oilfield
The Gudao oilfield, located in the Jiyang depression in the Bohaiwan Basin has reserves are 4 × 10 8 tons of oil and 47 × 10 8 m 3 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.
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.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).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.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.(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".

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

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

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.

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.8MPa 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.
where 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, is in the range (15.0, 22.9) for Track 132, Δθ is ~4° for Track 132 and 404, while is in the range (−180°, 180°) for any possible deformations.Hence the range for ∆ 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.

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

Figure 1 .
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 2 .
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 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.

Figure 4 .
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 5 .
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 6 .
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 7 .
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 8 .
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 9 .
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 B1 .
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 C1 .
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].

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 2 .
Best Fit Poroelastic disk reservoir model parameters and uncertainty.A: offset between InSAR and model estimated.B: mean value as InSAR reference.