A fusion model used in subsidence prediction in Taiwan

The Taiwan Water Resources Agency uses four techniques to monitor subsidence in Taiwan, namely data from leveling, global positioning system (GPS), multi-level compaction monitoring wells (MCMWs), and interferometry synthetic aperture radar (InSAR). Each data type has advantages and disadvantages and is suitable for different analysis tools. Only MCMW data provide compaction information at different depths in an aquifer system, thus they are adopted in this study. However, the cost of MCMW is high and the number of MCMW is relatively low. Leveling data are thus also adopted due to its high resolution and accuracy. MCMW data provide compaction information at different depths and the experimental data from the wells provide the physical properties. These data are suitable for a physical model. Leveling data have high monitoring density in spatial domain but lack in temporal domain due to the heavy field work. These data are suitable for a blackor grey-box model. Poroelastic theory, which is known to be more conscientious than Terzaghi’s consolidation theory, is adopted in this study with the use of MCMW data. Grey theory, which is a widely used grey-box model, is adopted in this study with the use of leveling data. A fusion technique is developed to combine the subsidence predicted results from poroelastic and grey models to obtain a spatially and temporally connected two-dimensional subsidence distribution. The fusion model is successfully applied to subsidence predictions in Changhua, Yunlin, Tainan, and Kaohsiung of Taiwan and obtains good results. A good subsidence model can help the government to make the accurate strategies for land and groundwater resource management.


Introduction
Subsidence is a worldwide hazard (Galloway et al., 1999;Chen and Rybczyk, 2005;Galloway and Burbey, 2011).The subsidence problem in Taiwan is suspected to be due to groundwater over pumping for aquaculture and irrigation along its coastal plain.The Taiwan Water Resources Agency (TWRA) uses four types of monitoring system to investigate land subsidence in Taiwan, namely data from leveling, global positioning system (GPS), multi-level compaction monitoring wells (MCMWs), and interferometry synthetic aperture radar (InSAR) (Hung et al., 2010).Leveling and MCMW data are adopted in this study due to their relatively long monitoring periods.MCMWs have several magnetic rings anchored at various depths.Variations are detected using the magnetic rings about once a month to measure the compaction in individual formations.The accuracy of the MCMW data is 0.1-0.5 cm.Leveling is a commonly used method for measuring height variations on a land surface by using backsight and foresight readings from leveling rods.The accuracy of leveling is 0.5-1.0cm.Although accuracy is high, measurement efficiency is low.The measurement period is about one to two years in Taiwan.
There are three types of model named white box, grey box, and black box models.White box models are purely theoretical, black box models have no model form, and grey box models are between white box and black box models.A physical model (white box model) avoids the incorrect behavior of a grey box model, whereas a grey box model avoids the assumptions and uncertainty of a physical model (white box model).Therefore, a technique that can fuse the results of physics-based and grey-box-based models would make the prediction of land subsidence more reasonable and reliable (Wang, 2015).Wang (2015) proposed a technique that can fuse two types of time series for subsidence prediction results in both the spatial and temporal domains and applied to subsidence fusion in Jhuoshuei River Alluvial Fan, Taiwan.The fusion method uses a slope criterion and a distance weighting factor to obtain a spatially and temporally connected twodimensional subsidence distribution by upgrading two onedimensional models.Wang et al. (2015) adopted this fusion technique to combine the physics-based nonlinear poroelastic model (NPM) and the grey-box-based grey system model (GSM) to evaluate the subsidence under the climate change effect in Tainan, Taiwan and got good results.The authors have also applied the fusion technique to predict the subsidence in Kaohsiung, Taiwan, which is not yet published.

Study area
The fusion technique used for subsidence prediction was successfully applied to Changhua and Yunlin counties and Tainan and Kaohsiung cities in Taiwan.In this study, only the results for Jhuoshuei River Alluvial Fan (including Changhua and Yunlin counties) and Tainan City are shown (Wang, 2015;Wang et al., 2015).There are less data in Kaohsiung and they are not discussed here.
Jhuoshuei River Alluvial Fan, located in central Taiwan (see Fig. 1), is the largest alluvial fan in Taiwan.Wu River and Beigang River flow through the north and south areas, respectively, and Jhuoshuei River flows across the middle of the alluvial fan.Changhua and Yunlin counties are located in the north and south divisions of Jhuoshuei River Alluvial Fan, respectively.The western boundary is the Taiwan Strait and the eastern boundary is Bagua and Douliu foothills, where are the main areas for groundwater recharge.
Tainan is located in southwestern Taiwan.The Tainan area is in the south central division of the Chianan Plain.Tainan is bounded by the Bajhang River in the north, the Erren River in the south, the Chiayi and Xinhua foothills in the east, and the Taiwan Strait in the west.The main geography is the Tainan Plain in the western part and some foothills and mountainous areas in the eastern part, as shown in Fig. 2. 2007.However, the compaction quantity at different depths can be different.

