Dimensional Upgrade Approach for Spatial-temporal Fusion of Trend Series in Subsidence Evaluation

Physical models and grey system models (GSMs) are commonly used to evaluate and predict physical behavior. A physical model avoids the incorrect trend series of a GSM, whereas a GSM avoids the assumptions and uncertainty of a physical model. A technique that combines the results of physical models and GSMs would make prediction more reasonable and reliable. This study proposes a fusion method for combining two trend series, calculated using two one-dimensional models, respectively, that uses a slope criterion and a distance weighting factor in the temporal and spatial domains. The independent one-dimensional evaluations are upgraded to a spatially and temporally connected two-dimensional distribution. The proposed technique was applied to a subsidence problem in Jhuoshuei River Alluvial Fan, Taiwan. The fusion results show dramatic decreases of subsidence quantity and rate compared to those estimated by the GSM. The subsidence behavior estimated using the proposed method is physically reasonable due to a convergent trend of subsidence under the assumption of constant discharge of groundwater. The technique proposed in this study can be used in fields that require a combination of two trend series from physical and nonphysical models.


Introduction
Subsidence is a worldwide hazard [1][2][3][4].Commonly used methods for studying land subsidence include data analysis [5][6][7][8][9][10], physical calculation [11][12][13][14][15] and numerical simulation [16][17][18][19][20][21][22][23][24][25].Data analysis analyzes the monitoring data, and prediction is mainly used for evaluating the type or trend of subsidence.The method used for prediction can be a black box model (e.g., regression analysis) or grey box model (e.g., grey system model (GSM)).Data analysis can avoid the assumptions in a physical model.However, the prediction results are easily affected by misleading monitoring data and the predicted series might deviate from physical behavior.Physical calculation uses experimental data, monitoring data of subsidence, or groundwater level data to determine subsidence using physics.The approach produces results that have physical support and are close to actual behavior.However, it is commonly used under a regular or simple conditions and uncertainty is introduced in the model assumptions and adopted parameters.Numerical simulation uses a physical theory to represent a complex situation.It is widely used in evaluating subsidence problems.Consolidation theory [26] and poroelastic theory [27] are commonly used physical theories for subsidence evaluation.A physical model avoids the incorrect behavior of a grey box model, whereas a grey box model avoids the assumptions and uncertainty of a physical 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.
Subsidence monitoring is commonly a point type system, with the data obtained in the temporal domain for each point [8].A one-dimensional model is thus commonly used for subsidence investigation due to its simplicity and convenience [28].However, uncertainty is embedded between monitoring points in the spatial domain, and the spatial information of subsidence is lacking between the data points in a one-dimensional model.A multi-dimensional model can be used but it is more complex and timeconsuming.A technique that can connect the one-dimensional points in the spatial domain would be useful.Previous studies on fusing different types of time series have focused on the temporal domain, with few considering the spatial domain [29][30][31][32][33][34][35].The present study thus proposes a technique that can fuse two types of time series for subsidence prediction results in both the spatial and temporal domains.
Poroelastic theory, proposed by Biot [27], is a rigorous physical theory for subsidence investigations [16][17][18][19][22][23][24][25].This theory uses a coupled system that can solve soil deformation and pore water pressure at the same time and obtain the dynamic interaction between the variables [36].Soil deformation can change the hydrogeological and mechanical properties and affect the evaluation of subsidence quantity [37][38][39].Wang and Hsu's model, which considers the deformation effect on porosity, hydraulic conductivity, and Young's modulus, are simpler and more complete than that of Kim and Parizek [37].The nonlinear poroelastic model (NPM) developed by Wang and Hsu [38] has been used in subsidence estimation in Taiwan with good results [24,25].The NPM is suitable for Taiwan subsidence problems and is thus adopted in this study to construct a physics-based subsidence model in Jhuoshuei River Alluvial Fan, Taiwan.
Grey system theory, proposed by Deng [40],is a grey box model that is used to study uncertain systems or systems with incomplete data.The GSM can be used to evaluate land subsidence for which the data set is limited [25].Wang et al. [25] combined the GSM and NPM by using the technique proposed in this study to investigate the land subsidence problem in Tainan, Taiwan.The present study proposes a fusion technique and shows the difference between the results with and without fusion in a subsidence evaluation of Jhuoshuei River Alluvial Fan, Taiwan.The fusion method combines the advantages of physics-based and grey-box-based models and fuses the one-dimensional predicted subsidence results in the spatial and temporal domains.The independent one-dimensional evaluations are then upgraded to a spatially and temporally correlative two-dimensional distribution.

