Land Subsidence in Chiayi, Taiwan, from Compaction Well, Leveling and ALOS/PALSAR: Aquaculture-Induced Relative Sea Level Rise

Abstract

Chiayi County is located in the largest alluvial plain of Taiwan with extensive aquaculture and rice farming sustained by water extracted from groundwater wells. Chiayi is a typical aquaculture area affected by land subsidence, yet such lands worldwide combine to provide nearly 90% of global aquaculture products, greatly reducing oceanic overfishing problems. This study uses precision leveling, multi-layer compaction monitoring well (MLCW) and spaceborne SAR interferometry (InSAR) to examine the cause and effect of land subsidence in Chiayi associated with groundwater extractions and changes. Heights at benchmarks in a leveling network are measured annually and soil compactions at 24–26 layers up to 300-m depths at 7 MLCWs are collected at one-month intervals. Over 2007–2011, 15 ALOS/PALSAR images are processed by the method of TCPInSAR to produce subsidence rates. All sensors show that land subsidence occur in most parts of Chiayi, with rates reaching 4.5 cm/year around its coast, a result of groundwater pumping from shallow to deep aquifers. MLCWs detect mm-accuracy seasonal soil compactions coinciding with groundwater level fluctuations and causing dynamic compactions. Compactions near Taiwan High Speed Rail may reduce the strength of the rail’s supporting columns to degrade its safety. The SAR images yield subsidence rates consistent with those from leveling and compaction wells after corrections for systematic errors by the leveling result. Subsidence in Chiayi’s coastal area leads to relative sea level rises at rates up to 15 times larger than the global eustatic sea level rising rate, a risk typical for world’s aquaculture-rich regions. At the fish pond-covered Budai Township, InSAR identifies subsidence spots not detected by leveling, providing crucial geo-information for a sustainable land management for aquaculture industry.

1. Introduction

Chiayi is a county located in the central coastal area of Taiwan. It is a part of the Chianan Plain, which is the largest plain in Taiwan. The Chianan Plain is the “rice warehouse” of Taiwan that supplies over 50% of rice, vegetables, fruits, flowers and aquaculture products to a population of 23.5 million [1]. As shown in Figure 1a, Chiayi is bordered with the Beigang River and the Choushui River alluvial fan to the north, the Beigang River and Tainan City to the south, Douliu hill to the east and the Taiwan Strait to the west. The Tropic of Cancer passes through Chiayi, and the highest mountain closest to the Tropic of Cancer, Mt. Jade (Yushan), is situated on Chiayi’s eastern border. Because of the increasingly economic development in Taiwan, water resources become stressful. In particular, for an agricultural county such as Chiayi, surface water can no longer meet the demand of water. As a result, heavy groundwater exploitations occur in most part of Chiayi. Figure 1b shows the densely distributed groundwater wells in Chiayi drilled to different depths. In the coastal area of Chiayi, the depths of groundwater wells can reach more than 250 m. Wells are also excavated along the Taiwan High Speed Rail (THSR), posing potential risks for its operation [2].

Chiayi is also a county with a large aquaculture industry. Figure 1a shows the distribution of fish farms in its coastal area. The groundwater wells in Chiayi’s coastal area (Figure 1b) are largely for supplying water to the fish farms seen in Figure 1a. In general, the depths of wells in this aquaculture zone are deeper than those in Chiayi’s inland area, where groundwater-based rice plantations are common. Dongshih, Budai and Yijhu are three major coastal townships with high aquaculture industries. Because the groundwater here has a constant temperature of about 20 °C, in cold weathers groundwater is used to adjust the temperature of fish ponds to avoid deaths of fish. As such, groundwater is not only a water resource but also a thermo-regulator for the aquaculture industry in Chiayi. According to the Ministry of Economy and the Water Resource Agency of Taiwan, from 1980 onward, a large portion of Taiwan’s coastal area has been transformed into farms for aquaculture, consuming a large amount of groundwater and in turn triggering coastal land subsidence. These aquaculture fisheries need clean water to grow fish, and this need is met by exploring groundwater around fish farms. As ocean fish sources become increasingly scarce, the need of aquaculture products has been growing to exacerbate groundwater extractions.

Since the completion of THSR in 2006, Chiayi has established several new industrial zones. Some large-scaled corporations in the industrial zones extract deep groundwater for multiple uses. The combined exploitations of shallow groundwater by farming and deep groundwater by industries make it increasingly difficult to track the overall use of groundwater in Chiayi. In 2015, a 6-month drought in Taiwan dramatically increased the need of groundwater, leading to a heavy groundwater pumping in Chiayi and worsened land subsidence here. Farming, industrialization and climate change have combined to prompt an urgent need for an effective and precise method to manage groundwater in Chiayi and to monitor land subsidence caused by groundwater extractions.

Land subsidence in Chiayi has long been monitored by geodetic methods largely based on precision leveling and GPS [3]. Geodetic sensors in Chiayi provide pointwise measurements of subsidence and the point density largely depends on the budget and the accessibility of monitoring locations. The current point density of geodetic sensors in Chiayi is not sufficient, requiring an augmented sensor to remove the deficiency in data density. An augmenting sensor is spaceborne SAR interferometry (InSAR), which has been significantly advanced over the past few decades due to the introductions of new processing methods and new satellite missions.