Subsidence data
There are only four MCMWs in Tainan due to the subsidence being relatively low.The depths of the MCMWs are all at 300 m.Since the wells were installed in different years, their monitoring periods are different.The compaction at the ground surface continuously increases at the Xiaying, Xuejia, and Chengda wells.The land surface at the Iian well shows an initial decrease and then recovery.The recovery behavior is due to an increase in groundwater level around this area and shows an elastic response behavior.Since this study focuses on subsidence evaluation, the subsidence of the Iian well is set to zero in the entire area estimation.

Nonlinear poroelastic model
The governing equations for the one-dimensional poroelastic model can be written as (Wang and Hsu, 2009a): where w is the vertical displacement; P is the change in pore water pressure; z is the vertical direction; t is time; a is called the final compressibility with a −1 = λ + 2µ, where  Following Wang and Hsu (2009b), the porosity (n), hydraulic conductivity (K), and Young's modulus (E) under the deformation effect can be written, respectively, as: where ε is the volume strain, and the subscripts 0 and 1 indicate before and after the deformation effect, respectively.The nonlinear parameters are embedded in the calculation during the modeling period.

Grey system model
Deng (1989) proposed grey system prediction theory.The GSM with one first-order variable is called the GM (1,1) model.The grey differential equation of the GM (1,1) model can be written as: where a and b are coefficients often evaluated using the least squares error method by fitting the theory and observation data and x (1) is a canonical series written as: where x (0) is the original data series (i.e., leveling data series in this study) and N is the data number.Substituting the evaluated a and b into Eq.( 5) and letting the initial condition be x (1) (1) = x (0) (1) yields: where x (1) (n + 1) is the predicted value for x (1) at the n + 1 point.Using back differential for x (1) yields: The predicted value x (0) (n+1) is thus obtained.Equation ( 8) is the basic equation for the GM(1,1) model.

Fusion method
A fusing method that can combine two types of trend series and obtain a spatial-temporal relative subsidence distribution  was proposed by the present author (Wang, 2015).The fusion concept is to use the results of NPM to re-evaluate the results of GSM under a slope criterion of the subsidence trend in the temporal domain and a weighting factor between the monitoring points in the spatial domain.
Considering a distance weighting factor between the locations of GSM i and NPM j , the mean slope of the subsidence trend for NPM j with respect to GSM i can be written as: where D ij is the distance from GSM i to NPM j , m is the number of monitoring points for NPM (i.e., the number of MCMWs), and SNPM j =1,m (t, D) is the mean slope of NPM j with respect to GSM i in the time interval dt.Assuming that the subsidence trends obtained from the NPM results and the initial subsidence quantities obtained from the GSM results are correct, the initial subsidence quantity of GSM is kept, but the slope of the subsidence trend obtained from the GSM is adjusted using the mean slope calculated from NPM.If the slope of the subsidence trend obtained from GSM i in the time interval dt is larger than the mean slope of NPM j in the same time interval, the GSM i slope is substituted with the mean NPM j slope.The slope criterion can be written as: where S GSM i and S NPM j represent the slopes of the subsidence trend for GSM i and NPM j , respectively.The criterion is used to evaluate the slope of the subsidence trend in each time interval for each leveling point with respect to all MCMWs by considering the distance weighting factor.Both the spatial and temporal results are combined in the fusion technique.The independent one-dimensional evaluations in the NPM and GSM are then upgraded to a spatially and temporally connected two-dimensional distribution.

Results for NPM
The MCMWs with maximum and minimum subsidence in Changhua and Yunlin, respectively, are chosen as examples for illustrating the fitting situation, as shown in Fig. 3.The Xinghua and Xizhou MCMWs are in Changhua County and the Yunchang and Jiaxing MCMWs are in Yunlin County.The subsidence fluctuation is more dominant in Yunlin County (Fig. 3c and d) than in Changhau County (Fig. 3a and b) due to the large fluctuation of groundwater level in Yunlin County.The fitting for Jiaxing well is poor due to the negative subsidence and the adopted compaction NPM, as shown in Fig. 3d.Negative subsidence indicates an uplift situation, which is not considered in the NPM in this study.From the fitting results, the quantity of discharge is 0 to 2.58 × 10 −9 m s −1 and that of loading is 0 to 1.07 × 10 4 N m −2 .Two and 15 MCMWs show zero discharge and zero loading, respectively.The main driving force for the land subsidence in Jhuoshuei River Alluvial Fan  is groundwater pumping, which is expressed as discharge (pumping rate in unit area).The influence of top loading is relatively small, occurring only at 11 MCMWs.The quantities of loading are also small, lower than atmospheric pressure (= 1.013 × 10 5 N m −2 ).The best fits for the compaction data from the wells in Tainan and the NPM are shown in Fig. 4. The fig- ures show that the fits are excellent at the Xiaying and Xuejia wells with correlation coefficients (R-squared) of about 0.9, good at the Iian well, but poor at the Chengda well due to both the hydraulic conductivity and Young's modulus being large.The additional loading is from 0 to 2.03 × 10 3 N m −2 , and the additional discharge is from pumping rate can be calculated by multiplying the additional discharge by the influence area.The quantities of additional pumping rate around the monitoring wells in Tainan are from −36 to 223 m 3 day −1 .However, the evaluated additional pumping rate is affected by the influence area and could be different in the field.

