Combination with precise leveling and PSInSAR observations to quantify pumping-induced land subsidence in central Taiwan

Choushui River Fluvial Plain (CRFP) is located in the western central Taiwan, where the geomaterials are composed of alluvial deposits. Because the CRFP area receives highly variable rainfall in wet and dry seasons, the groundwater becomes the main resource of residential water. The precise leveling monitoring from 1970s indicated that the coastal areas of CRFP had been threatened by serious pumping-induced land subsidence. On the basis of relatively accurate measurements of precise leveling measurements, we used cokriging technique to incorporate a number of InSAR images to quantify the surface deformation in CRFP. More specifically, the well-developed Persistent Scatterer InSAR (PSI) was employed to process 34 Envisat images (2005–2008) and the results of PSI was then used for improving the spatial resolution of data from precise leveling. The results of cokriging estimation indicate whether the rate or the area of the land subsidence slows down gradually from 2005 to 2008. The subsidence in the northern part of CRFP was influenced by the groundwater decline in aquifer III, and the southern part was influenced by groundwater decline in aquifer II and III. The cokriging estimation was also comparable with continuous GPS data, and their correlation coefficient is 0.9603 and the root mean square is 10.56 mm yr.


Introduction
Groundwater resource management is a notable issue due to growing population and rapid urban development, especially in areas where surface water is limited or shortage.Many cities, such as Shanghai, Las Vegas, and Bangkok had faced serious land subsidence because of unlimited groundwater withdrawal.To identify the land subsidence range effectively the geodetic techniques with different scales record the behaviors of surface deformation.Precise leveling, Global Positioning System (GPS) and Interferometric Synthetic Aperture Radar (InSAR) have been the most powerful techniques in geodesy (Galloway et al., 1998;Massonnet and Feigl, 1998;Galloway and Burbey, 2011;Hung et al., 2010).
In the central Taiwan, Choushui River fluvial plain (CRFP) faced the serious problem of land subsidence by poignant groundwater drawdown since 1970's.Central Geological Survey (CGS) has investigated several hydrologic and geologic studies.Thereafter, Water Resources Agency (WRA) have been collecting geodetic data such as precise leveling, GPS and compaction monitoring well (CMW) to monitoring the rate of land subsidence in CRFP (Central Geological Survey, 1999;Liu et al., 2001, Liu et al., 2004;Hwang et al., 2008).Hung et al. (2011) apply the Persistent Scatter InSAR (PSI), which was developed by Hooper et al. (2004, Hooper (2008), for identifying PS pixels and extracting more than 100 pixels km −2 in CRFP from 2006 to 2008.In this paper, the main aim is to integrate precise leveling data and PSI result with geostatistical method, variogram and cokriging, for retaining the accuracy of precise leveling and the spatial resolution of PSI and estimating the surface deformation, then comparing cokriging estimation with the difference of yearly groundwater level, at last verifying the cokriging estimation with the observed GPS data from 2005 to 2008 in CRFP.

Geological setting
Taiwan is located at the boundary between the Philippine Sea Plate (PSP) and Eurasian Plate (EUP), where the still ongoing collision started about 5-7 Ma.The PSP move toward the northwest to EUP at a rate of 82 mm yr −1 .Choushui River Fluvial Plain (CRFP) is one part of the coastal plain in the western Taiwan.This area covers about 2000 km 2 and the boundary is enclosed by Pakua tableland and Taolao hill to the east, Taiwan Strait to the west, Wu River to the north and Putzu River to the south.Choushui River is the longest river in Taiwan and flows cross Central Range, Western Foothill and Coastal Plain then in to Taiwan Strait.The range of the sedimentary thickness is from 750 to 3000 m and the grain size become finer from the east to the west.The materials of the sedimentary contain clay, fine sand, coarse sand, and gravel.According the drilling core in the depth of 300 m, the subsurface hydrogeology can be divided into eight overlapping sequences.There are four marine sequences (aquitard I-IV) and four non-marine sequences (aquifer I-IV).The right side of Fig. 1 shows four hydrogeological profiles in CRFP (Central Geological Survey, 1999;Liu et al., 2001;Liu et al., 2004).

Precise leveling
The land subsiding surveys determined by precise leveling in CRFP has recorded since 1975.These data acquired from Water Resource Agency (WRA) are based on the first or second order standard of leveling procedures proposed by Ministry of the Interior (MOI).The specifications demand that any loop misclosure should be below 3 mm √ K, where K is the distance between two neighboring benchmarks in km.In this study 298 precise leveling benchmarks were be used from 2005 to 2008 (Fig. 1).(Hooper et al., 2007;Hooper, 2008).The topographic effect was removed by the DEM from NASA's SRTM mission and the precise orbit parameters from Delft Institute for Earth-Oriented Space Research were used to correct the orbit errors.

Groundwater level
Since 1992 WRA installed hydrologic monitoring well evenly, and dug the depths ranging from 200 to 300 m in CRFP.According to the hydrogeology profiles (Fig. 1) CRFP can be divided into four aquifers in depth of 300 m.In this study we collect 98 monitoring well screens in the aquifer I, 97 screens in the aquifer II, 44 screens in the aquifer III, and 21 screens in the aquifer IV.The yearly average of groundwa-ter level from 2005 to 2008 was interpolated by using kriging method.

