The Relationship between Surface Displacement and Groundwater Level Change and Its Hydrogeological Implications in an Alluvial Fan: Case Study of the Choshui River, Taiwan

Balancing the demand of groundwater resources and the mitigation of land subsidence is particularly important, yet challenging, in populated alluvial fan areas. In this study, we combine multiple monitoring data derived from Multi-Temporal InSAR (MTI), GNSS (Global Navigation Satellite System), precise leveling, groundwater level, and compaction monitoring wells, in order to analyze the relationship between surface displacement and groundwater level change within the alluvial fan of the Choshui River in Taiwan. Our combined time-series analyses suggest, in a yearly time scale, that groundwater level increases with the vertical surface displacement when the effect of pore water pressure dominates. Conversely, this relationship is negative when the effect of water-mass loading predominates over pore water pressure. However, the correlation between the vertical surface displacement and the groundwater level change is consistently positive over the time scale of two decades. It is interpreted that the alluvial fan sequence in the subsurface is not fully elastic, and compaction is greater than rebound in this process. These findings were not well reported and discussed by previous studies because of insufficient monitoring data and analyses. Understanding the combined effect of groundwater level change and vertical surface displacement is very helpful for management of land subsidence and usage of groundwater resources. The spatial and temporal integration of multi-sensors can be applied to overcome the limitations associated with the single technique and provides further insights into land surface changes, particularly in highly populated alluvial fan areas.


Introduction
Alluvial fans are major reservoirs of groundwater in piedmont areas around the world. Formed by deposits of water-transported material, fan aquifers are a complex system of multiple aquifers and aquitards. The relationship between vertical surface displacement and groundwater change plays an important role in the management of groundwater resources. The withdrawal of groundwater from fan aquifers usually causes compaction of unconsolidated deposits and land subsidence in the time scale of decades due to the decrease of pore-water pressure (e.g., [1][2][3][4][5]). On the other hand, the hydrological mass loading by groundwater, free surface water, soil moisture, and snow is considered as an important factor for the seasonal surface displacement recorded in GNSS (Global Navigation System) time respectively. According to previous investigations [2,42,43], this alluvial fan is composed of unconsolidated sediments, mainly gravel, sand, and clay. The sediments originate from rock formations in the Western Foothills and the Central Range, including slate, shale, metamorphic quartzite, sandstone, and mudstone ( Figure 1). The thickness of the sediments varies from 750 to 3000 m, and the mean grain size decreases from east to west. Thus, the head of the alluvial fan mainly contains gravel and coarse sand, whereas its foot comprises soil and fine sand. The stratigraphic columns of the head and foot of the alluvial fan, which were derived from drilling data at the Pingting and Hefeng groundwater stations, are shown in Figure 2a. Due to the marine transgression and regression, frequent flooding, and channel migration of the Choshui River, the stratigraphic asset of this flood plain is very complex. Based on sedimentological data, the alluvial fan of the Choshui River can be divided into four aquitards and four aquifers at depths from 0 to 300 m. Figure 2b shows a conceptual hydrogeological model of the alluvial fan of the Choshui River. The four aquifers are essentially connected in the head of the alluvial fan [43]. Aquifer 2 (labeled as F2) is the thickest and largest among the four aquifers, with an average thickness of about 95 m. Previous studies demonstrated that most groundwater resources are pumped from Aquifers 2 and 3 [2,33,42,43]. The annual rainfall is approximately 1200-2000 mm in the area. A large portion of the precipitation occurs from May to September and provides abundant groundwater resources in the dry season.   Figure 1 of the alluvial fan of the Choshui River (modified from [42]). There are four aquitards and four aquifers at depths from 0 to 300 m, and Aquifer 2 (F2) is the thickest and largest among the four aquifers.