There are a number of methods that can be used to determine subsidence rates from SAR images, such as PSI [4,5,6], SBAS [7,8] and TCPInSAR [9,10]. In addition, there are many satellites, such as ERS, Envisat and ALOS/PALSAR, that supply SAR images for subsidence studies over alluvial plains, coastal zones and other surfaces [11,12,13,14]. In this paper, we will use the method of TCPInSAR [9] to determine land subsidence rates in Chiayi. TCPInSAR utilizes a least-squares-based method for multitemporal synthetic aperture radar interferometry that allows one to estimate deformations without the need of phase unwrapping. TCPInSAR has been successfully applied for mapping deformation associated with ground instability [10] and landslide movement [15].

This study will show how SAR images from Japan’s ALOS/PALSAR mission are used to determine land subsidence in Chiayi, and how the InSAR result is assessed by results from precision leveling, compaction monitoring well and GPS. Systematic errors in the ALOS/PALSAR result will be analyzed and reduced by fusion with precision leveling measurements, particularly over fish ponds and areas with large horizontal velocities. Chiayi is typical aquaculture area in the world to show how groundwater-sustained aquaculture can cause land subsidence and relative sea level rise [16,17,18], which in turn damages the aquaculture industry and affects the nearly 90% of global aquaculture products critical to sustaining world protein needs.

2. Geological Settings and Data

2.1. Geological Settings

Chiayi, as well as the Chianan Plain, has complicated stratum formations, consisting of multiple layers of gravel and sand that can lead to groundwater-induced soil compactions in both the coastal and inland areas of Chiayi. The coastal plain of Chiayi is a sedimentary wedge composed of unconsolidated gravel, sand, silt, and clay, along with variable amounts of shells. Figure 2a shows two profiles, along which geological drilling at key stations provide hydrogeological information around Chiayi. Figure 2b shows the hydrogeological settings from Dongshih Township to Kumho Township along the coast. Figure 2c shows the hydrogeological settings along a west-east profile from Dongshih Township to Kanjiao Township. These two profiles show clear distributions of unconsolidated deposits in Chiayi. Over the depth range of 0 to 250 m, there are four major layers of formations. The depth of the first layer is over the 0–60 m depths. Here the aquifer is primarily distributed near the Hsinwen and Kumho stations in the southern coastal area. The second layer is at depths between 60 m and 120 m. Here the aquifers are relatively thick around Dongshi and Cuogou in the north, gradually becoming thinner toward the south. The third layer is distributed between 120 m and 210 m with thicker aquifers around Dongshi and Anhe, much like the second layer. The fourth layer is below 210 m and here the soil is primarily composed of compressible, sludge blankets and silt layers. According to Hung and Wang (2016) [3], exploiting groundwater in these four layers in Chiayi has resulted in land subsidence.

2.2. Data from Precision Leveling and Compaction Wells

Figure 3 shows the distributions of the leveling network, continuous GPS stations and monitoring wells around Chiayi for monitoring land subsidence. The leveling network is bordered with the south bank of the Choushui River to the north, the north bank of the Beigang River to the south, Chiayi City to the east and the coastal region to the west. In the network, 63 primary leveling lines are connected to form 17 closed loops. The total mileage of the leveling routes is about 310 km. On average, two neighboring leveling benchmarks are 1.5 km apart and the leveling network covers an area of 1087 km². Over 2007–2011, we annually measured the heights of benchmarks in Figure 3. The heights were network-adjusted using the method described in Hwang and Hung (2008) [2]. Using such time-lapsed heights, we determined the rates of vertical displacements at the benchmarks.

To understand the mechanism of Chiayi’s land subsidence, magnetic-ring layered monitoring wells were used to identify the compressions in different layers below the surface, by devices called multi-layer compaction monitoring wells (MLCW). Figure 4 shows a sample picture of this device and how data are collected using the device. For each MLCW, at least 24 magnetic rings were anchored in different layers of specified depths through a borehole according to the soil compositions. The depth of each magnetic ring was measured using a magnetic probe and an invar tape at a 1-mm precision level. As shown in Figure 2, the soil composition of the stratum around Chiayi is very complex and the hydrogeological characteristics of the four layers are different. Compaction monitoring wells provide data for understanding the groundwater extraction depths.

There are 7 monitoring wells with a depth of 300 m installed in the Chiayi area for land subsidence monitoring (see Figure 3 for the locations of the MLCWs). In particular, at wells DSES, WLES, and BDES (Figure 3), the records exceed 10 years, allowing for examining the long-term, time-dependent aquifer-system compactions around the locations of these wells.

Table 1. Information about the four compaction monitoring wells in Chiayi.

Description	DSES	WLES	BDES	AHES
Time period	March 1997–August 2011	April 1998–August 2011	March 1997–August 2011	November 2004–August 2011
Total compaction (cm)	36.5	56.7	46.9	7.9
Number of magnetic rings	24	24	24	26
Depth range of major compaction (m)	47–282	0–300	10–275	0–185

summarizes the measuring times, depths of the primary compression layers and total compressions at DSES, WLES, BDES and AHES. The well AHES is about 1 km to the rail of THSR and has recorded data since 2004. The numbers of deployed magnetic rings at the 7 wells range from 24 to 26. Note that, at the writing of this paper (May 2017), the 7 wells shown in Figure 3 continue to generate data, but we use such data at the time span coinciding with those of the precision leveling and ALOS/PALSAR images (see Section 2.3), i.e., 2007–2011. The depth of each magnetic ring is measured every month. The depth variation between two successive magnetic rings can be used to determine the compression of the layer between the two rings. As a result, the compactions at individual layers can be determined to assess the contributions of compactions at these layers.