Study Area
Jhuoshuei River Alluvial Fan, located in central Taiwan (see Figure 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, which are the main areas for groundwater recharge.The Taiwan Water Resources Agency executed 93 hydrogeological drillings and setup 70 groundwater monitoring stations with 188 monitoring wells at various depths to investigate the groundwater resources in Jhuoshuei River Alluvial Fan [41].Each groundwater station includes one to five monitoring wells at different depths.The distance between each well is about 5 km.The sampling frequency for groundwater level monitoring is one hour.One of the hydrogeoloical profiles is shown in Figure 2. The hydrogeology formations can be separated into seven units within a depth of 300 m.They are referred to as aquifer one, aquitard one, aquifer two, aquitard two, aquifer three, aquitard three and aquifer four.The layer system is obvious in the mid and tail fans but not obvious in the top fan due to its composite material being mostly gravel.The aquifers thin out beneath the Taiwan Strait and stretch from the tail fan.Detailed hydrogeological descriptions can be found elsewhere [8,21].Jhuoshuei River Alluvial Fan contains a great amount of groundwater.However, there is little surface water in this area due to industrialization and urbanization.Groundwater is thus mainlyover-pumped for public, irrigation and industrial usage, which has induced subsidence.According to an investigation by the Taiwan Water Resource Agency, there are about 182,000 wells in Jhuoshuei River Alluvial Fan, with a total pumping quantity of around 2.16 billion m 3 /year [8,21].The maximum subsidence rate is 6.1 cm/year and the continuous subsidence area (subsidence rate larger than 3.0 cm/year) was about 309.1 km 2 in 2014 [42].Taiwan High Speed Rail passes across the subsidence areas, which is a safety concern [43].Subsidence in Jhuoshuei River Alluvial Fan is thus an important issue for the Taiwanese government.

Subsidence Monitoring
Four types of monitoring systems are used in Taiwan for subsidence investigation, namely leveling, Global Positioning System (GPS), multi-level compaction monitoring wells (MCMWs), and interferometry synthetic aperture radar (InSAR).Leveling and MCMW data are adopted in this study due to their relatively long monitoring periods.MCMW data provide compaction information at various depths, and the experimental data from the wells provide the physical properties.MCMW data are suitable for the study of NPMs.Leveling monitoring has been used for a long period for subsidence observation in Taiwan.The monitoring density is high and can compensate for the lack of MCMWs.The leveling data are suitable for GSMs.Detailed descriptions of the monitoring system can be found elsewhere [8].
There were 29 MCMWs and 708 leveling points in 2011 in Jhuoshuei River Alluvial Fan, Taiwan, as shown in Figure 1.The measurement frequencies are monthly for the MCMWs and yearly for leveling in this area.The land subsidence data from 2007 to 2011 are adopted, and thus the subsidence quantities shown in this study are relative to the level in 2007.Only    From the 26 monitoring wells, the subsidence quantity at the ground surface is the largest at Yuanchang well (= 0.252 m) and the smallest is at Jiaxing well (=0 m) in 2011 based on 2007.However, the compaction quantity at different depths can be different.Figure 3 shows the monitoring data for two of the MCMWs as anillustration.Figure 3a shows the compaction in the spatial domain at various periods at Xinghua well, which has the largest land subsidence quantity in Changhua County.The slope of the curve represents the compaction rate.Therefore, the main compacted formation in 2011 is at about 100-250 m, contributing 80% of total land subsidence.Figure 3b shows the compaction at various depths at Yuanchang well, which has the largest land subsidence quantity in Yunlin County.The compaction at different depths of Yuanchang well is more uniform than that of Xinghua well.The main compacted formation in 2011 was at about 50-280 m, contributing 89% of total land subsidence.

Nonlinear Poroelastic Model
The nonlinear poroelastic model (NPM) was developed by Wang and Hsu [38].The governing equations for the one-dimensional poroelastic model can be written as [44]: 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 λ and μ are Lame's constants (μ = E/2(1 + ν), and λ = Eν/(1 + ν) (1 − 2ν), where E is Young's modulus, and ν is Poisson's ratio); α is the dimensionless coefficient of the effective stress; Q −1 is the compressibility introduced by Biot; and κ is the Darcy conductivity, which is related to hydraulic conductivity K by K = γwκ, where γw is the unit weight of water (= 9810 N/m 3 ).
Following Wang and Hsu [38], 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
The GSM is based on the eigenvalue system and the model established based on the grey prediction theory.During the grey simulation, the accumulated generating operation is applied to the discrete and irregular initial number series to make the series generate an exponential rule, based on which the grey differential equation is established.The GSM can be established using data preprocessing and the grey differential equation.Generally, a method based on the time series model is used to process irregular data in order to smooth the series and increase prediction precision.
The leveling data are limited due to the sampling rate being yearly in Jhuoshuei River Alluvial Fan, Taiwan, and the measurement points being different in each year.Most of the leveling data have only four subsidence measurement points in the temporal domain.Therefore, the GSM is a suitable model for subsidence prediction with the use of leveling data.
Deng [45] proposed the 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: (1) (1) 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 Equation ( 5) and letting the initial condition be x (1) (1) = x (0) (1) yields: where  (1) ( 1) x n+ is the predicted value for x (1) at the n + 1 point.Using back differential for  (1)  x yields: The predicted value  (0) ( 1) x n+ is thus obtained.Equation ( 8) is the basic equation for the GM (1,1) model.The GSM developed by Hung et al. [46] is used to construct a GM (1,1) grey prediction model of subsidence in this study.Four data points of a leveling time series are input into the GSM model, and then the point for the next time step is predicted.Considering the time interval and subsidence quantity of the leveling data, the subsidence time series is evaluated by the GSM at each leveling point, as shown in Figure 1.After the subsidence time series for all the leveling points have been evaluated, the results of the GSM are obtained.

Model Settings
Land subsidence is mainly caused by surface loading and groundwater pumping.Considering both the mechanism of top loading and aquifer pumping, the conceptual model for a one-dimensional column is proposed [24].The quantities of top loading and bottom discharge are evaluated by fitting the results of the NPM to the monitoring data.The adopted parameters for the NPM are listed in Table 2.The property at each MCMW in this study is assumed to be homogeneous in the vertical direction with the mean value of the parameters.The mean porosity and hydraulic conductivity are obtained from triaxial permeability tests.The Young's modulus and Poisson's ratio are taken from Das [47].However, since some MCMWs do not have experimental data, mean values from the other MCMWs were used in the model.The settings of the NPM are similar to those in Wang et al. [24] except for the adopted parameters, and those of the GSM are the same as those in Wang et al. [25] except for the monitoring data.Detailed descriptions of the NPM and the GSM for subsidence evaluation can be found elsewhere [24,25].

Fusion Method
The subsidence predicted using the GSM shows a continuous increasing (or decreasing) trend due to the use of monitoring data from leveling, whereas that predicted using the NPM shows a convergence trend due to physical behavior.Based on the physical mechanism, the subsidence rate should decrease and tend to be stable if the driving force is constant.Therefore, the converging results of the NPM are reasonable, and thus the slopes in each time step of the convergent curve are adopted for the combination.The subsidence predicted using the GSM follows the monitoring data.The subsidence quantities in the initial period evaluated by the GSM should be reliable and are thus adopted for the combination.The number of leveling points is larger than the number of MCMWs in Jhuoshuei River Alluvial Fan, Taiwan.The fusion results thus can provide a better estimation than using each of the models.
The fusion concept is to use the results of the NPM to re-evaluate the results of the 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.Figure 4 shows a schematic map of the distribution of monitoring points (one for the GSM and three for the NPM).The distances from GSM1 to NPM1, NPM2, and NPM3 are d11, d12, and d13, respectively (i.e., dij is the distance from GSMi to NPMj).The slopes of the subsidence trend for GSMi and NPMj between time t and t − 1 can be written, respectively, as: where GSMi and NPMj are the evaluated subsidence values for the GSM at position i and the NPM at position j, respectively.Considering a distance weighting factor between the locations of GSMi and NPMj, the mean slope of the subsidence trend for NPMj with respect to GSMi can be written as: where is the mean slope of NPMj with respect to GSMi in the time interval dt and m is the number of monitoring points for the NPM (i.e., the number of MCMWs).A distance weighting between the points for the GSM and the NPM is considered in calculating the mean slope of the NPM.The influence decreases with increasing distance between the GSM and the NPM.
It is assumed that the subsidence trends obtained from the results of the NPM and the initial subsidence quantities obtained from the results of the GSM are correct.The initial subsidence quantity of the GSM is kept but the slope of the subsidence trend obtained from the GSM is adjusted using the slope calculated from the NPM.If the slope of the subsidence trend obtained from GSMi in the time interval dt is larger than the mean slope of NPMj in the same time interval, the slope of GSMi is substituted by the mean slope of NPMj.The slope criterion can be written as: The slope 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 the GSM are then upgraded to a spatially and temporally relative two-dimensional distribution.

Results for Nonlinear Poroelastic Model
The MCMWs with maximum and minimum subsidence in Changhua and Yunlin, respectively, are chosen as examples for illustrating the fitting situation, as shown in Figure 5.The Xinghua and Xizhou MCMWs are in Changhua County and the Yunchang and Jiaxing MCMWs are in Yunlin County.The land subsidence in these four MCMWs shows a continuous increasing trend and fluctuation, which is due to the recurrent dry and wet seasons.The subsidence fluctuation is more dominant in Yunlin County (Figure 5c,d) than in Changhau County (Figure 5a,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 Figure 5d.Negative subsidence indicates an uplift situation, which is not considered in the NPM in this study.From the data trend of Jiaxing well, the original land surface should be lower than the initial monitoring data in May 2008.The large fluctuation of the monitored land subsidence and the lack of the original land surface has led to the negative subsidence situation.If the original land surface can be found, the poor fitting situation can be improved.From the results of the NPM, the subsidence trend in the calibration period shows both convergence and divergence.However, the land subsidence approaches an asymptotic value for a long period due to the constant driving forces.Based on physical behavior, subsidence should converge under a constant driving force in the steady-state condition.
The verification results of the NPM are shown in Figure 5.The subsidence data for 2012-2014 are adopted for model verification.However, the data from March 2012 to March 2014 are missing due to an equipment operation error.From the verification results in Figure 5, the subsidence predicted using the NPM in Changhau County shows an overestimated pattern (Figure 5a,b) and that in Yunlin County shows an underestimated pattern (Figure 5c,d).The verification results are consistent with the monitoring results, which showed that land subsidence (groundwater level) in Changhua is decreasing (increasing) and that in Yunlin is still increasing (decreasing) [42].Therefore, the predicted subsidence in the following section might be overestimated in Changhua County and underestimated in Yunlin County.
The best fits for the subsidence data from the MCMWs and the numerical results from the NPM are listed in Table 3.The fitting R 2 values are good for the MCMWs except for the Jiaxing well.The reason for the poor fitting has been mentioned previously and is shown in Figure 5d.From the fitting results, the quantity of discharge is 0 to 2.58 × 10 −9 m/s 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) in Table 3.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 ).