Evolution of Land Subsidence
The alluvial fan of the Choshui River is the area that has suffered the most significant subsidence in Taiwan over the past 30 years [2]. In general, land subsidence may be natural and human induced. The main factors controlling natural land subsidence in the alluvial fan of the Choshui River are sediment compaction and the earthquake effect, such as soil liquefaction and co-seismic deformation [44]. On the other hand, human-induced land subsidence is related to groundwater overpumping and the loading of buildings. In the past 10,000 years, the natural land subsidence rate was about 0.6-2.1 mm/year based on borehole data [42,44]. With a growing population and rapid economic development, groundwater has been overpumped for residential, industrial, aquacultural, and agricultural use in Taiwan for decades. According to precise leveling data from 1995 to 2017, the land subsidence areas were located mainly along the shoreline before 2000, including Dacheng in the Changhua County, and Taishi in the Yunlin County (Figure 1). After 2000, the centers of land subsidence gradually shifted from the coastal area to the eastern inland area. There are three land subsidence centers located in Shihu, Erlin, and Shijou in the Changhua County, whereas the land subsidence centers in the Yunlin County are found at Baozhong, Tuku, Huwei, and Yuanchang ( Figure 1). Based on precise leveling data, the land subsidence rate had already decreased, but the maximum in 2017 still reached 3.5 cm/year and 6.7 cm/year in Shijou and Tuku, respectively. Furthermore, the safety concerns of the Taiwan High Speed Rail, which runs through the central land subsidence area in the Yunlin County, have received a lot of attention. The stratigraphic column of the head of the alluvial fan derived from the drilling data of the Pingting groundwater station shows sediment deposits of mainly gravel and coarse sand. The stratigraphic column derived from the drilling data of the Hefeng groundwater station at the foot of the alluvial fan shows sediment deposits of mainly soil and fine sand (both columns are modified from an internal project report by the Taiwan Water Resources Agency). The locations of these groundwater stations are shown in Figure 1 as blue squares. (b) The conceptual hydrogeological profile along the line aa' in Figure 1 of the alluvial fan of the Choshui River (modified from [42]). There are four aquitards and four aquifers at depths from 0 to 300 m, and Aquifer 2 (F2) is the thickest and largest among the four aquifers.

Evolution of Land Subsidence
The alluvial fan of the Choshui River is the area that has suffered the most significant subsidence in Taiwan over the past 30 years [2]. In general, land subsidence may be natural and human induced. The main factors controlling natural land subsidence in the alluvial fan of the Choshui River are sediment compaction and the earthquake effect, such as soil liquefaction and co-seismic deformation [44]. On the other hand, human-induced land subsidence is related to groundwater overpumping and the loading of buildings. In the past 10,000 years, the natural land subsidence rate was about 0.6-2.1 mm/year based on borehole data [42,44]. With a growing population and rapid economic development, groundwater has been overpumped for residential, industrial, aquacultural, and agricultural use in Taiwan for decades. According to precise leveling data from 1995 to 2017, the land subsidence areas were located mainly along the shoreline before 2000, including Dacheng in the Changhua County, and Taishi in the Yunlin County ( Figure 1). After 2000, the centers of land subsidence gradually shifted from the coastal area to the eastern inland area. There are three land subsidence centers located in Shihu, Erlin, and Shijou in the Changhua County, whereas the land subsidence centers in the Yunlin County are found at Baozhong, Tuku, Huwei, and Yuanchang ( Figure 1). Based on precise leveling data, the land subsidence rate had already decreased, but the maximum in 2017 still reached 3.5 cm/year and 6.7 cm/year in Shijou and Tuku, respectively. Furthermore, the safety concerns of the Taiwan High Speed Rail, which runs through the central land subsidence area in the Yunlin County, have received a lot of attention.

Data and Methods
In this study, we collected multiple datasets spanning a period of 20 years, including SAR (Synthetic Aperture Radar) images, precise leveling data, GPS data, groundwater, and compaction monitoring well records in order to analyze the vertical surface displacement of the alluvial fan of the Choshui River. The details of the datasets and the methods applied in this study are described below.