2.3. SAR Images from ALOS over Chiayi

To augment our ground-based land subsidence results, we obtained 15 SAR images from the ALOS/PALSAR mission. ALOS/PALSAR is a mission launched by the Japan Aerospace Exploration Agency (JAXA) and is equipped with an L-band SAR capable of measuring time-dependent surface elevations.

Table 2. Information about 15 ALOS/PALSAR images over Chiayi.

Orbit\Track\Frame	Ascending\447\470
Spatial Perpendicular Baselines boundary (m)	1100
Temporal baselines upper boundary (day)	900
Time Span	4 March 2007–15 March 2011
No. of interferogram pairs	29
TCPs	687,736
Area (Km2)	771
TCP point density (pixel/Km2)	890

summarizes the baselines of interferometric image pairs generated from the 15 images, which are along ALOS/PALSAR’s path 447 and frame 450 (ascending mode) over Taiwan and spanned from March, 2007 to March 2011. Each image covers a width of 70 km in the side-look direction and has a spatial resolution of approximately 30 m. The radar incidence angle is about 38.7°. To reduce the topographic phase residuals, we selected a total of 29 interferograms with the spatial and temporal baseline thresholds of 1100 m and 900 days.

3. Method for Processing ALOS/PALSAR Images for Subsidence Determination

3.1. Processing ALOS/PALSAR Images by TCPInSAR

This paper uses the method of TCPInSAR [9] to process the 15 ALOS/PALSAR images (Section 2.3) to estimate subsidence rates in Chiayi. Conventionally, high-quality interferograms are selected according to spatial, temporal and Doppler baselines. However, baseline might not be a reliable indicator of interferometric quality, especially in areas hit by droughts and wet seasons. Over such areas, interferograms may have high coherence between two successive drought seasons, even though their temporal baselines are longer than those between drought and wet seasons. Therefore, in this paper we adopt the interferometric coherence to select the basic observations of TCPInSAR, instead of baseline threshold. Because of the large volume of the interferometric pairs generated from 15 SAR images, the computation burden can be extremely high if we calculate a full resolution map for each interferogram. Therefore, we use an efficient procedure to select the needed interferograms as follows. First, we use the intensity-based quad tree method to select representative pixels. Then, we calculate the coherence at these pixels, which is the representing coherence that indicates the quality of the whole interferograms. Finally, interferograms are ranked and selected according to the representing coherence. In total, we selected 29 interferograms to balance the phase noise and pixel spatial density. The numerical procedure is summarized as follows:

(a) Selecting TCP candidate points

TCPInSAR uses the intensity image calculated by DInSAR to determine the standard deviation of the image offset for each pixel of each set of primary and secondary image pairs in the along-track direction and the perpendicular direction of a moving window [19]. Initially, a window containing 256 × 256 pixels was used to calculate the offset for each pixel. The size of the window was then changed from 4 × 4 to 64 × 64. Using iterations, we assigned an offset sequence to each candidate point. The following standard deviation of is calculated in this sequence:

$${\mathsf{\sigma }}_{\mathrm{x},\mathrm{i}}=\sqrt{\frac{3}{2\mathrm{N}}}·\frac{\sqrt{1-{\mathsf{\gamma }}^2}}{\mathsf{\pi }\mathsf{\gamma }}{\mathsf{\chi }}^{\frac{3}{2}}$$

(1)

where ${\mathsf{\sigma }}_{\mathrm{x},\mathrm{i}}$ is the standard deviation at the ith pixel in the xth image, N is the number of pixels in the moving window, χ is a oversampling parameter, which is a constant at the ith pixel, and γ is the cross-correlation between the image and window. If the standard deviation ${\mathsf{\sigma }}_{\mathrm{x},\mathrm{i}}$ is less than a pre-defined value (e.g., 0.1 in a pixel), the pixel is regarded as a potential TCP candidate point.

After the offset matrix O(l × m) corresponding to the image is formed, we obtained the offset peak value ${\mathrm{O}}_{\mathrm{c}}$ of the image by

$${\mathrm{O}}_{\mathrm{c}}=\mathrm{peak}({\mathrm{hist}}_{2\mathrm{D}}({\mathrm{O}}_{\mathrm{l}\times \mathrm{m}}))$$

(2)

where ${\mathrm{hist}}_{2\mathrm{D}}({\mathrm{O}}_{\mathrm{l}\times \mathrm{m}})$ is the 2D histogram statistics. Finally, for a given threshold value, A, if

$$|{\mathrm{O}}_{\mathrm{i},\mathrm{j}}-{\mathrm{O}}_{\mathrm{c}}|<\mathrm{A}$$

(3)

then pixel (i, j) is a TCP candidate.

(b) Finalizing TCP points

After the previous step, we then fitted a polynomial to the offset sequences. If the offset of a pixel deviates from the modeled value of the polynomial, the pixel is rejected. After the TCP candidate points of all the secondary images were calculated, their intersection is determined and the TCP candidate point is selected as the final TCP only if it appears in every secondary image.

(c) Determining range change rates along LOS

The differential phase at point p in the ith interferogram can be expressed as:

$${\mathsf{\varphi }}_{\mathrm{p}}^{\mathrm{i}}=\mathrm{W}\left\{{\mathsf{\varphi }}_{\mathrm{topo},\mathrm{p}}^{\mathrm{i}}+{\mathsf{\varphi }}_{\mathrm{def},\mathrm{p}}^{\mathrm{i}}+{\mathsf{\varphi }}_{\mathrm{atm},\mathrm{p}}^{\mathrm{i}}+{\mathsf{\varphi }}_{\mathrm{orb},\mathrm{p}}^{\mathrm{i}}+{\mathsf{\varphi }}_{\mathrm{N},\mathrm{p}}^{\mathrm{i}}\right\}$$

(4)

where $\mathrm{W}\{.\}$ is the SAR-derived phase, ${\mathsf{\varphi }}_{\mathrm{topo},\mathrm{p}}^{\mathrm{i}}$ is the contribution caused by topography, ${\mathsf{\varphi }}_{\mathrm{def},\mathrm{p}}^{\mathrm{i}}$ is the contribution by deformation at the pixel, ${\mathsf{\varphi }}_{\mathrm{atm},\mathrm{p}}^{\mathrm{i}}$ is the contribution due to atmospheric effect, ${\mathsf{\varphi }}_{\mathrm{orb},\mathrm{p}}^{\mathrm{i}}$ is the contribution due to orbit effect, and ${\mathsf{\varphi }}_{\mathrm{N},\mathrm{p}}^{\mathrm{i}}$ is the noise term.

The selected TCP points were used to form the Delaunay triangles, where wrapped phases were retrieved from the ALOS image pairs to form observation equations. Such equations were then used to form the normal equations using the method of least-squares residuals (LSQR) [20]. The normal matrix is a sparse matrix that can be inverted by a method specifically designed for such a matrix. The LSQR method also detected and removed outliers in the observations and iterated parameter estimations. The LSQR method simultaneously solved for the ambiguities in the wrapped phases and for the range change rates along the line-of-sight (LOS).

3.2. LOS-Projected Rates in the Up Direction

In order to compare the InSAR result with the leveling result, the LOS rates were then converted to vertical rates as follows. First, the relation among various velocity component is [11]:

$${\mathrm{V}}_{\mathrm{LOS}}=\sin (\mathsf{\theta })\cos (\mathsf{\phi })\mathrm{E}+\sin (\mathsf{\theta })\sin (\mathsf{\phi })\mathrm{N}+\cos (\mathsf{\theta })\mathrm{UP}$$

(5)

where ${\mathrm{V}}_{\mathrm{LOS}}$ is the deformation velocity in the LOS direction, $\mathrm{E},\mathrm{N}\mathrm{and}\mathrm{U}$ are the velocity components in the east, north and up directions, $\mathsf{\theta }$ is the incidence angle and $\mathsf{\phi }$ is the along-track azimuth of the ALOS satellite. In this paper, we assumed that the horizontal displacement rates in Chiayi are zero (E = N = 0 in Equation (5)), and then projected the InSAR-derived LOS rates along the up direction using

$${\mathrm{U}}_{\mathrm{INSAR}}=\frac{{\mathrm{V}}_{\mathrm{LOS}}}{\cos \mathsf{\theta }}$$

(6)

${\mathrm{U}}_{\mathrm{INSAR}}$ is called LOS projected rate in the up direction, which will be compared with vertical displacement rates from leveling. From Equations (5) and (6), the error associated with ${\mathrm{U}}_{\mathrm{INSAR}}$ is

$$\Delta {\mathrm{U}}_{\mathrm{H}}=\tan \mathsf{\theta }(\cos (\mathsf{\phi })\mathrm{E}+\sin (\mathsf{\phi })\mathrm{N})$$

(7)

4. Results

4.1. Inland and Coastal Subsidence from Leveling: Understanding the Mechanism by Compaction Well Data

First, we show the land subsidence results in Chiayi using the leveling and MLCW data described in Section 2.2. In addition to showing the state of subsidence, these sensor results will be used to remove the systematic errors from ALOS/PALSAR below (Section 5.2). Here we emphasize two issues: groundwater-induced dynamic compaction and coastal subsidence.

Figure 5 shows the rates of height change using the annually measured, least-squares adjusted heights at benchmarks from 2007 to 2011. The maximum rate is 4.5 cm/year in western Chiayi. The primary subsidence occurred in the coastal townships of Dongshih, Budai and Yijhu, where aquaculture is the major industry. The subsidence pattern in Figure 5 suggests that the aquaculture industry pumps groundwater from the aquifers in Figure 2 from the wells (Figure 1) to replenish water in fishing ponds, thereby creating land subsidence in this region. However, the depths at which groundwater is extracted in these villages cannot be estimated from the leveling-derived subsidence in Figure 5.

Figure 6a shows the cumulative compactions at the AHES MLCW in the inland area of Chiayi, associated with the variation of the groundwater level and rainfall. Figure 6a shows significant changes of groundwater level at depths 96 m and 164 m. The overall groundwater level has dropped since 2005. The groundwater level was the lowest in 2010 and 2011, corresponding to the years of least rainfall. As a result, the subsidence in 2010 and 2011 was the largest in the entire Chiayi. Figure 6a also suggests a high correlation between soil compaction and groundwater level change at the monthly time scales. During the lowering of groundwater level, an external force is applied to the soil, resulting in an effect similar to dynamic compaction [21,22], in which soil is compressed by continual, external forces. In the case of AHES, the seasonal and multi-year lowering of groundwater level increase soil stresses that are equivalent to external forces. The monitoring well at AHES provides evidence of dynamic compaction by groundwater lowering. The accuracy of the MLCW observations is about one mm, which cannot be achieved by continuous GPS observations at the daily intervals. Precision leveling is typically carried out once a year, so it cannot reveal the seasonal compactions seen by the multi-layer sensor at the AHES in Figure 6a. The mm-level accuracy (one single measurement, not time-average) is the unique feature of WLCW in subsidence monitoring that cannot be matched by other sensors, especially in detecting small dynamic compactions.