Results for Grey System Model
The predicted subsidence of the GSM was calculated time period by time period and point by point.The predicted subsidence distribution in 2039 is shown in Figure 6, as obtained by applying the interpolation method.The main land subsidence areas are located in the central areas of Yunlin and Changhua Counties.The pattern is similar to that of the monitoring data (not shown here) due to the nearly linear increment of the GSM.The maximum subsidence estimated using the GSM (2.31 m) occurs in Yunlin County.The maximum subsidence rate is 7.4 cm/year, which is larger than the rate for serious hazard potential of subsidence (5.0 cm/year) [48].
A relatively large subsidence for 2039 is obtained using the GSM due to the linear increasing trend of the monitoring data.Based on physical behavior, soil consolidation has a convergence trend if the exerted stress and/or pore water pressure changes are constant.A linear increasing subsidence trend is thus unreasonable under an unchanging driving force.Therefore, the predicted subsidence quantities shown in Figure 6 are overestimated.The fusion technique is thus applied to combine the two subsidence trend series from the NPM and the GSM in the spatial and temporal domains to obtain a reasonable and reliable estimation.

Subsidence Fusion Results
The slope criterion for the subsidence trend series (Equation ( 12)) is adopted to re-calculate the 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 Figure 7.Note that the subsidence distribution in Figure 7 is not an interpolated result but the fusion result.The pattern of the subsidence distribution is slightly different from that in Figure 6 but the quantity has dramatically decreased.The maximum subsidence decreases from 2.31 m to 0.59 m.The areas with subsidence quantities larger than 0.3 m are mainly located in the southwest 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 Figure 8.The locations of the two leveling points are marked in Figures 6 and 7.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 maximum subsidence is 0.59 m and the maximum subsidence rate is 1.0 cm/year for 2039.The subsidence behavior becomes reasonable for both the subsidence quantity and subsidence rate.Note that the subsidence trend for Fusion #2 in Figure 8 still increases with time.This fusion result indicates that the cumulative subsidence in Jhuoshuei River Alluvial Fan will increase in some areas after 2039 but will converge in the future, due to the application of physical theory.The subsidence predictions in this study include some assumptions and uncertainties.A one-dimensional NPM is adopted and lateral deformation is not considered in the estimation, which might have an influence on the predicted subsidence quantities.The adopted parameters were from laboratory tests and a previous study and thus might not represent the behavior in the field.The subsidence monitoring data may also have measurement errors.Further investigations are required for better estimation of land subsidence.Nevertheless, although the proposed subsidence models have their limitations and disadvantages, the fusion method can still work.The defects of the subsidence models do not affect the use of the proposed fusion technique.An improved fusion result can be obtained after the predicted subsidence results from the models being improved.