Methodology
The grid of 9867 active cells is set up in the study area and the size of each cell is 500 m × 500 m.Two geostatistical methods, variogram and cokriging, are adopted for analyzing and combining the precise leveling and PSI result.The FOR-TRAN codes gamv and cokb3d in GSLIB (Deutsch and Journel, 1998) make use of generating experimental variograms and cokriging respectively.

Variogram analysis
An experimental semivariogram, γ (h), is defined as half of the average squared difference between two values of an uniproc-iahs.net/372/77/2015/Proc.IAHS, 372, 77-82, 2015 variate approximately separated by vector h and can be written as: (1) where N(h) is the number of pairs; z(u i ) and z(u i +h) are the value of precise leveling or PSI at the location u i and u i + h of the pair i; and the vector h can be specified with particular directions and lag distance.

Cokriging interpolation
The cokriging estimate is a linear combination of both primary and secondary data value and minimizes the variance of the estimation error by exploiting the cross-correlation between two variables (Isaaks and Srivastava, 1989).The precise leveling is set to the primary variable, L(u a 1 ), the PSI results is sets to the secondary variable, P (u a 2 ), and the cokring estimation, Z COK (u), can be written as: where λ a 1 are the weights applied to the n 1 precise leveling data, u a 1 are the n 1 locations of precise leveling benchmarks, λ a 2 are the weights applied to the n 2 PSI results, u a 2 are the n 2 locations of PSI pixel.

Cokriging estimation and groundwater level
The cokriging estimations in the three intervals display on first column of Fig.In general the subsidence in the northern part of CRFP was influenced by groundwater decline in aquifer III, and the southern part was influenced by groundwater decline in aquifer II and III.In 2005-2006 the land subsidence is most serious however the groundwater decline is slightest.It reveals that the elastic compaction would not be a major factor to induce the land subsidence during this interval.In 2005-2007 and 2005-2008, the regional groundwater decline corresponds the same location of land subsidence.Especially in aquifer III the depth is about 175-300 m, and the regional decline of groundwater level distributes in CRFP obviously.It may reflect the pumping-water location or behavior in aquifer III.The geo-material near eastern CRFP contains the gravel higher than 90 % and the permeability of the gravel is between 1.0 × 10 −4 and 9.9 × 10 −4 (Central Geological Survey, 1999).It gave a good path for groundwater flow and supply.Therefor even the groundwater level near Pakua Tableland draws down almost 15 m, the surface deformation reveals slight.
The five continuous CPS stations in CRFP were acquired from WRA.Two stations, Hsi-Kang (HK) and Hsin-Hsing (HH) recorded the vertical displacement from 2005 to 2008 and the others stations, Tu-Ku (TK), Lin-Nei (LN) and Ko-Tso (KT) recorded the vertical displacement from 2006 to 2008. Figure 4 shows the direct comparisons of cokriking estimation with GPS data.The correlation coefficient between the cokriging estimation and GPS data within the same cell is 0.9603 and the RMS is 10.56 mm yr −1 .Thus the cokriging estimation is a reliable data for detecting the land subsidence in this study.
To sum above results, the mechanism of land subsidence in CRFP is quite complication.The models not only include elastic and inelastic compaction behaviors but also should consider the characteristics of different aquifers.

Figure 1 .
Figure 1.The study area.The black circle is the leveling benchmarks.The red triangle is continuous GPS station.The black lines indicate the location of hydrogeological profile.The right side of figure shows the hydrogeological profiles.The hydrogeological profiles were modified from Central Geological Survey (1999).

Figure 2 .
Figure 2. The left side of the figure is the vertical deformation of precise leveling and PSI between 2005 and 2008.The cold color means the land subsidence.The right side of the figure from top to bottom are precise leveling variogram, PSI variogram and cross-variogram between precise leveling and PSI.Every variogram has 6 directions and a Gaussian model.

Figure 3 .
Figure3.First column is the cokriging estimation in the three intervals(2005-2006, 2005-2007 and 2005-2008).The rest columns are the difference of yearly average of groundwater level in four aquifers.
from Track 232, Frame 3123 and 17 Envisat interferogram form Track 232, Frame 3141 were used to estimate the displacement of the radar line-of-sight (LOS).At the same period 298 precise leveling data recorded the vertical deformation in CRFP.The left side of Fig. 2 shows the vertical displacement rate of precise leveling and PSI from 2005 to 2008.The cold color represents the land subsidence relative to the reference point and the warm color mean surface uplift.The most subsidence rate of PSI result almost reaches 60 mm yr −1 and the most subsidence rate of precise leveling data reaches 80 mm yr −1 in the center of southern CRFP.The lag distance of 1 km was used to analyze the semivariograms and cross semiveriogram in 6 horizontal directions (right side of Fig.2).The red triangle is north direction, the green square is N-30 • E direction, the blue star is N-60 • E, the yellow diamond is N-90 • E, the purple circle is N-60 • W, the gray hexagon is N-30 • W, and the black solid line is the omnidirectional semivariogram fitted by Gaussian model.

Figure 4 .
Figure 4.The comparison of cokriging estimation with GPS data.The code name of GPS station is written in text.
Due to the area of CRFP covers two SAR image modes (Track 232, Frame 3123 and Track 232, Frame 3141), 34 Envisat radar images (16 images from Track 232, Frame 3123; 18 images from Track 232, Frame: 3141) from European 3. The negative values of the rate (cold