Figure 6b suggest that the soil at the inland station AHES is primarily composed of coarse sand and clay, with the major compression occurring mainly at layers between 0 and 185 m. The major compression is related to the declines of groundwater level at the depths between 96 and 164 m (Figure 6b).

Figure 7a shows the cumulative compactions at three coastal MLCWs (BDES, WLES and DSES), and changes of groundwater level at BDES, and rainfalls at BDES. In addition to linear trends, the compactions also show seasonal oscillations coinciding with groundwater changes. At BDES, the overall groundwater level has dropped since 2005 and reached the lowest levels in 2010 and 2011, during which the rainfalls were also the least (see the second plot in Figure 7a; the spikes of rainfall are caused by typhoons and extreme monsoonal rains). The groundwater level at the 201-m depth experiences the most rapid level change. In fact, the groundwater levels here have dropped to a level below the past maximum drawdown (or the past maximum effective stress). However, the groundwater level here has started to rise again from 2011 onward. The MLCW data show a close relation between the land subsidence depth and groundwater level: a soil layer is compressed as the groundwater level drops, and the layer is expanded as the groundwater level rises. The compression and expansion of the layer increased and decreased the surface elevations. Although the surface elevations fluctuate with the groundwater levels, the elevations show a continuous subsiding tendency due to continual groundwater loss.

Figure 7b shows the cumulative compactions and drilling data at WLES, which is a coastal station. While the major compaction at AHES (inland station) is steep and occurs at depths 0–200 m, the major compaction at WLES (coastal station) varies gently from a depth of 300 m to the surface. Groundwater level variations at WLES are representative of those in the coastal area of Chiayi, where groundwater extractions occur at both shallow and deep aquifers. As indicated in Figure 7b, the soil in Chiayi’s coastal region is primarily composed of fine sand and clay. The entire stratum is filled with highly compressible soil. In particular, the layer at the depths of 200–300 m has the largest compression.

Finally, we compare land subsidence rates from precision leveling, MLCW and continuous GPS at the coastal station BDES, as shown in Figure 8. The GPS-derived vertical displacements are based on the daily solutions by the software Bernese 5.0 using relative positioning to KMNM (near the mainland China), which is an IGS station. Figure 8 indicates that the vertical displacement rates over 2007–2011 from these three sensors are consistent to about 1 cm/year. Typically, groundwater level declines in the dry season and rises in the wet season. The trends of the vertical displacements from GPS and monitoring well are all negative, indicating land subsidence. When groundwater level rose (Figure 8), GPS and MLCW records show elevation rebounds, which are not seen in the leveling record. The rising elevations are associated with rising groundwater levels in the wet season and during short raining times. This association of elevation and groundwater level is seen at both the coastal station (BDES, Figure 7a) and the inland station (AHES, Figure 6a). The coincident movements of groundwater level and elevation suggest short-term elastic rebounds of the soil at these stations, which can be used to partially alleviate land subsidence problems.

The above analysis on the state of Chiayi subsidence is primarily based on the measurements by precision leveling and magnetic ring sensors, which deliver results at scattered sites. If the budget is a concern, only a low site density is allowed, making it possible to miss areas with land subsidence. This spatial resolution problem may be resolved by using satellite-based methods such as InSAR. In Section 3, we will show how InSAR makes up the spatial density deficiency of leveling and MLCW by providing densely covered subsidence measurements in Chiayi, particularly in the coastal area with fish ponds.

4.2. Chiayi Subsidence from InSAR

Using TCPInSAR in Chiayi, we identified 687,736 TCP points (pixels) (Figure 9), covering an area of approximately 771 km². This is equivalent to a point density of 890 pixel/km². In particular, the spatial density of TCP over Chiayi’s coastal area is very high compared to those from precision leveling and monitoring wells, except near the townships of Budai and Dongshih. The two coastal townships are covered with fish ponds, where only interferograms on sufficiently wide paved roads can be produced. In the two townships, our leveling benchmarks were installed along the roads bordering the fish ponds, where land subsidence is evident (Figure 5). The roads containing the benchmarks are also covered with TCPs. Furthermore, the TCP point density around the two townships is much higher than the point density of leveling benchmarks.

In this paper, the method of Kriging [23] was then used to interpolate the InSAR-derived rates ${\mathrm{U}}_{\mathrm{INSAR}}$ on a grid (Figure 10; see also Section 3.2) for comparison with the vertical displacement rates on benchmarks (Section 5.1). Note that However, as indicated in Figure 1a, the horizontal velocities in Chiayi are significantly larger than zero, so this assumption can introduce errors in the InSAR result. In Figure 10, we show the errors in ${\mathrm{U}}_{\mathrm{InSAR}}$ due to neglecting the horizontal velocities at the GPS stations in Figure 1a. The largest error is 8.3 mm/year in eastern Chiayi (within our study area). In the coastal area, the errors are between 0.1–4.6 mm/year. Because of the sparse point density of GPS (Figure 1a), currently a reliable velocity field from GPS for Chiayi is not achievable.