Conclusions
This study proposed a technique that fuses two trend series in the spatial and temporal domains using a slope criterion and a distance weighting factor and then applied it to investigate subsidence in Jhuoshuei River Alluvial Fan, Taiwan.The linear increasing trend of subsidence predicted using a grey-box-based GSM was corrected to a convergence trend with a physics-based poroelastic model.The fusion results show dramatic decreases of subsidence quantity and subsidence rate.The maximum subsidence quantity decreased from 2.31 to 0.59 m and the subsidence rate decreased from 7.4 to 1.0 cm/year for 2039.The subsidence behavior became physically reasonable for both the subsidence quantity and subsidence rate 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 one-dimensional evaluations to a spatially and temporally connected two-dimensional distribution.The proposed method can also improve the subsidence distribution prediction when there are few subsidence monitoring wells.The proposed fusion method can be used for areas that require a combination of trend series in the spatial and temporal domains.

Figure 1 .
Figure 1.Study area of Jhuoshuei River Alluvial Fan, Taiwan, and distribution of leveling points and multi-level compaction monitoring wells.A-A' is the line representing the hydrogeological profile shown in Figure 2.

Figure 2 .
Figure 2. One hydrogeological profile of Jhuoshuei River Alluvial Fan, Taiwan.The profile trace is shown in Figure 1.