Data Acquisitions
We used multiple SAR images derived from different satellites, including ERS-1/2, ENVISAT, ALOS-1, and Sentinel-1A datasets (Table 1), for analyzing vertical surface displacement. ERS-1/2, ENVISAT, and Sentinel-1A are C-band satellites with a wavelength of about 5.6 cm. The maximum detectable velocity between neighboring PS pixels of ERS-1/2 and ENVISAT is about 14.6 cm/year, whereas that of Sentinel-1A is about 42.6 cm/year. This is because Sentinel-1A has a shorter revisit interval of 12 days. ALOS-1 is an L-band satellite with a wavelength of about 23.6 cm and a maximum detectable velocity of about 46.8 cm/year between neighboring PS pixels. The raw Synthetic-aperture radar (SAR) images derived from ERS-1/2, ENVISAT, and ALOS-1 were processed into the SLC (single look complex) images using ROI_PAC software developed by JPL/Caltech. Interferograms were derived with the SLC images using the Doris software [45]. The interferograms of the SLC images from Sentinel-1A were generated by the SNAP software, which was developed by the European Space Agency (ESA). The PSI procedure of all interferograms was processed with StaMPS/MTI (Stanford Method for Persistent Scatterers/Multi-Temporal InSAR). The distribution of the in situ data used in this study is shown in Figure 3. Precise leveling is one of the most accurate methods for measuring vertical surface displacement in the scale of several millimeters along leveling routes [46]. There are 807 precise leveling benchmarks collected from the Water Resources Agency (WRA) in the alluvial fan of the Choshui River. The earliest precise leveling record in the area dates back to 1992. In the network, the total length of the precise leveling routes is about 970 km with an average benchmark spacing of about 1.5 km. The specifications of the leveling survey demand that the maximum misclosure between forward and backward runs in a section is below 3 mm √ K (K being the distance between two neighboring benchmarks in km) [47]. The leveling measurements followed the standards of the Taiwan Vertical Datum 2001 (TWVD2001), which is determined by the mean sea level on 1 January 2001 of Keelung Harbor tide gauge in northern Taiwan [48,49]. By subtracting two adjacent precise leveling measurements, the land elevation change of every independent benchmark can be obtained. However, the measurement frequency of precise leveling survey was about one year, which is much lower than the GNSS and the Multi-Temporal InSAR (MTI) methods. The GNSS indicates several satellite positioning systems, including American GPS, Russian GLONASS, China's BeiDou, and European Galileo. In this study, the processed GNSS signals were only from the American GPS satellites over two decades for consistency as available positioning data in the early years were GPS data. There are 41 continuous GNSS monitoring stations in the alluvial fan of the Choshui River, which provide regular and frequent observations of three-dimensional Remote Sens. 2020, 12, 3315 6 of 20 surface displacement. Three GNSS stations (PINT, HNES, and YSLL) and nearby monitoring systems were selected for further analysis of vertical surface displacement, groundwater level change, and sediment compaction in different geological conditions of the alluvial fan. The GPS data received from the GNSS stations are provided by the GPS Lab of Academia Sinica, Taiwan [50]. These data were processed into the kinematic positioning of daily solution based on the static method by the GIPSY-OASIS II software [51] in the International Terrestrial Reference Frame (ITRF) 2008. Precise ephemerides provided by the International Global Navigation Satellite System Services (IGS) are employed to reduce satellite orbit and clock errors during the data processing [52]. The precision of GPS data in the vertical direction is generally about two to three times larger than that in the horizontal direction. The accuracy of GPS data in the vertical component can be better than 1 cm on a daily solution of the static mode [53,54]. Groundwater monitoring wells were gradually established and have been maintained throughout the alluvial fan since 1992, and they form a groundwater monitoring network for collecting complete and comprehensive groundwater information. There are 93 distributed groundwater stations in the network consisting of 226 monitoring wells. The monitoring wells reach different depths ranging from 24 to 306 m for continuously recording the groundwater level of different aquifers. The groundwater levels of the wells have been recorded automatically hourly since 1997. Multi-level compaction monitoring wells reach a depth up to 300 m and they are divided into 25 levels.
The device provides high-accuracy measurement (about 1-5 mm) of the aquifer compaction at different depths. By measuring the depth of each magnetic ring every month, the relative compaction between two rings can be determined. Then, the main compaction depth and range can be located within the monitoring wells.