The InSAR result (Figure 10) shows a prominent subsidence pattern in the coastal area, but a rather minor pattern in the inland area. However, the magnitudes of the InSAR-derived vertical displacement rates are different from those from precision leveling (Figure 5), especially in the coastal area; the differences will be explained in Section 5. Like the leveling result, the subsidence rates from InSAR decrease progressively from the coastal area to the inland area. The maximum subsidence rate in the coastal area is 4.0 cm/year and the pattern of subsidence is consistent with that from leveling (Figure 5). In conclusion, our InSAR result (Figure 10) is similar to that from leveling (Figure 5), but shows more detailed subsidence signatures not seen by leveling, especially in the coastal zone covered with fish farms.

5. Discussion

5.1. Assessing the InSAR Result and Identifying Systematic Errors

5.1.1. Comparison between the InSAR and Leveling Results

Despite the similar patterns seen in the InSAR and leveling results (Figure 10 vs. Figure 5), there exist substantial differences between the two results. Here we quantify their differences and identify the weakness of InSAR in Chiayi as follows. Figure 11 compares the vertical displacement rates derived from precision leveling and InSAR at 144 benchmarks. Note that here the InSAR result refers to the LOS projected rates in the up direction (Equation (6)), which have an opposite sign with that of subsidence rate from leveling. The differences between the two fields range from 1 to 2 cm/year (absolute values). The root mean squared (RMS) value of the differences is 0.89 cm, and about 74% (107 points) of the differences fall within 1 cm. For a more in-depth assessment of the InSAR result, we divided the areas with large differences into three blocks (A–C, Figure 11). The average differences between InSAR and leveling in Blocks A and B are about 1.0–1.5 cm/year. The large differences in Block A are along the bank of a river bordering Chiayi County and Yunlin County. The cause of the large differences in Block A may be the low coherence over the river bank. Block B is mostly covered by fish ponds around the Yijhu Township. As stated earlier, around fish farms, TCP points (Section 3.1) can exist only over sufficiently large paved roads and objects that can serve as stable scatterers. As such, the number of TCP points in Block B is relatively low (Figure 10) and the overall accuracy from TCPInSAR can be degraded. However, despite the dense fish ponds in areas south of Block A and in areas south and southwest of Block B, the differences are small, except at three benchmarks in the Budai Township.

Block C is on the eastern side of THSR. The differences here (Figure 11) are between 1–2 cm. These consistently large differences are partially caused by the assumption of zero horizontal displacement. According to Figure 1a and Yen et al. (2008) [24], significant horizontal displacements exist in the eastern part of Chiayi due to diastrophism. Because the incident angle of ALOS satellite is 38.7° and the horizontal displacement is set to zero, the LOS-projected vertical rates can be severely biased (see Equations (5) and (6)). In addition, large differences tend to occur in areas of small subsidence: if one draws a NE-SW line passing through the center of Chiayi, the areas of small subsidence are roughly located on the eastern side of this line.

5.1.2. InSAR’s Systematic Errors in Chiayi

Our assessment of the InSAR result in the previous section shows that, although the 15 ALOS/PALSAR images and the sophisticated TCPInSAR technique can generate convincing rates of land subsidence in Chiayi, some of the rates can differ from those from leveling by up to 2 cm/year (Figure 10). Here we list potential factors that can degrade the accuracies of the rates from our InSAR processing. The purpose of presenting these factors is to provide a lesson for ALOS/PALSAR applications in a place such as Chiayi that is covered by rice and fish farms.

(1) Sensitivity of L-band radar

The ALOS radar has a wavelength of 23.6 cm that enables a stronger penetrating power to vegetation compared to radar with a shorter wavelength, e.g., a C-band radar. Although ALOS/PALSAR is more likely to detect deformations in areas with vegetation, it is not sensitive to small deformations (<one cm/year). Specifically, notable interferograms can be acquired only if the deformation exceeds half of the radar wavelength. In addition, over wetlands covered with vegetation, interferograms from an L-band radar such as ALOS show better coherence than those obtained using a C-band radar such as Radarsat-1 [25].

(2) The quality of InSAR images

We used 15 SAR images spanning 4 years to determine vertical displacement rates, mostly at the cm accuracy level. We expect to improve the vertical accuracy given a longer time span and a large number of SAR images. Despite the short time span and low number of images, the present result is important in demonstrating the potential of ALOS/PALSAR in detecting large land subsidence over an area covered mostly by fish ponds (Chiayi’s coastal zone), as well as small subsidence over areas covered with vegetation and other attributes (inland zone). Use of a longer time span can also increase the vertical accuracy from InSAR. Moreover, more quality assessments can be made to select best SAR images. Another problem is caused by the time-varying vertical rates: The InSAR images span over March 2007–March 2011, while the leveling data span over September 2007–August 2011. The mean vertical rates over these two different time spans can be different to affect the quality assessment of the ALOS result. Given a longer time span for InSAR images, it is more likely to achieve optimizations in best image selection and best quality assessment.

(3) Effect of un-modeled horizontal displacement

The incident angle of ALOS/PALSAR is larger than that of a C-band satellite such as ERS and Envisat, making ALOS/PALSAR more susceptible to errors in vertical rate than the latter, especially when significant horizontal displacements exist in the study area. If horizontal velocities are known, they should be used when converting LOS velocities to vertical velocities. As stated earlier, the horizontal velocities in Chiayi are assumed to be zero and this assumption has led to errors in the ALOS-derived vertical rates. To reduce this error, GPS-derived horizontal velocities in a dense network is needed.