3039
Jianyang and Jinhu wells, which were installed in December 1994 at a depth of 200 m.Three recently installed wells (in August 2009), namely Dongguang, Yiwu and Zhengmin wells, are not adopted in this study due to their short monitoring period.Zhengmin well has the largest monitoring depth (=330 m).

Figure 3 .
Figure 3. Monitoring data from MCMWs at (a) Xinghua and (b) Yuanchang wells in spatial domain at various periods.

Figure 4 .
Figure 4. Schematic map of distribution of monitoring points.Grey system models (GSM)i and NPMj indicate monitoring points for grey system and poroelastic models, respectively.d11, d12, and d13 are distances from GSM1 to NPM1, NPM2, and NPM3, respectively.

Figure 6 .
Figure 6.Subsidence distribution in Jhuoshuei River Alluvial Fan in 2039 obtained using GSM.Triangle symbols numbered #1 and #2 are locations of demonstration points used for fusion results in Figure 8.

Figure 7 .
Figure 7. Subsidence distribution with fusion in Jhuoshuei River Alluvial Fan in 2039 obtained using NPM and GSM.Triangle symbols numbered #1 and #2 are locations of demonstration points used for fusion results in Figure 8.

Figure 8 .
Figure 8.Comparison of results obtained with and without fusion for two leveling points marked in Figure 7.
26 MCMWs and 452 leveling points are used in this study due to a lack of monitoring data.The basic information for the 26 MCMWs is listed in Table 1.The monitoring depth is 200 to 330 m.The earliest installed monitoring wells are Haifeng, Lunfeng,

Table 1 .
Basic information for multi-level compaction monitoring wells (MCMWs) adopted in this study.
* Stations not adopted in this study due to short monitoring period.
[47]draulic conductivity value of 8.17 × 10 −8 m/s and porosity value of 0.449 were calculated their means for other MCMWs.bValues of Young's modulus and Poisson's ratio were taken from Das[47].

Table 3 .
Calibration results for NPM at MCMWs.