Subsidence prediction
The slope criterion for the subsidence trend series is adopted to re-calculate results of the GSM by considering the mean slope calculated from the NPM with a distance weighting factor.The fusion results for subsidence distribution in Jhuoshuei River Alluvial Fan for 2039 are shown in Fig. 5.The maximum subsidence decreases from 2.31 to 0.59 m.The areas with subsidence quantities larger than 0.3 m are mainly located in the south-west areas of Changhua County and south-central areas of Yunlin County.Comparisons of the results obtained with and without fusion for two leveling points are shown in Fig. 6.The locations of the two leveling points are marked in Fig. 5.The subsidence predicted using the GSM for the leveling points shows a linear increasing trend.The subsidence quantities for points #1 and #2 are 0.71 and 1.23 m, respectively, for 2039.After using the fusion technique, both subsidence curves show convergence trends.The subsidence quantities are 0.55 and 0.53 m, respectively, and the differences are 0.16 and 0.70 m, respectively.The decreases vary with leveling point.The subsidence behavior becomes reasonable for both the subsidence quantity and subsidence rate.The subsidence distributions obtained by combining the NPM and GMS in the years 2020 and 2039 in Tainan are shown in Fig. 7.The results show that the main areas with large subsidence move from Xuejia and Beimen Districts to Yanshui and Xiaying Districts, as evidenced by the large set-tlement increment at the Xiaying well.The pattern in the subsidence areas obtained by combining the NPM and GSM is similar to that obtained using only the GSM.However, the largest subsidence areas are different, and the subsidence quantity is much lower for the former due to the constraint of the physics-based NPM.The similar pattern might be due to the limited number of MCMWs (four) in Tainan, and thus there is relatively low influence exhibited in the combination.However, with more MCMWs, the combination could produce better results.

Conclusion
A physics-based NPM and a grey-box-based GSM were used to investigate land subsidence in Changuha, Yunlin, and Tainan, Taiwan.A technique was proposed to fuse the subsidence results from NPM and GSM in the spatial and temporal domains using a slope criterion and a distance weighting factor.The linear increasing trend of subsidence predicted using a grey-box-based GSM was corrected to a convergence trend with a physics-based NPM.The fusion results show dramatic decreases of subsidence quantity and subsidence rate.The subsidence behavior became physically reasonable due to the convergence trend of subsidence under the assumption of constant discharge of groundwater and top loading.The proposed fusion technique upgrades the independent onedimensional evaluations to a spatially and temporally connected two-dimensional distribution.
The parameters in the NPM were taken from laboratory experiments and the data were insufficient.Field or more laboratory experiments will provide the typical values and improve the modeling results.There are few MCMWs installed in Taiwan.More available monitoring wells will provide better predicted subsidence results from the combination of the NPM and GSM.The recent NPM assumed a soil column with uncertainty embedded for subsidence evaluation.Layer system or complex concept could be used to develop a detail model for subsidence evaluation.Therefore, there are uncertainties embedded in the model.More data and further investigations can improve the estimations.

Figure 2 .
Figure 2. River system and geographic map with distribution of districts and monitoring points in Tainan City (From Wang et al., 2015).

Figure 6 .
Figure 6.Comparison of results obtained with and without fusion for two leveling points marked in Fig. 5(From Wang, 2015).

− 1 .Figure 7 .
Figure 7. Subsidence distribution in Tainan predicted using combination of the NPM and GSM for the target years (a) 2020 and (b) 2039 (From Wang et al., 2015).

Wu River Beigang River Jhuoshuei River Yunlin County Changhua County Bagua foothills Douliu foothills Taiwan Strait A A′
Jiaxing Figure 1.Study area of Jhuoshuei River Alluvial Fan, Taiwan, and distribution of leveling points and MCMWs (From Wang, 2015).