(4) Limited TCP point density around fish farms

Only sparse TCP points, mainly on paved roads and man-made objects, were identified over Chiayi’s fish farm area. If a paved road is too narrow to reflect ALOS/PALSAR’ radar signals, no TCP can be generated. The limited TCP point density has reduced its capacity in monitoring coastal subsidence and relative sea level rise. A solution will be to deploy simple reflection corners in places where major subsidence has occurred but only few TCP points are found. See also the example of Budai in Section 5.3.

(5) The effects of atmosphere and ionosphere

The atmospheric delay is one of the most important limitations in InSAR measurements of deformation [26,27]. Taiwan is situated in a subtropical area surrounded by oceans (Figure 1), and is subject to monsoonal rains and typhoons. The yearly mean temperature, relative humidity and cloud-cover percentage are about 23.8 °C, 74.1% and 79%, respectively. These combine to create large atmospheric delays for ALOS/PALSAR and maybe other satellites. In addition, the L-band of ALOS/PALSAR is more prone to ionospheric effects and more sensitive to moisture than a C or X-band radar, which has a shorter wavelength and can suffer from stronger scattering by rain drops or aerosols. A solution will be to use a mesoscale meteorological model to account for atmospheric delays over Taiwan.

5.2. Correcting the InSAR Result Using Rates from Leveling

Removing the systematic errors shown in Section 5.1 may require a considerable time and effort. Here we use the leveling-derived vertical displacement rates to make a quick fix of the errors, following the approach of Hung et al. (2011) [28]. First, the differences between the vertical displacement rates from InSAR (Figure 5 and Figure 10, Section 4.1) and leveling at the leveling benchmarks (Figure 3) were computed. The differences were then used to form corrections on a grid using the Kriging method of interpolation [23,28]. Finally, the gridded corrections were added to the gridded vertical displacement rates from InSAR (Figure 10) to form a grid of combined rates (Figure 12). In this way, the InSAR-derived rates are corrected by the leveling result, and in theory contain less horizontal motion-induced biases (Figure 10) compared to the original InSAR result.

Figure 13 compares the rates from the corrected InSAR field and from leveling at the leveling benchmarks along the three profiles (see Figure 12 for their locations). The InSAR rates at the benchmarks were interpolated from the grid in Figure 12. Compared to the original InSAR rates (Figure 10), the corrected rates contain much less systematic errors (differences in Figure 12 vs. differences along the three profiles in Figure 13). In addition, the corrected InSAR field in Figure 12 shows more subsidence details than those in the leveling-only field. The improved detail is due to the fact that InSAR delivered 687,736 TCP points, compared to 144 leveling benchmarks (Section 4.2). The corrected InSAR field shows the overall subsidence trend in Chiayi: subsidence is large in the coastal area and decreases toward the southeast. This trend cannot be correct without the correction by the leveling result. The corrected field can also identify potential locations of subsidence that are missed by the current leveling result and this is a subject to explore.

Here we present an analysis of the merged (corrected InSAR) field along the three profiles in Figure 13 to address the subsidence problem in Chiayi. Figure 13a shows the vertical displacement rates from the merged field along profile A-A’ in eastern Chiayi, roughly parallel to the rail of THSR. The subsidence rates along this profile decrease gradually from the south to the north. At the distances of 4–8 km (counted from the starting point) along A-A’ near AHES, a bowl-shaped subsidence pattern occurs. The subsidence here is also indicated by the compactions recorded at AHES (a MLCW in Figure 7b). The compaction records at AHES show major compressions occurring at shallow depths. These shallow compactions may reduce the strength of the columns supporting the rail of THSR [2]. Figure 13b shows the vertical displacement rates along B-B’s in western Chiayi. Along B-B’, the subsidence rates increase from about 2.5 cm/year to 4.5 cm/year (near 11 km), and then remain at the 3.5 cm/year level (14 km and beyond). This pattern of rate variation is consistent with that indicated by the four coastal compaction wells, DSES, WLES, BDES and NSES (Figure 3; their records are not fully shown in this paper). That is, the largest subsidence occurs at the midpoint between the townships of Dongshih and Budai, with a maximum rate of 4.5 cm/year. Finally, Figure 13c shows the vertical displacement rates along a west-east profile (C-C’) in central Chiayi. The rates along C-C’ suggest that the overall subsidence rate is the largest in the coastal area and the rate decreases towards the inland area. Again, the pattern of subsidence rates seen along C-C’ is associated with the large groundwater loss and the compressible soil formation in Chiayi’s coastal area, compared to the subsidence pattern in Chiayi’s inland area that suffers from less groundwater pumping and compressible soils.

The three profiles in Figure 13 show representative, subsidence-affected locations in Chiayi. With the hydrogeological data (Figure 2) and the groundwater level and MLCW records (Figure 6, Figure 7 and Figure 8), the subsidence rates along the profiles can be used to show the cause and consequence of land subsidence in Chiayi, providing useful information for safe operations of THSR and for planning a sustainable land development here.

5.3. A Detailed Aquaculture-Induced Subsidence over Fish Farm from the Corrected InSAR