Persistent Scatterer SAR Interferometry (PSInSAR)
In this study, we applied the PSInSAR method for overcoming the limitations of Differential SAR Interferometry (DInSAR) [19,20]. DInSAR can detect surface displacement with centimeter-level accuracy [23,55,56]. However, DInSAR is ineffective for monitoring a long time period of surface displacement because of spatial and temporal decorrelation, such as baseline errors, DEM residual effects, and atmospheric disturbances [24,[57][58][59]. If the value of a pixel is dominated by a stable and strong scatterer (such as a building, a bridge, or a rock body), which is brighter than the background

Persistent Scatterer SAR Interferometry (PSInSAR)
In this study, we applied the PSInSAR method for overcoming the limitations of Differential SAR Interferometry (DInSAR) [19,20]. DInSAR can detect surface displacement with centimeter-level accuracy [23,55,56]. However, DInSAR is ineffective for monitoring a long time period of surface displacement because of spatial and temporal decorrelation, such as baseline errors, DEM residual effects, and atmospheric disturbances [24,[57][58][59]. If the value of a pixel is dominated by a stable and strong scatterer (such as a building, a bridge, or a rock body), which is brighter than the background scatterers, the underlying displacement signal may be extracted from the backscattered phase information (Figure 4). This kind of pixel is defined as a persistent scatterer (PS) pixel. Various methods and thresholds have been developed to select stable and strong scatterers in interferograms [19][20][21]60]. PSInSAR is a method for identifying stable pixel scatterers based on coherence or amplitude from a series of interferograms to detect a long time period of surface displacement. First, a master image is selected based on temporal interval, perpendicular baseline, and Doppler centroid frequency baseline. By creating interferograms between master and slaves, and removing the topography of DEM, the result of DInSAR can be obtained. The PS candidate pixels are initially selected based on amplitude and phase stability. Then, a PS probability based on a combination of amplitude and estimated phase stability is applied to determine a PS pixel. With 3-D phase unwrapping, the displacement along the radar line of sight (LOS) direction can be measured from the phase difference of persistent scatterers in SAR images covering a long period of time.
Remote Sens. 2020, 12, x FOR PEER REVIEW 8 of 23 frequency baseline. By creating interferograms between master and slaves, and removing the topography of DEM, the result of DInSAR can be obtained. The PS candidate pixels are initially selected based on amplitude and phase stability. Then, a PS probability based on a combination of amplitude and estimated phase stability is applied to determine a PS pixel. With 3-D phase unwrapping, the displacement along the radar line of sight (LOS) direction can be measured from the phase difference of persistent scatterers in SAR images covering a long period of time.

Multi-Reference Points for PSInSAR
Through interferograms formation, phase stability estimation, PS selection, and displacement calculation, a long time period of surface displacement can be obtained. Furthermore, contrary to previous studies, which only used one reference point to adjust the PSInSAR result, here, we attempted to use all GPS data in the same time period covered by the SAR images as reference points. Based on the data fusion method [35,36], which was applied for improving the land subsidence monitoring resolution, the PSInSAR result could be integrated with in situ monitoring data. The workflow of data processing is shown in Figure 5. By projecting the three-dimensional vectors of GPS to the LOS direction of the satellite, the differences along LOS direction between the GPS data and the nearby average PS pixels can be obtained in a radius of about 150 m. Then, interpolating the differences by using the IDW method of smoothing, the calibrated value of PS pixels can be calculated. Small horizontal surface displacement of the alluvial fan of the Choshui River [40,41] is still considered in this study. The influence of the horizontal effect along the LOS direction is subtracted for accurate assessment of the land subsidence results derived from the PSInSAR method. Finally, land subsidence can be derived by projecting the residual value along LOS to the vertical direction with the incident angle of each satellite.

Multi-Reference Points for PSInSAR
Through interferograms formation, phase stability estimation, PS selection, and displacement calculation, a long time period of surface displacement can be obtained. Furthermore, contrary to previous studies, which only used one reference point to adjust the PSInSAR result, here, we attempted to use all GPS data in the same time period covered by the SAR images as reference points. Based on the data fusion method [35,36], which was applied for improving the land subsidence monitoring resolution, the PSInSAR result could be integrated with in situ monitoring data. The workflow of data processing is shown in Figure 5. By projecting the three-dimensional vectors of GPS to the LOS direction of the satellite, the differences along LOS direction between the GPS data and the nearby average PS pixels can be obtained in a radius of about 150 m. Then, interpolating the differences by using the IDW method of smoothing, the calibrated value of PS pixels can be calculated. Small horizontal surface displacement of the alluvial fan of the Choshui River [40,41] is still considered in this study. The influence of the horizontal effect along the LOS direction is subtracted for accurate assessment of the land subsidence results derived from the PSInSAR method. Finally, land subsidence

Results
The variation of land subsidence in the study area from 1995 to 2017 was calculated using the ERS-1/2 (1996)(1997)(1998)(1999), ENVISAT (2006ENVISAT ( -2008, ALOS-1 (2006-2011), and Sentinel-1A (2014-2017) SAR images, as well as in situ monitoring data, including precise leveling  and GPS (1995GPS ( -2017. The temporal trends of the vertical surface displacement were compared with the groundwater level change and the compaction of sediments at different depths. The vertical surface displacements in the past two decades can be characterized and analyzed to reveal possible relationships among the observed phenomena.

Spatial Variation of Land Subsidence
The spatial patterns of the land subsidence based on PSInSAR generally agree with those derived from the in situ monitoring data (precise leveling and GPS) ( Figure 6). However, the land subsidence rates in the coastal zone derived from the ERS-1/2 images were underestimated in comparison with the data derived from the precise leveling, which were interpolated by the inverse distance weighting (IDW) method. The average subsidence rate in the coastal area in Dacheng derived from ERS-1/2 is only about 43 ± 3.6 mm/year, whereas the measurement obtained from precise leveling nearby exceeds 200 mm/year. This difference may be attributed to the phase unwrapping errors on a lower density of scatterers and limitation of the maximum rate detectable with the ERS-1/2 dataset, which is around 146 mm/year under certain assumptions [24,38]. The average subsidence rate in the main subsidence area in Erlin of the Changhua County derived from ENVISAT-based PSInSAR results is about 71 ± 5.4 mm/year. The average subsidence rate derived from ALOS-based PSInSAR results in the main subsidence area in Tuku of the Yunlin County is greater than 63 ± 4.1 mm/year. According to our analysis of the Sentinel-1A images from 2014 to 2017, the recent conditions of land subsidence can be obtained. The results show that the main subsidence area is in Tuku of the Yunlin County, where the average subsidence rate is about 50 ± 3.4 mm/year. Table 2 shows comparisons between the subsidence rates derived from PSInSAR and precise leveling in the main subsidence area. Through comparisons of the results from different C-band satellite images (ERS-1/2, ENVISAT, and Sentinel-1A), the result from Sentinel-1A provides better spatial and temporal information due to its shorter satellite revisit frequency (12 days). In general, the land subsidence rate is gradually decreasing with time in the alluvial fan of the Choshui River.

Subsidence Velocity along the Taiwan High Speed Rail
The land subsidence rate along the Taiwan High Speed Rail can be acquired by interpolating precise leveling and PSInSAR data in different time periods. Due to the construction of the Taiwan

Subsidence Velocity along the Taiwan High Speed Rail
The land subsidence rate along the Taiwan

Time Series of Surface Displacement and Groundwater Level
To further understand the characteristics of the alluvial fan, three GNSS stations were selected with available nearby precise leveling benchmarks. The PS pixels in a radius of 1.5 km from the GNSS stations were used for displacement analysis. The GPS, leveling, and PS data were calculated to present the time series of vertical surface displacement (Figure 8). The temporal trends of displacement derived from ENVISAT and Sentinel-1A images are more consistent with those derived

Time Series of Surface Displacement and Groundwater Level
To further understand the characteristics of the alluvial fan, three GNSS stations were selected with available nearby precise leveling benchmarks. The PS pixels in a radius of 1.5 km from the GNSS stations were used for displacement analysis. The GPS, leveling, and PS data were calculated to present the time series of vertical surface displacement (Figure 8). The temporal trends of displacement derived from ENVISAT and Sentinel-1A images are more consistent with those derived from GPS and precise leveling. On the head of the alluvial fan, around the PINT station, the time series of the vertical surface displacement shows a very small amount of land subsidence of 0.46 cm/year from 1997 to 2017 (Figure 8a)

Calculated Compaction of the Aquifer System
The compaction conditions of the aquifer system at different depths from the inland area to the coastal area can be calculated from the records obtained from the multi-level compaction monitoring wells near the PINT, HNES, and YSLL GNSS stations (Figure 9). These three monitoring wells provide data for the periods of 2007-2017, 2005-2017, and 1996-2017. To compare the compaction condition during different time periods, the data from each monitoring well were divided into two time periods. The compaction ratio was calculated by dividing the compaction of a time period with the compaction during the entire time period. Based on the numbers of the compaction ratio in different time periods, the compaction condition gradually slowed down. On the head of the alluvial fan, the records from the monitoring well near PINT GNSS station show very small amounts of compaction and expansion (Figure 9a). The sediment compaction between the head and foot of the alluvial fan is presented in Figure 9b. Compaction is reduced at depths of about 50, 150, and 250 m, which indicates the existence of an aquitard. The aquifers and aquitards can be identified in the compaction profiles where the amount of compaction is large for the aquifers and small for the aquitards. The amount of consolidation can be inferred by local geological conditions. The amount of the groundwater level change is almost the same at different depths (Figure 8b), but the compaction mainly occurs in the depth range of 150-250 m (Aquifer 3). We infer that the water level change at deep depths resulted in more sediment compaction than the water level change at shallow depths. This is likely related to the overall mass loading. According to the monitoring data near the YSLL GNSS station at the foot of the alluvial fan, compaction mainly occurs at depths below 250 m (Figure 9c). Sediment compaction continued although the groundwater level became steady and even increased after 2004 (Figure 8c). Therefore, the effect of overpumping on land subsidence may continue for many years with the subsidence rate decreasing gradually.
with the rise of the groundwater level from 2007 to 2017 (Figure 8c). Thus, the timing of land subsidence is likely related to and triggered by the groundwater level change.

Discussion
The surface displacement of a highly populated alluvial fan area was investigated by integrating multiple monitoring datasets both spatially and temporally. The integrated results were used to obtain correlations between groundwater level change and vertical surface displacement for the

Discussion
The surface displacement of a highly populated alluvial fan area was investigated by integrating multiple monitoring datasets both spatially and temporally. The integrated results were used to obtain correlations between groundwater level change and vertical surface displacement for the entire alluvial fan of the Choshui River. Based on the correlations, a seasonal interaction model of the groundwater level change and the surface displacement is proposed in this study.

Model of the Interaction between Vertical Surface Displacement and Groundwater Level Change
According to previous studies, the relationship between groundwater level change and vertical surface displacement depends on the pore water pressure and the water mass loading [62][63][64]. Here, we used groundwater and GPS data to explore this interaction in areas of the alluvial fan of the Choshui River. We collected more complete monitoring data covering a longer time period than the data presented in Chiang et al. [64] in order to analyze the combined effects of the pore water pressure and the water mass loading on the vertical surface displacement. We applied detrending on the dataset with linear regression to acquire the seasonal variation, and used the centered moving average with a one-month window size to reduce fluctuation. Thus, we obtained a seasonal time series of the vertical surface displacement and the groundwater level change (Figure 10). It shows both negative and positive relationships between the vertical surface displacement and the groundwater level change on the head of the alluvial fan near the PINT GNSS station and at the foot of the alluvial fan near the YSLL GNSS station. The outliers of the detrending time series were deleted by using the interquartile range method, which defines the outliers as data points 1.5 interquartile away from the median. The relationships between the two datasets were calculated by Pearson's correlation coefficient per year in each site, and the mean correlation coefficients of each site were then derived. On the head of the alluvial fan, the correlation coefficient between the groundwater level change and the vertical surface displacement derived from the GPS data is negative. For instance, the correlation coefficient is −0.25 between the Pingting groundwater station and the PINT GNSS station (Figure 10a), which are near the Pingting Tableland at the head of the alluvial fan. On the other hand, a positive correlation coefficient of 0.66 is obtained at the foot of the alluvial fan, based on the Hefeng groundwater station and the YSLL GNSS station (Figure 10b).
Remote Sens. 2020, 12, x FOR PEER REVIEW 16 of 23 entire alluvial fan of the Choshui River. Based on the correlations, a seasonal interaction model of the groundwater level change and the surface displacement is proposed in this study.

Model of the Interaction between Vertical Surface Displacement and Groundwater Level Change
According to previous studies, the relationship between groundwater level change and vertical surface displacement depends on the pore water pressure and the water mass loading [62][63][64]. Here, we used groundwater and GPS data to explore this interaction in areas of the alluvial fan of the Choshui River. We collected more complete monitoring data covering a longer time period than the data presented in Chiang et al. [64] in order to analyze the combined effects of the pore water pressure and the water mass loading on the vertical surface displacement. We applied detrending on the dataset with linear regression to acquire the seasonal variation, and used the centered moving average with a one-month window size to reduce fluctuation. Thus, we obtained a seasonal time series of the vertical surface displacement and the groundwater level change (Figure 10). It shows both negative and positive relationships between the vertical surface displacement and the groundwater level change on the head of the alluvial fan near the PINT GNSS station and at the foot of the alluvial fan near the YSLL GNSS station. The outliers of the detrending time series were deleted by using the interquartile range method, which defines the outliers as data points 1.5 interquartile away from the median. The relationships between the two datasets were calculated by Pearson's correlation coefficient per year in each site, and the mean correlation coefficients of each site were then derived. On the head of the alluvial fan, the correlation coefficient between the groundwater level change and the vertical surface displacement derived from the GPS data is negative. For instance, the correlation coefficient is −0.25 between the Pingting groundwater station and the PINT GNSS station (Figure 10a), which are near the Pingting Tableland at the head of the alluvial fan. On the other hand, a positive correlation coefficient of 0.66 is obtained at the foot of the alluvial fan, based on the Hefeng groundwater station and the YSLL GNSS station (Figure 10b). Under the complex hydrogeological conditions, the positive and negative relationships between the groundwater level change and the vertical surface displacement can be explained by the effects of pore water pressure and water mass loading. The idea is that with the groundwater level rising, Under the complex hydrogeological conditions, the positive and negative relationships between the groundwater level change and the vertical surface displacement can be explained by the effects of pore water pressure and water mass loading. The idea is that with the groundwater level rising, the increased pore water pressure causes expansion of the sediment formation, whereas the increased water mass loading leads to compaction of the sediment formation. The overall amount of expansion and compaction is reflected in the elevation change of the land surface. Thus, a conceptual model is proposed to explain the observed phenomenon ( Figure 11). Chiang et al. [64] proposed a similar model for interpreting the phenomenon. Given larger storage coefficients of aquifers, the vertical surface displacement is dominated by the effect of water mass loading. The amount of compaction of the bottom sediment due to the increase of water mass loading is much greater than the amount of expansion of the aquifer caused by the increase of pore water pressure (Figure 11a). The correlation between the vertical surface displacement and the groundwater level change is negative, because the effect of water mass loading predominated over the effect of pore water pressure. Thus, the decrease of the groundwater level is supposed to reflect the decrease of the water mass loading, which should result in an uplift of land surface. On the other hand, when the storage coefficients of the aquifers are smaller, the vertical surface displacement is controlled mainly by the pore water pressure. The amount of expansion of the aquifer caused by the effect of pore water pressure is greater than the amount of compaction of the bottom sediment due to the effect of water mass loading (Figure 11b). In the condition dominated by pore water pressure, the correlation between the vertical surface displacement and the groundwater level change is positive. This indicates that the increase of the groundwater level is accompanied by the uplift of the land surface.
Remote Sens. 2020, 12, x FOR PEER REVIEW 17 of 23 the increased pore water pressure causes expansion of the sediment formation, whereas the increased water mass loading leads to compaction of the sediment formation. The overall amount of expansion and compaction is reflected in the elevation change of the land surface. Thus, a conceptual model is proposed to explain the observed phenomenon ( Figure 11). Chiang et al. [64] proposed a similar model for interpreting the phenomenon. Given larger storage coefficients of aquifers, the vertical surface displacement is dominated by the effect of water mass loading. The amount of compaction of the bottom sediment due to the increase of water mass loading is much greater than the amount of expansion of the aquifer caused by the increase of pore water pressure (Figure 11a). The correlation between the vertical surface displacement and the groundwater level change is negative, because the effect of water mass loading predominated over the effect of pore water pressure. Thus, the decrease of the groundwater level is supposed to reflect the decrease of the water mass loading, which should result in an uplift of land surface. On the other hand, when the storage coefficients of the aquifers are smaller, the vertical surface displacement is controlled mainly by the pore water pressure. The amount of expansion of the aquifer caused by the effect of pore water pressure is greater than the amount of compaction of the bottom sediment due to the effect of water mass loading (Figure 11b). In the condition dominated by pore water pressure, the correlation between the vertical surface displacement and the groundwater level change is positive. This indicates that the increase of the groundwater level is accompanied by the uplift of the land surface.

Hydrogeological Implication of the Interaction Model for the Alluvial Fan of the Choshui River
Groundwater is an important resource in alluvial fan areas, where the population density is usually high. The positive or negative interactions of the groundwater level change and the surface displacement are important for water resources management in alluvial fan areas. The detailed spatial correlation between the groundwater level change and the surface displacement is also critical for properly judging and utilizing groundwater resources. The density of the PS pixels derived from Sentinel-1A images by the MTI technique is much higher than the density derived from the GNSS stations. The time series data derived from the PS pixels are applied to represent the vertical surface displacement around the groundwater stations, which have no GNSS station nearby. The PS pixels are selected within a 1.5 km radius of the groundwater monitoring wells for calculating the correlation coefficient. The spatial distribution of the correlation coefficient between the vertical surface displacement and the groundwater level change is shown in Figure 12. One groundwater station may record the groundwater level change at different aquifer depths. We used the records of

Hydrogeological Implication of the Interaction Model for the Alluvial Fan of the Choshui River
Groundwater is an important resource in alluvial fan areas, where the population density is usually high. The positive or negative interactions of the groundwater level change and the surface displacement are important for water resources management in alluvial fan areas. The detailed spatial correlation between the groundwater level change and the surface displacement is also critical for properly judging and utilizing groundwater resources. The density of the PS pixels derived from Sentinel-1A images by the MTI technique is much higher than the density derived from the GNSS stations. The time series data derived from the PS pixels are applied to represent the vertical surface displacement around the groundwater stations, which have no GNSS station nearby. The PS pixels are selected within a 1.5 km radius of the groundwater monitoring wells for calculating the correlation coefficient. The spatial distribution of the correlation coefficient between the vertical surface displacement and the groundwater level change is shown in Figure 12. One groundwater station may record the groundwater level change at different aquifer depths. We used the records of the groundwater level change mostly at the depths of around 150-200 m, where high-quality seasonal and low-frequency signals are usually observed. Figure 12. Calculated spatial distribution of the correlation coefficient between the vertical surface displacement derived from GPS and PSInSAR and the groundwater level change derived from monitoring wells. Due to local geological conditions, areas with larger values of storage coefficients, such as the head of the alluvial fan, the mountain area, and area near the Choshui River, show negative relationships between the vertical surface displacement and the groundwater level change. Squares represent the vertical surface displacement derived from GPS data, whereas circles represent the vertical surface displacement based on PS pixels, which were selected in a radius of about 1.5 km around the groundwater stations. The background colors indicate land elevation with four color ramps from 0-20, 20-45, 45-120, and 120-200 m.

Conclusions
The management of groundwater resources is a critical issue in highly populated alluvial fan areas. Due to overpumping, serious land subsidence may occur, causing damages in the traffic infrastructure and buildings. In the alluvial fan of the Choshui River in central Taiwan, the spatial variation of land subsidence from 1995 to 2017 is revealed based on SAR images, precise leveling, and GPS data. Our combined analysis indicates that this land subsidence was likely triggered by groundwater level change and the sediment consolidation varied with the local geological sediment conditions. In order to interpret the combined effects of the pore water pressure and the water mass loading on the vertical surface displacement, a seasonal interaction model of the groundwater level change and the vertical surface displacement is proposed here. In areas with small values of the storage coefficient of the aquifers, such as at the foot of the alluvial fan, pore water pressure plays a Figure 12. Calculated spatial distribution of the correlation coefficient between the vertical surface displacement derived from GPS and PSInSAR and the groundwater level change derived from monitoring wells. Due to local geological conditions, areas with larger values of storage coefficients, such as the head of the alluvial fan, the mountain area, and area near the Choshui River, show negative relationships between the vertical surface displacement and the groundwater level change. Squares represent the vertical surface displacement derived from GPS data, whereas circles represent the vertical surface displacement based on PS pixels, which were selected in a radius of about 1.5 km around the groundwater stations. The background colors indicate land elevation with four color ramps from 0-20, 20-45, 45-120, and 120-200 m.
In most areas with a smaller storage coefficient of about 0.0019 [64], indeed the correlation between the vertical surface displacement and the groundwater level change is positive. The results show that the negative correlation between the vertical surface displacement and the groundwater level change appears in areas such as the top of the alluvial fan, the mountain area, and the area close to the Choshui River. In these areas, the storage coefficient is larger, about 0.15 [64], due to the geological conditions. Still, few nearby locations show the opposite relationship between the vertical surface displacement and the groundwater level change, which may be caused by local geologic variations. Land subsidence has occurred due to the groundwater withdrawal, both at the head or the foot of the alluvial fan, in the time scale of two decades. At the yearly time scale, correlations between vertical surface displacement and groundwater change are positive in some cases and negative in others. On the head of the alluvial fan, this relationship is positive in the time scale of two decades, yet it is negative in the yearly time scale. Consequently, compaction exceeds rebound because the alluvial fan sequence in the subsurface is not fully elastic. Thus, due to the groundwater level drop, the land subsided consistently at the time scale of two decades. In contrast, the correlation between the vertical surface displacement and the groundwater level change in the yearly time scale remained negative in areas with larger storage coefficients. From the viewpoint of the utilization and sustainability of water resources in an alluvial fan, pumping groundwater from areas with large storage coefficients might be a better option, not only to satisfy the demand for water resources but also to mitigate land subsidence.

Conclusions
The management of groundwater resources is a critical issue in highly populated alluvial fan areas. Due to overpumping, serious land subsidence may occur, causing damages in the traffic infrastructure and buildings. In the alluvial fan of the Choshui River in central Taiwan, the spatial variation of land subsidence from 1995 to 2017 is revealed based on SAR images, precise leveling, and GPS data. Our combined analysis indicates that this land subsidence was likely triggered by groundwater level change and the sediment consolidation varied with the local geological sediment conditions. In order to interpret the combined effects of the pore water pressure and the water mass loading on the vertical surface displacement, a seasonal interaction model of the groundwater level change and the vertical surface displacement is proposed here. In areas with small values of the storage coefficient of the aquifers, such as at the foot of the alluvial fan, pore water pressure plays a predominant role in determining the vertical surface displacement, resulting in a positive relationship between the vertical surface displacement and the groundwater level change. On the other hand, in areas where the storage coefficient of the aquifers is large, such as at the head of the alluvial fan, the mountain area, and the area near the Choshui River, the vertical surface displacement is dominated by the effect of the water mass loading. The water mass loading results in a negative relationship between the vertical surface displacement and the groundwater level change. However, in the time scale of two decades, the correlation between the vertical surface displacement and the groundwater level change is consistently positive in all areas. This may have been that the amount of compaction exceeds rebound since the alluvial fan sequence in the subsurface is not fully elastic. In this study, we combined monitoring data from various sources, not only to overcome the limitations associated with the single technique but also to provide the spatial and temporal patterns of the aquifer-system response and gain insight into the surface behavior. This study suggests that pumping groundwater in areas with large storage coefficients within an alluvial fan may be a better option for balancing the demand of groundwater resources and mitigating land subsidence.