Here we show the advantage of InSAR over leveling in detecting detailed land subsidence over an area covered with fish ponds. Figure 14 shows the vertical rates from the corrected InSAR (Section 5.2) and from leveling (Section 2.2) over Chiayi’s Budai Township. Budai’s major industry is fish farming, and it is best known for its milkfish and oysters. In Figure 14, there are only 21 leveling benchmarks, compared to thousands of TCP points here. Thus, InSAR has clear advantages over leveling in increasing point density and in showing areal extent of subsidence. The vertical rates from InSAR are over paved roads, buildings, and any natural and man-structures that have resulted in TCPs using our data selection procedure in Section 3. In Section 5.1, we have shown the differences between the rates from InSAR and leveling. We now show more detail of the differences in Figure 14. In general, the differences in rate in Figure 14 are within 1 cm/year. In Figure 14, there are several regions with subsidence rates close to 5 cm/year. Without the InSAR result in Figure 14, these subsidence-affected regions cannot be identified. Furthermore, our InSAR result can provide information for the Budai government to take measures to prevent further subsidence and to manage the spots that are undergoing serious subsidence. The example here shows the value of InSAR in detecting detailed subsidence over an aquaculture area such as Budai.

5.4. Implication for Relative Sea Level Rise

The analysis in Section 5.2 and Section 5.3 indicates large land subsidence rates in Chiayi’s coastal area, which are clearly caused by groundwater pumping for the aquaculture industry (Section 4.2). Chiayi’s coastal subsidence raises the issue of relative sea level rise [16,17,18,28]. In this paper, the maximum coastal subsidence rate in Chiayi is 4.5 cm/year (Figure 5 and Figure 12), which is about 15 times of the global eustatic sea level rising rate of 2–3 mm/year [18]. In Figure 15, we show recent tidal records at the Dongshi and Wengang tide gauges in Chiayi; their locations are shown as DSTG and WGTG in Figure 2. The rates of sea level change at the two stations are 43 and 26 mm/year, respectively, which far exceed the global rate. The two rates from the tidal records are consistent with those from leveling, MLCW and InSAR (Figure 5 and Figure 15, with opposite signs). This implies that ALOS/PALSAR can detect relative sea level rise in coastal Chiayi, despite the unfavorable conditions for a successful InSAR detection of land subsidence here (see the systematic errors in Section 5.1).

The direct risk of Chiayi’s coastal subsidence is poor drainage that can easily lead to flooding during rainstorms or even minor rain events. Additional risks are salinization of coastal soils and aquifers. Salinization will gradually deplete organic matters from the affected soils to destroy the ecological system here, rendering permanent waste lands. Coastal aquifers, if managed properly and without salinization, can supply long-term drinkable water and the needs for industries and other vital functions, but salinized aquifers will cease to provide all such functions.

Another risk of coastal land subsidence, unique to Chiayi, is arsenic pollution. In the 1950s, arsenic pollutions have caused blackfoot disease in southwestern Taiwan, including Chiayi [29]. A recent study [30] showed that arsenic can move from shallow aquifers to deep aquifers through deformed aquitards induced by land subsidence. As such, deep aquifers with previously little arsenic can be contaminated by arsenic and can no longer supply clean groundwater for drink and for industrial uses.

Despite these risks, aquaculture is still a highly efficient food production method and has clear environmental benefits over other forms of animal food production [16]. Furthermore, land-based aquaculture has become increasingly important for a stable supply of protein for countries that see a decreasing oceanic supply of fish and the need to protect their oceanic ecology systems from over-fishing. In addition to Taiwan, many countries, particularly those in Asia, have developed extensive aquaculture industries. According to FAO (2013) [31], countries in Asia supplies nearly 89% of world’s aquaculture products. It is reasonable that these countries may have extracted groundwater for their aquaculture industries, and suffered from relative sea level rise seen in Chiayi’s coastal area [32,33,34]. Monitoring coastal land subsidence is the first step to prevent the risks of relative sea level rise.

6. Conclusions

This paper shows the benefit of using multiple sensors to detect and to correct the InSAR-derived land subsidence in Chiayi, a region rich in fish and rice farms sustained by extensive groundwater extractions. Precision leveling is a highly accurate method, but is costly and renders only sparse data points. A MLCW is even more costly than leveling, but it can detect seasonal and secular compactions at different depths, and at different aquitards and aquifers, supplying unique data sets to understand the nature of land subsidence in connection with groundwater extractions. InSAR is fast, less costly and potent in supplying dense points (pixels) of subsidence observations. However, InSAR-derived subsidence rates can be prone to many systematic errors; in the case of Chiayi, we fixed such errors by the leveling result. In summary, all these sensors detect land subsidence in most parts of Chiayi, and the maximum subsidence rate is 4.5 cm/year near the coast. The major cause of the subsidence is groundwater pumping at different aquifer depths, as revealed by MLCW observations.

The dense fish farms in Chiayi’s coastal area were converted from rice farms in the 1980s. The essential difference between fish farming and rice farming in terms of groundwater use is the former pumps groundwater progressively from shallow to greater depths, resulting in larger and widely spread land subsidence. In fact, Chiayi’s aquaculture has led to a relative sea level rise at a rate that is one order of magnitude larger than the global rate of sea level rise. This lesson of land subsidence in Chiayi is important for any aquaculture-rich coastal area in the world that can suffer from land subsidence problems similar to that in Chiayi. Monitoring potential land subsidence in an aquaculture-rich region is the first step toward a sustainable, healthy aquaculture industry, which has been recognized by United Nations’ Food and Agriculture Organization (FAO) as a critical method to supplement oceanic fishing for supplying world’s protein needs.